HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS		QUICK SEARCH:   [advanced] Author: Keyword(s): Year:  Vol:  Page:

This Article Abstract Figures Only Full Text (PDF) Free Alert me when this article is cited Alert me if a correction is posted Services Related articles in JEQ Similar articles in this journal Similar articles in ISI Web of Science Similar articles in PubMed Alert me to new issues of the journal Download to citation manager Citing Articles Citing Articles via HighWire Citing Articles via ISI Web of Science (8) Citing Articles via Google Scholar Google Scholar Articles by Larsbo, M. Articles by Jarvis, N. Search for Related Content PubMed PubMed Citation Articles by Larsbo, M. Articles by Jarvis, N. Agricola Articles by Larsbo, M. Articles by Jarvis, N. Related Collections Vadose Zone Processes and Chemical Transport Pesticides Dual Porosity/Permeability Models Solute Transport Models

TECHNICAL REPORTS

Vadose Zone Processes and Chemical Transport

Simulating Solute Transport in a Structured Field Soil

Uncertainty in Parameter Identification and Predictions

Department of Soil Sciences, SLU, Box 7014, 750 07 Uppsala, Sweden

^* Corresponding author (Mats.Larsbo{at}mv.slu.se)

Received for publication April 13, 2004.

	ABSTRACT

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
Dual-permeability models have been developed to account for the significant effects of macropore flow on contaminant transport, but their use is hampered by difficulties in estimating the additional parameters required. Therefore, our objective was to evaluate data requirements for parameter identification for predictive modeling with the dual-permeability model MACRO. Two different approaches were compared: sequential uncertainty fitting (SUFI) and generalized likelihood uncertainty estimation (GLUE). We investigated six parameters controlling macropore flow and pesticide sorption and degradation, applying MACRO to a comprehensive field data set of bromide andbentazone [3-isopropyl-1H-2,1,3-benzothiadiazin-4(3H)-one-2,2dioxide] transport in a structured soil. The GLUE analyses of parameter conditioning for different combinations of observations showed that both resident and flux concentrations were needed to obtain highly conditioned and unbiased parameters and that observations of tracer transport generally improved the conditioning of macropore flow parameters. The GLUE "behavioral" parameter sets covered wider parameter ranges than the SUFI posterior uncertainty domains. Nevertheless, estimation uncertainty ranges defined by the 5th and 95th percentiles were similar and many simulations randomly sampled from the SUFI posterior uncertainty domains had negative model efficiencies (minimum of –3.2). This is because parameter correlations are neglected in SUFI and the posterior uncertainty domains were not always determined correctly. For the same reasons, uncertainty ranges for predictions of bentazone losses through drainflow for good agricultural practice in southern Sweden were 27% larger for SUFI compared with GLUE. Although SUFI proved to be an efficient parameter estimation tool, GLUE seems better suited as a method of uncertainty estimation for predictions.

Abbreviations: EF, model efficiency • GLUE, generalized likelihood uncertainty estimation • SUFI, sequential uncertainty fitting

	INTRODUCTION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
THERE IS TODAY a widespread acceptance that preferential flow of pesticides is an important process in structured soils that contributes to ground water and surface water contamination (Flury, 1996; Jarvis, 2002). Dual-permeability solute transport models account for preferential flow by including a separate flow domain describing rapid non-equilibrium flow in soil macropores with first-order mass exchange between the two domains (Feyen et al., 1998; imnek et al., 2003). In this respect they are more complex than models based solely on the convection–dispersion equation since they require additional parameters to describe macropore flow. One difficulty with these models has been the lack of reliable procedures for estimating these parameters (Forum for the Coordination of Pesticide Fate Models and Their Use, 1995; imnek et al., 2003), since some are either difficult or impossible to measure. A calibration procedure is, therefore, required. Intuitively, the greater degree of complexity of dual-permeability models would imply an increase in the data requirements for efficient calibration (Jarvis, 1999, 2001), but no attempts have yet been made to quantify this.

Inverse modeling has become a popular calibration method since in some respects it limits the subjectivity in the process and because it provides quantitative estimates of parameter uncertainty. Recently, attempts have been made with varying success to use inverse modeling to parameterize preferential flow models on data from column experiments (Schwartz et al., 2000; Kätterer et al., 2001; Roulier and Jarvis, 2003). They all found that the parameter determining the mass exchange between flow domains was difficult to identify. This was attributed to parameter insensitivity (Kätterer et al., 2001; Roulier and Jarvis, 2003), a sparse parameter sampling scheme (Roulier and Jarvis, 2003), and failure of the first-order mass transfer concept used in the model to handle different time scales (Schwartz et al., 2000). Roulier and Jarvis (2003) also suggested that higher time resolution data for both flux and resident concentrations might facilitate parameter identification. To date, no attempts have been made to parameterize dual-permeability solute transport models by inverse modeling using field data.

Even though calibration has become a cornerstone in pesticide fate modeling (Dubus et al., 2002) it may not always be possible to find an optimal parameter set that best describes the observations. Multiple combinations of input parameters will often provide equally good fits to experimental data (Beven and Binley, 1992). This is often referred to as "equifinality" (Beven, 1993). Equifinality occurs, for example, if different time periods are simulated well by different combinations of parameters or if the parameters are correlated. Furthermore, in the case of multiple data sets a multi-objective goal function must be used. A characteristic of a multi-objective problem is that the solution will not, in general, be unique (Yapo et al., 1998). The maximization of an individual goal function often leads to the deterioration of others, which means that many combinations of individual goal function values can result in the same multi-objective goal function value (often referred to as pareto-optimality). An increasing awareness of these problems, in addition to the well-known uncertainties and errors in field measurements and the effects of spatial heterogeneity (Dubus et al., 2003), may lead to the abandoning of the traditional quest for one optimal parameterization of a model in favor of parameters conditioned on observations of a system where the parameter uncertainties are made explicit. However, the selection and implementation of methods designed to account for uncertainties involves a number of subjective choices including the selection of parameters to be included in the analysis, the distribution function of input parameters, the way correlations are handled, and the sampling scheme used (Dubus et al., 2003). Furthermore, the parameter uncertainties must in some way be transferred into uncertainties in model predictions to be of any value in environmental risk assessment.

