Testing a river basin model with sensitivity analysis and autocalibration for an agricultural catchment in SW Finland

Modeling tools are needed to assess (i) the amounts of loading from agricultural sources to water bodies as well as (ii) the alternative management options in varying climatic conditions. These days, the implementation of Water Framework Directive (WFD) has put totally new requirements also for modeling approaches. The physically based models are commonly not operational and thus the usability of these models is restricted for a few selected catchments. But the rewarding feature of these process-based models is an option to study the effect of protection measures on a catchment scale and, up to a certain point, a possibility to upscale the results. In this study, the parameterization of the SWAT model was developed in terms of discharge dynamics and nutrient loads, and a sensitivity analysis regarding discharge and sediment concentration was made. The SWAT modeling exercise was carried out for a 2nd order catchment (Yläneenjoki, 233 km2) of the Eurajoki river basin in southwestern Finland. The Yläneenjoki catchment has been intensively monitored during the last 14 years. Hence, there was enough background information available for both parameter setup and calibration. In addition to load estimates, SWAT also offers possibility to assess the effects of various agricultural management actions like fertilization, tillage practices, choice of cultivated plants, buffer strips, sedimentation ponds and constructed wetlands (CWs) on loading. Moreover, information on local agricultural practices and the implemented and planned protective measures were readily available thanks to aware farmers and active authorities. Here, we studied how CWs can reduce the nutrient load at the outlet of the Yläneenjoki river basin. The results suggested that sensitivity analysis and autocalibration tools incorporated in the model are useful by pointing out the most influential parameters, and that flow dynamics and annual loading values can be modeled with reasonable accuracy with SWAT. Sensitivity analysis thus showed the parameters which should be known better in order to result in more realistic parameter values. Moreover, the scenario runs for CWs made with SWAT revealed the high demand of land area for this protective measure to be substantially effective.


Introduction
Due to the obligations set by the EU Water Framework Directive (WFD), the environmental authorities are in need for information on the effects of different management options aimed at improved water quality.In addition to the knowledge obtained from the Finnish field-scale experiments on the effects of agricultural practices and protection measures on loading (e.g.Turtola and Kemppainen 1998, Koskiaho et al. 2003, Puustinen et al. 2005, Uusi-Kämppä 2005), mathematical modeling tools are needed to generalize the effects in varying environmental conditions on the catchment scale.Models like SOIL/SOILN (Johnsson et al. 1987), GLEAMS (Knisel 1993) and ICECREAM (Tattari et al. 2001) have been used to assess phosphorus (P) and nitrogen (N) losses from agricultural land in Finland (Granlund et al. 2000, Knisel and Turtola 2000, Tattari et al. 2001).The catchment scale INCA-N model has also been applied for selected river basins (Rankinen et al. 2004).These models have, however, limitations partly in terms of catchment-scale evaluations of loading and the effects of management actions.For these purposes, the SWAT model (Soil and Water Assessment Tool (Arnold et al. 1998, Neitsch et al. 2005)) offers an attractive alternative.Unlike many other models, SWAT has an open-source code and since there are numerous SWAT users around the world, the model is regularly updated and the developers offer support for users.The SWAT model was chosen for this study since it simulates all relevant variables, suspended sediment as well as P and N loading on catchment scale.SWAT has also a potential to estimate the effects of a broad range of cultivation practices like reduced tillage and water protection measures such as buffer zones and CWs.SWAT has already been found to have potential with respect to the future requirements set by WFD in Scotland (Dilks et al. 2003) and in Finland (Bärlund et al. 2007).In cold winter conditions SWAT has been applied in Canada (Lévesque et al. 2008).In Finland, the SWAT model has been previously applied to the rivers basins of Vantaanjoki (Grizzetti et al. 2003) and Yläneenjoki (Bärlund et al. 2007, Bärlund and Kirkkala 2008, Tattari et al. 2008).
The objectives of this study are twofold: first to determine the most influential parameters which result in a substantially reduced parameter set for model calibration and to apply the SWAT autocalibration tool to improve model fit.Second objective was to test and evaluate the usefulness of the SWAT model (SWAT2005) for assessing the load estimates as well as to test the model's ability to estimate the effect of CWs in an agriculturally dominated catchment.

