Environ. Sci. Technol. 2007, 41, 7444-7450
Including Spatial Variability in Monte Carlo Simulations of Pesticide Leaching BERTRAND LETERME,† M A R N I K V A N C L O O S T E R , * ,‡ TON VAN DER LINDEN,§ AALDRIK TIKTAK,| AND MARK D. A. ROUNSEVELL† Department of Geography, Universite´ Catholique de Louvain, Place Louis Pasteur 3, B-1348 Louvain-la-Neuve, Belgium, Department of Environmental Sciences and Land Use Planning, Ge´nie Rural, Universite´ Catholique de Louvain, Croix du Sud 2, BP 2, B-1348 Louvain-la-Neuve, Belgium, National Institute for Public Health and the Environment, PO Box 1, 3720 BA Bilthoven, The Netherlands, and Netherlands Environmental Assessment Agency (MNP), PO Box 303, 3720 AH Bilthoven, The Netherlands
A methodology is developed to quantify the uncertainty in a pesticide leaching assessment arising from the spatial variability of non-georeferenced parameters. A Monte Carlo analysis of atrazine leaching is performed in the Dyle river catchment (Belgium) with pesticide half-life (DT50) and topsoil organic matter (OM) content as uncertain input parameters. Atrazine DT50 is taken as a non-georeferenced parameter, so that DT50 values sampled from the input distribution are randomly allocated in the study area for every simulation. Organic matter content is a georeferenced parameter, so that a fixed uncertainty distribution is given at each location. Spatially variable DT50 values are found to have a significant influence on the amount of simulated leaching. In the stochastic simulation, concentrations exist above the regulatory level of 0.1 µg L-1, but virtually no leaching occurs in the deterministic simulation. It is axiomatic that substance parameters (DT50, sorption coefficient, etc.) are spatially variable, but pesticide registration procedures currently ignore this fact. Including this spatial variability in future registration policies would have significant consequences on the amount and pattern of leaching simulated, especially if risk assessments are implemented in a spatially distributed way.
Introduction Leaching modeling of plant protection products (PPPs) within the soil-crop continuum is often used in the assessment of the risk of groundwater contamination by surface applied PPPs. For example, the EU pesticides directive EU/91/414 (1) stipulates that modeled groundwater concentrations of PPPs may not exceed 0.1 µg L-1 for a single product. Modeling procedures for groundwater risk assessment are standardized * Corresponding author. E-mail: marnik.vanclooster@ uclouvain.be. Tel.: 32 10 47 37 10. Fax: 32 10 47 38 33. † Department of Geography, Universite ´ Catholique de Louvain. ‡ Department of Environmental Sciences and Land Use Planning, Universite´ Catholique de Louvain. § National Institute for Public Health and the Environment. | Netherlands Environmental Assessment Agency. 7444
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 21, 2007
by the FOCUS (FOrum for the Coordination of pesticide fate models and their USe) working groups following a step by step, tiered approach (2). At the European level, leaching models are used in combination with a limited number of standard scenarios to make conservative predictions of PPP concentrations in groundwater across the European agricultural areas. Four pesticide leaching models are currently used in these predictions: PEARL (3), MACRO (4), PELMO (5), and PRZM (6). At the level of individual member states, a variety of assessment methods are applied (e.g., 7). In The Netherlands, for example, assessments are completed on the basis of the GeoPEARL model (8, 9), which couples PEARL with spatial data derived from standard Geographical Information Systems. Tiktak et al. (10) demonstrated that a similar approach could be adopted at the pan-European level. This would take account of soil and climate variability and thus allow the identification of potential high or low risk areas. The harmonization of pan-European risk assessment methods is currently in progress (11). Unfortunately, uncertainty is often insufficiently addressed in risk assessment. Dubus et al. (12) described the different sources of uncertainty in pesticide fate modeling. Uncertainty is present in the primary data (physical, chemical, and environmental properties), arising from (i) the spatial and temporal variability of environmental variables, (ii) fieldbased sampling procedures, and (iii) laboratory analysis. There is also uncertainty in the estimation of model input parameters. This includes the direct estimation of, e.g., pesticide half-life (DT50) from laboratory data, even though harmonization of the boundary conditions and settings can limit this source of uncertainty (13), and the use of indirect methods such as pedotransfer functions and/or the interpolation of spatially referenced variables (12, 14). The structure and the numerical solution of the model provide additional sources of uncertainty (14, 15). Dubus et al. (12) note that typical probabilistic approaches (e.g., Monte Carlo analysis) often ignore the latter two sources of uncertainty. Understanding the consequences of uncertainty is needed to improve risk assessment as a decision-support tool (14). This paper addresses uncertainty in spatially distributed PPP leaching modeling. A distinction is made between georeferenced and non-georeferenced parameters. For a spatially distributed simulation, georeferenced parameters are defined as parameters that were already distributed in deterministic simulations (e.g., soil properties), whereas non-georeferenced parameters are those that were considered to be spatially constant (e.g., constant pesticide properties). This approach is analogous to the distinction between uncertainty and variability in second-order Monte Carlo simulations (16, 17). This study differs from second-order Monte Carlo simulations in using spatially distributed parameters. At a given location, georeferenced parameters are uncertain, whereas nongeoreferenced parameters display spatial variability. Notice that uncertainty refers to limited knowledge of the system, which can, in principle, be reduced by adding more data or increasing our understanding of the system, whereas variability is inherent uncertainty, e.g., spatial variability of organic matter, and cannot be reduced. Studies examining the effect of spatial variability of PPP leaching by stochastic simulation have suggested that spatial variability can significantly affect the leached fractions of PPPs (18, 19). In general, stochastic simulations are likely to generate extreme events that are not captured within an “average”, deterministic simulation. A number of modeling studies have used Monte Carlo techniques to include the spatial variability of soil and/or pesticide properties in risk 10.1021/es0714639 CCC: $37.00
2007 American Chemical Society Published on Web 09/21/2007
FIGURE 1. Parametrization of OM content and texture for the simulation profiles. assessment methods. Earlier work focused mainly on the variability of soil properties within soil units, whereas pesticide properties were kept constant (20-23). Recently, attempts were made to include non-georeferenced variables, such as the pesticide half-life (61). Monte Carlo simulation was also applied to non-georeferenced parameters. Zacharias et al. (24) developed a stochastic framework at the field scale in which multiple realizations of a hypothetical field were undertaken with the same statistical distribution of soil properties. The Monte Carlo simulation, however, considered all soil and pesticide parameters as non-georeferenced. This means that the outputs could be examined only after aggregation at the fieldscale. Lindahl et al. (25) combined georeferenced field management parameters (such as pesticide application dose) with a range of non-georeferenced parameters to study the contamination of surface water. They found that simulated values of pesticide loads reached concentrations similar to measurements mainly because of small but critical summer flows captured in the stochastic approach. The objective of this study is to determine whether leaching at the regional scale is significantly affected by the inclusion of the spatial variability of non-georeferenced parameters. This objective has been tackled by refining, applying, and evaluating a Monte Carlo approach that samples both georeferenced and non-georeferenced parameters in a spatially distributed assessment of PPP leaching. The innovative part of this approach is that the inclusion of non-georeferenced parameters is performed in a spatially distributed way (no spatial aggregation). A stochastic simulation of atrazine (2-chloro-4-ethylamino-6-isopropylamino1,3,5-triazine) leaching was performed for a well-documented case study with soil organic matter content (OM; georeferenced parameter) and pesticide half-life (DT50; non-georeferenced) as the uncertain input parameters. The issues of reproducibility and truncation levels of the input probability density functions were also addressed. Finally, although data availability limits direct validation, the results are discussed with respect to information about groundwater contamination in the study area.
Materials and Methods Catchment Description. The study was undertaken in the part of the Dyle catchment located in the Walloon region (580 km2) of Belgium. The Brusselian aquifer is located in tertiary sands overlain by a quaternary loess layer (0-5 m). The aquifer is of primary importance for the supply of
drinking water and has been classified as a “vulnerable area” following the implementation of the European Union Nitrate Directive (26). Groundwater depth varies from 0 to more than 30 m. The climate is temperate maritime, with a mean annual temperature of 9.7 °C and mean summer and winter temperatures of 16.4 and 3.0 °C, respectively. Average annual precipitation is 804 mm and no significant variation is found within the catchment. Typical land uses are urban (generally located in the valleys), pasture and forests (found on valley slopes), and arable (loamy soils on the plateau). The simulation of atrazine leaching was restricted to arable land, which covers approximately 46% of the study area (27) and for which the dominant soil type can be classified as a luvisol. Model Description. GeoPEARL couples the PEARL model (3) with spatial data derived from standard Geographical Information Systems. PEARL is a one-dimensional, dynamic, multilayered model of the fate of pesticides and relevant transformation products in the soil-plant system. The model is linked to the soil water atmosphere plant model (SWAP, 28) for hydrology. Soil water flow is described with the Richards’ equation and reference evapotranspiration is calculated with the FAO modified Penman-Monteith approach (29). Pesticide transport is simulated with the convection-dispersion equation. Sorption onto the soil solid phase is described with a Freundlich isotherm. The degradation of pesticides is described with a first-order rate equation and a number of reduction factors, which account for the influence of temperature, soil moisture, and soil depth (30). Parametrization. Annual (spring) applications of atrazine were simulated for silage maize cropping. The atrazine coefficient of sorption on organic matter (Kom) was assumed constant to be (74 L kg-1; confidential data). The simulation period was fixed as 1980-2002, including six initialization years. Annual average concentrations of atrazine at 1 m depth were extracted from the results by dividing the annual cumulative leached amount of pesticide with the annual cumulative water drainage. Climatic variability was ignored, as test simulations in the catchment showed it to be unimportant (results not shown). Over larger areas, climatic conditions become more important (e.g., 31). Figure 1 shows a flow chart of the parametrization of organic matter content and texture for the stochastic simulation. Soil properties were interpolated from soil profiles and a digital soil map (scale 1:20 000). The Aardewerk database (32) for Belgium consists of more than 10 000 soil profile descriptions (texture, organic matter, pH, etc.) for the different horizons down to 1.5 m. Three-hundred ninetyVOL. 41, NO. 21, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
7445
three profiles with arable land use were available within the arable part of the study area. A buffer zone of 4 km around the study area was included because neighboring soil profiles could provide valuable information. Pedotransfer functions (PTFs) were used to estimate soil input parameters not available in Aardewerk. Soil dry bulk density was calculated with the PTF of Bollen et al. (33) and 21 parameters of the water retention equations (35) were derived with the PTFs of Wo¨sten et al. (34). Organic matter content of the surface horizon was interpolated using ordinary kriging. Topsoil OM has changed since the 1950s (when the data were collected) and so values were adjusted using a correction factor available for each soil association and land use (36). The resulting map is shown in Figure A1 (in the Supporting Information). Finally, an empirical relationship was used to determine organic matter content with depth (37)
OM(z) ) OMb + (OM0 - OMb)exp(-kz)
(1)
where z is the depth (cm); OM0 and OMb are the organic matter contents in the top horizon (25 cm) and at the bottom of the soil profile (150 cm), respectively; and k is a study area specific constant. The average profile of organic matter content and the mathematical fit are given in Figure A2 (in the Supporting Information). Interpolation of texture fractions in the surface horizons was on the basis of the Bayesian maximum entropy (BME) approach, which allows the inclusion of hard (accurate) and soft (vague) data in a spatial estimation context (38, 39). Hard data were the texture fractions of Aardewerk profiles, and soft data consisted of the texture class given by the 1:20 000 soil map at every estimation location. A variant of the regular BME algorithm using a Monte Carlo procedure, called BME/ MC (40), was used because it takes into account the fundamental constraints on the textural fractions (they sum to one and belong to the (0, 1) interval). The 1:20 000 soil map was then used to define average soil profiles with depth. The Aardewerk profiles were associated with map units using a “class matching” procedure (41). The following rules were used to determine the groups of soil series needing an average profile with depth: • Texture class (from the Belgian texture triangle; 42) • Profile development or not • Presence/absence of a sandy substrate Some groups were split into further subgroups if these had a high number (>30) of Aardewerk profiles available to parametrize the subgroups. This led to the definition of 16 groups of soil series. For each group, the available Aardewerk profiles were used to calculate average values of the texture fractions and depth to a sandy substrate for the deeper horizons (from 25 to 150 cm). Variable horizon thicknesses were accounted for by discretizing the soil profiles. Because of the relatively small variability in OM and texture over the study area (e.g., silty soils cover more than 90% of the study area), it was possible to reduce the number of plots from 26 606 to 100 using k-means cluster analysis. Test simulations proved that the results would not be significantly affected by this clustering (results not shown). This substantially reduced the run time of the Monte Carlo analysis. The variables used in the cluster analysis were (i) the soil series group (ensuring the correct retention of these groups in the cluster analysis), (ii) organic matter content, (iii) silt fraction, and (iv) clay fraction, in the top horizon (r ) 0.0344 between silt and clay fractions). Uncertainty Analysis. The uncertainty analysis was based on a Monte Carlo (MC) approach, modified to account for the spatial variability of non-georeferenced parameters. Selection of the Parameters for the Monte Carlo Procedure. To evaluate the methodology, we included only 7446
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 21, 2007
FIGURE 2. Cumulative density functions and lognormal fits of atrazine DT50 values. Squares refer to site specific data of Pussemier et al. (48) and circles to data from the Dutch registration dossier (49). two parameters in the uncertainty analysis (one georeferenced and one non-georeferenced). The choice of the parameters was on the basis of a previous sensitivity analysis performed on PESTLA (a predecessor of PEARL) (43), which found the Freundlich exponent, Kom, and DT50 to be the most sensitive pesticide parameters, and organic matter content and bulk density the most sensitive soil parameters. Thus, the analysis presented here focused, a priori, on the uncertainty of atrazine DT50 and organic matter content, i.e., a non-georeferenced and a georeferenced parameter, respectively. Restricting the analysis to two parameters allowed us to decrease the simulation computation cost. Definition of the PDFs. Organic matter content is a georeferenced parameter, because the interpolation using ordinary kriging led to a spatially distributed characterization of organic matter content. Following the cluster analysis, 100 values of organic matter content were assigned to the corresponding clusters. For each cluster, the uncertainty distributions around these values were obtained using block kriging. In other words, the ordinary kriging (together with the class matching procedure) allowed the clusters delineation, whereas block kriging performed on these clusters was necessary to evaluate the uncertainty around the central values. DT50 is a non-georeferenced parameter because no georeferenced data are available for this parameter (i.e., initial parametrization of GeoPEARL assumed a spatially constant value). However, it is widely recognized that pesticide DT50 values can show a large variability at the field or catchment scale (44-46). Coquet et al. (47) argued that whenever available, site-specific data should be preferred over databases to limit bias in pesticide leaching risk assessments at the catchment scale. Therefore, data from Pussemier et al. (48) were used to assess the distribution of DT50 values. These data consist of 33 laboratory measurements of the atrazine degradation rate in loam or sandy loam soils; most located in the Dyle catchment. DT50 of less than 10 days was found in more than 60% of the soil samples. The site-specific DT50 measurements of Pussemier et al. (48) are much lower than the Dutch registration dossier data (see Figure 2). The rapid dissipation observed by Pussemier et al. (48) was linked to the adaptation of microbial communities resulting from repeated pretreatments with atrazine (intensive maize cropping). Such observations are excluded from the Dutch registration dossier if the effect is known to occur because these data are considered as inappropriate in a registration context.
Parameter Sampling. Organic matter content and DT50 input parameters were generated using a Latin hypercube sampling (LHS) procedure. LHS uses a probabilistic sampling scheme that provides better coverage of the input distribution and has been shown to be more efficient than the random sampling in MC analysis (50). A single GeoPEARL simulation with the MC procedure corresponds to 100 PEARL runs over the whole study area. The total number of model runs is therefore equal to N × 100, where N is the number of MC simulations. For each unique cluster, the normal distribution of the georeferenced parameter (organic matter content) was sampled N times. Furthermore, for each of the N simulations, 100 values from the lognormal distribution of the non-georeferenced parameter (DT50) were drawn with LHS and randomly allocated to the 100 unique clusters. Thus a major difference between the georeferenced and non-georeferenced parameters lies in the way their PDFs were sampled in the MC analysis. Truncation levels were arbitrarily set to the first and 99th percentiles of the distributions used for organic matter content, and to the fifth and 95th percentiles of the distribution used for atrazine DT50. Automatic Running of GeoPEARL. A Matlab (version 7.0) routine was created for the generation of LHS input parameters and these replaced the default values in the GeoPEARL input files. The N MC simulations (i.e., N × 100 model runs) were processed on a grid of 256 central processing units (CPUs) using batch files. Stability of the results was tested by iteratively increasing the number of samples, because the accuracy of the MC technique is inversely proportional to the square root of the number of runs N (14). N was successively fixed at 10, 50, 100, 150, and 200. Analysis of the Model Outputs. An indicator of atrazine leaching potential was defined as the 80th percentile of the annual average concentrations at a 1 m depth (for a total of 17 simulation years). For a given location in the catchment, this corresponds to the 80th percentile of leaching because of variations in weather conditions (2). This indicator could be mapped for any of the N MC simulations. Indication of leaching risk at the catchment scale was chosen as (i) the 80th percentile in space (51) and (ii) the percentage of the area exceeding the regulatory limit of 0.1 µg L-1. Minimum, median, and maximum values of these two indicators were reported against the increase in the number N of simulations and were compared to the results of a deterministic GeoPEARL assessment with average OM and DT50 values. Reproducibility and the Effect of Truncation. Reproducibility should be examined when undertaking uncertainty analysis. The study of reproducibility of MC simulations looks at the influence that the generation of a random sample with a particular seed number may have on the overall outcome of MC modeling (52). The robustness of the MC analysis was assessed for N ) 50 by running the simulations with six different seed numbers in the LHS method. The assessment of stability was restricted to median values of model outputs, as the minimum and maximum values would be expected to vary in an unpredictable way. A major issue concerning the influence of modeler subjectivity is the truncation level of the input PDFs (57). The MC analysis for N ) 50 was repeated with truncation set to the first and 99th percentiles of the lognormal distribution for atrazine DT50. Truncation was unchanged for organic matter content.
Results and Discussion Deterministic Simulation. The 80th percentile of the annual average concentrations at 1 m depth was lower than 10-6 µg L-1 for all the study area using average values for OM and DT50, i.e., several orders of magnitude below standard
FIGURE 3. (a) Eightieth percentile of the annual average concentrations at 1 m depth for the deterministic GeoPEARL assessment. (b) Maximum simulated values of the 80th percentile of the annual average concentrations of atrazine (µg L-1) at 1 m (Monte Carlo analysis; N ) 100 simulations). White areas are nonarable land. detection limits (Figure 3a). This result is explained by the low mean DT50 value of 8.97 days, which was derived from the distribution fitted to the data of Pussemier et al. (48, Figure 2). Stochastic Simulation. Table 1 summarizes the results of the Monte Carlo simulations. The two leaching indicators are significantly different than in the deterministic simulation. The median 80th percentile in space is always around 0.0030 µg L-1, about 2 orders of magnitude lower than the regulatory limit of 0.1 µg L-1. The median value of this indicator appears to be stable with increasing N. Furthermore, the range VOL. 41, NO. 21, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
7447
TABLE 1. Results of the Monte Carlo Analysis N 80th percentile in space (µg L-1) part of the catchment with > 0.1 µg L-1 (%)
minimum median maximum minimum median maximum
between the minimum and maximum values tends to increase with a higher number of simulations (N). This could be expected because as N increases, so does the probability that the (random) spatial allocation of DT50 values to the 100 plots will produce extreme cases (e.g., high DT50 values located in areas with low organic matter content, leading to higher vulnerability). Figure 3b shows the maximum simulated values of the 80th percentile of the annual average concentrations of atrazine (µg L-1) at a 1 m depth. The location of high concentrations is similar in the stochastic and deterministic simulations (Figure 3a), although the stochastic simulations give much larger absolute values. This suggests that the spatial pattern of leaching risk does not depend on the adoption of a probabilistic approach. The pattern of risk appears to be negatively correlated with organic matter content in the surface horizon (Figure A1 in the Supporting Information), reflecting the influence of sorption on atrazine fate. The 80th percentile in space remains below the regulatory limit of 0.1 µg L-1 (cf. Table 1), and the percentage of the study area exceeding 0.1 µg L-1 ranges from 5.1 to 12.1% for N ) 100 (see Figure 4 and Table 1). The median value of the area above 0.1 µg L-1 is also found to be relatively stable for all N. The results suggest that the median indicator values can be estimated with confidence using N ) 50 ( ) 25 × number of parameters). The simulation of atrazine leaching was strongly affected by the inclusion of spatial variability in DT50. The uncertainty in organic matter content (a georeferenced parameter) had a negligible influence on the results because of a lower coefficient of variation (3.9% on average compared to 117% for DT50). This does not seem to support the choice of organic matter content as the most important georeferenced parameter in the present study. This conclusion is probably specific to the study area, where the limited variability in
FIGURE 4. Cumulative density functions of the results of the Monte Carlo simulations with N ) 100. The vertical dash-dotted line indicates the regulatory limit of 0.1 µg L-1. The x-axis shows the 80th percentile of the annual average concentrations at 1 m depth for a total of 17 simulation years. The y-axis starts at 0.8. 7448
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 21, 2007
10
50
100
150
200
0.0014 0.0028 0.0046 6.8 9.2 10.7
0.0015 0.0031 0.0092 6.7 8.9 10.6
0.0010 0.0033 0.0073 5.1 8.6 12.1
0.0011 0.0033 0.0097 6.2 8.7 11.9
0.0010 0.0034 0.0102 5.7 8.6 12.3
organic matter content leads to low block kriging variances. It should also be noted that the conclusions derived from the uncertainty analysis pertain to only one substance, namely atrazine. For other compounds, the results could be different. Furthermore, the approach is essentially a parameter uncertainty analysis (although uncertainty due to interpolation is also addressed for organic matter content). Model conceptual errors, such as the absence of preferential flow, have not been taken into account. This reflects the difficulty in separating the different sources of error in dynamic nonlinear cases (53). Therefore, the conclusions presented here are valid only in the absence of significant preferential flow processes at the regional scale. The deterministic GeoPEARL assessment produced virtually no leaching, mainly because of the very low average DT50 value extracted from the field data of Pussemier et al. (48). However, significant leaching was simulated in the MC analysis using a lognormal distribution of DT50. Thus, the inclusion of spatial variability of any non-georeferenced parameter is likely to modify the conclusions drawn from simulations performed with spatially distributed models. Taking this variability into consideration could provide a better evaluation of the risk of contamination (54), especially if future risk assessments are to be implemented in a spatially distributed way. Geisler et al. (55) have advocated that the spatial variability of substance properties should be considered in regulatory procedures. However, the huge variability found even at the field scale (e.g., 56) suggests that substance properties should be considered as non-georeferenced parameters in most applications. Reproducibility and Truncation. Using different seed numbers did not appear to modify the examined outputs in a significant way, especially considering the median values (cf. Table A1 in the Supporting Information). The 80th percentile in space was stable around 0.0031 and the percentage of the study area above 0.1 µg L-1 remained between 8.5 and 9%. These results suggest that the MC analysis was robust with respect to the randomness of the seed number. Table A1 (in the Supporting Information) also gives the results of the MC analysis using a lower level of truncation in the sampling of the DT50 PDF. When the first and 99th percentiles of the lognormal distribution were chosen, the sampled values of atrazine DT50 reached a maximum at about 78 days, whereas the maximum was at about 41 days for the initial configuration (truncation at the fifth and 95th percentiles). This led to a significant increase in the amount of simulated pesticide leaching. The 80th percentile in space doubled and the area above 0.1 µg L-1 increased to 12%. This finding illustrates the importance of the subjective choices made by the modeler in typical MC analyses, which introduce variability in the probabilistic modeling. Beulke et al. (57) argued that user subjectivity biases truncation, type and parametrization of the distributions, correlation between parameters, methods for sampling, and size of random samples. In the present study, truncation of the input PDFs clearly had an impact on the results, but it is rather difficult to define truncation thresholds in an objective way. However, performing the Monte Carlo analysis with a different level
of truncation for atrazine DT50 provided information on the impact of this choice on the simulation outputs. Field Samples versus Databases. Different studies (47, 58) recommended field-sample measurements rather than databases to parametrize the distribution of DT50; this was found to be important here. Using DT50 data from the Dutch registration dossier led to concentrations above 0.1 µg L-1 for more than 80% of the study area, in all MC simulations with N ) 50 (results not shown). However, the contamination of groundwater by atrazine is believed to be mainly from non-agricultural sources (59, 60). Thus, the DT50 of Pussemier et al. (48) produced more realistic simulations, i.e., rarely exceeding 0.1 µg L-1 for the agricultural sources. This supports the appropriateness of field-sample measurements when deriving pesticide parameters. It should be noted, however, that model error was not assessed in this study (among other sources of uncertainty), e.g., GeoPEARL does not account for preferential flow. A good match of modeled and field observations would, therefore, not be expected. Future work will seek to add more sensitive parameters to the analysis and to assess whether the results are still relevant in the context of pesticide registration. There is also a need to evaluate the interaction between several nongeoreferenced parameters in MC simulations.
Acknowledgments This work was supported by the Fonds National de la Recherche Scientifique (FNRS-Belgium).
Supporting Information Available Two figures presenting the parametrization of organic matter content, and a table reporting the results of the simulations concerning the reproducibility of the MC analysis and the effect of truncation (PDF). This material is available free of charge via the Internet at http://pubs.acs.org.
Literature Cited (1) European Commission. Council Directive 91/414/EEC of 15 July concerning the placing of plant protection products on the market. Official J. Eur. Communities 1991, 230. (2) FOCUS Groundwater Scenarios in the EU Review of Active Substances; Reference Sanco/321/2000 rev.2; European Commision: Brussels, Belgium, 2000. (3) Tiktak, A.; van den Berg, F.; Boesten, J. J. T. I.; Leistra, M.; van der Linden, A. M. A.; van Kraalingen, D. Pesticide Emission Assessment at Regional and Local Scales: User Manual of FOCUS Pearl Version 1.1.1; RIVM Report 711401008, Alterra Report 28; National Institute for Public Health and Environment: Bilthoven, The Netherlands, 2000. (4) Larsbo, M.; Jarvis, N. J. MACRO 5.0. A model of water flow and solute transport in macroporous soil. Technical description. Studies in the Biogeophysical Environment; Emergo 2003:6; Department of Soil Sciences, Swedish University of Agicultural Sciences: Uppsala, Sweden. (5) Klein, M. PELMO: Pesticide Leaching Model, User Manual Version 2.01; Fraunhofer-Institut fur Umweltchemie und Okotoxikogie: Schmallenberg, Germany, 1995. (6) Carsel, R. F.; Imhoff, J. C.; Hummel, P. R.; Cheplick, J. M.; Donigian, A. S., Jr. PRZM-3, Model for Predicting Pesticide and Nitrogen Fate in the Crop Root and Unsaturated Soil Zones: Users Manual for Release 3.0; U.S. Environmental Protection Agency: Athens, GA, 1998. (7) Van der Linden, A. M. A.; Boesten, J. J. T. I.; Cornelese, A. A.; Kruijne, R.; Leistra, M.; Linders, J. B. H. J.; Pol, J. W.; Tiktak, A.; Verschoor, A. J. The New Decision Tree for the Evaluation of Pesticide Leaching from Soils; RIVM report 601450019/2004; National Institute for Public Health and the Environment: Bilthoven, The Netherlands, 2004. (8) Tiktak, A.; de Nie, D. S.; van der Linden, T.; Kruijne, R. Modelling the leaching and drainage of pesticides in The Netherlands: the GeoPEARL model. Agronomie 2002, 22, 373-387. (9) Tiktak, A.; van der Linden, A. M. A.; Boesten, J. J. T. I. The GeoPEARL Model. Model Description, Applications and Manual; RIVM report 716601007; National Institute for Public Health and the Environment: Bilthoven, the Netherlands, 2003.
(10) Tiktak, A.; de Nie, D. S.; Pin ˜ eros, Garcet, J. D.; Jones, A.; Vanclooster, M. Assessment of the pesticide leaching risk at the Pan-European level. The EuroPEARL approach. J. Hydrol. 2004, 289, 222-238. (11) Report of the FOCUS Groundwater WorkGroup; EC Document Reference Sanco/???/2007, version 1; European Commission: Brussels, Belgium, in preparation. (12) Dubus, I. G.; Brown, C. D.; Beulke, S. Sources of uncertainty in pesticide fate modelling, Sci. Total Environ. 2003, 317, 53-72. (13) FOCUS. Guidance Document on Estimating Persistence and Degradation Kinetics from Environmental Fate Studies on Pesticides in EU Registration. Report of the FOCUS Work Group on Degradation Kinetics, EC Document Reference Sanco/10058/ 2005 version 2.0; European Commission: Brussels, Belgium, 2006. (14) Brown, J. D.; Heuvelink, G. B. M. Assessing uncertainty propagation through physically based models of soil water flow and solute transport. In The Encyclopedia of Hydrological Sciences; Anderson, M.G., Ed.; John Wiley & Sons: New York, 2005; pp 1181-1195. (15) Addiscott, T. M.; Smith, J.; Bradbury, N. Critical evaluation of models and their parameters, J. Environ. Qual. 1995, 24, 803807. (16) EUFRAM. Introducing probabilistic methods in to the ecological risk assessment of pesticides; EUFRAM report Vol. 1, version 6; European Commission: Brussels, Belgium, May 11, 2005. (17) Wu, F.-C.; Tsang, Y.-P. Second-order Monte Carlo uncertainty/ variability analysis using correlated model parameters: application to salmonid embryo survival risk assessment. Ecol. Modell. 2004, 177, 393-414. (18) Jury, W. A.; Gruber, J. A stochastic analysis of the influence of soil and climatic variability on the estimate of pesticide groundwater pollution potential. Water Resour. Res. 1989, 25, 2465-2474. (19) van der Zee, S. E. A. T. M.; Boesten, J. J. T. I. Effects of soil heterogeneity on pesticide leaching to groundwater. Water Resour. Res. 1991, 27, 3051-3063. (20) Carsel, R. F.; Parrish, R. S.; Jones, R. L.; Hansen, J. L.; Lamb, R. L. Characterizing the uncertainty of pesticide leaching in agricultural soils. J. Contam. Hydrol. 1988, 2, 111-124. (21) Petach, M. C.; Wagenet, R. J.; DeGloria, S. D. Regional water flow and pesticide leaching using simulations with spatially distributed data. Geoderma 1991, 48, 245-269. (22) Foussereau, X.; Hornsby, A. G.; Brown, R. B. Accounting for variability within map units when linking a pesticide fate model to soil survey. Geoderma 1993, 60, 257-276. (23) Zhang, H.; Haan, C. T.; Nofziger, D. L. An approach to estimating uncertainties in modeling transport of solutes through soils. J. Contam. Hydrol. 1993, 12, 35-50. (24) Zacharias, S.; Heatwole, C. D.; Persaud, N.; Bruggeman, A. C.; Kumar, D.; Smith, C. N. Stochastic simulation of field-scale pesticide transport using Opus and GLEAMS. J. Environ. Qual. 1999, 28, 411-423. (25) Lindahl, A. M. L.; Kreuger, J.; Stenstrom, J.; Gardenas, A. I.; Alavi, G.; Roulier, S.; Jarvis, N. J. Stochastic modeling of diffuse pesticide loss from a small agricultural catchment. J. Environ. Qual. 2005, 34, 1174-1185. (26) European Commission. Council Directive 91/676/EEC of 12 December concerning the Protection of Waters against Pollution caused by Nitrates from Agricultural Sources. Official J. Eur. Communities L 1991, 375. (27) SIGEC. Agricultural land use data from the 1999 agricultural census; Syste`me Inte´gre´ de Gestion Et de Controˆle, Ministe`re de la Re´gion Wallonne, Direction Ge´ne´rale de l’Agriculture: Wallone, Belgium, 1999. (28) Van Dam, J. C. Field-Scale Water Flow and Solute Transport. SWAP Model Concepts, Parameter Estimation, and Case Studies. Ph.D. thesis, Wageningen Agricultural University, Wageningen, The Netherlands, 2000. (29) Allen, R. G.; Pereira, L. S.; Raes, D.; Smith, M. Guidelines for Computing Crop Water Requirements; FAO Irrigation and Drainage paper; Food and Agriculture Organization: Rome, 1998; Vol. 56. (30) Boesten, J. J. T. I.; van der Linden, A. M. A. Modelling the influence of sorption and transformation on pesticide leaching and persistence. J. Environ. Qual. 1991, 20, 425-435. (31) Bleecker, M.; DeGloria, S. D.; Hutson, J. L.; Bryant, R. B.; Wagenet, R. J. Mapping atrazine leaching potential with integrated environmental databases and simulation models. J. Soil Water Conserv. 1995, 50, 388-394. (32) Van Orshoven, J.; Vandenbroucke, D. Guide de l’utilisateur de AARDEWERK. Base de donne´es de profils pe´dologiques; IRSIAVOL. 41, NO. 21, 2007 / ENVIRONMENTAL SCIENCE & TECHNOLOGY
9
7449
(33)
(34)
(35) (36)
(37) (38) (39) (40) (41) (42) (43) (44) (45) (46) (47)
(48) (49)
COBIS-Instituut voor Land en Waterbeheer, Katholieke Universiteit Leuven: Leuven, Belgium, 1993. Bollen, M. J. S.; Bekhuis, F. H. W. M.; Reiling, R.; Scheper, E. To a Spatial View of the Vulnerability of the Soil and Groundwater. Part 1: The Vulnerability for Loading with Heavy Metals, Pesticides, Groundwater Discharge and Ammonia Deposition; RIVM report 711901012 (in Dutch); National Institute for Public Health and the Environment: Bilthoven, the Netherlands, 1995. Wo¨sten, J. H. M.; Veerman, G. J.; de Groot, W. J. M.; Stolte, J. Waterretentie en doorlatendheidskarakteristieken van boven en ondergronden in Nederland; Alterra report 153; Alterra: de Staringreeks, The Netherlands, 2001. van Genuchten, M. T. A closed-form equation for predicting the hydraulic conductivity of unsaturated soils. Soil Sci. Soc. Am. J. 1980, 44, 892-898. van Wesemael, B.; Lettens, S.; Roelandt, C.; Van, Orshoven, J. Changes in soil carbon stocks from 1960 to 2000 in the main Belgian cropland areas, Biotechnol., Agron., Soc. Environ. 2004, 8, 133-139. Sleutel, S.; De Neve, S.; Hofman, G. Estimates of carbon stock changes in Belgian cropland. Soil Use Manage. 2003, 19, 166171. Christakos, G. A Bayesian/maximum-entropy view to the spatial estimation problem. Math. Geol. 1990, 22, 763-777. Christakos, G. Modern Spatiotemporal Geostatistics; Oxford University Press: New York, 2000. Bogaert, P.; D’ Or, D. Estimating soil properties from thematic soil maps: the Bayesian Maximum Entropy approach. Soil Sci. Soc. Am. J. 2002, 66, 1492-1500. Van Orshoven, J. Assessing hydrodynamic land qualities from soil survey data. Ph.D. thesis, K.U. Leuven, Leuven, Belgium, 1993. Tavernier, R.; Mare´chal, R. Carte des associations de sols de la Belgique. Pe´dologie 1958, 8, 134-182. Dubus, I. G.; Brown, C. D.; Beulke, S. Sensitivity analyses for four pesticide leaching models. Pest Manage. Sci. 2003, 59, 962982. Walker, A.; Brown, P. A. Measurement and prediction of chlorsulfuron persistence in soil. Bull. Environ. Contam. Toxicol. 1983, 30, 365-372. Smith, C. N.; Parrish, R. S.; Carsel, R. F. Estimating sample requirements for field evaluations of pesticide leaching. Environ. Toxicol. Chem. 1987, 6, 343-357. Charnay, M.-P.; Tuis, S.; Coquet, Y.; Barriuso, E. Spatial variability in 14C-herbicide degradation in surface and subsurface soils. Pest Manage. Sci. 2005, 61, 845-855. Coquet, Y.; Hadjar, D.; Gilliot, J.-M.; Charnay, M.-P.; Moeys, J.; Dufour, A.; Beaudoin, N. Biases in the spatial estimation of pesticide loss to groundwater. Agron. Sustainable Dev. 2005, 25, 465-472. Pussemier, L.; Goux, S.; Vanderheyden, V.; Debongnie, P.; Tressnie, I.; Foucart, G. Rapid dissipation of atrazine in soils taken from various maize fields. Weed Res. 1997, 37, 171-179. Dorgelo, F. O. Eindrapportage CTBase 2005; Board for the Authorisation of Pesticides: Wageningen, The Netherlands, 2006; in Dutch.
7450
9
ENVIRONMENTAL SCIENCE & TECHNOLOGY / VOL. 41, NO. 21, 2007
(50) Helton, J. C.; Davis, F. J. Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems. Reliab. Eng. Syst. Saf. 2003, 81, 23-69. (51) Vanclooster, M.; Pin ˜ eros Garcet, J. D.; Boesten, J. J. T. I.; van den Berg, F.; Leistra, M.; Smelt, J. H.; Jarvis, N. J.; Roulier, S.; Burauel, P.; Vereecken, H.; Wolters, A.; Linneman, V.; Fernandez, E.; Trevisan, M.; Capri, E.; Padovani, L.; Klein, M.; Tiktak, A.; van der Linden, A. M. A.; de Nie, D. S.; Bidoglio, G.; Baouroui, F.; Jones, A.; Armstrong, A. Effective Approaches for Asessing the Predicted Environmental Concentrations of Pesticides (APECOP) Final Report; European Commision: Brussels, Belgium, 2004. (52) Dubus, I. G.; Janssen, P. H. M. Issues of replicability in Monte Carlo modeling: a case study with a pesticide leaching model. Environ. Toxicol. Chem. 2003, 22, 3081-3087. (53) Beven, K. On the concept of model structural error. Water Sci. Technol. 2005, 52, 167-175. (54) Soulas, G.; Lagacherie, B. Modelling of microbial degradation of pesticides in soils. Biol. Fertil. Soils 2001, 33, 551-557. (55) Geisler, G.; Hellweg, S.; Liechti, S.; Hungerbuhler, K. Variability assessment of groundwater exposure to pesticides and its consideration in life-cycle assessment. Environ. Sci. Technol. 2004, 38, 4457-4464. (56) Elabd, H.; Jury, W. A.; Cliath, M. M. Spatial variability of pesticide adsorption parameters. Environ. Sci. Technol. 1986, 20, 256260. (57) Beulke, S.; Brown, C. D.; Dubus, I. G.; Galicia, H.; Jarvis, N. J.; Schaefer, D.; Trevisan, M. User-subjectivity in Monte Carlo modelling of pesticide exposure. In Pesticide Behaviour in Soils, Water and Air; SCI, University of Warwick: Coventry, U.K., 2006. (58) Muller, K.; Smith, R. E.; James, T. K.; Holland, P. T.; Rahman, A. Spatial variability of atrazine dissipation in an allophanic soil. Pest Manage. Sci. 2003, 59, 893-903. (59) Pre´vention de la Pollution des eaux par les Pesticides en Re´gion Wallonne. Partie 3: Recherche des Causes de cas de Contamination des eaux Potabilisables en Re´gion Wallonne; Convention SPGE-Phytofar-SPF Sante´ Publique, Se´curite´ de la Chane alimentaire et Environnement, Final Report 2004; Centre d’Etude et de Recherches Ve´te´rinaires et Agrochimiques: Brussels, Belgium, 2004. (60) Leterme, B.; Vanclooster, M.; Rounsevell, M. A. D.; Bogaert, P. Discriminating between point and non-point sources of atrazine contamination of a sandy aquifer. Sci. Total Environ. 2006, 362, 124-142. (61) Van der Linden, A. M. A.; Tiktak, A.; Heuvelink, G. B. M.; Leijnse, A. Spatial uncertainty analysis of pesticide leaching using a metamodel of PEARL. In Proceedings of Accuracy 2006 7th International Symposium on Spatial Accuracy Assessment in Natural Resources and Environmental Sciences; Caetano, M., Painho, M., Eds.; Instituto Geogra´fico Portugueˆs: Lisboa, Portugal; p 359-366.
Received for review June 16, 2007. Accepted July 31, 2007. ES0714639
