VZJ Journal of Natural Resources and Life Sciences Education
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 23 January 2008
Published in Vadose Zone J 7:10-21 (2008)
DOI: 10.2136/vzj2006.0186
© 2008 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by van der Keur, P.
Right arrow Articles by Refsgaard, J. C.
Right arrow Search for Related Content
PubMed
Right arrow Articles by van der Keur, P.
Right arrow Articles by Refsgaard, J. C.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by van der Keur, P.
Right arrow Articles by Refsgaard, J. C.
Related Collections
Right arrow Water Quality
Right arrow Nutrients

ORIGINAL RESEARCH

Uncertainty in Simulation of Nitrate Leaching at Field and Catchment Scale within the Odense River Basin

P. van der Keura,*, J. R. Hansena, S. Hansenb and J. C. Refsgaarda

a Geological Survey of Denmark and Greenland, Østervoldgade 10, DK-1350 Copenhagen K, Denmark
b Dep. of Agricultural Sciences, Univ. of Copenhagen, Højbakkegaard Allé 9, DK-2630 Taastrup, Denmark
* Corresponding author (pke{at}geus.dk).

All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.


Received 22 December 2006.



   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 Case Study Location
 Methodology
 Results and Discussion
 Conclusions
 REFERENCES
 
There is a clear need to conduct uncertainty analyses with respect to nitrate leaching at a range of spatial scales, from the field to the catchment scale, dependent on the spatial scale for which an assessment is relevant. At the scale of the Odense catchment in Denmark, an uncertainty assessment most likely leads to an unacceptably high number of simulations as a result of the large number of combinations of distributed data available. In this study we present a framework for an uncertainty assessment for model parameters for the unsaturated component in a linked rootzone and groundwater model and exemplified for the Odense River catchment. The proposed framework consists of (i) a simplification of the model, (ii) identification of uncertain parameters, (iii) generation of model parameter sets by means of Latin hypercube sampling (LHS), (iv) rootzone model simulation of realizations under (iii), (v) uncertainty analysis based on the standardized rank regression coefficient index, (vi) ranking of simulated nitrate leaching, and (vii) simulation of selected parameter sets sampled from the output distribution under (vi). The results show a high sensitivity of the van Genuchten soil water release characteristics parameters for percolation and leaching. The amount of applied slurry has a relatively high sensitivity for leached nitrate from the rootzone. The ranking nitrate leaching derived from the LHS simulations for a 25-yr period is preserved for longer 55-yr simulations. This implies that the proposed framework is applicable to the Odense catchment and likely has general applicability. In this way, uncertainty assessments are feasible at the catchment scale.

Abbreviations: GIS, geographic information system • IQR, interquartile range • LHS, Latin hypercube sampling • pdf, probability density function • SOM, soil organic matter • VIF, variance inflation factor • WFD, Water Framework Directive


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
Case Study Location
 Methodology
 Results and Discussion
 Conclusions
 REFERENCES
 
Nonpoint-source pollution from agriculture is a widespread problem in Europe and North America. This typically concerns compounds like nitrogen from fertilizers and slurry, phosphorus, and pesticides, among other sources. High levels of nitrogen in groundwater reservoirs pose a threat to the availability of clean drinking water, particularly in countries that are almost entirely dependent on groundwater reservoirs, such as Denmark. In Europe, the Water Framework Directive (WFD, 2000/60/EC) has prompted the need to develop methodologies and tools for implementing the WFD and integrating uncertainties associated with this process (Refsgaard et al., 2005). Research has usually focused on transport and fate of such compounds at the spatial scale of plots and fields and temporal scales of months and a few years. Policymakers, however, are interested in long-term effects of agricultural practices as part of longer-term policies and at scales larger than the plot and field scale (Mulla and McBratney, 2003). Physically based rootzone models for predicting nitrate leaching from the unsaturated zone to groundwater reservoirs to assess water quality are quite useful and have been routinely applied in leaching vulnerability assessments in water resources management for many decades. Mapping of groundwater vulnerability may be derived from an index-based method, such as DRASTIC (Aller et al., 1987), or regression types of models, such as the recently developed logistic regression model (Nolan et al., 2002, 1997). The main advantage of a process-based approach such as that used here, as opposed to an index- or regression-oriented method, is the potential of scenario analyses under temporally and spatially changing conditions, extrapolating to future changes in, for example, climatic regimes.

Vachaud and Chen (2002a,b) assessed the uncertainty related to nitrogen leaching at various scales using physically based models for a small (200-km2) catchment in southern France, focusing on the importance of within soil class variability of transport parameters. Wriedt and Rode (2005) demonstrated the uncertainty in modeling nitrate transport and turnover as a result of complex and spatially distributed physical and chemical interactions for a 20-km2 lowland catchment system in a northern Germany lowland. At the national scale, Kroeze et al. (2003) and de Vries et al. (2003) focused their studies on uncertainties in the fate of nitrogen by an assessment of an overall nitrogen budget experimentally and by means of a model, respectively, for the Netherlands. For a wide range of spatial scales, Quinn (2004) advocated a scale-appropriate modeling technique to best understand nitrate losses at that particular scale. This was demonstrated for a 1400-km2 catchment in the United Kingdom by combining knowledge from simulation, geographic information system (GIS) databases, and expert judgment. Hansen et al. (1999) performed an uncertainty analysis for nitrate leaching from the rootzone using the DAISY model (Hansen et al., 1991) for several land management types in the Karup Basin in the western part of Denmark. Thorsen et al. (2001) followed up on the studies by Refsgaard et al. (1999) and Hansen et al. (1999) on catchment-scale modeling of nitrate leaching for the Karup Basin, applying the MIKE-SHE/DAISY (Styczen and Storm, 1993) code to the Karup Basin using Monte Carlo techniques for assessing model input parameter uncertainty for the same five key parameters as in Hansen et al. (1999). These parameters—precipitation, soil hydraulic properties, soil organic matter content, slurry composition and depth of the oxidized zone (redox interface) in the aquifer—are known to be important for water balance, nitrate leaching, and transformation. Thorsen et al. (2001) concluded that uncertainty in simulated nitrate flux concentrations depends on the considered spatial and temporal scale; that is, uncertainty in simulated flux concentrations from the rootzone is much larger than simulated concentrations at the catchment scale. Uncertainty in model structure was not considered, although known to be important and usually considered much larger than the first (e.g., Carrera and Neuman, 1986; Refsgaard et al., 2006).