The SWAT model
The SWAT model (Soil and Water Assessment Tool) is a continuous time model that operates on a daily time step at catchment scale (Arnold et al. 1998, Neitsch et al. 2005).It can be used to simulate water and nutrient cycles in agriculturally dominated large catchments.The catchment is generally partitioned into a number of sub-basins based on the threshold area which defines the minimum drainage area required to form the origin of a stream.The smallest unit of discretization is a unique combination of soil and land use overlay referred to as a hydrologic response unit (HRU).SWAT is a partly process-based model and partly a distributed model, including many empirical relationships.The water quantity processes simulated by SWAT include precipitation, evapotranspiration, surface run-off and lateral subsurface flow, ground water flow and river flow.Water quality processes are calculated with various well-known equations.For example, erosion caused by rainfall and runoff is computed with the Modified Universal Soil Loss Equation (MUSLE) (Williams 1975).In terms of P, the primary mechanism of soluble fraction movement in the soil is by diffusion.Organic and mineral P attached to soil particles may be transported by surface runoff to the main channel (Neitsch et al. 2005).For channel flow simulation, SWAT uses Manning's equation coupled with variable storage or Muskingum routing method.Interactions and relationships of the QUAL2E model (Brown and Barnwell 1987) are used as in-stream water quality processes in SWAT.The model has been widely used but also further developed in Europe (e.g.Eckhardt et al. 2002, Krysanova et al. 1999, van Griensven and Meixner 2003).

The experimental catchment
The Yläneenjoki catchment, 233 km 2 in area, is located on the coastal plains of south-western Finland, thus the landscape ranges in altitude from 50 to 100 m above the sea level (Fig 1 .).The soils in the river valley are mainly clay, whereas moraine and peat soils dominate elsewhere in the catchment (Fig. 1).Forests and natural wetlands cover 65% of the catchment the rest being agricultural (34%) and urban (1%) areas.Agriculture in the Yläneenjoki catchment consists of mainly cereal production (Fig. 1) and poultry husbandry and in this part of Finland agriculture can be considered as intensive for Finnish standards.Subsurface drainage is applied for major part of the agricultural soils in Yläneenjoki catchment.The recommended drain depth is circa 1.2 m and drain spacing around 20 m depending on hydraulic conductivity and drainage demands.All farmers in the Yläneenjoki area follow the recommended fertilizer levels which were derived from studies of the Finnish Agri-Environmental Programme (Palva at al. 2001).Both winter-time vegetation cover and reduced tillage have become more common during the recent years.Long-term (1961Long-term ( -1990) ) average annual precipitation is 630 mm of which approximately 11% falls as snow (given as the maximum water equivalent of the snow cover, assuming no sublimation) (Hyvärinen et al. 1995).The average monthly temperature for the period November to March ranges between -0.5 and -6.5 ºC.Average discharge in the Yläneenjoki main channel is 2.1 m 3 s -1 (Mattila et al. 2001), which equates to an annual water yield of 242 mm (1980)(1981)(1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990).The highest discharges typically occur in spring and late autumn.The portion of groundwater flow is not measured but according to typical annual water balances groundwater accounts for less than 20% of annual rainfall.Erosion rates in Yläneenjoki area vary owing to local natural conditions and management practices.Agricultural field plots and cathments dominated by agriculture produce higher erosion rates than forested areas.The clayey soil of the south-western coastal plains is more susceptible to erosion than inland predominantly sandy and till soils (Tattari and Rekolainen, 2006).
The regular monitoring of water quality of the river Yläneenjoki has been started as early as 1970s.The nutrient load into the Lake Pyhäjärvi via the river has been calculated based on monthly averages from manual water sampling results (ca 12 samples per year) and daily water flow records obtained by site-specific head-discharge relation at one point, Vanhakartano (Kauppila and Koskiaho, 2003).Both the monitoring time-series and calculated loadings are used to compare model results with observations in order to evaluate model performance in this study.

