Influence of spatial interpolation methods for climate variables on the simulation of discharge and nitrate fate with SWAT

For ecohydrological modeling climate variables are needed on subbasin basis. Since they usually originate from point measurements spatial interpolation is required during preprocessing. Different interpolation methods yield data of varying quality, which can strongly influence modeling results. Four interpolation methods to be compared were selected: nearest neighbour, inverse distance, ordinary kriging, and kriging with external drift (Goovaerts, 1997). This study presents three strategies to evaluate the influence of the interpolation method on the modeling results of discharge and nitrate load in the river in a mesoscale river catchment (∼1000 km2) using the Soil and Water Assessment Tool (SWAT, Neitsch et al., 2005) model: I. Automated calibration of the model with a mixed climate data set and consecutive application of the four interpolated data sets. II. Consecutive automated calibration of the model with each of the four climate data sets. III. Random generation of 1000 model parameter sets and consecutive application of the four interpolated climate data sets on each of the 1000 realisations, evaluating the number of realisations above a certain quality criterion threshold. Results show that strategies I and II are not suitable for evaluation of the quality of the interpolated data. Strategy III however proves a significant influence of the interpolation method on nitrate modeling. A rank order from the simplest to the most sophisticated method is visible, with kriging with external drift (KED) outperforming all others. Responsible for this behaviour is the variable temperature, which benefits most from more sophisticated methods and at the same time Correspondence to: S. van der Heijden (vdheijden@iww.uni-hannover.de) is the main driving force for the nitrate cycle. The missing influence of the interpolation methods on discharge modeling is explained by a much higher measuring network density for precipitation than for all other climate variables.


Introduction
Climate data as input for ecohydrological modeling are usually measured at individual points.For the modeling of larger catchments, these data are needed on subbasin or raster basis, consistently covering the whole study area.Thus, spatial interpolation is necessary to upscale the data from point to area.A large number of different interpolation methods has been proposed over time (for a recent review see e.g.Tveito et al., 2006).From the early, simple nearest neighbour technique or Thiessen polygon method (Thiessen, 1911) via the inverse distance weighting (e.g.Shepard, 1968), interpolation methods evolved to more sophisticated geostatistical methods like kriging in all its variants (Goovaerts, 1997), which are becoming more and more popular.Many other techniques e.g. based on splines (Hutchinson, 1998a,b) or genetic algorithms (Demyanov et al., 1998;Huang et al., 1998) have also been applied recently.
With such abundance of methods to select from the question arises which method a modeler should apply for his particular case.Many studies have been carried out comparing the predictive performance of interpolation methods via cross validation (e.g.Dubois et al., 1998;Vicente-Serrano et al., 2003;Stahl et al., 2006;Hofstra et al., 2008), often on a monthly or annual basis.However, to our knowledge only few studies so far examined the influence of different interpolation methods on the modeling results of rainfall-runoff models on a daily or sub-daily basis concerning hydrological validation (e.g.Kneis and Heistermann, 2009), and no studies concerning the validation of interpolation methods relating to nitrogen loads.This study focuses on four interpolation techniques, namely nearest neighbour, inverse distance weighting, ordinary kriging, and kriging with external drift, and aims to evaluate their influence on daily discharge and monthly nitrate simulations with the ecohydrological model SWAT, when all necessary climate input data is interpolated with each of the four methods.Three strategies of model calibration and application have been developed in order to evaluate the influence of the interpolation techniques.

Study area and data
The investigation area is the upper Leine river catchment with a size of approximately 1000 km 2 (see Fig. 1).It is part of the larger Aller river catchment, which covers most of the south-eastern part of the federal state of Lower Saxony in Northern Germany.The average annual precipitation is about 700 mm and the mean annual temperature is around 8 • C. The elevation is between 533 m on the southern edge of the catchment and 115 m at the gauge Leineturm, which is the outlet of the catchment.
The geology in the study area is dominated by the Mesozoic sedimentary rock strata of the Muschelkalk and the Buntsandstein west and east of the Leine trough.The Keuper, which is found in the subsurface of the actual trough, is mostly covered by loess of several meters thickness.Also, the broad Leine floodplains are filled with the fluvial sediments gravel and coarse sand.Under these climatic and geologic conditions Cambisols, Luvisols and gley soils developed as predominant soils, as well as Rendzinas and Rankers on the more elevated terrain.
The main landuse is rain-fed agriculture, predominantly in the flatter low-lying regions of the Leine valley, while the more sloped terrain is dominated by semi-natural deciduous and mixed forest.The main settlement in the area is the city of Göttingen.
For this study, the investigation area is divided into 25 subcatchments (Fig. 1).Climate data is available from the climate and precipitation network of the German Weather Service (DWD), with a total of 94 climate stations and 378 precipitation gauges for the whole Aller catchment.Only about 10 climate and 40 precipitation stations are near or within the actual study area, but spatial interpolation was done for the larger catchment.Out of the total number of stations about 50 climate stations and 250 precipitation gauges were recording at any one time step, which were available for interpolation.Six climate variables are needed by the selected numerical model and thus were selected for interpolation.These are daily values for the sum of precipitation, minimum and maximum temperature, relative humidity, wind speed, and sunshine duration (which is converted into solar radiation after interpolation).
The digital elevation model is the freely available SRTM grid (USGS, 2004) with 90 m cell size, land use information is taken from the CORINE Land Cover data set (Umweltbundesamt, 2004).Daily discharge measurements are available from five gauges with four of them on the Leine river itself and one on the tributary Garte.Out of the nine available water quality sampling sites with biweekly or monthly measurements only three correspond to a discharge gage and could be used in the study (see Fig. 1).Discharge and water quality data were obtained from the Lower Saxony Water Management, Coastal Defence and Nature Conservation Agency (NLWKN).Due to the large spans of time between the irregular samplings it was impossible to examine nitrate on a daily basis.Since no dependence between nitrate concentration and discharge value was found, measured nitrate concentration values were simply interpolated linearly and converted into monthly nitrate load values.