The objectives of this study were to evaluate the data requirements for efficient parameter identification in a solute transport model accounting for macropore flow and to compare two conceptually different calibration and uncertainty estimation methodologies, sequential uncertainty fitting (SUFI) (Abbaspour et al., 1997) and generalized likelihood uncertainty estimation (GLUE) (Beven and Binley, 1992). The model MACRO 5.0 (Larsbo and Jarvis, 2003) was used to simulate water flow and transport of bromide and the weakly sorbed herbicide bentazone in a drained field during a one-year period. Parameters controlling macropore flow and pesticide sorption and degradation were determined through inverse modeling by applying the iterative parameter estimation procedure SUFI to a comprehensive data set of soil water content, drainflow, and resident and flux concentrations of bromide and bentazone. The GLUE procedure was also applied to the same field observations to enable a comparison between SUFI and GLUE. Moreover, GLUE was used to examine the significance of different groups of field observations for effective conditioning of selected model parameters. The field experiment did not represent "good agricultural practice" (GAP) since bentazone was applied on bare soil in autumn at nearly double the maximum recommended dose. The results from the SUFI and GLUE analyses were, therefore, translated into predictions of maximum concentrations in drainflow and total loss of bentazone through field drains for a typical GAP scenario for bentazone in Sweden. This was accomplished both by a straight-forward Monte Carlo approach based on SUFI calibrations and by the GLUE approach. The significance of the availability of different groups of field observations for uncertainty in predictions was also analyzed.

	MATERIALS AND METHODS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
The Field Site and Measured Data
The model simulations presented here are based on field experiments at Lanna, Sweden (58°21' N, 13°08' E; Fig. 1) . The soil is a well-structured silty clay (Table 1), classified as a Typic Eutrochrept (USDA), with three tile drains installed at a 13.5-m spacing and 1-m depth, draining a plot 0.4 ha in size (Fig. 1). Lanna is situated on a flat plain, with the slope at the field site less than 1%. A simultaneous dose of potassium bromide (44.4 kg Br^– ha^–1) and the weakly sorbed herbicide bentazone (2.51 kg ha^–1) was applied by a tractor-mounted pesticide sprayer on 18 Oct 1994. The crop grown during the summer of 1995 was spring-sown rape (Brassica napus L.). A comprehensive field data set of flows and storages of water, bromide, and bentazone is available for model testing and calibration (Table 2). A more detailed description of the experiment is reported in Larsson and Jarvis (1999).

View larger version (30K): [in this window] [in a new window]   Fig. 1. Location map and layout of the drainage system in the experimental field at Lanna.

The division between flow domains is given by a water potential, _b (m), and the corresponding saturated water content, _b (m³ m^–3), and hydraulic conductivity, K_b (m s^–1), in the micropores. Water flow in the micropores is governed by Richards' equation:

[1]

where C = / (m^–1) is the differential water capacity, (m³ m^–3) is the volumetric micropore water content, (m) is the soil water pressure head, t (s) is time, z (m) is depth, K (m s^–1) is the unsaturated hydraulic conductivity, and S_d, S_r, and S_w (s^–1) are source–sink terms accounting for drainage, root water uptake, and water exchange with macropores, respectively. In the micropores, the water retention curve () is given by the van Genuchten (1980) function whereas the hydraulic conductivity function K() is given by Mualem's (1976) model. Water flow in the macropores, q_ma (m s^–1), is described by a modified kinematic wave approach (Germann and Beven, 1985), where the macropores are assumed to drain by gravity only. The hydraulic conductivity in the macropores, K_ma (m s^–1), is expressed as a power function of the degree of saturation in the macropores, S_ma:

[2]

where K_s (m s^–1) is the saturated conductivity of the total pore system and n* (unitless) is a "kinematic" exponent reflecting macropore size distribution and tortuosity.

Lateral water flow from macropores to micropores is described as a first-order approximation to the diffusion equation:

[3]

where d (m) is an effective diffusion pathlength related to aggregate size, D_w (m² s^–1) is an effective water diffusivity, and _w (unitless) is a scaling factor introduced to match the approximate and exact solutions to the diffusion problem (Gerke and van Genuchten, 1993). Water flow can occur in the reverse direction if the micropores are saturated. Here any excess water is instantaneously transferred to the macropores. Water uptake by roots can take place from both flow domains, but is preferentially extracted from the macropores.

Two different kinds of drainage systems are considered: (i) a primary drainage system located in the soil profile, and (ii) a secondary drainage system surrounding the field. In both cases, flux rates from saturated layers above the drainage depth are predicted using seepage potential theory for layered soils (Leeds-Harrison et al., 1986).

Solute transport in the micropores is calculated using the convection–dispersion equation with source–sink terms (kg m^–3 s^–1) representing mass exchange between flow domains, U_e, crop uptake, U_c, degradation, U_d, losses to field drains, U_s, and losses due to regional ground water flow, U_g:

[4]

where c_mi (kg m^–3) is the solute concentration in the liquid phase, s (kg m^–3) is the sorbed concentration in the solid phase, f (unitless) is the mass fraction of the solid material in contact with water in the macropore domain, (kg m^–3) is the soil bulk density, _m (m³ m^–3) is the micropore water content, accounting for an inaccessible soil volume due to anion exclusion, q (m s^–1) is the water flow rate, and D (m² s^–1) is the dispersion coefficient, calculated as the sum of an effective diffusion coefficient and a dispersion term. Solute transport in the macropores is assumed to be dominated by convection.

The mass transfer term, U_e, accounts for both diffusion and convective flow:

[5]

where D_e (m² s^–1) is an effective diffusion coefficient, c_ma (kg m^–3) is the solute concentration in the liquid phase in macropores, and c' (kg m^–3) indicates either the solute concentration in macropores or in "accessible water" in the micropores, depending on the direction of water flow, S_w. The solute concentration in the water routed into the macropores at the soil surface is calculated assuming instantaneous equilibrium in a thin surface layer or mixing depth, z_d (m).