To capture the complexity and nonlinearity in processes controlling water flow and solute transport, physically based hydrological modeling in integral water management is required (Blöschl and Sivapalan, 1995). A wide range of approaches for performing uncertainty analyses exist, including methods based on Taylor expansion (Iman and Helton, 1985), differential analysis (Rogers et al., 1985), and error propagation (Janssen et al., 1990) (refer to Saltelli et al., 2004, for more information). All of these methods have advantages and disadvantages. It is beyond the scope of this paper to assess their possible use in the present context; see Christiaens (2001) for a review. The Latin hypercube sampling (LHS) technique adopted here is an efficient way of dealing with output uncertainty of models with a distributed, physically based character and many model parameters (Iman et al., 1980) for generating model parameter distributions. This technique uses a stratified sampling approach, contrary to Monte Carlo, which uses random sampling. It thus requires far less samples than needed for pure Monte Carlo. Christiaens and Feyen (2001, 2000) used LHS for an uncertainty analysis on hydraulic properties and propagation through the distributed MIKE-SHE model (Refsgaard and Storm, 1995). As mentioned above, Hansen et al. (1999) and Thorsen et al. (2001) also used LHS for parameter distribution propagation in the MIKE-SHE/DAISY model.

A severe limitation for all mentioned uncertainty assessment methods is that only parameter uncertainty is considered, while the conceptual model structure is assumed to be correctly captured. Uncertainty assessments for a distributed modeling approach using a combined rootzone and groundwater model are hampered by the sheer number of Monte Carlo–generated simulation runs when grid-scale distributed data on soil properties, climate, upper- and lower-boundary conditions, and land management are used. The rootzone component in a distributed groundwater model can be parameterized for a large number of combinations of distributed data, each combination representing a geographical location (grid) in the catchment. Output distributions from rootzone model simulations must then be propagated through the catchment-scale groundwater model, potentially leading to prohibitively long and complicated uncertainty analyses. Therefore, there is a clear need to study whether such distributed systems can be simplified, rendering uncertainty analyses feasible.

The need for an uncertainty assessment in integrated water management at the river basin scale to support decision makers in implementing the WFD in Europe led to the EU-FP5 HarmoniRiB project (Refsgaard et al., 2005). This work is part of the Odense case study within the HarmoniRiB project. It puts forth a novel way to analyze uncertainties involved in assessing policy measures for reducing nitrate pollution due to agricultural activities and its impact on groundwater resources and is part of a broader aim to consider balancing uncertainty in effects for the overall model-predicted improvement of water quality and uncertainty in economic consequences. This is highly relevant for implementing the WFD. To our knowledge, no prior attempts have been made to investigate whether simplifying complex hydrological models can help support required uncertainty analysis in integrated water management. We propose a framework to facilitate uncertainty assessments for modeling nitrate leaching within the context of integrated water management at the basin scale using a combined rootzone–groundwater model. The focus is on uncertainty in the rootzone domain as a first step toward a combined rootzone–groundwater modeling uncertainty assessment approach.

The objectives of this study are (i) to develop a framework for facilitating an uncertainty assessment for simulating nitrate leaching at the field and catchment scale with focus on the rootzone component, and (ii) to apply the proposed framework to the Odense River catchment.


   Case Study Location
TOP
 ABSTRACT
 INTRODUCTION
 Case Study Location
Methodology
 Results and Discussion
 Conclusions
 REFERENCES
 
The Odense River Basin (Fig. 1 ), on the island of Funen in Denmark, encompasses 622 km2 and is part of the Odense Fjord Basin (1,046 km2). The Odense River Basin is the largest basin in the Odense Fjord basin, and the Odense River catchment is taken as a representative part of the Odense River Basin for the analyses in this work.


Figure 1
View larger version (24K):
[in this window]
[in a new window]

  		FIG. 1. The Odense River Basin, Denmark.

 
The Odense floodplain was formed by a meltwater river during the last glacial period. Clayey soil types are slightly dominant and encompass approximately 51% of the basin, while sandy soil types cover approximately 49% (Fyns County, 2003). The glacial sediment soils of Funen are particularly well suited to agricultural crops. By far the most important contributor to the fjord, affected by impacts of nutrients, is the Odense River Basin. Monitoring performed by the Funen County Council since 1976 shows that the objectives of the Regional Plan for the catchment have still not been obtained for the Odense fjord, although diffuse, nonpoint, nitrogen input has been reduced by 25 to 35%. Water Framework Directive typologies, reference conditions, and future objectives for rivers, streams, lakes, groundwater, and transitional waters have been established for the Odense River Basin by Funen County and is thoroughly described in Fyns County (2003). The amount and quality of existing data for the Odense catchment is very good, and many modeling studies have been performed in the past. Basin modeling has been performed with a focus both on groundwater conditions and on runoff and nutrient loads. Several nutrient models have been tested for this basin under the EUROHARP project (www.euroharp.org). As part of the National Aquatic Environmental preparatory studies in Denmark, a comprehensive model was established for simulation of nitrate leaching from the Odense Fjord Basin (Nielsen et al., 2004).


   Methodology
TOP
 ABSTRACT
 INTRODUCTION
 Case Study Location
 Methodology
Results and Discussion
 Conclusions
 REFERENCES
 
