Published online 8 March 2006
Published in Vadose Zone J 5:222-233 (2006)
DOI: 10.2136/vzj2005.0016
© 2006 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
SPECIAL SECTION: FROM FIELD- TO LANDSCAPE-SCALE VADOSE ZONE PROCESSES
On the Value of Local Measurements for Prediction of Pesticide Transport at the Field Scale
Keith Bevena,*,
Danrong Zhangb and
André Mermoudc
a Environmental Science/Lancaster Environment Centre, Lancaster Univ., Lancaster LA1 4YQ, UK
b Key Lab. of Water Resources Development, Hohai Univ., Nanjing, 210098, China
c Institute of Environmental Science and Technology, Swiss Federal Institute of Technology, ISTE/HYDRAM, ENAC, EPFL, 1015 Lausanne, Switzerland
* Corresponding author (k.beven{at}lancaster.ac.uk)
Received 31 January 2005.
 |
ABSTRACT
|
|---|
Pesticide transport through the soil profile at the field scale is notoriously difficult to predict because of a lack of appropriate field-scale models and model parameters that take adequate account of heterogeneity in local flow rates and parameters. It is generally impossible to measure all the parameters required to describe that heterogeneity. A way of approximating pesticide transport at the field scale, given uncertainty in the representation of both local characteristics and field-scale distributions, is presented. The methodology allows prior estimates of field-scale distribution parameters to be conditioned on uncertain column-scale measurements. It is shown how this conditioning can drastically constrain steady-state flow predictions of atrazine (2-chloro-4-ethylamino-6-isopropylamino-1,3,5-triazine) transport to groundwater for a site in the Rhone Valley, Switzerland. Some of the simplifying assumptions of the current analysis, in particular, to allow for transient simulations, can be relaxed in the future as more computer power becomes available. The methodology, however, should remain valid.
Abbreviations: BTC, breakthrough curve • CDE, convection–dispersion equation • CEC, cation exchange capacity • GLUE, Generalized Likelihood Uncertainty Estimation • OC, organic carbon
|
INTRODUCTION
|
|---|
LABORATORY undisturbed column experiments are frequently used to quantify local parameters since they can mimic satisfactorily the factors affecting local transport in the field. It has been shown by many studies that the analytical solution of variants of the convection–dispersion equation (CDE) can effectively describe the solute transport in soil columns (e.g., Comegna et al., 2001; Zurmühl, 1998; Zhang et al., 2005). Practical environment management, however, requires predictions at larger field and catchment scales. Because column experiments are costly and time consuming, only a limited number of experiments can be performed before making larger scale predictions. Even in research plots, it is difficult to obtain sufficient experimental information about the heterogeneity of the transport parameters in the model upon which to base the field-scale predictions. It follows that any predictions at the scales of interest for practical management will be necessarily associated with significant uncertainty. In this respect, it is important to take full advantage of the information provided by the parameter values derived from large, undisturbed column experiments to constrain the predictive uncertainties at the larger scale.
Following Bresler and Dagan (1981, 1983), Amoozegar-Fard et al. (1982), Beven (1993) and others, a field can be considered to be made up of a distribution of parallel vertical columns, with a common lower boundary condition (here a specified depth of the groundwater table). It is assumed that the columns are large enough to minimize the effects of local intercolumn interactions through preferential or other flow pathways, so that they can be treated as independent. Furthermore, the field is considered to be large enough to include several variogram ranges of local variability or that the nugget variance is a very large proportion of the total variance (e.g., van Wesenbeeck and Kachanoski, 1991). Consequently, it is possible to deal directly with distributions of parameters rather than worrying about the spatial covariation.
The methodology followed here is an extension of the original "value of a single measurement" paper of Beven (1993). That analysis showed how field-scale transport could be conditioned by a single column experiment. It also suggested how adding further column measurements might further constrain the uncertainty in field-scale predictions. It was assumed, however, that the parameter values for the single measurement site were known precisely (i.e., by optimization) and were restricted only to the dispersion of a conservative solute. This study extends the methodology of Beven (1993) by: (i) considering reactive substances such as pesticides as well as water flow and conservative transport, (ii) testing the additional conditioning that might be provided by multiple column experiments, and (iii) allowing for uncertainty in the estimation of parameter values for each column experiment. The approach differs from a forward uncertainty analysis based on prior distributions of parameters as a way of representing within-field variability (e.g., Beulke et al., 2004) by adding conditioning on the local field data and allowing for uncertainty in the calibration step. Some of the simplifying assumptions of the current analysis, in particular, that of a quasi-steady average annual flow rate in the profile, can be relaxed in the future as more computer power becomes available. The methodology, however, should remain valid.
|
The Study Site
|
|---|
This study aims to predict pesticide transport at the field scale in a particular field in the Rhône Valley between Martigny and Charrat in Valais, Switzerland (latitude: 46°7'12'' N, longitude: 7°7'35'' E). The field is part of the flat alluvial valley bottom at 450 m above sea level and is under intensive agricultural use. The soils are typical fluvisols (Entisols) of high hydraulic conductivity with relatively low clay and organic matter, and the site is in a zone of high vulnerability to groundwater contamination. The pH of the soil lies within the range 7.8 to 8.5. The organic carbon (OC) content varies from 5.3 to 14.2 mg g–1. Cation exchange capacity (CEC) varies from 4.6 to 11.6 cmole kg–1. Sand, silt and clay contents are within the ranges of 14 to 93, 6 to 72, and 1 to 16%, respectively. Higher OC content, CEC, and clay content appear in deeper samples taken below depths of 70 cm (Zhang, 2003).
Four large undisturbed soil columns, each about 53 cm long and 491 cm2 in cross section, were collected from the field and subjected to long time scale breakthrough curve (BTC) experiments using KBr and atrazine to trace water and pesticide movement (Zhang, 2003). These experiments were performed at a quasi-steady flow rate, at or near saturation. Atrazine is a widely used herbicide to control broadleaf and many other weeds. It has been found in the groundwater of the Rhone flood plain near the sites where the soil samples were collected. The results of fitting transport models to these BTC experiments were reported in Zhang et al. (2005). The results of this analysis suggested that the pesticide transport could be described adequately in these soils by an equilibrium convection-dispersion model with only four parameters (v, D, R, and µ; see Table 1), although some discrepancies appear in describing the tail of the BTCs.
|
Representing Field-Scale Transport
|
|---|
The representative field was assumed to be made up of a collection of vertical columns having the same cross-sectional area as the experimental columns (491 cm2). Thus, a field of 10 ha can be represented by a distribution of 2 036 660 columns. In this approximate analysis it is necessary to assume that interconnections between the columns and spatial correlations in the parameters could be neglected. The latter is not a significant limitation for fields in which the field is large relative to the range of any spatial correlation. Two sources of information about the variability and distributions of the local column parameters can be used: information from the literature and information obtained by fitting observed BTCs for the experimental columns taken from the field. The methodology proposed here aims to combine these sources of information in a consistent way, using the literature estimates as prior distributions and the measurement derived values to condition posterior field distributions, taking account of the uncertainty in both sources of information.
Lacking a large number of measurements in the field, it is also necessary to estimate the form and scale of variation of the distribution of local-scale values for each parameter from information found in the literature. This will allow the range of possible field-scale distributions of any parameter to be defined in a way that is consistent with the chosen maximum and minimum expected values. The most commonly used distributional forms include the uniform, normal, and lognormal distribution, with different ways of defining the scale of variation. Once the form and scale of variation of the distribution are chosen, it was assumed that the field might contain all the potential parameter distributions that lie within the maximum and minimum values for that parameter (Fig. 1
). For normal and lognormal distributions it was also assumed that the distributions have the same CV as the mean varies. Information about the potential SD or CV of the distributions for different parameters can be found in the literature. Having defined a range for a parameter, and an expected field-scale CV for that parameter, the range within which the mean of the field-scale distribution can vary is limited by the requirement that most of the local values (e.g., >95%) of the field-scale distributions should be within the feasible range defined by the maximum and minimum values. The distributions can then be distinguished from one another by their means (
, where x is one of the parameters required).
View larger version (13K):
[in this window]
[in a new window]
|
Fig. 1. An illustration of the field-scale parameter distribution, with changing mean but fixed standard deviation, moving across the feasible parameter range [x0, xn].
|
|
A major assumption of the current analysis is that long-term recharge to the groundwater can be represented as a quasi-steady unsaturated flow. While it is known that at depths >1 m or more this can be a reasonable assumption in many soils, the short-time scale variability in hydrology of the upper part of the profile, as well as the effects of heterogeneities and preferential flows, will have an effect on pollutant transport to depth. It has been shown elsewhere that some of this variability can be accounted for by treating time-variable inputs in terms of an equivalent steady-state cumulative mass flux. However, this will not always be the case, particularly where preferential flows are important. In this paper we present a methodology that would also be applicable to transient and preferential flows (given an acceptable model of the latter). Relaxing the quasi-steady-state assumption would undoubtedly require considerably more computer time, but would not change the principles of the analysis.
|
Specification of Prior Distributions of Column Parameters for the Field Application
|
|---|
The analysis requires specifying prior distributions for the four parameters defined in Table 1. It is known that these parameters can be dependent on the measurement scale. Ideally, in making the assumption that the field can be represented as a collection of vertical columns equivalent to the undisturbed soil columns sampled from the site, the prior distributions should reflect the measurement scale of those columns (noting, however, that at the scale of the columns used here, the measurements will reflect the artificial lateral boundaries resulting from sealing the sides of the columns). However, there is no theory available as to how parameter values change (jointly) with scale, nor is there much empirical evidence upon which to condition parameters by soil type, for example. Moving from column to field scale is a process of integration and will result in the expected reduced variance of the field-scale average as the scale increases since the column fluxes can be treated additively, while the conditioning processes in the current methodology will impose some constraints on the uncertainty of the posterior predictions. In what follows, we consider only the uncertainty in these averaged predictions after integration to the scale of the complete field of 10 ha.
The prior distributions for each parameter were estimated as follows:
Mean Pore Water Velocity
The pore water velocity, v (m s–1), is dependent on the spatial variability of infiltration rates and soil properties. Local infiltration into the soil is a very complex process, varying with space and time. However, ignoring surface runoff on this relatively flat field, the multi-year average annual recharge rate at the base of the profile is simply derived by subtracting the multi-year average annual evapotranspiration from average precipitation during the same period.
The average annual precipitation in the experimental area is 759 mm, and the estimate of annual actual evapotranspiration 741 mm. Accordingly, the long-term average recharge rate is only 18 mm yr–1 (i.e., 0.00205 mm h–1). The long-term average volumetric soil moisture content, 
, is assumed to be 0.2 m3 m–3. The mobile pore water fraction, m, is assumed to be uniformly distributed, in the range 0.5 to . Such a range includes most of the fractional mobile water values given in the literature (van Genuchten and Wierenga, 1976; Gaudet et al., 1977; Khan and Jury, 1990; Jacobsen et al., 1992; Beven et al., 1993; Kamra et al., 2001). The pore water velocity is then calculated by dividing the estimate of the recharge rate, q (m s–1), by the mobile water content, m (m3 m–3). The resulting range of the long-term average pore water velocities extends from 0.01025 to 0.0205 mm h–1; a uniform distribution over this range is assumed because no other information is available to estimate the distribution.
Dispersion Coefficient
The dispersion coefficient, D (m2 s–1), depends strongly on the average pore water velocity. It is, therefore, simpler to estimate the dispersivity, 
(m), which is less dependent on the velocity. Beven et al. (1993) summarized the values of dispersion parameters obtained from different experiments; recently published dispersion parameters values are given in Table 2. Most of the values vary from a few millimeters to a few decimeters. Comegna et al. (2001), Dubus and Brown (2002), and Thomasson and Wierenga (2002) suggested that should be lognormally distributed. Thus, the field dispersivity was assumed to be lognormally distributed. Taking into account the data summarized in Table 2, the range of was set from 0.05 to 285 cm, and the SD of the logarithmic transformation of values was assumed to be 0.5 in light of the values published by Beven et al. (1993), Comegna et al. (2001), Kamra et al. (2001), Schwartz et al. (2000), and Thomasson and Wierenga (2002). The dispersion coefficient D was then calculated as the product of and v.
Retardation Coefficient
The retardation factor, R, was calculated from estimated values of the distribution coefficient, Kd (m3 kg–1); soil dry bulk density, 
(kg m–3), and mobile water content; m (Table 1). The atrazine distribution coefficient is dependent on the soil OC content. Many R, Kd, or OC data can be obtained from published batch, column, and field experiments; some recently published Kd values for different soils and different OC contents are presented in Table 3. Most studies suggest that Kd values are normally distributed (e.g., Rao et al., 1986). This is confirmed by the results of the small sample batch experiments for this field site presented in Zhang (2003), which revealed that the atrazine Kd is normally distributed, with a mean of 0.329 cm3 g–1 and a SD of 0.134 cm3 g–1, leading to a CV of 0.4. Therefore, it was assumed that the field Kd values are normally distributed. The feasible range was considered to be situated between 0.1 and 3.1 cm3 g–1, and the CV was set equal to 0.4. Subsequently, R was calculated according to the sampled values of Kd, m, and .
Degradation Coefficient
The first-order decay coefficient, µ (s–1), was calculated from estimated values of the first-order degradation rate constant, k1 (s–1), under the assumption that degradation rates in both liquid and solid phases are identical, which results in the relation µ = Rk1 (Table 1). Pesticide degradation in soil is affected by many factors, including biological activity, soil properties, and environmental and experimental conditions. Thus, the degradation rate constants might show a high variability. Table 4 presents recently published atrazine degradation data. Although some researchers reported lognormal distribution of degradation rate constants at the regional scale (Brejda et al., 2000; Dubus and Brown, 2002), the batch experiments presented in Zhang (2003) showed that for the site studied here the atrazine k1 values are normally distributed with a mean of 0.0027 d–1 and CV of 0.3. Thus, it was assumed that the field k1 values are also normally distributed, and that the CV of all the potential field-scale distributions is 0.30. The limits of the k1 feasible range were fixed to 0.001 and 0.116 d–1, corresponding to half-lives of 6 and 700 d, respectively. A value of µ is then calculated from the calculated value of R (given Kd, m, and ) and the sampled value of k1.
View this table:
[in this window]
[in a new window]
|
Table 4. Summary of atrazine degradation parameters from recently published batch, column, and field experiments (k1 = 0.693/t1/2).
|
|
Correlation among the Parameters
It is commonly the case that, in fitting BTCs to experimental tracer data, high correlations among parameter values are found in the vicinity of the optimum. This is correlation in the functioning of the model representation of transport in fitting the data rather than correlation in the occurrence of parameter values that might be locally representative in the field. However, given that parameter values are required to represent every column in the field it seems reasonable to assume that parameter correlation might be found in fitting the data for all columns and should, therefore, be part of defining a multivariate distribution for the joint expectation of parameter values in forming parameter sets at the local (column) scale in contributing to field-scale transport. Unfortunately, it is not quite that simple because optimization of parameters does not give an adequate representation of the different parameter sets that might provide behavioral simulations of column-scale transport, as shown by Zhang et al. (2005) in applying the Generalized Likelihood Uncertainty Estimation (GLUE) methodology to parameter identification at the column scale. Generalized Likelihood Uncertainty Estimation is a methodology based on Monte Carlo simulation that allows many different parameter sets to be evaluated in fitting a set of observed data. Those that fit well (i.e., are behavioral) are given high likelihood values; those that do not provide acceptable fits are rejected as no-behavioral, given likelihoods of zero, and not used in prediction (e.g., Beven and Freer, 2001; Beven, 2001, 2005). In addition, any correlation between parameters in model identification is not often reported. Thus, in what follows, we assumed that the local-scale parameter values can be sampled independently, with the effects of any parameter interactions being introduced by the conditioning resulting from the multiple column scale behavioral parameter sets available from fitting the observed BTCs.
Specification of the Input Boundary Conditions for CXTFIT
In the area of the selected representative field, the groundwater level is located between 1 and 2 m below the soil surface. Consequently, the column length was set at a standard 1.5 m, and the atrazine field-scale transport was predicted to this depth, assuming a steady-state flow rate and constant effective mobile pore volume.
The long-term average recharge rate was specified as explained above. Pesticides are generally applied by farmers in a short time period. Therefore, a dirac delta (pulse) input was used to model the application of 10 mg of atrazine to the top of each individual column (491 cm2) at time zero (i.e., about 2 kg ha–1, the dose usually applied by farmers). Initial resident concentrations in the profile were assumed to be negligible.
|
Implementation of Integration at the Field Scale for Prior Estimates of Atrazine Fluxes to Groundwater
|
|---|
The feasible range of a given parameter x (as chosen from the sampled parameters , , Kd, and k1) is first defined based on the values reported from column or field experiments in the literature as outlined above. The feasible range is the widest reasonable range that can be assigned to the parameter. The range is divided into n discrete increments, as shown in Fig. 1, with the lower limit denoted x0, and the upper limit xn. The parameter value at the ith discrete increment is xi. In the following text, we will use i to indicate the location of a parameter value. If several parameters are considered, the subscripts i1, i2, ... will be assigned to the parameter increments. The subscript j, representing the jth discrete increment of x, is used to indicate the location of the mean value, xj.
The number of combinations of parameters in a finite discretization of the four-dimensional parameter space is potentially very large. For example, the relatively coarse discretization into 21 values used in what follows produces 214 = 194 481 parameter sets. A local-scale BTC at the chosen lower boundary depth is calculated for each set. Given the assumed mean, variance, and shape of the distribution associated with each parameter, each BTC is attributed a weight according to the location of the parameter set in the four-dimensional parameter space. Numerical integration is made by a linear weighted sum of the local BTCs–fluxes to get a cumulative curve at the field scale. The procedure is then repeated for all feasible parameter distributions. Only distributions for which 95% of the values fall within the feasible range for each parameter are considered.
An example is given below to show how the weights of the local BTCs are calculated. First, a one-dimensional case is treated. Let us assume a parameter x to be normally distributed with a mean at xj and standard deviation 
x. The feasible range of x is split into n discrete increments, xi (i from 0 to n) of size 
x. The weight of xi is given by the discrete form of the density function, determined by the variance of the distribution and its location in the distribution, l(xi):
 | [1] |