Pesticide degradation, U_d, follows first-order kinetics and is in this study assumed to proceed at the same rate in both liquid and solid phases in both flow domains. The degradation rate coefficient, µ (s^–1), is adjusted for soil temperature by a modified form of the Arrhenius equation (Boesten and van der Linden, 1991) and soil moisture by a modified form of Walker's function (Walker, 1974).

Equilibrium sorption partitioning is calculated using the Freundlich isotherm. Although MACRO 5.0 allows for kinetic sorption, an instantaneous equilibrium between the liquid phase and the sorbed phase was assumed.

Solute loss to field drainage systems, U_s, is calculated assuming complete lateral mixing of solutes within a flow domain for each soil layer. Solute lost in lateral shallow ground water flow, U_g, is calculated for each saturated soil layer using a retention time concept (Larsbo and Jarvis, 2003).

Initial Conditions, Boundary Conditions, and Driving Data
The bottom boundary condition was chosen to allow simulation of a fluctuating water table in the soil profile. If the base of the profile is saturated, a no-flow condition is applied, which allows the water table to rise. When the soil dries out and the bottom layer in the profile becomes unsaturated, a zero potential condition is applied, which causes water to flow upward into the profile. Initial water contents were taken from measurements on replicated core samples (Table 2) at the start of the simulation period. Initial concentrations of both bromide and bentazone were set to zero based on analyses of these samples.

Hourly data of precipitation corrected for wind drift (mm) and daily mean short wave radiation (W m^–2), relative humidity (unitless), wind speed (m s^–1), and air temperature (°C) were used as driving data for the model. The accumulated precipitation for the simulation period 12 Oct. 1994 to 28 Nov. 1995 was 767 mm.

Parameterization
Six key parameters in the model were selected for calibration using the SUFI global search algorithm and for the GLUE analysis (Table 3). The Lanna soil is a heavy clay in which considerable macropore flow should be expected (Larsson and Jarvis, 1999). The diffusion pathlength, d, and the saturated micropore hydraulic conductivity, K_b, are important parameters determining the strength of macropore flow (Dubus and Brown, 2002). The diffusion pathlength is impossible to measure directly and was, therefore, included in the calibration procedure both for the topsoil (0–30 cm depth), d_top, and the subsoil (30–175 cm), d_sub. The saturated micropore hydraulic conductivity was calibrated for all depths except the top 1 cm where measurements were available (Jarvis and Messing, 1995). Bentazone degradation in Lanna soil followed first-order kinetics and rate coefficients have been derived in laboratory incubations (Bergström et al., 1994). However, measurements of degradation and sorption in laboratory batch experiments are often difficult to extrapolate to field conditions because sample treatment and storage may affect soil biological processes and physical conditions (Beulke et al., 2000; Boesten, 2000). Moreover, sensitivity analyses for the MACRO model show that the degradation rate and sorption coefficients are among the most sensitive parameters for pesticide leaching (Dubus and Brown, 2002). Therefore, both parameters were included in the calibration. The degradation rate coefficients for the topsoil, µ_top, and the subsoil, µ_sub, were calibrated separately. The sorption constant was calculated as the product of the organic carbon fraction in the different soil horizons (Table 1) and the soil organic carbon partition coefficient, K_oc (cm³ g^–1), which was treated as a calibration parameter.

View this table: [in this window] [in a new window]   Table 5. Hydraulic parameters treated as constants for different soil horizons.

[6]

where w_i is the weight given to each data set, m is the number of data sets, n is the number of observations in each group, O_ij and P_ij are the observed and simulated values, and _i is the average of the observations for each group. The weights are constrained by:

[7]

The same two equations with m equal to one were used to calculate EF values for individual groups of observations. If all observed and simulated values are identical, EF will be equal to one, while a negative value indicates a poor fit, meaning that the average value of the observations would be a better estimator than the model simulations. The strata associated with the best values of the goal function after each iteration define the new reduced uncertainty domains. A critical tolerance, T_crit, determines the strata that are removed between iterations and a stopping rule determines when the iterative procedure is stopped. The uncertainty domains when the stopping rule is violated define the posterior uncertainty domains.

We assigned equal weights to each group of observations (Table 2). The term T_crit was defined in absolute terms as the maximum EF_tot value, EF_tot,max, minus 0.2. All strata with zero "hits" at T_crit were removed before the subsequent iteration. Following Roulier and Jarvis (2003), we decided to first calibrate parameters that influence soil hydrology and tracer transport (K_b, d_top, and d_sub) simultaneously against the observed water contents, drainflow, and resident and flux concentrations of bromide. In a second step, parameters controlling sorption and degradation of bentazone (K_oc, µ_top, and µ_sub) were calibrated against observations of resident and flux concentrations of bentazone, retaining the parameter values obtained in the first calibration step. All parameter uncertainty domains were divided into six strata in both calibration steps. The initial uncertainty domain specified for each calibrated parameter can be found in Table 3.

We transferred the uncertainties in parameter values to uncertainties in estimations by running simulations with parameter values sampled randomly from the posterior uncertainty domains. Uncertainty ranges defined by the difference between the 95th and 5th percentiles of the cumulative distributions of the simulated output variables were compared with the observations. If some of the observations fall outside these uncertainty ranges, it means that the estimations do not "honor" all the observations (Abbaspour et al., 1999). However, in many cases this may not be meaningful since both the observations and the driving data for the simulations may be subject to error. A compromise between honoring the data and achieving reasonably low uncertainty in predictions is often desirable, which introduces an element of subjectivity into the procedure. The SUFI procedure results in parameters that are conditioned on the observations. Conditioned parameters have an advantage over fitted parameters whenever Monte Carlo type simulations are to be run (Abbaspour et al., 1999).