Background GIS data
For the SWAT simulations the available data on elevation, land use and soil types had to be aggregated.The resolution (5 m) of the digital elevation model (DEM) proved to be inaccurate for low-lying areas for successful set-up of the Yläneenjoki catchment.Hence, we used a modified DEM where the main channels of the catchment were somewhat deepened to emphasize the actual routes of water.The agricultural crops and their location were available from Information Centre of the Ministry of Agriculture and Forestry in Finland (TIKE).According to this database 100 different crop types are grown in the Yläneenjoki region.To build up a reasonable SWAT set-up this data was divided into 5 crop classes (autumn cereals, spring cereals, root crops, grasses and other crops) (Fig. 1).All the rest land area was classified as forest.Soil type data was based on soil textural information of the Geological Survey of Finland in which soil type information covers southern Finland in a scale of 1:100 000 and it is available in 25 × 25 meter cells.Our SWAT project had 4 soil types namely clay, moraine, silt and peat (Fig. 1).

The modeling approach, sensitivity analysis and autocalibration
The SWAT project was build for the Yläneenjoki catchment making use of the earlier SWAT application, for example, the parameters of Bärlund's (2007) earlier project were set as initial values.In contrast to the earlier work, more detailed GIS based land use map was used and in addition, the simulation years differed from the earlier set-up.The present SWAT project resulted in 29 sub-catchments.In the project, threshold values were used to distinguish different land use and soil types within each sub-catchment.For example, if more than 1% of a sub-catchment is under grass and these areas were divided on clay and silt soils both soils types representing more than 10% of the sub-catchment area, this would result in two HRUs (grass-clay and grass-silt) within this sub-catchment.In all, this approach resulted in 257 HRUs in the project.
The sensitivity analysis method implemented in SWAT is called Latin Hypercube One-factor-At-a-Time (LH-OAT) (Morris 1991).The procedure of LH sampling is simple.First, the distribution of each selected parameter is subdivided into m ranges, each with a probability of occurrence equal to 1/m.In the second phase, random values of the parameters are generated.Then the model is run m times with the random combinations of the parameters.In the LH-OAT sensitivity analysis, the OAT (change each parameter sequentially at a time) method is repeated for each point sampled in the LH sampling.Sensitivity is expressed by a dimensionless index I, which is calculated as the ratio between the relative change in model output and the relative change of a parameter.The sensitivity analysis was performed for a set of 32 parameters that have been found to be most influential for river discharge and sediment transport with many earlier modeling studies (Krysanova andArnold 2008, Bärlund et al. 2007).The ranges of variation of these parameters are based on earlier studies and expert knowledge.For some parameters a relative change was preferred to absolute values (Table 1).The sensitivity classes are shown in Table 2.In the latest version of SWAT, autocalibration is implemented by using the Shuffled Complex Evolution (SCE-UA) algorithm that optimizes an objective function (SSQ= sum of the squares of the residuals) by systematically searching the entire parameter space (global optimization) ( van Griensven and Meixner 2006).Daily averages of discharge were calibrated against the corresponding values determined from the observations made at the Vanhakartano measurement station.Here, the autocalibration tool of SWAT was utilized.Nash-Sutcliffe coefficients were calculated for the calibration results.Years 1995-1999 were chosen as the calibration period.In terms of hydrology, the years differed quite a lot from each other.

Scenario runs for constructed wetlands
The calibration period 1995-1999 was suitable for testing of the scenarios (Table 3) because almost every CW in the area has been established in 2000s.For the scenario runs, reduction performance of each CW was determined according to their dimensioning as presented in Puustinen et al. (2007): where TP ret = Total phosphorus reduction per year (% of the input loading) TN ret = Total nitrogen reduction per year (% of the input loading) A rel = CW-to-watershed area ratio (%) The effects of CWs were assessed by comparing the material fluxes produced by scenarios 1, 2 and 3 with those produced by 0-scenario at the Vanhakartano measurement point.