Interpolation methods
Four established spatial interpolation methods were selected for this study.The two simple methods of nearest neighbour (abbr.NN) and inverse distance weighting (IDW), as well as the more sophisticated methods ordinary kriging (OK) and kriging with external drift (KED).Let {z(u i ),i = 1,...,n} be the set of a measured climate variable at n locations u i .The four interpolation methods are used to estimate the climate variables at unsampled locations u.

Nearest neighbour (NN)
The simplest approach for this task is the Thiessen polygon method.Each unsampled location u is assigned the value of the nearest observation u i .
This method produces unrealistic maps with sharp boundaries between polygons of equal value, but preserves the variance of the observed data.

Inverse distance weighting (IDW)
To avoid sudden jumps in values between neighbouring locations, the variable z can be estimated taking several surrounding observations into account.The estimation is a weighted linear combination, with the weights λ i being inversely proportional to the square distance between sampled location u i and the point of estimation u: (2)

Ordinary kriging (OK)
Kriging techniques are generalized least-squares regression algorithms.Kriging considers the spatial variability based on the semivariogram.The experimental semivariogram γ (h) can be calculated for every time step, but an average experimental semivariogram over all time steps has been found to be sufficient (Haberlandt, 2007).The average semivariogram is obtained by weighting the daily semivariograms with the variance for each day.Then a theoretical semivariogram model is fitted to the average experimental semivariogram.Different theoretical semivariograms were used for different climate data, with the spherical semivariogram used the most.
When the theoretical semivariogram is fitted, it is used to determine weights λ OK i (u) in a way to ensure unbiasedness of the estimator, E {Z OK (u) − Z(u)} = 0, and to minimize the estimation variance, Var{Z OK (u) − Z(u)}.Like with IDW, these weights are needed to estimate the unknown value at an unsampled location u as a linear combination of neighbouring observations:

Kriging with external drift (KED)
While the three former methods are univariate, the fourth method KED is bivariate and uses secondary information to interpolate the climate data.Additionally to the spatial dependence of a variable, the linear relation to another variable can be processed.The basic assumption is that the expected value of the target variable Z(u) is a linear function of an additional variable Y (u): The KED estimate is the same as for OK (see Eq. 3), but the weights are determined as the solution of a system of linear equations (the kriging system) under the employ of the secondary information y(u).For further details the reader is referred to geostatistical textbooks (e.g.Isaaks and Srivastava, 1989;Goovaerts, 1997).

Interpolation software
The six climate variables were interpolated using the program SAINT (Haberlandt, unpublished), which allows the user to chose from several interpolation methods, including the four selected NN, IDW, OK, and KED.In addition to interpolation, SAINT also allows to compare the prediction quality of the methods by crossvalidation.Interpolation was done over the whole Aller river catchment on a grid with 900×900 m 2 cell size (10×10 cells of the SRTM DEM).The program interpolates one time series for every cell, and subsequently aggregates all the cells within a given subbasin to a mean areal time series for the subbasin.Thus, for every method, the spatial interpolation returns one time series for each of the 25 subbasins.This procedure yielded four sets of equally interpolated climate data, one for each interpolation method.Variable is HRU-specific, but all HRUs were calibrated at the same time keeping the ratio to each other constant, thus it represents only a single degree of freedom during calibration.
For OK and KED semivariograms were calculated for each climate variable with a FORTRAN program based on the Geostatistical Software Library (GSLIB, Deutsch and Journel, 1992) using all available stations during the time period from 1960 to 2006.Interpolation was then carried out employing those semivariograms.Since many climate variables correlate with terrain elevation, this was taken as secondary information for KED.