Generalized Likelihood Uncertainty Estimation
The GLUE procedure (Beven and Binley, 1992) was used in this study as an alternative to the SUFI/Monte Carlo approach for uncertainty estimation with the MACRO 5.0 model. The GLUE methodology explicitly recognizes the underlying limitations of environmental models by accepting that all models and measurements are to some extent in error. Consequently, we cannot expect to find one optimal unique parameter set for a specific model application by some calibration procedure. There may be many combinations of parameters that, according to some measure of goodness-of-fit, a goal function, or likelihood function, represent the observations equally well. This is referred to as "equifinality" of parameter sets (Beven, 1993). Equifinality does not mean that the estimations are identical, only that they fit the calibration data equally well. In contrast to most calibration procedures, the objective of GLUE is not to identify optimal values or uncertainty domains of specific parameters. The GLUE procedure is only concerned with evaluating the "likelihood" of combinations of parameters as simulators of the observations. Likelihood is here used in a broad sense, meaning a specified measure of how well the outcome of a model and parameter set describes the observations. Not all parameter sets will be acceptable simulators of the observations. The values of the likelihood function can subsequently be used as weights in predictive simulations. A threshold value defining acceptable, or "behavioral" (Spear and Hornberger, 1980), parameter sets must then be chosen. All "nonbehavioral" parameter sets are discarded by assigning them zero weight. The GLUE procedure does not provide any information on interactions between parameters, but these are implicitly reflected in the likelihood values.

The basic requirements of the GLUE procedure are that distributions of parameters and maximum and minimum values are specified for all parameters considered in the analysis and that some measure of goodness-of-fit, a likelihood function, is defined. In addition, a procedure for using likelihood weights in the estimation of uncertainty in model predictions must be defined when GLUE is used predictively. To enable a comparison with SUFI, we used EF_tot as the likelihood function and maximum and minimum values corresponding to the initial uncertainty domains used in the SUFI calibrations (Table 3). A uniform distribution within the domain was chosen for all parameters, in accordance with the procedure outlined by Beven and Binley (1992). Since MACRO 5.0 cannot simulate two solutes at the same time, two rounds of 30000 Monte Carlo simulations with the same parameter sets were run, one for bromide and one for bentazone. The parameter sets were generated using latin hypercube sampling from the uncertainty domains. The results from each simulation were compared with the observations using Eq. [6], with each group assigned equal weight.

The groups of data were combined in four different ways to evaluate with GLUE the significance of data availability for conditioning of parameters related to macropore flow and pesticide transport: (i) all observations (All), (ii) only soil water contents and resident concentrations in the soil profile (Res), (iii) only drainflow and flux concentrations of bromide and bentazone (Flux), (iv) soil water contents, drainflow, and both resident and flux concentrations of bentazone (NoTracer).

In this study, all parameter sets with EF_tot values within 0.2 of EF_tot,max were considered behavioral. This definition of the threshold for behavioral simulations was chosen because it allows comparisons of parameter conditioning between different data groups with different EF_tot,max values. In the same way as for SUFI, 5th and 95th percentiles of model outputs for the behavioral simulations were compared with the observations. Clearly, these uncertainty limits depend on the choice of goal function and on the threshold value defining behavioral parameter sets.

Even though the GLUE procedure focuses on parameter sets, information on an individual parameter can be obtained from its cumulative likelihood distribution. Those parameters showing more deviation from the initial uniform distribution have been more conditioned by the procedure. A steep slope indicates that the measure of likelihood (EF) is sensitive to that parameter. However, a visual inspection of sensitivity will depend on the parameter range, which makes comparisons between parameters difficult. This type of parameter conditioning differs from the SUFI approach in that parameter values only have validity within their parameter set. However, if desired, uncertainty ranges for individual parameters can be identified from percentiles of the cumulative likelihood distributions. The cumulative likelihood distributions were calculated using only the behavioral simulations with EF_tot values rescaled so that they varied between zero and 0.2. The degree of parameter conditioning is influenced by the criterion for behavioral parameter sets. Since the number of behavioral parameter sets decreases as the criterion gets stricter, the range of behavioral parameter values is also likely to decrease and the degree of parameter conditioning increase. This effect was investigated by reducing the threshold for behavioral parameter sets to EF_tot,max minus 0.4.

Predictions
The experiment at Lanna did not reflect normal agricultural practice since bentazone was applied on bare soil in autumn at nearly twice the recommended dose. In Sweden, bentazone is commonly used in spring to control weeds in field peas (Pisum sativum L.). The knowledge of parameter uncertainty gained by applying SUFI and GLUE to the experimental data was, therefore, translated into predictions of environmental risks by re-parameterizing the model to represent good agricultural practice in southwest Sweden. We simulated an application of the maximum recommended dose (1.305 kg ha^–1) of bentazone (Swedish Board of Agriculture, personal communication, 2004) on 2 May. Crop parameters for peas were taken from Forum for the Coordination of Pesticide Fate Models and Their Use (2001).

The SUFI procedure does not provide any information about the distributions of parameter values within the posterior uncertainty domains. Hence, it seems reasonable to transfer the distributions from the initial uncertainty domains onto the posterior uncertainty domains. Using latin hypercube sampling from the posterior uncertainty domains, 235 parameter combinations were generated for predictive simulations. Since we used uniform distributions throughout, all predictions received the same weight. A distribution of any simulated output variable can then be constructed from which standard statistical measures can be calculated.

The GLUE procedure, on the other hand, implicitly provides information on the parameter distributions from the likelihood values assigned to each of the behavioral parameter sets. These were rescaled as before to remove the effect of different absolute values of EF_tot for different groups of outputs and then used as weights for the predictions. Again, distributions of the predictions were calculated.

We tested the effects of the availability of different groups of observations on predictions using the same data groups as in the analysis of parameter conditioning. To limit the number of predictive simulations, 235 parameter sets were randomly sampled from the behavioral parameter sets when necessary. All parameter sets were used when the number of behavioral simulations was less than 235.