Framework for Uncertainty Assessments
In this section, a framework for uncertainty assessments with respect to simulating nitrate leaching is presented and is then applied to the Odense catchment. The framework consists of the following consecutive steps:

  1. Reducing model complexity. The combined rootzone–groundwater model for a basin must be simplified while still able to capture the essential controls of the system. This needs to be checked by comparison to the original model, for example, by calculating the Nash–Sutcliffe efficiency (Nash and Sutcliffe, 1970).
  2. Identifying model parameters that most likely contribute to uncertainty in the simulated output of interest. Output in this case is simulated leaching from the rootzone to groundwater.
  3. Generating model parameter sets from the parameters identified in Step 2. Usually, a stratified sampling technique by means of LHS is most effective.
  4. Simulating selected output by Monte Carlo. The parameter sets generated under Step 3 are input to the rootzone model and running for a sufficient period of time.
  5. Performing an uncertainty analysis with respect to the output of interest and justification of identified uncertain model parameters under Step 2.
  6. Ranking of model output and corresponding parameters sets. Selecting percentile values of output so as to cover the entire output distribution.
  7. Running the rootzone model for selected parameter sets corresponding to selected percentile values for an extended period to improve model results. Consistency check to ensure that ranking of output for an extended period conforms to the results under Step 6.

Reduction of Model Complexity (Step 1)
The DAISY Model
The DAISY model is a one-dimensional soil–plant–atmosphere system (Hansen et al., 1991; Svendsen et al., 1995; Abrahamsen and Hansen, 2000) designed to simulate water balance, solute balance, and crop production in agroecosystems. Soil water dynamics includes water flow described by the Richards equation in the soil matrix as well as macropore flow (Christiansen et al., 2004). Furthermore, it includes water uptake by plants and pipe drainage. The solute balance model simulates transport, sorption, transformation processes, and plant uptake by solving the convection–dispersion equation. Special emphasis is put on nitrogen dynamics in agroecosystems. Mineralization immobilization, nitrification and denitrification, sorption of ammonium, uptake of nitrate and ammonium, and leaching of nitrate and ammonium are simulated. A schematic representation is shown in Fig. 2 .


Figure 2
View larger version (54K):
[in this window]
[in a new window]

  		FIG. 2. The DAISY model. (SVAT = Soil-Vegetation-Atmosphere Transfer scheme.)

 
The atmosphere and the groundwater constitute the boundaries of the considered system. Considered fluxes at the surface are precipitation and irrigation (gains) and evapotranspiration and surface runoff (losses). Considered fluxes at the lower boundary of the system are deep percolation (loss), capillary rise (gain), and if the considered system contains artificial drain, also drain flow (loss). The model considers the recently added fresh organic matter with slow turnover rate, the living soil microbial biomass, and old native soil organic matter (SOM), subdivided in SOM1 and SOM2, with a high turnover rate.

The DAISY model has been linked to a fully distributed catchment model by Styczen and Storm (1993) for including the groundwater component and making groundwater simulations with DAISY input possible. A full coupling between DAISY and MIKE-SHE and the significance of feedback from the groundwater to the rootzone model is reported by Christiansen et al. (2004) but is not operational in the present study.

Modifying the DAISY Component of the Original Distributed Model
Under Step 1, an original model developed by Nielsen et al. (2004) as a comprehensive tool for assessing the national plan for the aquatic environment based on the MIKE-SHE/DAISY model code is modified to make an uncertainty analysis viable. In the original model, the different information sources on climate, soil type, land management, and boundary conditions are integrated in a GIS containing a large number of polygons. The Nielsen et al. (2004) model was set up for the Odense Fjord catchment and built on the Danish National Water Resource Model (Henriksen et al., 1997, 2003). The groundwater domain is divided into nine geological layers, and the grid size is 500 x 500 m. A digital elevation model was developed for this grid size from 25 x 25-m elevation data. Climate data included 27 10 x 10-km corrected precipitation grids and 20 x 20-km temperature grids, while global radiation data from one climate station was assumed to represent the whole catchment. Soil data for A (0–30 cm), B (30–80 cm) and C (80–350 cm) horizons was based on various Danish databases. The dominating soils consist of sandy loam (74%), coarse sand (21%), clay (1.5%), and peat (3.5%). Hydraulic parameters have been estimated from the textural composition of these soils using the HYPRES pedotransfer function (Wösten et al., 1999). Land use for the catchment is estimated from a combination of national databases and further described in Fyns County (2003). In the original model, the DAISY model was set up for a large number of soil columns (2094) derived from combinations of climate zones (grids), soil type, groundwater level, and land management. In addition, within one type of land management, crop rotations were permutated; that is, the sequence was changed to assure that all crop combinations occur for all years simulated. Model results of permutated crop rotations were then averaged for further computations. These combinations led to a total of 5577 DAISY time series for input into the MIKE-SHE module. For an uncertainty study based on Monte Carlo simulations, this large number of DAISY columns clearly must be reduced dramatically since each column requires several hundreds of model parameterizations, rendering Monte Carlo–type modeling too heavy a task. Modification—that is, simplification—of the original model will basically lead to another model, which may behave differently compared with the full-scale model. Justification for this is discussed further below. Hence, to reduce the complexity of the model, the modeled area is reduced. In addition, the number of combinations of climate zones, soil types, land management, and lower boundary must be reduced substantially. This is accomplished and justified by the following:

  1. Reducing the simulated area of the entire Odense Fjord Basin (1046 km2) to the Odense River catchment (622 km2). As mentioned above, the Odense River drains the major part of the Odense Fjord Basin, and the Odense River is therefore by far the most important contributor of water to the Odense Fjord. An uncertainty analysis for the Odense River Basin is thus likely to be representative of the Odense Fjord Basin.
  2. Reducing the number of rainfall stations from 27 climate grid zones to one rain gauge within the area. This rain gauge has measured the same amount of yearly average precipitation for the period 1992 to 2001 as a weighted average of the 27 10 · 10-km grids. The spatial variation in measured grid precipitation across the catchment is ~25%. Uncertainties in rainfall may be spatially correlated when estimated from observations that are spatially correlated. Storm et al. (1988) indicated on the basis of geostatistical analyses that the length scale for precipitation under Danish conditions is on the order of magnitude of 10 km. Our study adopts the approach by Thorsen et al.(2001), where precipitation is assumed to be approximately spatially constant at the length scale of the catchment (~25 km). This imposes limitations as to the interpretation of water flow and leaching processes at spatial scales smaller than the catchment scale—that is, field level—but is not expected to severely impact the lumped discharge response at the monitoring station used in this study. Hansen et al. (1999) varied precipitation in the Karup Basin by adding a random error (zero mean and a standard deviation of 50%) to measured precipitation series and found that precipitation contributed least to the uncertainty associated to water balance, mineralization, and plant N uptake. Thorsen et al. (2001) for the same Karup Basin introduced the same uncertainty to the precipitation time series as in Hansen et al. (1999) and found that this affected the simulated yearly water balance but not the simulated nitrogen balance, suggesting that the timing of the percolated water controlled by the soil hydraulic properties is more important for the simulated nitrogen loads than the total annual amounts of percolation. In this study, it is assumed that with long simulation periods (25 and 55 yr) with the year-to-year variation, the uncertainty in precipitation is sufficiently taken into account.
  3. Reducing the soil texture classes in the Odense River catchment from four to one. The four classes included were sandy loam (JB6), coarse sand (JB4), peat (JB11), and glacial clay (JB8), covering 74, 21, 3.5, and 1.5%, respectively. This reduction is justified because both JB4 and JB6 are classified as sandy loam according to the USDA classification system, and differences in simulated water percolation and nitrate leaching from the rootzone are likely to be small (together, JB4 and JB6 represent 95% of the soil types within the catchment). The soil in the area is therefore represented by a JB6, sandy loam soil.
  4. Reducing the number of land management types. The number of crop rotations is reduced to one for each group of cattle and pig farms, plant production, low vegetation (permanent grass), and forest. The selected crop rotations were visually inspected and seemed to capture most of the variation in total leached nitrate from the rootzone, as well as the total amount of leached nitrate as simulated by the Nielsen et al. (2004) model. The cropping patterns related to cattle and pig farms and plant production have to be permutated to ensure that each crop is equally represented for each climate year. Hence, the total number of crop rotations amounts to 15, that is, four permutations for the cattle, pigs, and plant production land management types, and only one rotation for deciduous and coniferous forest, as well as for permanent grass.
  5. The most dominant lower boundary of the Odense river catchment is a varying groundwater level between 0.75 and 1.5 m below surface constituting 90% and an average of 1.25 m (Nielsen et al., 2004). Therefore, the DAISY lower boundary option of drain pipes at 1 m depth is applied for the sandy loam soil type.