The SWAT model
In this study the influence of the different interpolation methods on modeling results of the Soil and Water Assessment Tool (Neitsch et al., 2005) is examined.The SWAT model is a semi-distributed, process-oriented model capable of simulating runoff, nutrient and other agricultural chemicals dynamics as well as sediment yield in large complex watersheds with varying soils, land use, and management conditions.The simulated water balance comprises interception, evapotranspiration, snow melt, surface runoff, infiltration, soil percolation, lateral flow, groundwater flow, and river routing.Evapotranspiration is calculated here after Penman-Monteith, snow melt with the degree day method, infiltration based on the SCS curve number method, runoff transformation using a surface runoff lag method and flood routing is calculated with the variable storage method.Simulations are carried out on the basis of hydrotopes (Hydrologic Response Units, HRUs), which are characterized by unique combinations of land use, soil, and topography in each subbasin.Calculated flows are aggregated at the subbasin outlet and routed through the stream network.
To facilitate the use of GIS-based data for model setup, the ArcSWAT interface (Winchell et al., 2008) has been created to link the SWAT model and ArcGIS 9.For this study ArcGIS 9.2 and ArcSWAT 2.0.0 were employed, which includes the SWAT model in its 2005 version.
Altogether 20 parameters were used for automated calibration using the PEST program (Doherty, 2004) for strategies I and II, as described below.Seventeen were responsible for calibration of the hydrology, while the remaining three governed the nitrate calibration (see Table 1).Modeling time period was from 1979 till 2005.The first two years were taken as warm-up period for the model.For discharge, calibration was from 1981 to 1986, validation from 1987 to 1991.For nitrate, since only monthly values were available, calibration was from 1981 to 1995 and validation from 1996 to 2005.

Three strategies to evaluate the influence
Three strategies were developed to examine the influence of the four different interpolation methods on runoff and nitrate simulation results: Strategy I : in this simple approach, the SWAT model was calibrated only one time, using PEST and one selected set of interpolated climate data from mixed interpolation methods.For each climate variable for the mixed set the best performing interpolation method during cross-validation was selected.Then, all four interpolation climate data sets were subsequently applied to the calibrated model and the modeling results were compared running the SWAT model for both the calibration and the validation periods.The climate data for calibration was selected from the interpolated data sets via crossvalidation.Data resulting from the method which performed best in crossvalidation considering correlation and RMSE were chosen for each climate variable.This was KED for temperature, humidity, and wind speed, OK for precipitation, and IDW for sunshine duration.One possible problem with this strategy is that it might be biased towards one or more of those interpolation methods used for calibration.
Strategy II : to circumvent a biased calibration the second approach was designed as fourfold calibration of the model, separately with each of the four interpolated climate data sets.Again automated calibration with the PEST tool was employed.Results of those four calibrated models were compared, again for simulations using both the calibration and validation periods.Individual calibration however might cover the differences in modeling results arising from the different interpolation methods.
Strategy III : to overcome the possible problem of compensation of poor climate data by separate calibration the third strategy was developed.No calibration is necessary for this approach.Of the 20 model parameters used for calibration in strategies I and II, 1000 random parameter sets were generated, uniformly distributed between the minimum and maximum values allowed by SWAT.These 1000 parameter sets were applied to each of the four interpolated climate data sets.As evaluation the number of all realisations above a selected threshold in Nash-Sutcliffe-Efficiency (determined over a period of 10 years from 1981 to 1990) for both discharge and nitrate load was counted for each interpolation method.

Results and discussion
In Table 2 some exemplary results of the spatial interpolation are shown.The numbers represent average long-term values of precipitation sum and maximum daily temperature for the Upper Leine catchment.The presented numbers already show clearly the variation in the different interpolation methods, while the differences can be much more pronounced for smaller regions or on individual days.
Simulation results were generally good with the calibrated SWAT model.Figure 2 exemplarily shows two graphs, one for daily discharge and one for monthly nitrate load at gauge Leineturm.

Strategy I
Strategy I did not show any significant differences between interpolation methods, neither for discharge nor for nitrate simulation.Figure 3 shows results for gauges Leineturm and Reckershausen, the other gauges are similar.Nash-Sutcliffe efficiency (NSE) is higher for the discharge and nitrate gauges with larger catchment areas, as it is usually the case, but there is no conclusive difference between interpolation methods at each gauge.This is surprising, since it was expected that this approach would perform better for interpolated climate data used for calibration.A reason might be the random influence of automated calibration.Not only the choice of the climate data set for calibration, but also the PEST optimization start parameters have a strong effect on the final calibrated model, and thus the performance of the modeling results.This means that strategy I is not suitable to detect any differences in modeling results arising from the interpolation method for the climate data.