Two target output variables were considered to reflect acute and chronic toxicity: the maximum concentration of bentazone in drainflow at an hourly time resolution and the accumulated loss of bentazone through drainage during two years. The same driving data as in the calibration simulations were used.

	RESULTS AND DISCUSSION

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 
Overall Model Performance
The EF values calculated using different groups of observations are presented in Table 6 for the optimal SUFI parameter values, the GLUE simulation with the largest EF_tot, and the GLUE simulation with the largest EF value for each data group. The model simulated the field observations fairly well (EF_tot = 0.308 for the optimal SUFI simulation and EF_tot = 0.241 for the best GLUE simulation). For the sake of brevity, we only show visual comparisons of the measurements and GLUE simulations for flux (Fig. 2 and 3) and resident concentrations (Fig. 4 and 5) .

View larger version (28K): [in this window] [in a new window]   Fig. 2. Measured flux concentrations of bromide (triangles) compared with the 5th and 95th percentiles of the generalized likelihood uncertainty estimation (GLUE) behavioral simulations obtained using all observations (gray lines) (n = 235) and the GLUE simulations with the largest overall model efficiency (black line). The percentile curves include values from different simulations.

View larger version (25K): [in this window] [in a new window]   Fig. 3. Measured flux concentrations of bentazone (triangles) compared with the 5th and 95th percentiles of the generalized likelihood uncertainty estimation (GLUE) behavioral simulations obtained using all observations (gray lines) (n = 235) and the GLUE simulations with the largest overall model efficiency (black line). The percentile curves include values from different simulations.

View larger version (29K): [in this window] [in a new window]   Fig. 4. Resident concentration profiles of bromide on five sampling occasions. Measured values (triangles) compared with the 5th and 95th percentiles of the generalized likelihood uncertainty estimation (GLUE) behavioral simulations obtained using all observations (gray lines) (n = 235) and the GLUE simulation with the largest overall model efficiency (crosses). Note the different scales on the x axes.

View larger version (18K): [in this window] [in a new window]   Fig. 5. Resident concentration profiles of bentazone on three sampling occasions. Measured values (triangles) compared with the 5th and 95th percentiles of the generalized likelihood uncertainty estimation (GLUE) behavioral simulations obtained using all observations (gray lines) (n = 235) and the GLUE simulation with the largest overall model efficiency (crosses). Note the different scales on the x axes.

The EF_tot was larger for the optimal SUFI simulation than the best GLUE simulation, which indicates that SUFI is an efficient parameter estimation tool and that the number of GLUE simulations was too small to fully cover the parameter space. The largest EF value for each data group was in most cases much larger than the corresponding EF values both for the optimal SUFI parameter set and the best GLUE simulation, demonstrating that the model cannot simultaneously match all observations. If optimization is performed on a subset of data it will be at the expense of the data not included. This reflects errors in the measurements, in the model process descriptions, and in the parameterization.

Sequential Uncertainty Fitting Results
The posterior uncertainty domains from the SUFI calibration are presented in Table 3 together with the optimal values. The optimal value for K_b is of the same order of magnitude as the measured value at the soil surface (0.13 mm h^–1). A smaller value of K_b in the subsoil could be expected due to the increasing clay content with depth (Table 1). On the other hand, crust formation is likely to decrease K_b at the soil surface in clay soils (Messing and Jarvis, 1993). The optimal value obtained for K_oc is consistent with the weak or negligible sorption of bentazone in clay soils reported by Gaston et al. (1996) and Scorza Junior (2002). The optimal value of the degradation rate coefficient of bentazone in the Lanna topsoil is similar to the values determined in the laboratory incubation experiments (0.052 and 0.059 d^–1 in two samples; Bergström et al., 1994). Degradation could not be detected in the subsoil samples in the incubation experiments (Bergström et al., 1994), whereas SUFI suggests a half-life of 33 d based on the field data. This confirms the difficulties that are sometimes encountered when extrapolating laboratory measurements to the field (Beulke et al., 2000, Boesten, 2000). The optimal values for d_top and d_sub are typical for structured soils exhibiting strong macropore flow (Kätterer et al., 2001; Roulier and Jarvis, 2003). The fact that the optimized parameters seem physically sound and similar to measured data where available gives us confidence both in the model and in the parameter identification procedure.

All initial parameter uncertainty domains were reduced in the SUFI procedure, except for d_sub, indicating that the initial uncertainty domains were large enough. On the other hand, SUFI, as applied here, decreased the uncertainty domains a little too much for some parameters (e.g., d_top, µ_top) since some GLUE parameter sets with parameter values lying outside the SUFI posterior uncertainty domains had EF_tot values larger than the threshold value defining behavioral parameter sets (Fig. 6) . Moreover, it seems reasonable that the part of the uncertainty domain including the smallest values of d_sub should have been removed since no GLUE parameter sets sampled from this region had EF_tot values close to the threshold (Fig. 6). These differences may have arisen partly because SUFI was applied as a two-step procedure. A denser sampling scheme in SUFI would probably increase the chances of finding the correct uncertainty domains. As the SUFI sampling scheme gets denser, the difference between stratified sampling and Monte Carlo sampling vanishes. Since the GLUE simulations were sampled from uniform distributions, a selection of GLUE simulations with parameter values that fall inside all SUFI posterior uncertainty domains (Fig. 6) can be used as a representative sample of simulations from the SUFI posterior uncertainty domains. Of the 30000 GLUE simulations, 148 satisfy this criterion. Many of these parameter sets are clearly poor simulators of the observations (Fig. 6). The minimum EF_tot value for this subset of simulations was –3.2. There are two main reasons for this. First, as noted above, the posterior uncertainty domains may not always have been determined correctly, which allowed sampling from parts of the parameter space characterized by only small EF_tot values. Second, the parameters are more or less correlated. This can be seen in Table 7, which shows the correlation coefficients for the parameter values used in the behavioral GLUE simulations. There are probably additional, more complex correlations between more than two parameters at a time that cannot easily be captured by a correlation matrix. Neglecting correlations between parameters can result in substantial bias in model estimations (Smith et al., 1992).

View larger version (75K): [in this window] [in a new window]   Fig. 6. Model efficiencies (EFtot) for all generalized likelihood uncertainty estimation (GLUE) simulations (gray) and for the GLUE simulations with parameter values that lie within all sequential uncertainty fitting (SUFI) posterior uncertainty domains (black). The EFtot values were calculated using all observations. The horizontal lines indicate the threshold for behavioral simulations. The term Kb is the saturated micropore hydraulic conductivity, Koc is the soil organic carbon partition coefficient, dtop and dsub are the diffusion pathlengths in the topsoil and subsoil, respectively, and µtop and µsub are the degradation rate coefficients in the topsoil and subsoil, respectively.