In total, the number of DAISY soil columns are in this way reduced from 5577 to 42.

Setup of Modified DAISY
The simplified DAISY setup is derived from the original model setup as described in Nielsen et al. (2004) and outlined in the previous section. In this study, the simplified model is not calibrated against crop yield as it is assumed that the overall dynamics of the simplified model is preserved (see "Comparison of Simplified and Original Models" under "Results and Discussion").

Two sets of simulations are considered. The first set comprises simulations run from 1 Jan. 1989 to 31 Dec. 2013, where the climatic forcings of the years 1992 to 2002 are repeated for 2003 to 2013, in total, 25 yr. The first 10 yr are run ("warming up") to allow the organic matter pools reach a "natural" state to avoid strong dependence of their initial state on output variables (Styczen et al., 2004). Yearly accumulated nitrate leaching after the warming-up period from the rootzone is recorded for the uncertainty analysis.

The second set of simulations is derived from the first by ranking the outcome of accumulated nitrate leaching from the rootzone from the 25-yr period (15 yr for accumulated nitrate output). Then, 1st, 5th, 10th, 25th, 50th, 75th, 90th, 95th, and 99th percentile runs are selected, covering the whole output distribution for simulation of solute transport over a period of 25 yr. Next, the parameter sets selected corresponding to the output percentiles are run for the extended period of 55 yr by repeating three times the meteorological input from 1992 to 2002 beyond 2013. As argued previously, the longer period most likely ensures an improved approximation of soil organic matter pools equilibrium (Styczen et al., 2004). In addition, the longer simulation period facilitates the modeling process for model-specific reasons. The time series for those selected runs for water percolation and nitrate leaching from the rootzone can then be further propagated into the groundwater model. It is important to note that the output distribution from the 25-yr simulation series is not proportionally sampled according to the distribution densities. Percentile values are selected from the distribution to ensure that the entire output distribution is covered and to facilitate a simple procedure for selecting from the distribution.

Because of the input requirements to the groundwater model, analyses of uncertainty in DAISY time series for both periods are restricted to water percolation and nitrate leaching. The overall setup for DAISY is summarized in Table 1 .


View this table:
[in this window]
[in a new window]

  		TABLE 1. Overview of the DAISY setup. Right-hand column lists which parameters are varied in Monte Carlo analysis.

 
Identifying and Quantifying Key Sources of Uncertainty for the DAISY Model (Step 2)
DAISY model parameters were preselected on the basis of literature and are summarized in Table 2 . The contribution of the selected parameters to uncertainty in simulated leaching is analyzed below in "Uncertainty Analysis (Step 5)." Statistics on soil hydraulic properties (Table 2 properties 1–4) for the van Genuchten–Mualem retention–hydraulic conductivity functions are taken from Meyer et al. (1997). Of the included properties, soil hydraulic properties (Table 2, properties 1–4) relate to transport of nitrate from the rootzone, whereas slurry amount (Table 2, property 5) relates to turnover processes. Root depth is known to contribute to uncertainty related to soil water balance (e.g., Seyfried and Wilcox, 2006) and is included in the uncertainty analysis.


View this table:
[in this window]
[in a new window]

  		TABLE 2. Parameter uncertainties that are assumed most important. Properties 1–4 are described in the section "Soil Hydraulic Properties."

 
The following sections elaborate on the quantification of parameter uncertainties listed in Tables 3–5GoGo , with a focus on DAISY simulation. Properties 1 through 6 (Table 2) are further explored for use in the uncertainty analysis.


View this table:
[in this window]
[in a new window]

  		TABLE 3. Hydraulic parameters statistics for soil class sandy loam (Meyer et al., 1997).

 

View this table:
[in this window]
[in a new window]

  		TABLE 4. Hydraulic parameters statistics for A, B, and C horizons for JB6 (Styczen et al., 2004).

 