where j is the mean of the jth normal distribution.
The rescaled weights are
 | [2] |
Thus, a parameter value that is closer to the mean value is given a higher weight. For other distributions the weights are calculated according to their respective probability density function.
Assuming independence of the prior distributions of parameter values (lacking better information about the covariations as noted above) for a multidimensional parameter space, weights for each parameter are determined separately, and the weight of one parameter set is calculated by multiplying the weights of its four components. Thus, the weight associated to a parameter set (mi1, i2, Kdi3, k1i4) is:
 | [3] |
with the four parameters being at the i1th, i2th, i3th, i4th discrete increments, respectively.
The form and variance of the distributions of the four parameters were defined as specified above. The weight is now determined by the location of the parameter set in the joint distribution. The location is dependent on the values of f(mi1), f(i2), f(Kdi3), and f(k1i4), that is, on the means of the distributions of m, , Kd, and k1. For a single parameter x, if the field-scale distribution is a normal distribution, the location of the mean, j, should satisfy the following condition:
 | [4] |
where x is either the given constant standard deviation, or j CV if the coefficient of variation is given. This condition will guarantee that more than 95% of any field-scale distribution will be within the range [x0, xn].
Following the same approach as Beven (1993), the means of the feasible field-scale parameter distributions are assumed to be uniformly distributed, since no prior information about these means would normally be available. Again, the probability function is assigned to the means as weights, and the discrete form is used to describe the uniform distribution. If altogether m(m < n) field-scale distributions are found satisfying the above condition, the weights of the distributions with a mean j will be
 | [5] |