View this table: [in this window] [in a new window]   Table 7. Correlation coefficients for overall model efficiencies for different parameters. Only the behavioral simulations using all observations are considered.

The 5th and 95th percentiles for bromide and bentazone concentration in drainflow are shown together with the GLUE simulation with the largest EF_tot and the observations in Fig. 2 and 3. Many of the field observations, especially the peak flux concentrations, fall outside the estimated uncertainty range. Figures 4 and 5 show depth profiles of resident bromide and bentazone, respectively. Again, some of the observations fall outside the estimated uncertainty ranges. Although not shown here, the uncertainty ranges for SUFI are of the same magnitude as those for GLUE both for flux and resident concentrations of bromide and bentazone.

Cumulative likelihood distributions for each parameter conditioned on different combinations of data are shown in Fig. 7 together with the uniform initial distributions. Each group of observations results in a different number of behavioral simulations ranging from 67 when only flux measurements are used to 2668 when only storage measurements are used. All parameters evaluated by All were highly conditioned by the evaluation procedure except for d_top, which was poorly conditioned, regardless of which groups of observations were used. All parameters except µ_top were poorly conditioned by Res compared with All. Parameter uncertainty ranges defined as the difference between the 5th and 95th percentiles of the cumulative distributions were on average 45% larger for Res. This should perhaps be expected for K_b and d considering the lack of sensitivity of macropore flow parameters to soil water contents and resident concentrations (Jarvis, 1999). However, the failure to condition bentazone sorption and the subsoil degradation using Res is also striking and perhaps more surprising. This is probably a result of the infrequent sampling in time of Res (Table 2) and the large uncertainty in these observations due to the spatial variability in the field (Larsson and Jarvis, 1999). All parameters, except d_top and K_oc, were highly conditioned by Flux. The parameter uncertainty ranges were on average 24% larger compared with conditioning by All. However, it should be noted that the cumulative likelihood distributions for K_b, d_sub, and µ_top differ considerably for conditioning by Flux and All. The reason why K_oc was poorly conditioned by Flux whereas µ_top and µ_sub were highly conditioned is not clear to us. The parameter uncertainty ranges for NoTracer conditioning were on average only 22% larger compared with All. A marked difference between the conditioning of K_b and d_sub by NoTracer and All shows that the bromide observations contain information on solute transport that is not included in the bentazone observations. This implies either that some process or mechanism that influences bromide and bentazone transport differently is not properly accounted for in the model or that correlations exist between the transport parameters and sorption and degradation parameters. This latter explanation is supported by the difference between the NoTracer and All conditioning of the pesticide properties. It is well known that sorption and macropore flow can have similar effects on the movement of pesticides that are difficult to distinguish without tracer data (Gaber et al., 1995; Jarvis et al., 1995). In our case, larger values of K_b for NoTracer compared with All generate more water flow in the micropores, the effects of which may be compensated by larger values of K_oc. This hypothesis is supported by a correlation coefficient of 0.34 between K_b and K_oc for the NoTracer behavioral simulations.

View larger version (37K): [in this window] [in a new window]   Fig. 7. Cumulative likelihood distributions for parameter conditioning with generalized likelihood uncertainty estimation (GLUE) using different combinations of observations. The number of behavioral simulations included in each combination of groups is denoted by n. The term Kb is the saturated micropore hydraulic conductivity, Koc is the soil organic carbon partition coefficient, dtop and dsub are the diffusion pathlengths in the topsoil and subsoil, respectively, and µtop and µsub are the degradation rate coefficients in the topsoil and subsoil, respectively.

Scenario Predictions for Good Agricultural Practice
Figure 8 shows the cumulative distributions of the predicted maximum bentazone concentrations and the accumulated losses of bentazone through drainflow for the SUFI and GLUE approaches. The 5th, 50th, and 95th percentiles are presented in Table 8. Using All, the uncertainty ranges are 35 and 21% smaller for GLUE predictions compared with SUFI for the maximum concentrations and accumulated losses, respectively. The results of the GLUE predictions using different combinations of observations clearly show that All and NoTracer give the smallest uncertainty ranges. The uncertainty ranges for Flux are large compared with All and NoTracer even though the degree of parameter conditioning was almost as high for Flux as for NoTracer. This is probably an effect of the bias in the cumulative likelihood distributions of K_b, d_sub, and µ_top (Fig. 7) compared with All.

View larger version (32K): [in this window] [in a new window]   Fig. 8. Comparison of predictions of maximum concentrations of bentazone in drainflow and accumulated losses of bentazone through drainflow for sequential uncertainty fitting (SUFI) and generalized likelihood uncertainty estimation (GLUE) using different combinations of observations. The number of predictive simulations is denoted by n.

	CONCLUSIONS

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

The GLUE analysis showed that observations of soil water contents, drainflow, and both flux and resident concentrations for both tracer and pesticide gave the highest degree of parameter conditioning. However, not all selected parameters were highly conditioned by the observations even with such a comprehensive data set as Lanna. The diffusion pathlength in the topsoil was especially difficult to identify with GLUE. The parameters were generally poorly conditioned when only data on soil water contents and resident concentrations were available. This was attributed to the lack of sensitivity of macropore flow parameters to these variables, infrequent sampling, and large spatial variability in the observations. Data on only drainflow and flux concentrations gave highly conditioned parameters except for K_oc and the diffusion pathlength in the topsoil, but apparently biased estimates for others (e.g., degradation rate coefficient in the topsoil). Marked differences in conditioning of transport, sorption, and degradation parameters were found when bromide data were excluded. This shows that the bromide observations contain information on solute transport that is not included in the bentazone observations and that parameter correlations can result in incorrect estimates when tracer data are lacking.