View this table:
[in this window]
[in a new window]

  		TABLE 5. Soil hydraulic and N-related parameter statistics after Monte Carlo sampling for JB6 and crop rotation for four cattle farm (cs1–4) and three pig farm (ps1–3) slurry amounts and one plant production (pp) fertilizer amount applied for different crops at different times in tonnes wet weight per hectare. Interquartile range (IQR) is 75th–25th percentile.

 
Soil Hydraulic Properties
For the original model (Nielsen et al., 2004), hydraulic properties were estimated through the HYPRES pedotransfer functions (Wösten et al., 1999). In the present work, the Danish classification JB6 soil type is reclassified to the USDA soil system in order to apply the soil hydraulic property statistics derived for the USDA system. The JB6 soil class is equivalent to a USDA sandy loam, and the soil hydraulic parameter statistics—probability density functions (pdf) and correlations between parameters from Meyer et al. (1997)—for sandy loam are applied. Meyer et al. (1997) derived these distribution functions by resampling the original distributions of Carsel and Parrish (1988) based on soil pedotransfer functions from Rawls and Brakensiek (1985) and also included correlations between parameters.

DAISY requires a parameterization of soil hydraulic functions. In this case the van Genuchten (1980) soil water release characteristics curve Se(h) and the unsaturated hydraulic conductivity function K(Se) (Mualem, 1976) function are applied:

Formula 1[1]

Formula 2[2]
where n and {alpha} are van Genuchten fitting parameters, and m is assumed to equal 1 – 1/n, h denotes soil water potential, {theta} is soil water content, and {theta}s and {theta}r are saturated and residual water content, respectively. The statistic properties of the hydraulic van Genuchten parameters n and {alpha}, as well as {theta}s and {theta}r, are taken from Meyer et al. (1997) for USDA soil texture classes, including parameter correlations, and are listed in Table 3.

The actual values of hydraulic parameters for the A, B, and C horizons are taken from recommendations by Styczen et al. (2004) for a JB6 soil (Table 4), and the statistical properties are scaled relative to the mean according to Meyer et al. (1997) for a sandy loam as shown in Table 3. Thus, the mean values in Table 3 for a sandy loam were rescaled to fit the JB6 parameter values in Table 4 and the statistical properties accordingly. Residual water content ({theta}r) was set to 0.0 by Styczen et al. (2004) and therefore not included in Table 4.

The clay content is not varied since clay content for water percolation is explicitly taken into account by varying the hydraulic properties and does not contribute to the main uncertainties for simulated mineralization (Hansen et al., 1999). Because root depth is known to be important for the rootzone water balance (Seyfried and Wilcox, 2006), this parameter is included in the present analyses. Its variation is estimated by expert judgment (Table 5). Although hydraulic properties vary slightly across the soil profile (Table 4), to keep the number of varied parameters as low as possible for the LHS, the coefficients n, {alpha}, and {theta}s for horizons A, B and C are first averaged and subsequently sampled for the Monte Carlo–based uncertainty analysis. For the same reason, the hydraulic conductivity Ksat for A and B horizons are (lognormally) averaged while Ksat for the C horizon is kept. Clearly, the variation of the soil hydraulic properties reflects the variability in texture within one soil class and does not account for the heterogeneity at the basin length scale, that is, the systematic uncertainty of hydraulic parameters at the spatial scale of the Odense basin. This is related to the question as to what extent the pedotransfer function that translates soil physical properties to hydraulic parameters is representative at the catchment scale (Wösten et al., 2001). The implications of this are discussed further under "Results and Discussion."

Soil Organic Matter and Nitrogen Content
Nitrogen content and dry matter in cattle and pig slurry, as sampled by Hansen et al. (1999) and Thorsen et al. (2001), are substituted with the amount of slurry since little is known about the variability of slurry composition. For initialization of the DAISY mineralization model, SOM2 is important; uncertainty associated with this parameter is substantial, and an assessment often relies on expert judgment (Hansen et al., 1999; Thorsen et al., 2001). The initial SOM allocated to SOM2 is taken from Nielsen et al. (2004). An overview of all included parameters and their statistical properties is provided in Table 5. A measure for the parameter statistics is the coefficient of variation (CV), which can be defined as the interquartile range (IQR)—that is, the difference between the 25th and 75th values—normalized by the median (e.g., Christiaens and Feyen, 2001), and is included in Table 5.

Latin Hypercube Sampling (Step 3)
Step 3 of the framework involves analyzing whether the uncertainty in selected model parameters has significance for simulated water percolation and nitrate leaching. The uncertainty analysis of DAISY model parameters with respect to accumulated nitrate leaching from the rootzone for the simulated period of 25 yr was conducted by sampling from the model parameter distributions using the effective LHS technique (McKay et al., 1979; Iman et al., 1980). The LHS technique is a stratified sampling rather than a random sampling in the traditional Monte Carlo technique and was also used by Hansen et al. (1999) and Thorsen et al. (2001) for the DAISY model. The UNCSAM (Janssen et al., 1992) statistical software package was used for the LHS technique and for postprocessing of the data for the uncertainty analysis. Table 5 shows all the parameters varied for each land management type. The hydraulic parameters, vG{alpha}, vGn, {theta}s, Ks-u, Ks-l (Tables 3 and 4), and root depth, the latter assumed normally distributed (Dubus et al., 2000), are varied for each management type. In addition, four cattle slurry (cs1–4), three pig slurry (ps1–3) amounts, and one plant fertilizer amount (wet weight) are varied. Data on the distribution type for fertilizer amounts are not known, and therefore a normal distribution is assumed. Hence, the number of model parameters considered in the uncertainty analysis is 10 (6 + 4), 9 (6 + 3), and 7 (6 + 1) for cattle farms, pig farms, and plant production, respectively. For the coniferous and deciduous forests and for permanent grass, only hydraulic parameters are varied; that is, six parameters. Correlation between hydraulic parameters as listed in Table 3 have been taken into account in the LHS sampling. Correlation between slurry amount parameters is not known and not included here. The sample size for LHS has been the subject of studies for quite some time (e.g., Stein, 1987). Recently, Hansen et al. (1999) argued that 25 LHS runs are sufficient for the six parameters they varied for the Karup catchment in Denmark. Iman and Helton (1985) argued that four-thirds the number of parameters is sufficient, whereas Christiaens (2001) recommends between two- and fivefold the number of parameters involved. In the present study, we adopted a safe number of 100 runs. The statistical properties of the LHS sampled parameters are listed in Table 5.