Sensitivity analysis
Tables 4a and 4b summarizes the sensitivity ranking for the performance of discharge and sediment concentration.In terms of discharge, none of the parameters got the ranking "very high sensitivity".
The most influential parameters representing the sensitivity class "high sensitivity" were connected to equations describing processes such as ground water flow (parameter GWQMN), evaporation and surface and subsurface runoff (SOL_z), snow melt (TIMP), soil evaporation (ESCO), soil water dynamic ( SOL_AWC), surface and subsurface runoff (CN2 and SOL_K).Detailed descriptions of these parameters are given in Table 1.Surface and groundwater flow, snow melt, and evaporation are all highly relevant processes when simulating discharge in Finnish hydrological conditions.Our ranges of these parameters differ slightly from the ones presented by Holvoet et al. (2005), where the most influential parameters for hydrology were in descending order: CN2, RCHRG_DP, SURLAG, GWQMN, ESCO, SOL_z, ALPHA_BF (in our study the ranking order was: GWQMN, SOL_z, TIMP, ESCO, SOL_AWC, CN2 and SOL_K).The differences in parameter order between these two studies can be almost totally explained by the different parameter ranges and by different catchment characteristics.For example, Holvoet et al. (2005) applied the relative change of ± 50% while our choice of variation was only ± 10-20% resulting in In terms of sediment concentration, four parameters, SPCON and SURLAG representing general watershed properties and CH_COV and CH_EROD representing main channel characteristics, respectively, were ranked in the sensitivity class "very high sensitivity".SPCON is a coefficient in sediment transport equation and is defined by the user.SURLAG is surface runoff lag coefficient and also defined by the user.The factors CH_COV and CH_EROD are conceptually similar to the soil erodibility factors used in the USLE equation and in theory they can be measured.For the Yläneenjoki area, these data were not available.The parameters ranked as "high sensitivity" in the case of discharge were also here equally ranked.The results of the sensitivity study were further utilized in the calibration procedure as discussed in the following.

Calibration
Daily and monthly averages of flow were calibrated against the corresponding values determined from the observations made at Vanhakartano.The calibration was done with the parameters which were ranked as most influential ones in the sensitivity analysis study.Autocalibration results of 6 simulation runs suggested the best parameter values including their good range in terms of discharge (Table 5).For example, SURLAG decreased from the initial value 4 to 0.4 indicating that the surface runoff lag is very short in the Yläneenjoki area.For parameter SOL_AWC the analysis gave quite a high value which at the same time indicated that soil water holding capacity needed to be greater.Nash-Sutcliffe (NS) coefficients were calculated for the calibration results.NS coefficients were done for monthly values because for daily values it is very difficult to achieve a good fit even with intensive calibration efforts.The best coefficient value was 0.7 for the monthly discharge.Figure 2 shows the clear improvement obtained by the calibration process in the fit between simulated and observed daily flow.Particularly in spring and autumn the amendment was clear: not only the simulated peaks were much closer to the observed, but also the autumnal low flow period appeared much more realistic.In mid-winter and summer the results remained, even after calibration, relatively weak.There seems to be substantial snow melt in January and heavy summer rains, during July, that are missed by the model.As for annual average flow, the simulated values were generally lower than the measured.
In spite of the visual improvements in flow dynamics achieved by the autocalibration process, the fit between the simulated and observed daily flow remained poor when it was assessed by NS coefficients.For phosphorus and erosion, the compatibility between simulations and the calculated loads based on measurements were examined only with annual values (Fig. 3).In the first year (1995) of calibration suspended sediment loading was highly overestimated and total phosphorus load, on the contrary, underestimated.One reason for this could be that the warm-up period of the modeling set-up was only one year.On the other hand, the load estimates based on manual grab sampling inherently lack reliability because long periods between the sampling occasions remain unknown.Our recent unpublished study shows that e.g.seasonal erosion can be considerably underestimated (19-72%) when based on manual sampling estimates.The years 1996-1999 show better fit exclusive of year 1999 for total P load.