The largest EF value for each of the data groups was in most cases much larger than the corresponding EF value obtained from the best overall GLUE simulation, reflecting errors in the measurements, in the model process descriptions, and in the parameterization. Without comprehensive data sets there would be no way of knowing that such problems existed, and little possibility to reduce predictive uncertainty and error by improving the model and its parameterization. Since different data groups are simulated with a varying degree of accuracy, it might seem reasonable to put more weight in the calibration process on the data to be predicted (e.g., flux concentrations). However, this approach is risky since it is impossible to foresee the effects on extrapolations of errors in the model and in the parameterization.

The SUFI procedure proved to be an efficient parameter estimation method, reducing all initial uncertainty domains except for the diffusion pathlength in the subsoil. Even though SUFI decreased the initial parameter uncertainty domains significantly, many parameter combinations sampled from within these domains resulted in poor simulations. This is because SUFI did not determine the posterior uncertainty domains correctly for all parameters and does not account for correlations between parameters. It does not seem reasonable to make predictions using parameter sets that we know are poor at simulating the system. Indeed, compared with GLUE, random sampling from the SUFI posterior uncertainty domains resulted in 54 and 17% larger prediction uncertainty ranges for maximum bentazone concentrations in drainflow and total losses of bentazone through drainflow, respectively, for a typical scenario assuming good agricultural practice.

	ACKNOWLEDGMENTS

 
This work was funded by the EU 5th framework project APECOP ("Effective approaches for assessing the predicted environmental concentrations of pesticides. A proposal supporting the harmonized registration of pesticides in Europe", QLK4-CT1999-01238). The authors thank Dr. Martin Larsson for making the field data available to us and Fredrik Stenemo for supplying a preliminary parameterization for Lanna (both at the Department of Soil Sciences, SLU, Uppsala). This paper has benefited from discussions with Dr. Stephanie Roulier and Anna Lindahl at the Department of Soil Sciences, SLU, Uppsala, and with Klas Hansson at the Department of Hydrology, Uppsala University.

	REFERENCES

TOP ABSTRACT INTRODUCTION MATERIALS AND METHODS RESULTS AND DISCUSSION CONCLUSIONS REFERENCES

 