Strategy II
As strategy I before, the second approach neither showed any significant differences between interpolation methods (Fig. 4, gauges not shown are similar).Again NSE for both discharge and nitrate varies between gauges, but not between interpolation methods.Although the influence of a biased calibration has been removed, this outcome is plausible, since separate calibration might compensate for differences resulting from the climate data sets.Furthermore, the influence of the PEST optimization start parameters is very strong, as it turned out in a follow-up investigation.Strategy II was carried out several times using different PEST start parameters, and each time the result was different, and no realization clearly favoured any of the interpolation methods in modeling results.

Strategy III
The results of strategy III are presented in Fig. 5  interpolation method on discharge simulation.Although this strategy is not biased towards any method and no calibration is carried out, all discharge gauges show a similar number of "good" realizations with a Nash-Sutcliffe efficiency above the chosen threshold for each method.In agreement with the previous methods, higher numbers of good realizations are found for gauges with larger catchment areas.This holds true for all tested threshold values.Nitrate simulation on the other hand shows an extreme and clear influence of the interpolation method.All three gauges consistently reveal a rising number of good realizations from the simplest to the most sophisticated interpolation method, with KED clearly performing best, often with about twice as much good realizations than the second-best method.This again holds true for all tested threshold values.

Conclusions and discussions
Although strategies I and II did not show any significant difference between interpolation methods, strategy III proves that there is a strong influence at least on nitrate simulation.KED clearly outperforms all other evaluated interpolation methods, and a well-defined rank order from the simplest to the most sophisticated method is visible.The question arises, why this only shows in the modeling results for nitrate, and not for discharge.Most likely the measuring network density plays a major role in this.Previous studies proved an influence of the network density on modeling results (e.g.Bárdossy and Das, 2008).The main driving force for discharge is precipitation, while nitrate dynamics strongly depend on plant growth, with temperature as its main driving force.The precipitation measuring network for the Aller river catchment has a four times higher density than the measuring network for the other climate variables, which even raises to six times when taking the average number of "non-missing" stations into account.It is selfevident that the higher the network density, the smaller the differences between interpolation methods.Further analysis showed that indeed temperature makes up for the largest part of the differences in performance between NN and KED in strategy III.Since temperature also is the climate variable benefitting most from more sophisticated interpolation methods, especially KED with elevation as secondary information, as shown by cross-validation, the large difference in performance between interpolation methods regarding nitrate simulation in strategy III can be explained.To test the influence of measuring network density, the analysis could be extended considering different densities of the precipitation network.
One major problem with the first two strategies is the random influence of the automated calibration.The optimization starting parameters for PEST have a strong influence on the final calibrated model parameters, which of course influences which interpolated climate data performs best.
Even when one climate data set is used for calibration, as in strategy II, it does in no way guarantee that this data set performs best with the calibrated model.Subsequent trials revealed that different optimization starting parameters yield totally different calibrated model parameters and different rank orders of performance of the interpolation methods.These findings allow two major conclusions: firstly, any strategies involving automated calibration in the way presented in this paper (one random set of optimization starting parameters, minimum and maximum parameter values as allowed by the model) are not suitable to evaluate the performance of interpolation methods.Even a combination of strategies I and II did not give any other results (fourfold calibration with each climate data set, application of all four data sets to all four calibrated models in turn).Secondly, if for a model application a single automated calibration as shown is used, it is very unlikely that the modeler will get the best out of his climate input data, even when testing several interpolation methods.
Automated calibration could become suitable for evaluation of the interpolation methods with some alterations and improvements of the evaluation strategies.For one, the permitted parameter space could be narrowed down to physically meaningful limits for the examined study area instead of employing the maximum limits set by the model.Furthermore, a Monte Carlo-like automated calibration with varying optimization starting parameters seems necessary to eliminate the random influence of the starting parameters.A Monte Carlo-like automated calibration might also reveal other aspects of the influence of the interpolation methods, e.g. if a more robust parameter calibration would be possible using climate data from more sophisticated interpolation methods.Other aspects for future research are to analyze how well the findings from this study are transferable to other investigation areas and/or other models.Also, if differences arise, it is worthwhile to invesigate which climate variables contribute most to the differences, when interpolated with different methods.This may lead to conclusions about whether one interpolation method is more suitable for a specific climate variable than another, with different methods performing best for different variables.A useful further study would be to establish a guideline for modelers how to interpolate the climate input data for their models in order to receive the best performance and the most robust model parameters, while minimizing the amount of time needed for calibration.

Table 1 .
SWAT Parameters used for automated calibration.
*The study area was divided into three hydrological landscapes for groundwater, R: from spring to gauge Reckershausen, G: from Reckershausen to Göttingen, L: from Göttingen to Leineturm.

Table 2 .
Exemplary results of the spatial interpolation of the climate variables.