Simulation of Latin Hypercube Sampling–Generated DAISY Realizations (Step 4)
Data flow and file input/output for accumulated data for water and solute transport related output variables was automated through a DAISY–SENSAN interface. The SENSAN package developed for sensitivity and uncertainty analyses is part of PEST (Doherty et al., 1995) and controls data flow from the DAISY model to output files that subsequently are processed by auxiliary software for further processing by UNCSAM.

Uncertainty Analysis (Step 5)
The results of the analysis under the framework's Step 5 as described here are used to confirm the significance of the preselected parameters under Step 2 and crucial for proceeding to Step 6.

The variation in simulated output Y(t) as a result of the varying parameter input X(t) by means of the LHS simulations is captured by a linear regression model that fits Y(t) and X(t). This procedure is described further in Janssen et al. (1992) and Christiaens (2001). The regression coefficient (SRRC) was used by performing a regression on rank-transformed data rather than raw standardized data since this procedure resulted in higher R2 values. Moreover, the SRRC index has been shown to be one of the most robust and reliable means of assessing model sensitivity (Dubus and Janssen, 2003; Saltelli and Marivoet, 1990) and allows a comparison of the relative contribution of each input parameter in the prediction of the model (Hamby, 1994). The SRRC index method has been applied successfully in pesticide leaching vulnerability studies (Dubus and Brown, 2002; Dubus et al., 2000) but has not yet been applied to nitrate leaching studies. The linear regression model for SRRC can be written as (Janssen et al., 1992)

Formula 3[3]
where Yranked,s(o) and Xranked,s(i) are rank transformed and standardized output and input, respectively, and {varepsilon} is the regression model error. The higher the SRRC value, the more influence on model predictions this parameter has. When applying linear regression analysis, accuracy and numerical problems can occur if the dependent variables X(i) show too high correlations (collinearity problem). Snee (1983) proposed the variance inflation factor (VIF) as a measure for collinearity; VIF is further discussed in Janssen et al. (1992). In the case of no mutual correlation, VIF = 1, and according to Snee (1983), linear regression can be applied for VIF < 10. In this study, the VIF is always less than 10 (Table 6 ).


View this table:
[in this window]
[in a new window]

  		TABLE 6. R2 and colinearity for different management types.

 
Ranking DAISY Output and Simulate Extended Period (Steps 6 and 7)
In Steps 6 and 7 of the framework, simulated accumulated nitrate leaching from the rootzone for the 25-yr period (not including the 10-yr warming-up period) is ranked, and the LHS-generated DAISY parameter sets corresponding to the 1st, 5th, 10th, 25th, 50th, 75th, 90th, 95th, and 99th percentile leached nitrate are prepared for the extended 55-yr simulation period. The consistency between leached output from the selected percentiles from the 25-yr period is then checked with the ranked output based on the 55-yr series. This means that the 99th percentile parameter set should lead to the highest simulated nitrate leaching for the long period, and vice versa. If so, a strongly reduced number of parameter sets would suffice while preserving a realistic range of rootzone model output for further groundwater modeling.


   Results and Discussion
TOP
 ABSTRACT
 INTRODUCTION
 Case Study Location
 Methodology
 Results and Discussion
Conclusions
 REFERENCES
 
Comparing Simplified and Original Models (Step 1)
The hydrograph of the simplified MIKE-SHE/DAISY model against observed values at Odense River station 45-26 compares well with the original model of Nielsen et al. (2004). The hydrograph of the simplified model as shown in Fig. 3 is the 50th percentile (median) ranked DAISY output for all land management types propagated through the MIKE-SHE model, resulting in simulated discharge at a monitoring station for the Odense River (45-26) for a period of 4 out of a total of 25 yr. The Nash–Sutcliffe model efficiency coefficient E (Nash and Sutcliffe, 1970) is used to assess the explained variance of the simplified model compared with the original model with respect to simulated discharge. The coefficient E is defined as

Formula 4[4]
In Eq. [4], Qtsim and Qsimmean are discharge of simplified model at time t, and averaged over entire time period T, respectively; Qtorig is discharge of the original model. The Nash–Sutcliffe coefficient E for the period 1990 to 2000 is 0.83, implying a rather good correspondence between the two models with respect to discharge. The correlation coefficient (r) for the same period is 0.97. The daily discharge simulated for the period 1990 to 2000 at the monitoring station for the simplified vs. the original is shown in Fig. 4 . The simplified model was not calibrated against harvest since the main objective is an uncertainty analysis, where focus is on relative differences in simulated nitrate leaching as a result of varying DAISY parameters under equal boundary conditions, which justifies a relaxation toward calibration against harvest of dry matter and nitrogen. Also, the (calibrated) parameterization of Nielsen et al. (2004) has been adopted to facilitate comparison with the original model. Independent calibration could have changed this, making direct comparison between the original and simplified model as shown in Fig. 3 and 4 more difficult. Clearly, for environmental assessment, using the DAISY model as embedded in the original MIKE-SHE/DAISY model, calibration against harvest data as performed in Nielsen et al. (2004) is critically important.


Figure 3
View larger version (21K):
[in this window]
[in a new window]

  		FIG. 3. Hydrograph of modified MIKE-SHE/DAISY model (50th percentile ranked time series) against original model (Nielsen et al., 2004) and observed data at Odense River station 45-26.

 

Figure 4
View larger version (22K):
[in this window]
[in a new window]

  		FIG. 4. Simulated daily discharge at Odense River station 45-26 for the period 1 Jan. 1990 to 31 Dec. 2000 of simplified model vs. original model (Nielsen et al., 2004).

 
Uncertainty in Simulated Nitrate Leaching and Water Percolation (Steps 2–4)
The results of identification of uncertain parameters and subsequent simulations based on LHS-generated DAISY parameter sets are presented in this section. Simulated nitrate leaching and water percolation from the rootzone for the six land management types and four permutations (p1, p2, p3, and p4) are shown in Fig. 5 and 6 , respectively. The CV, defined as the IQR normalized by the median, is also indicated (Table 5).