• Abbaspour, K.C., M.A. Sonnleitner, and R. Schulin. 1999. Uncertainty in estimation of soil hydraulic parameters by inverse modeling: Example lysimeter experiments. Soil Sci. Soc. Am. J. 63:501–509.[Abstract/Free Full Text]
• Abbaspour, K.C., M.T. van Genuchten, R. Schulin, and E. Schläppi. 1997. A sequential uncertainty domain inverse procedure for estimating subsurface flow and transport parameters. Water Resour. Res. 33:1879–1892.
• Armstrong, A.C., A.M. Portwood, P.B. Leeds-Harrison, G.L. Harris, and J.A. Catt. 1996. The validation of pesticide leaching models. Pestic. Sci. 48:47–55.
• Bergström, L., N.J. Jarvis, and J. Stenström. 1994. Pesticide leaching data to validate simulation models for registration purposes. J. Environ. Sci. Health. Part A Environ. Sci. Eng. A29:1073–1104.
• Beulke, S., I.G. Dubus, C.D. Brown, and B. Gottesburen. 2000. Simulation of pesticide persistence in the field on basis of laboratory data—A review. J. Environ. Qual. 29:1371–1379.[Abstract/Free Full Text]
• Beven, K.J. 1993. Prophecy, reality and uncertainty in distributed hydrological modelling. Adv. Water Resour. 16:41–51.
• Beven, K., and A. Binley. 1992. The future of distributed models: Model calibration and uncertainty prediction. Hydrol. Processes 6:279–298.
• Boesten, J.J.T.I. 2000. From laboratory to field: Uses and limitations of pesticide behaviour models for the soil/plant system. Weed Res. 40:123–138.
• Boesten, J.J.T.I., and A.M.A. van der Linden. 1991. Modeling the influence of sorption and transformation on pesticide leaching and persistence. J. Environ. Qual. 20:425–435.[Abstract/Free Full Text]
• Dubus, I.G., S. Beulke, and C.D. Brown. 2002. Calibration of pesticide leaching models: Critical review and guidance for reporting. Pest Manage. Sci. 58:745–758.
• Dubus, I.G., and C.D. Brown. 2002. Sensitivity and first-step analysis for the preferential flow model MACRO. J. Environ. Qual. 31:227–240.[Abstract/Free Full Text]
• Dubus, I.G., C.D. Brown, and S. Beulke. 2003. Sources of uncertainty in pesticide fate modelling. Sci. Total Environ. 317:53–72.[Medline]
• Feyen, J., D. Jacques, A. Timmermann, and J. Vanderborght. 1998. Modelling water flow and solute transport in heterogeneous soils: A review of recent approaches. J. Agric. Eng. Res. 70:231–256.
• Flury, M. 1996. Experimental evidence of transport of pesticides through field soils—A review. J. Environ. Qual. 25:25–45.
• Forum for the Coordination of Pesticide Fate Models and Their Use. 1995. The final report of the work of the regulatory modelling work group of FOCUS. EU Doc. 4952/VI/95. EU Commision, Brussels.
• Forum for the Coordination of Pesticide Fate Models and Their Use. 2000. FOCUS groundwater scenarios in the EU plant protection review process. Report of the FOCUS Groundwater Scenarios Workgroup. EC Document Reference Sanco/321/2000. EU Commision, Brussels.
• Forum for the Coordination of Pesticide Fate Models and Their Use. 2001. FOCUS surface water scenarios in the EU evaluation process under 91/414/EEC. Report of the FOCUS Working Group on Surface Water Scenarios. EC Document Reference Sanco/4802/2001-rev.1.221. EU Commision, Brussels.
• Gaber, H.M., W.P. Inskeep, S.D. Comfort, and J.M. Wraith. 1995. Nonequilibrium transport of atrazine through large intact soil cores. Soil Sci. Soc. Am. J. 59:60–67.[Abstract/Free Full Text]
• Gaston, L.A., M.A. Locke, S.C. Wagner, R.M. Zablotowics, and K.N. Reddy. 1996. Sorption of bentazon and degradation products in two Mississippi soils. Weed Sci. 44:678–682.
• Gerke, H.H., and M.Th. van Genuchten. 1993. Evaluation of a first-order water transfer term for variably saturated dual-porosity flow models. Water Resour. Res. 29:1225–1238.
• Germann, P.F., and K. Beven. 1985. Kinematic wave approximation to infiltration into soils with sorbing macropores. Water Resour. Res. 21:990–996.
• Jarvis, N.J. 1999. Using preferential flow models for management purposes. p. 521–535. In J. Feyen and K. Wiyo (ed.) Modelling of transport processes in soils at various scales in time and space. Wageningen Pers, Wageningen, the Netherlands.
• Jarvis, N.J. 2001. Modelling pesticide environmental fate: Process understanding and knowledge gaps. p. 11–18. In A. Walker (ed.) Pesticide behaviour in soils and water. BCPC Symp. Proc. 78. The British Crop Protection Council, Farnham, Surrey, UK.
• Jarvis, N.J. 2002. Macropore and preferential flow. p. 1005–1013. In J. Plimmer (ed.) The encyclopedia of agrochemicals. Vol. 3. John Wiley & Sons, New York.
• Jarvis, N.J., M. Larsson, P. Fogg, and A.D. Carter. 1995. Validation of the dual-porosity model ‘MACRO’ for assessing pesticide fate and mobility in soils. p. 161–170. In A. Walker et al. (ed.) Pesticide movement to water. BCPC Monogr. 62. The British Crop Protection Council, Farnham, Surrey, UK.
• Jarvis, N.J., and I. Messing. 1995. Near-saturated hydraulic conductivity in soils of contrasting texture measured by tension infiltrometers. Soil Sci. Soc. Am. J. 59:27–34.[Abstract/Free Full Text]
• Kätterer, T., B. Schmied, K.C. Abbaspour, and R. Schulin. 2001. Single- and dual-porosity modelling of multiple tracer transport through soil columns: Effects of initial moisture and mode of application. Eur. J. Soil Sci. 52:25–36.
• Larsbo, M., and N. Jarvis. 2003. MACRO 5.0. A model of water flow and solute transport in macroporous soil. Technical description. Emergo 2003:6. Swedish Univ. of Agric. Sci., Dep. of Soil Sci., Uppsala, Sweden.
• Larsson, M.H., and N.J. Jarvis. 1999. Evaluation of a dual-porosity model to predict field-scale solute transport in a macroporous soil. J. Hydrol. (Amsterdam) 215:153–171.
• Leeds-Harrison, P.B., C.J.P. Shipway, N.J. Jarvis, and E.G. Youngs. 1986. The influence of soil macroporosity on water retention, transmission and drainage in a clay soil. Soil Use Manage. 2:47–50.
• Loague, K., and R.E. Green. 1991. Statistical and graphical methods for evaluating solute transport models: Overview and application. J. Contam. Hydrol. 7:51–73.
• Messing, I., and N. Jarvis. 1993. Temporal variation in the hydraulic conductivity of a tilled clay soil as measured by tension infiltrometers. J. Soil Sci. 44:11–24.
• Mualem, Y. 1976. A new model for predicting the hydraulic conductivity of unsaturated porous media. Water Resour. Res. 12:513–522.
• Roulier, S., and N. Jarvis. 2003. Modeling macropore flow effects on pesticide leaching: Inverse parameter estimation using microlysimeters. J. Environ. Qual. 32:2341–2353.[Abstract/Free Full Text]
• Schwartz, R.C., A.S.R. Juo, and K.J. McInnes. 2000. Estimating parameters for a dual-porosity model to describe non-equilibrium, reactive transport in a fine-textured soil. J. Hydrol. (Amsterdam) 229:149–167.
• Scorza Junior, R.C. 2002. Pesticide leaching in macroporous clay soils: Field experiment and modeling, Ph.D. diss. Wageningen Univ., Wageningen, the Netherlands.
• imnek, J., N.J. Jarvis, M. Th. van Genuchten, and A. Gärdenäs. 2003. Review and comparison of models for describing non-equilibrium and preferential flow and transport in the vadose zone. J. Hydrol. (Amsterdam) 272:14–35.
• Smith, A.E., P.B. Ryan, and J.S. Evans. 1992. The effect of neglecting correlations when propagating uncertainty and estimating the population distribution of risk. Risk Anal. 12:467–474.[Medline]
• Spear, R.C., and G.M. Hornberger. 1980. Eutrophication in peel inlet: II. Identification of critical uncertainties via generalized sensitivity analysis. Water Res. 14:43–49.
• Van Genuchten, M.Th. 1980. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 44:892–898.[Abstract/Free Full Text]
• Walker, A. 1974. A simulation model for prediction of herbicide persistence. J. Environ. Qual. 3:396–401.[Abstract/Free Full Text]
• Wösten, J.H.M., A. Lilly, A. Nemes, and C. Le Bas. 1998. Using existing soil data to derive hydraulic parameters for simulation models in environmental studies and in land use planning. Rep. 156. DLO-Staring Centre, Wageningen, the Netherlands.
• Yapo, P.O., H.V. Gupta, and S. Sorooshian. 1998. Multi-objective global optimization for hydrologic models. J. Hydrol. (Amsterdam) 204:83–97.
• Yates, S.R., M.Th. van Genuchten, A.W. Warrick, and F.J. Leij. 1992. Analysis of measured, predicted, and estimated hydraulic conductivity using the RETC computer program. Soil Sci. Soc. Am. J. 56:347–354.[Abstract/Free Full Text]

Related articles in JEQ:

This Issue in Journal of Environmental Quality

JEQ 2005 34: 403-407. [Full Text]  

This article has been cited by other articles:

				A. Boivin, J. Simunek, M. Schiavon, and M. Th. van Genuchten Comparison of Pesticide Transport Processes in Three Tile-Drained Field Soils Using HYDRUS-2D Vadose Zone J., June 21, 2006; 5(3): 838 - 849. [Abstract] [Full Text] [PDF]

This Article Abstract