Scenario runs for constructed wetlands
In a well-dimensioned (A rel =5%) CW, annual reductions of TSS and TP may be even 70% (Koskiaho et al. 2003).The dimensioning of the CWs established in early 2000s in Yläneenjoki basin was however not as generous.Hence, the load reduction achieved with scenario 1 expectedly remained rather low (1.3% for both TP and TN).With scenario 2 TP loading was reduced by 3.7% and TN loading by 3.9%, i.e. still not very much.The results suggest that even high increase of the number of CWs does not lead to substantial load reductions if their dimensioning is inadequate.Meanwhile, when the dimensioning of the 10 CWs was changed according to scenario 3 (Table 3), TP loading was reduced by 17% and TN loading by 18%.In practice, however, this would mean that even more than 1% of the Yläneenjoki basin should be converted to CWs, i.e. very large land areas should be available.The most realistic and cost-effective approach is probably to try to concentrate the CWs in such parts of the catchment, where the above area is not very large and input  1.1999 1.2.1999 1.3.1999 1.4.1999 1.5.1999 1.6.1999 1.7.1999 1.8.1999 1.9.1999 1.10.1999 1.11.1999 1.12.1999Discharge  1.1999 1.2.1999 1.3.1999 1.4.1999 1.5.1999 1.6.1999 1.7.1999 1.8.1999 1.9.1999 1.10.1999 1.11.1999 1.12.1999Discharge  concentrations are high (high field-%, steep slopes, high number of farms with animal husbandry).

Conclusions
The sensitivity analysis proved to be useful not only by pointing out the most important parameters for discharge and sediment concentration calibration but at the same time also demonstrated which parameters should be better known or measured in the catchment for future modeling approaches.The most influential parameters for discharge were 1) Threshold water depth in the shallow aquifer for return flow to occur 2) Soil depth 3) Snow pack temperature lag factor 4) Soil evaporation compensation factor 5) Available soil water capacity 6) Initial SCS II value and 7) Saturated hydraulic conductivity.In terms of sediment concentration, four parameters, namely 1) Linear re-entrainment parameter for channel sediment routing 2) Channel cover factor and 3) Channel erodibility factor and 4) Surface runoff lag coefficient, were the most influential ones.Narrowing the range of parameter values improved the outcome.The parameterization of the model to achieve satisfactory calibration results in terms of flow and sediment dynamics proved to be a laborious task.Achieving satisfactory model fit on daily basis is quite challenging for a long calibration period with varying flow patterns.The sensitivity analysis, however, gave a good starting point for the calibration, which greatly improved the fit between simulated and observed daily flow particularly in spring and autumn.For some parameters e.g.available soil water capacity (SOL_AWC) the autocalibration produced somewhat unrealistic values.As for annual average flow, the simulated values were generally lower than the measured.SWAT simulated quite well the sediment and phosphorus concentrations along the main channel giving, according to the monitoring experiences, higher values for the upper reaches in the catchment and lower values for the measurement points closer to the lake Pyhäjärvi.
Invaluable information was available by cooperative local farmers and authorities.They were willing to give us material not only about the agricultural management practices and protective measures implemented in the Yläneenjoki catchment so far, but also about the measures planned for the future to protect the Lake Pyhäjärvi.For example, there were seven CWs already established in the area and many more in planning phase.In terms of CWs the results suggest that quite large areas are needed if substantial reductions in agricultural loading are to be aimed for.In order to be efficient, a CW has to have high CW-to-watershed area ratio and this requirement of dimensioning tends to increase the total land area needed.If dimensioning of single CWs is inadequate, not even high increase in their number seems to significantly contribute to water pollution control.Hence, in reality, CWs should be seen as a good supplementary part of comprehensive water protection planning rather as the epoch-making solution for the problem.

Fig. 1 .
Fig. 1.Location of the study area, and topographic, agricultural land use and soil maps of the Yläneenjoki catchment.

Fig. 2 .
Fig. 2. Measured and simulated discharge for Yläneenjoki outlet (Vanhakartano) with original manually calibrated parameter set (upper figure) and with the autocalibrated parameter set (lower figure).

Table 1 .
Model parameters and their ranges used in the sensitivity analysis.
* = relative percent change

Table 4 .
Sensitivity ranking for discharge (a) and sediment concentration (b) parameters.

Table 5 .
Autocalibration results of 6 simulations made for a selected set of parameters with SWAT.