Figure 5
View larger version (20K):
[in this window]
[in a new window]

  		FIG. 5. Box-plot diagram of simulated accumulated leached nitrate per year (25-yr period) for cattle and pig farms and plant production for permutations P1–P4, coniferous (C) and deciduous (D) forest, and grass (G). The upper and lower levels of the boxes denote 75th and 25th percentile values, respectively; the horizontal line in the boxes is the median value. The upper and lower values are maximum and minimum values, respectively. The upper value is the interquartile range/median (CV).

 

Figure 6
View larger version (18K):
[in this window]
[in a new window]

  		FIG. 6. Box-plot diagram of simulated accumulated water percolation per year for cattle and pig farms and plant production for permutations P1–P4, coniferous (C)- and deciduous (D) forest, and grass (G). The upper and lower levels of the boxes denote 75% and 25% percentile values, respectively; the horizontal line in the boxes is the median value. The upper and lower values are maximum and minimum values, respectively. The lower value is the interquartile range/median (CV).

 
It is clear from Fig. 5 and 6 that the CV for nitrate leaching is five to eight times larger than the CV of water percolation and that differences for nitrate leaching between different land management types are larger than between permutations of the same crop rotation. For both water percolation and nitrate leaching, differences due to permutations are not significant. For water percolation, there is little difference between land management types and permutations, except for the two forest types. This is in accordance with the findings of Hansen et al. (1999) for the management types in the Karup catchment; however, the CV values for leaching are generally higher than reported by Hansen et al. (1999) and Thorsen et al. (2001), whereas percolation is on the same order of magnitude. Variation in simulated nitrate leaching can indeed be very high; Kyllmar et al. (2005) under Swedish conditions and Hansen et al. (1999) for the Karup catchment found similar ranges as shown in Fig. 5. The modest and appreciable differences between land management types with respect to simulated water percolation and nitrate leaching from the rootzone respectively indicate, as expected, that simulated leaching contributes most to the overall uncertainty. The rather high CV values for leaching confirm this and underline the need to take into account uncertainty in rootzone model parameters at the field scale. To arrive at an uncertainty assessment at the scale of the Odense catchment, the uncertain DAISY time series for water percolation and nitrate leaching must be distributed over the catchment following the simplified representation. This was illustrated for the 50th percentile and shown in Fig. 3 and 4.

Parameter Uncertainty Analysis (Step 5)
Table 6 shows the overall R2 as well as a measure for collinearity, the VIF. The overall R2 for fitting the output variable "leached nitrate" to a linear combination of hydraulic parameters and nitrogen input shown in Table 7 ranges from 0.78 to 0.90. The VIF is close to 1 for all management types, implying near uncorrelated parameter values in the regression analysis. As mentioned above, the largest VIF must not exceed 10 for applying a linear regression model. For the measure of uncertainty, the SRRC (Eq. [3]) is used, as explained in the section "Simulation of Latin Hypercube Sampling–Generated DAISY Realizations (Step 4)" above.


View this table:
[in this window]
[in a new window]

  		TABLE 7. Standardized rank regression coefficient (SRRC) and ranking (in parentheses) for water percolation for different management types for 25-yr simulation.

 
Tables 7 and 8 provide the SRRC values for accumulated water percolation and nitrate leaching from the rootzone for a period of 15 yr. The higher (absolute) SRRC values, and not exceeding 1, the more sensitive the parameter is relative to the others and thus contributing to a higher extend to the overall uncertainty. Note that the SRRC value is not the percentage variance explained. The UNCSAM package enables a t test for calculated SRRC values, and numbers with asterisks in Table 7 and 8 indicate a SRRC value significantly different from 0. The results in Table 7 and 8 show that the hydraulic parameters {theta}s and n are very sensitive to both water percolation and nitrate leaching from the rootzone. Hydraulic conductivity for the upper and lower zones is ranked lower in this respect but is significant in most cases. Similar results for percolated water were obtained by Dubus et al. (2000) using the MACRO model in a pesticide leaching vulnerability study. The low and almost always insignificant values for root length suggest that this parameter has negligible influence on both percolation and leaching under the conditions here. The insensitivity of the rooting depth was unexpected because other studies have shown that this parameter greatly affects the simulated percolation and leaching. However, it is suspected that a deep water table was simulated with the applied DAISY setup, making the water uptake and denitrification less subject to the rooting depth. Clearly, medium- to high-ranked and significant SRRC values indicate that the amount of applied slurry has importance for crops and cattle and pig farms with respect to nitrate leaching and none for water percolation. It may be argued that more potentially uncertain model parameters should have been included in the LHS analysis. However, several studies have indicated that the parameters included in the present study are among the most important in this respect.


View this table:
[in this window]
[in a new window]

  		TABLE 8. Standardized rank regression coefficient (SRRC) and ranking (in parentheses) for nitrate leaching for different management types for 25-yr simulation.

 
Simulations Based on 55-Year Series (Steps 6 and 7)
The output statistics for the extended period of 55 yr for the 1st, 5th, 10th, 25th, 50th, 75th, 90th, 95th, and 99th percentile runs are listed in Tables 9 and 10 . As for the average accumulated yearly nitrate leaching, it is clear from Fig. 5 and Tables 9 and 10 that 25th, 50th (median), and 75th percentiles are generally higher for the 55-yr simulations compared with the 25-yr, except for pig farms that have the same order of magnitude. The maximum and minimum accumulated yearly water percolation from the rootzone for the long 55-yr simulation are also higher than the 97.5th and 2.5th percentiles for the shorter 25-yr simulations in Fig. 5 and Tables 9 and 10, whereas ranges, calculated as the difference between maximum and minimum values, are comparable to each other. Importantly, the ranking of nitrate leaching derived from the LHS simulations for 25 yr is preserved for the longer 55-yr simulations (Tables 9 and 10). This means that DAISY parameter sets resulting in high percentile accumulated leaching for the 25-yr simulation also result in the high accumulated value for the 55-yr simulation, and vice versa for the lower percentile and the lowest value. This implies that it is sufficient to sample rootzone parameter sets according to the output distribution from the short 25-yr simulations and propagate these sets to the groundwater model for extended simulations and subsequent uncertainty analyses.