Assuming the distributions of the parameters are independent from each other, the weight of a four-dimensional joint distribution (
mj1, 
j2, 
dj3, 
j4) is the product of the weights of the four single distributions.
For the specified field-scale distribution (mj1, j2, dj3, j4), the integral field-scale BTC can be calculated as the weighted sum of the local-scale BTCs. Its associated weight is obtained by multiplying the weights of the four parameter distributions given by Eq. [3].
In this way, the field distributions prior to any conditioning on the experimental column results and their respective cumulative curves are of equal probability, subject to the assumptions made about the feasible range of each parameter and the form of the prior distribution within that range. In principle, some of these assumptions can easily be relaxed, for example by assuming a distribution of the variances of the local parameter distributions to allow for uncertainty in the estimation of the variance, but at the expense of a greater computational burden.
Given prior likelihoods for each of the field-scale distributions, at each time a vector of predicted BTC concentrations over all field-scale distributions and an associated vector of the field likelihood weights are assembled. The concentration–likelihood pairs are sorted on the concentrations for that time step. The sorted likelihoods are summed up to get a cumulative distribution for the associated concentration. This is repeated for each time step and the prediction limits are then plotted.
The 50, 90, and 99% quantiles of atrazine concentration predicted at the lower boundary of the profile (1.5-m depth) at each time step for a period of 100 yr are plotted in Fig. 2
. The lower limits, such as the 5 and 10% quantiles, have not been plotted because they are so low that they cannot be distinguished from the horizontal axis. The highest concentration appears about 8 yr after the atrazine application. For the 99% quantiles, which corresponds to extreme cases, the highest value is around 0.3 µg L–1, slightly higher than the usual limit admitted for a given pesticide, namely 0.1 µg L–1.
|
Introduction of the Parameter Values Obtained from the Column Experiments
|
|---|
The ranges specified for the parameters must be large enough to cover all possible values that might appear in the field. They are based on results published in the literature, and thus, they cover different soils, experimental techniques, and environmental conditions (Tables 2, 3, and 4). However, for the one specific field that is of interest, local parameter values should not present such wide variations and feasible field-scale distributions should not be attributed equal probability. Some of them should be more likely to appear than others. In this section, the weights will be updated, given the available experimental results from the large column experiments.
In this field, BTCs were determined for four undisturbed soil column experiments using both KBr tracer and atrazine. A large number of behavioral parameter sets was obtained for each experiment with the GLUE methodology (Zhang et al., 2005). Each parameter set is associated with a likelihood weight (that is set to zero for a parameter set considered evaluated as nonbehavioral). The behavioral parameter sets and their associated weights can be used to condition the weights of the parameter distributions. Note that any covariation of effective parameter values in producing a behavioral parameter set in calibration will introduce some covariation into the posterior field-scale distribution weights as a result of this conditioning process.
The marginal distributions of the Kd values of the behavioral parameter sets obtained from experiments performed on Columns A, B, C, and D are shown in Fig. 3
as an example. Given that all the values derived from experiments are located within the lowest quarter of the feasible range, low values of Kd and field-scale Kd distributions with low means will be given larger weights.
View larger version (15K):
[in this window]
[in a new window]
|
Fig. 3. Feasible ranges for Kd (0.1–3.1 cm3 g–1) and for field-scale Kd distributions (0.55–1.8 cm3 g–1) with a CV of 0.4 and distributions of the Kd values of the behavioral parameter sets obtained from Column Experiments A, B, C, and D (the field-scale normal distributions with means of 0.55 and 1.8 and CV of 0.4 are also shown). The bar shows the prior range of the feasible field-scale means under the assumptions of the analysis before conditioning on the column data. Kd = distribution coefficient.
|
|
The conditioning process for the field-scale distributions is performed as follows. Within the GLUE methodology, those parameter sets that gave the best fit to the observed BTCs on a column have the highest likelihood weight; those that were considered as giving unacceptable predictions were given a likelihood of zero. Taking Column A as an example, 1868 behavioral parameter sets were obtained for Column A (Zhang et al., 2005). Each parameter set consists of four parameter values (m, , d, 1) and is associated with a GLUE likelihood weight w4(m, , d, 1). For a specific field-scale distribution (mj1, j2, dj3, 1j4), the more that the behavioral parameter sets are located around the mean, the larger the weights of these sets and the higher the likelihood weight this distribution will have. If none of the behavioral parameter sets lies within a given potential distribution, the likelihood will be set to zero. Thus, the weight of the distribution can be calculated as
 | [6] |