View this table:
[in this window]
[in a new window]

  		TABLE 9. Output statistics for accumulated nitrate leaching and water percolation from rootzone for 55 yr simulation for two management types including four permutations (p1–p4). The first column is the percentile corresponding to the ranking of 100 runs for 25 yr.

 

View this table:
[in this window]
[in a new window]

  		TABLE 10. Output statistics for accumulated nitrate leaching and water percolation from rootzone for 55-yr simulation for five management types, including four permutations (p1–p4) for plant production. The first column is the percentile corresponding to the ranking of 100 runs for 25 yr.

 
Limitations of Framework Approach
There are several limitations to the proposed framework. The first relates to the question of how uncertainty at the catchment scale is approached. Contrary to the method used by Thorsen et al. (2001), where uncertainty at the catchment scale was estimated by not allowing spatial variation between model grid and a smaller catchment scale variance was used, this framework uses a larger field-scale variance for hydraulic and land management parameters. The field-scale parameter uncertainty is then propagated through the rootzone model and subsequently distributed over the catchment for groundwater modeling to arrive at an aggregated catchment scale uncertainty. Hence, focus is on uncertainty of input variables to the DAISY model, which is estimated from within soil class variability for a sandy loam, whereas uncertainty of land management parameters is estimated from expert judgment. Not included in this approach is a systematic deviation from the mean value for model parameters as estimated from the class pedotransfer function used here (Meyer et al., 1997). This implies that representativity of the soil class at the catchment scale due to differences in geological setting or soil structure may lead to a systematic deviation of hydraulic behavior at the catchment scale. This limitation boils down to model structure uncertainty (Refsgaard et al., 2006) on the applied pedotransfer function and can be addressed by considering multiple pedotransfer function models in estimation of hydraulic and transport parameters (unpublished data). Another related limitation of the approach is connected to the rootzone modeling itself and model structure uncertainty; that is, to what extent the DAISY model process dynamics regarding water flow, solute transport, and nitrate turnover are correctly captured. This too requires, in principle, incorporation of multiple process descriptions (model structure), such as a GLUE approach (Beven and Binley, 1992). It is, however, outside the scope of the present study to do such an analysis since the same process description for the DAISY model is used in both the original and simplified model. Another limitation of the approach is whether a simplification as done here can be justified. This is addressed in the section "Modification of the DAISY Component of the Original Distributed Model," above, for the Odense catchment and it seems to apply generally to catchments with a similar degree of homogeneity. If soil types cannot be grouped as demonstrated here and stronger gradients in meteorological conditions exist, a simplification may not apply. Finally, the preselection of parameters is assumed to contribute significantly to output uncertainty. It may be argued that other model parameters should have been included as well. However, a rigorous systematic sensitivity analysis for the model system used does not exist, and the relevance of the selected parameters is well documented.


   Conclusions
TOP
 ABSTRACT
 INTRODUCTION
 Case Study Location
 Methodology
 Results and Discussion
 Conclusions
REFERENCES
 
We have proposed a framework for uncertainty assessment with respect to nitrate leaching at the catchment scale and applied it to the Odense catchment. The framework is based on deriving uncertain output from a rootzone model at the field scale for further propagation to a groundwater model at the catchment scale, exemplified here by the Odense River catchment. DAISY model parameters that contribute most to the uncertainty in simulated nitrate leaching from the rootzone have been identified by computing the SRRC value on ranked accumulated leaching output for a 25-yr period. Parameters related to the van Genuchten soil water release characteristics and the Mualem unsaturated conductivity function are important, as is the amount of applied slurry representing the N content in slurry. It was found that crop sequence (permutations) for 25-yr simulations has little impact on accumulated percolation, with CV, calculated as the IQR/median, about 10, 7.5, and 9 for cattle farms, pig farms, and plant production, respectively. Permutations affect nitrate leaching to a higher extent than percolation, with CV ranges from 65 to 81 for cattle farms and 42 to 53 for pig farms. The CV for leached nitrate shows a range from ~15 for forest to 80 for cattle farms. The ranking of nitrate leaching derived from the LHS simulations for the shorter 25-yr period are, with few exceptions, preserved for the longer 55-yr simulations. Limitations of the framework are related to (i) the extent to which the modeling system can be simplified under different conditions, (ii) the preselection of parameters assumed to contribute significantly to model output, and (iii) the uncertainties that are not included in this approach, mainly, model structure uncertainty. The proposed framework applies to the Odense catchment and has a varying degree of general applicability depending on the properties of the catchment. This work addresses the problem of uncertainty assessment of hydrologic modeling at the field and catchment scales using numerical hydrologic models that can potentially describe the substantial complexity of a catchment with regard to both the conceptual model of the catchment and process descriptions. The framework provides an easily applicable procedure for selecting rootzone model output time series covering the entire output distribution for further propagation to a catchment-scale groundwater model. Future work should focus on how to capture uncertainty associated with model structure at the catchment scale.


   ACKNOWLEDGMENTS
 
This paper builds on results obtained from the EU-FP5 project Harmonised Techniques and Representative River Basin Data for assessment and Use of Uncertainty Information in Integrated Water Management (HarmoniRiB) under contract EVK-CT-2002-00109.


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 Case Study Location
 Methodology
 Results and Discussion
 Conclusions
 REFERENCES
 





This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Right arrow reprints & permissions
Citing Articles
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by van der Keur, P.
Right arrow Articles by Refsgaard, J. C.
Right arrow Search for Related Content
PubMed
Right arrow Articles by van der Keur, P.
Right arrow Articles by Refsgaard, J. C.
GeoRef
Right arrow GeoRef Citation
Agricola
Right arrow Articles by van der Keur, P.
Right arrow Articles by Refsgaard, J. C.
Related Collections
Right arrow Water Quality
Right arrow Nutrients

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Soil Science Society of America Journal
Journal of Plant Registrations Journal of
Environmental Quality		 The Plant Genome