with p the number of the behavioral parameters sets obtained within the GLUE methodology, w4(mk, k, Kdk, k1k) the weights of the parameter sets, and
 | [7] |
in which
where f() is the density function of each parameter and k is the kth behavioral parameter set.
The weight of the joint parameter distributions conditioned by the experiment results are then combined with the prior likelihood weights of the distributions to get posterior weights for the field-scale distributions. The weights of the field-scale distributions are updated as
 | [8] |
where n1 to n4 are the numbers of distributions considered for each of the four parameters and some of the weights in the denominator may be zero. In this way, the measured parameter sets that have a high likelihood (according to the GLUE evaluation) and that are close to the mean of the prior distribution sample will result in a high contribution to the joint posterior probability for that parameter distribution. Correlation between the parameters, introduced by the conditioning on the behavioral column parameter sets, will be reflected in these weights when used in prediction.
The posterior field-scale distributions of the parameter sets consist of the same distributions sampled for the prior estimates, but with a modified likelihood (some which might be zero if the distribution contains none of the uncertain local-scale parameter sets for one or more of the columns). Therefore, the predictions of each field-scale distribution of parameter sets will not change (since the inputs to the model are assumed known). It is only the likelihood associated with each field-scale distribution that will change under the assumptions made here. At the prior stage each field-scale distribution was equally likely, but after conditioning on the parameter sets determined from the columns in the field, some of the posterior likelihoods will be relatively large and some will be very small. Thus, the predicted distribution of field-scale pesticide leaching will be different.
Here again, the fluxes/concentrations and associated weights are calculated, and the pairs are sorted for each time step. Afterward, different prediction limits can be calculated with the obtained cumulative distribution of atrazine concentrations. For the other columns, the experimental results can be introduced following the same procedure. The 50 and 90% quantiles of atrazine concentration predictions at the lower boundary of the profile, obtained after individual conditioning with the data of Columns A, B, C and D, are shown in Fig. 4
. Comparison with Fig. 2, shows that the 50% quantile predictions are not very different, with or without conditioning, but values corresponding to the upper 90% quantile are reduced by the conditioning. That means that the effect of conditioning is to diminish extreme values because experiment results will restrict the feasible field-scale parameter distributions.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 4. The (a) 50% and (b) 90% quantiles of atrazine concentration predictions at the lower boundary of the profile (1.5-m depth) with field-scale parameter distribution conditioned individually by behavioral parameter sets of Columns A, B, C, and D.
|
|
In a second step, the column results were taken into account consecutively. After conditioning on Column A parameter sets, the results were used for a new conditioning with a second set of behavioral models obtained from fitting the observations of Column B. Each set of behavioral simulations for the experimental data (keeping with the GLUE assessments of parameter uncertainty) can be assessed as a potential sample in the conditioning process of the previous section, and the likelihood weights associated with each field-scale distribution modified to reflect the added information.
Then results from Column C are introduced, and lastly those of Column D. Note again that during the whole process, the predictions of each field-scale distribution of parameter sets do not change. Only the set of likelihood weights associated with each field-scale distribution changes. Thus, the distribution of pesticide field-scale leaching will be modified as the distribution of posterior likelihoods changes with added information.
The 50 and 90% quantiles of atrazine concentration predictions at the lower boundary of the profile (1.5-m depth) with field-scale parameter distribution conditioned consecutively by behavioral parameter sets of Columns A, B, C, and D are shown in Fig. 5
. Comparing with Fig. 4, the 90% quantiles are not very different, but values of the 50% quantile are increased by conditioning on the additional column behavioral parameter sets. That means that the difference between the 50 and 90% quantiles is further decreased.
View larger version (28K):
[in this window]
[in a new window]
|
Fig. 5. The (a) 50% and (b) 90% quantiles of atrazine concentration predictions at the lower boundary of the profile (1.5-m depth) with field-scale parameter distribution conditioned consecutively by behavioral parameter sets of Columns A, B, C, and D.
|
|
|
Summary of the Conditioning Procedure for Field-Scale Transport
|
|---|
The methodology outlined above is an extension of that proposed by Beven (1993) to reactive transport and multivariate conditioning. The whole methodology can be summarized as follows:- Choose appropriate appropriate prior distributions–ranges of variation for all the column-scale model parameters (here m, , Kd, and k1).
- Reevaluate the maximum ranges of each parameter by looking at the extremes of the column parameter sets (as determined in GLUE) and the assumed range of variation in Step 1 (the range must include all the behavioral parameter values from the GLUE analysis of the transport model fitted to the column experiments).
- Create equally probable prior field-scale distributions by moving the specified field-scale distributions across the parameter ranges for all combinations of parameters. Calculate the field-scale BTC or mass fluxes for each distribution.
- Using the likelihoods for the behavioral parameter sets for the first Column A experiment (as determined in the GLUE calibration), evaluate the joint likelihoods of the different parameter values in every behavioral "measurement" set being included in each of the field-scale distributions. This gives a one-dimensional array of likelihoods that can be summed over all the behavioral parameter sets for this column to give a posterior likelihood for the different field-scale distributions. Recalculate the posterior field-scale BTCs or mass fluxes by rescaling the likelihood weight of each field-scale distribution to obtain a unit cumulative posterior distribution.
- Repeat Step 4 consecutively for the other available columns (here B, C, and D), using the posterior weights for the field-scale distributions of the previous step as a prior distribution for adding the information associated with the behavioral parameter sets of the next column.
Figure 6
shows the 99% quantiles of atrazine concentration predictions at the lower boundary of the profile, with the field-scale parameter distribution unconditioned, conditioned with data of Column A, and conditioned with data of Columns A, B, C, and D. Some field-scale parameter distributions are restricted by the conditioning, so that extreme concentration values are decreased and the predictions of the potential risk of atrazine groundwater contamination is reduced.
View larger version (13K):
[in this window]
[in a new window]
|
Fig. 6. The 99% quantiles of atrazine concentration predictions at the lower boundary of the profile (1.5-m depth) with field-scale parameter distribution (a) unconditioned, (b) conditioned with data of Column A, and (c) conditioned with data of Columns A, B, C, and D. Note changes in concentration scale.
|
|
|
CONCLUSIONS
|
|---|
Using the proposed methodology, different prediction limits, such as the 50, 90, and 99% quantiles of atrazine concentration, provide a description of what might happen in the field, as conditioning on the limited available data. Such results will be helpful for decision-makers to estimate the risk of possible impacts caused by application of chemicals.
It was shown that the behavioral parameter sets obtained from column experiments can effectively condition the field-scale parameter distributions while reflecting the uncertainty associated with fitting the BTC experiments. Conditioning with the behavioral parameter sets obtained with the GLUE methodology from the undisturbed columns causes a significant decrease of the 90 and 99% quantiles. Consequently, and subject to the assumptions made, the long-term risk of groundwater contamination by atrazine was found to be very low.
There are some important limitations of the analysis as performed above. One of the strongest assumptions is related to the adopted upper boundary condition that considers a steady supply of water at a rate estimated from the annual average recharge. The resulting very low flow rate means that the atrazine takes a very long time to reach the groundwater table, with predicted BTCs that reflect the consequent opportunity for degradation. In reality, at least locally in the field, chemicals might be transported quickly to the groundwater by periods of rapid recharge through the vadose zone following heavy rains or preferential flows. In addition, while Zhang (2003) has shown that the first-order sorption–desorption process and degradation representations used are valid at this site, this may not be the case for other soils. The undisturbed column experiments used in conditioning also involved long time periods and higher flow rates. The potential for preferential flows under wet conditions will be reflected in the column experiments, but a simple assumption of constant diffusivity, retardation, and degradation coefficients may not reflect the nonlinearities of transport under variable flow rates in natural conditions. It would be interesting to extend the proposed methodology to fully dynamic simulations that include the possible effects of preferential flow during short rainy periods with uncertain inputs. Such an extension would require more data, parameter information, and (a lot more) computer time, but the methodology to be applied would be essentially the same.
|
ACKNOWLEDGMENTS
|
|---|
Funding for this study was provided by the Swiss National Science Foundation (Contract 2100-055793) and is gratefully acknowledged. The continuing development of the GLUE methodology is supported in part by Grant NER/L/S/2001/00658 from the UK Natural Environment Research Council. Comments from two anonymous referees helped improve the presentation of the paper.
|
REFERENCES
|
|---|
- Abu-Zreig, M., R.P. Rudra, W.T. Dickinson, and L.J. Evans. 1999. Effect of surfactants on sorption of atrazine by soil. J. Contam. Hydrol. 36:249–263.
- Amoozegar-Fard, A., D.R. Nielsen, and A.W. Warrick. 1982. Soil solute concentration distributions for spatially varying pore water velocities and apparent diffusion coefficients. Soil Sci. Soc. Am. J. 46:3–9.
- An, Y., Y.-C. Jin, and T. Viraraghavan. 1999. Solute transport in heterogeneous fields. J. Hydrol. Eng. ASCE. 4(3):293–296.
- ARS. 2001. The ARS Pesticide Properties Database [Online]. Available at http://www.ars.usda.gov/Services/docs.htm?docid=6433 (accessed 22 Aug. 2005, verified 19 Dec. 2005). USDA-ARS, Southeastern Watershed Research Lab, Tifton, GA.
- Beulke, S., C.D. Brown, I.G. Dubus, C.J. Fryer, and A. Walker. 2004. Evaluation of probabilistic modelling approaches against data on leaching of isoproturon through undisturbed lysimeters. Ecol. Modell. 179:131–144.[CrossRef]
- Beven, K.J. 1993. Estimating transport parameters at the grid scale: On the value of a single measurement. J. Hydrol. (Amsterdam) 143:109–124.
- Beven, K.J. 2001. Rainfall-runoff modelling—The primer. Wiley, Chichester, UK.
- Beven, K.J., and J. Freer. 2001. Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems. J. Hydrol. (Amsterdam) 249:11–29.
- Beven, K.J., D.E. Henderson, and A.D. Reeves. 1993. Dispersion parameters for undisturbed partially saturated soil. J. Hydrol. (Amsterdam) 143:19–44.
- Brejda, J.J., T.B. Moorman, J.L. Smith, D.L. Karlen, D.L. Allan, and T.H. Dao. 2000. Distribution and variability of surface soil properties at a regional scale. Soil Sci. Soc. Am. J. 64:974–982.[Abstract/Free Full Text]
- Bresler, E., and G. Dagan. 1981. Convective and pore scale dispersive solute transport in unsaturated heterogeneous fields. Water Resour. Res. 17:1683–1689.
- Bresler, E., and G. Dagan. 1983. Unsaturated flow in spatially variable fields. 3. Solute transport models and their application to two fields. Water Resour. Res. 19:429–435.
- Comegna, V., A. Coppola, and A. Sommella. 2001. Effectiveness of equilibrium and physical non-equilibrium approaches for interpreting solute transport through undisturbed soil columns. J. Contam. Hydrol. 50:121–138.[Medline]
- Dubus, I.G., and C.D. Brown. 2002. Sensitivity and first-step uncertainty analyses for the preferential flow model MACRO. J. Environ. Qual. 31:227–240.[Abstract/Free Full Text]
- Ersahin, S., R.I. Papendick, J.L. Smith, C.K. Keller, and V.S. Manoranjan. 2002. Macropore transport of bromide as influenced by soil structure differences. Geoderma 108:207–223.[CrossRef][ISI]
- Gaber, H.M., W.P. Inskeep, S.D. Comfort, and J.M. Wraith. 1995. Nonequilibrium transport of atrazine through large intact cores. Soil Sci. Soc. Am. J. 59:60–67.[Abstract/Free Full Text]
- Gaudet, J.P., H. Jégat, G. Vachaud, and P.J. Wierenga. 1977. Solute transfer, with exchange between mobile and stagnant water through unsaturated sand. Soil Sci. Soc. Am. J. 41:665–671.[Abstract/Free Full Text]
- Grundl, T., and G. Small. 1993. Mineral contributions to atrazine and alachlor sorption in soil mixtures of variable organic carbon and clay content. J. Contam. Hydrol. 14:117–128.
- Hornsby, A.G., D. Wauchope, and R. Herner (ed.) 1996. Pesticide properties in the environment. Springer, New York.
- Jacobsen, O.H., F.J. Leij, and M.Th. van Genuchten. 1992. Parameter determination for chloride and tritium transport in undisturbed lysimeters during steady flow. Nord. Hydrol. 23:89–104.
- Jacques, D., J. Simunek, A. Timmerman, and J. Feyen. 2002. Calibration of Richards' and convection-dispersion equations to field-scale water flow and solute transport under rainfall conditions. J. Hydrol. (Amsterdam) 259:15–31.
- Kamra, S.K., B. Lennartz, M.Th. van Genuchten, and P. Widmoser. 2001. Evaluating non-equilibrium solute transport in small soil columns. J. Contam. Hydrol. 48:189–212.[Medline]
- Khan, A.U.-H., and W.A. Jury. 1990. A laboratory study of the dispersion scale effect in column outflow experiments. J. Contam. Hydrol. 5:119–131.[CrossRef]
- Kruger, E.L., P.J. Rice, J.C. Anhalt, T.A. Anderson, and J.R. Coats. 1997. Comparative fates of atrazine and deethylatrazine in sterile and nonsterile soils. J. Environ. Qual. 26:95–101.[ISI]
- Ma, L., and H.M. Selim. 1997. Evaluation of nonequilibrium models for predicting atrazine transport in soils. Soil Sci. Soc. Am. J. 61:1299–1307.[Abstract/Free Full Text]
- Mbuya, O.S., P. Nkedi-Kizza, and K.J. Boote. 2001. Fate of atrazine in sandy soil cropped with sorghum. J. Environ. Qual. 30:70–77.
- Momii, K., Y. Hiroshiro, K. Jinn, and R. Berndtsson. 1997. Reactive solute transport with a variable selectivity coefficient in an undisturbed soil column. Soil Sci. Soc. Am. J. 61:1539–1546.[Abstract/Free Full Text]
- Moreau, C., and C. Mouvet. 1997. Sorption and desorption of atrazine, deethylatrazine, and hydroxyatrazine by soil and aquifer solids. J. Environ. Qual. 26:416–424.
- Pang, L., M.E. Close, J.P.C. Watt, and K.W. Vincent. 2000. Simulation of picloram, atrazine, and simazine leaching through two New Zealand soils and into groundwater using HYDRUS-2D. J. Contam. Hydrol. 44:19–46.[CrossRef]
- Qiao, X., L. Ma, and H.E. Hummel. 1996. Persistence of atrazine and occurrence of its primary metabolites in three soils. J. Agric. Food Chem. 44:2846–2848.[CrossRef]
- Rao, P.S.C., K.S.V. Edvardsson, L.T. Ou, R.E. Jessup, P. Nkedi-Kizza, and A.G. Hornsby. 1986. Spatial variability of pesticide sorption and degradation parameters. p. 101–115 In W.Y. Garner et al. (ed.) Evaluation of pesticides in ground water. American Chemical Society, Washington, DC.
- Schulin, R., A. Papritz, H. Flühler, and H.M. Selim. 1991. Parameter estimation for simulation of binary homovalent cation transport in aggregated cation transport in aggregated soils at variable ionic strength. J. Contam. Hydrol. 7:1–19.
- 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.
- Short, T.E., and C.G. Enfield. 1988. Large laboratory column study of the transport and degradation of atrazine, carbofuran, and diuron in soils. p. 139–150. In P. Jamet (ed.) Methodological aspects of the study of pesticide behaviour in soil/Aspects methodologiques de l'etude du comportement des pesticides dans le sol. INRA, Versailles, France.
- Spark, K.M., and R.S. Swift. 2002. Effect of soil composition and dissolved organic matter on pesticide sorption. Sci. Total Environ. 298:147–161.[Medline]
- Spongberg, A.L., and G. Lou. 2000. Adsorption of atrazine and metolachlor in three soils from Blue Creek wetlands, Watervill, Ohio. Sci. Soils 5:1–9.
- Thomasson, M.J., and P.J. Wierenga. 2002. Spatial variability of the effective retardation factor in an unsaturated field soil. J. Hydrol. (Amsterdam) 272:213–225.
- Toride, N., F.J. Leij, and M.Th. van Genuchten. 1995. The CXTFIT Code for estimating transport parameters from laboratory or field tracer experiments. Version 2.0. Res. Rep. 137. U.S. Salinity Lab., Riverside, CA.
- van Genuchten, M.Th., and P.J. Wierenga. 1976. Mass transfer studies in sorbing porous media: I. Analytical solutions. Soil Sci. Soc. Am. J. 40:473–480.[Abstract/Free Full Text]
- van Wesenbeeck, I.J., and R.G. Kachanoski. 1991. Spatial scale dependence of in situ solute transport. Soil Sci. Soc. Am. J. 55:3–7.[Abstract/Free Full Text]
- Vogeler, I., D.R. Scotter, B.E. Clothier, and R.W. Tillman. 1998. Anion transport through intact soil columns during intermittent unsaturated flow. Soil Tillage Res. 45:147–160.[CrossRef]
- Wauchope, R.D., T.M. Buttler, A.G. Hornsby, P.W.M. Augustijn-Beckers, and J.P. Burt. 1992. The SCS/ARS/CES pesticide properties database for environmental decision making. Rev. Environ. Contam. Toxicol. 123:1–164.[ISI][Medline]
- Zhang, D. 2003. Measurement scale effects on the determination of sorption and degradation parameters for modelling chemical transport in the soil. Thesis 2768. Swiss Federal Institute of Technology (EPFL), Lausanne, Switzerland.
- Zhang, D., K.J. Beven, and A. Mermoud. 2005. Goodness of fit and fitness for purpose: Model calibration and prediction uncertainty for pesticide transport in soils. Adv. Water Resour. (in press).
- Zurmühl, T. 1998. Capability of convection-dispersion transport models to predict transient water and solute movement in undisturbed soil columns. J. Contam. Hydrol. 30:101–128.
This article has been cited by other articles:
|
|
|
|
 
J. A. Vrugt, P. H. Stauffer, Th. Wohling, B. A. Robinson, and V. V. Vesselinov
Inverse Modeling of Subsurface Flow and Transport Properties: A Review with New Developments
Vadose Zone J.,
May 27, 2008;
7(2):
843 - 864.
[Abstract]
[Full Text]
[PDF]
|
|
|
|
|
|
|
D. L. Corwin, J. Hopmans, and G. H. de Rooij
From Field- to Landscape-Scale Vadose Zone Processes: Scale Issues, Modeling, and Monitoring
Vadose Zone J.,
March 8, 2006;
5(1):
129 - 139.
[Abstract]
[Full Text]
[PDF]
|
|
|
| This Article |
|
| |
TOP
ABSTRACT
INTRODUCTION
C/
