Journal of Environmental Quality 31:175-187 (2002)
© 2002 American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America
TECHNICAL REPORT
Heavy Metals in the Environment
Assessment of Uncertainty and Risk in Modeling Regional Heavy-Metal Accumulation in Agricultural Soils
A. Keller*,a,
K. C. Abbaspourb and
R. Schulina
a Swiss Federal Institute of Technology (ETH Zürich), Institute of Terrestrial Ecology, Grabenstr. 3, CH-8952 Schlieren, Switzerland
b Swiss Federal Institute for Environmental Science and Technology (EAWAG), Ueberlandstr. 133, P.O Box 611, CH-8600 Dübendorf, Switzerland
* Corresponding author (armin.keller{at}ito.umnw.ethz.ch)
Received for publication October 30, 2000.
 |
ABSTRACT
|
|---|
Present agricultural land use and atmospheric deposition may lead to heavy-metal accumulation rates in soils that may violate soil quality standards in the future. To undertake suitable preventive measures against heavy-metal enrichment, flux balances in agroecosystems and their uncertainties have to be assessed. For this reason we developed an empirical stochastic model, PROTERRA-S, that considers heavy-metal inputs through agricultural management as well as outputs by crop removal and leaching on a regional scale. In this manuscript we describe application of PROTERRA-S to the Sundgau region in Switzerland. Considering uncertainty in informational and natural variability, large variations of the aggregated regional cadmium and zinc balances were found, with standard deviations that were of the same order of magnitude as their average values. Uncertainty in the simulated net zinc flux originated mainly from uncertainty in the zinc concentrations of manure and crops and from uncertainty in atmospheric deposition of zinc. For cadmium, the main contribution to the total uncertainty came from uncertainty in crop concentration, regression functions to estimate Freundlich parameters, atmospheric deposition, and from spatial variation of soil pH and cation exchange capacity (CEC). For both zinc and cadmium, informational uncertainty in input data were large, indicating that significant uncertainty reduction could be achieved by additional data collection campaigns. A monetary risk value for the regional zinc accumulation rate in Sundgau was calculated to be on the order of 22 million Euro.
Abbreviations: CEC, cation exchange capacity RTU, root of uncertainty SRC, standardized regression coefficient UIF, uncertainty impact factor
 |
INTRODUCTION
|
|---|
THE SUSTAINABILITY of current land use in agroecosystems can be assessed with respect to heavy-metal accumulation in soils by balancing their input and output fluxes. Many studies point to the fact that, at the present rate, agricultural land use may not be sustainable because of heavy-metal accumulations. Some examples include studies by von Steiger and Baccini (1990), Reiner et al. (1996), and Moolenaar and Lexmond (1998) at the field and farm scales, and Andersson (1992) and Schütze and Nagel (1998) at the regional and national scales.
To provide meaningful information for taking suitable preventive measures against heavy-metal accumulation in soils, the uncertainty and reliability of the inputs and the outputs have to be known (von Steiger et al., 1998). There are two major sources of uncertainty in a modeling study: the uncertainty caused by parameters and the uncertainty due to model error. To quantify the model error in modeling heavy-metal balances is demanding, as it requires long-term monitoring data. Instead, we focus on the uncertainty caused by input parameters in this work. Parameter uncertainties in heavy-metal balances for soils with agricultural land use may be "natural" (i.e., arising from spatial and temporal variation that is inherent in natural processes) or "informational" (i.e., arising from measurement errors or incomplete or unreliable data). In addition, a common source of informational uncertainty in heavy-metal balance models is the uncertainty introduced by regression functions to derive model input parameters. Consequently, realistic modeling requires simulation procedures that account for these data characteristics (Abbaspour et al., 1998).
Aside from the limited amount of data that is usually available for heavy-metal balance analyses, the informational uncertainty may also arise as a result of scale differences. Soil and agricultural data necessary for heavy-metal balance studies often refer to different spatial scales. Metal output fluxes by leaching, erosion, or crop uptake are generally estimated at a plot or field scale. Other data, in particular those concerning agricultural management (e.g., livestock and crop production), are usually only available at a regional scale (Moolenaar, 1998). These data cannot easily be assigned to geographic coordinates, contrary to basic soil information that is obtained from soil maps in a GIS.
Various models have been developed to assess heavy-metal balances of agricultural soils. Few attempts were made to address the sensitivity of model results based on parameter uncertainty (Boekhold and van der Zee, 1991; Palm, 1994; Moolenaar and Lexmond, 1998). Tiktak et al. (1999) and de Vries and Bakker (1998) determined the contribution of uncertain model input parameters to the total uncertainty of the model results. However, these models did not consider metal inputs from land management and their uncertainties.
Uncertainty analysis includes the following tasks: (i) determination of the uncertainty in the model input parameters, (ii) propagation of parameter uncertainty through the model taking account of cross correlations between uncertain parameters, (iii) analysis of the contribution of each uncertain component to the total uncertainty of the results, and (iv) analyses of alternative strategies to reduce the uncertainties and associated risks of inadequate management decisions. In general, while risk imparted by informational uncertainties must be reduced by further detailed investigations, the uncertainties caused by natural variation must be maintained for a realistic representation.
To address these tasks in the assessment of heavy-metal balances at the regional scale for agricultural land, we developed the empirical stochastic model, PROTERRA-S. The model is an extended version of PROTERRA (von Steiger and Obrist, 1993) that was designed for assessing phosphorus, cadmium, zinc, and copper balances of agricultural land over regions of some 100 km2. The basic units of the balances are land-use systems, defined by livestock composition and cropping systems. Metal balances on the level of the land-use systems can be disaggregated to the level of individual crop types or aggregated to the regional level. To facilitate the performance of PROTERRA-S, the model could be linked to regional databases, such as agricultural statistics and soil information systems. The agricultural statistics of the land-use systems are coupled with the soil characteristics by a stochastic approach, considering each possible combination of soil type and fertilization management, both of which are related to livestock production and cultivated crops. Uncertain model input parameters can be treated as random variables, quantified with probability distributions, and propagated to model outputs using a Latin hypercube sampling scheme with correlated variables.
In this paper we demonstrate the applicability of the model and show how it deals with data and parameter uncertainty, illustrated by the regional aggregated zinc and cadmium balances for a case study.
 |
THEORETICAL BACKGROUND
|
|---|
Basic Equations
In the following, we briefly summarize the basic model equations. A detailed description of the PROTERRA-S model can be found in Keller et al. (2001). Heavy-metal balances calculated by PROTERRA-S are based on the following mass conservation equation proposed by van der Zee et al. (1990) and Boekhold and van der Zee (1991) for the plow layer of an agricultural soil:
 | [1] |
where M (g ha-1) is the total heavy-metal content, t is the time (yr), IAtm (g ha-1 yr-1) is the metal input flux by atmospheric deposition, IAgr (g ha-1 yr-1) is the sum of the metal input fluxes by agricultural activities, and kC (yr-1) and kL (yr-1) are the crop-uptake rate coefficient and the leaching-rate coefficient, respectively. The exponents m and n, are empirical parameters. The total metal content M is the sum of metal sorbed to the solid phase and metal dissolved in the soil solution:
 | [2] |
where
is the dry soil bulk density (kg m-3), z is the depth of the plow layer (m),
is the volumetric water content (L m-3), and ct (mg kg-1) and cs (mg L-1) are the sorbed and dissolved metal concentrations of the soil, respectively. To convert M from mg m2 to g ha-1 a factor of 10 was used in Eq. [2]. PROTERRA-S considers agricultural metal inputs resulting from application of animal manure (IMan), sewage sludge (ISe), commercial fertilizer (ICF), and pesticides (IPes), all expressed in g ha-1 yr-1:
 | [3] |
Metal inputs with fertilizer applications are linked to the following phosphorus (P) balance:
 | [4] |
Here, OP,Crop (g ha-1 yr-1) is the P export with crop harvest or grazing, and IP,tol (g ha-1 yr-1) represents an empirically determined rate of P surpluses in agroecosystems. IP,Man, IP,Se, and IP,CF are the P inputs with manure, sewage sludge, and commercial fertilizers (g ha-1 yr-1), respectively. OP,Crop, IP,tol, IP,Man, and IP,Se are obtained from agricultural statistics and guidelines, nutrient balances on experimental farms, and from bookkeeping of wastewater treatment plants. IP,CF is then calculated from these data according to Eq. [4], assuming that farmers usually follow the fertilization strategy recommended by the Swiss guidelines for water protection (Federal Office of Environment, Forests and Landscape, 1994). After the P balance is determined, metal inputs from fertilizer application are estimated from the P fluxes and from the ratios between metal and P concentrations of the fertilizers.
Assuming a linear relationship between the metal concentration in crop type i, cCrop,i (mg kg-1), and its concentration in soil, ct (see, for comparison, Gerritse et al., 1983 and Kuboi et al., 1986), the uptake-rate coefficient kC (yr-1) of crop type i in Eq. [1] is given by:
 | [5] |
where Yi (kg ha-1 yr-1) is the yield (dry matter) of crop type i, and the exponent m in Eq. [1] is set equal to one. A linear relationship between metal concentration in the soil and in the crop does not hold in general. However, there is no significant difference between the assumption of linearity and nonlinearity at metal concentrations at or near background as in our case, so the linear assumption was used. Leaching of heavy metals from the plow layer into the subsoil, OL (g ha-1 yr-1), is determined assuming equilibrium partitioning of metals between the adsorbed and dissolved phases described by Freundlich isotherms. As the dissolved amount of heavy metal in soil can usually be neglected, the leaching-rate coefficient is approximately given by:
 | [6] |
where qw (L m-2 yr-1) is the (Darcy) water flux at the lower boundary of the plow layer, and KF (L kg-1) and n (unitless) are the Freundlich parameters. These parameters are related to soil properties (van der Zee and van Riemsdijk, 1987; Buchter et al., 1989) and are derived from regression functions explained below.
Application of the Model
Metal balances are calculated in two steps. First, flux balances of heavy metals are simulated for individual land-use systems within the region of interest. Second, the input and output fluxes of these land-use systems are aggregated by their area fractions to estimate the average flux balance for the region. The land-use systems are the basic units for the areal balances, which are categories defined by the user of PROTERRA-S based on the available statistical databases. We classified the land-use systems by farm type (i.e., predominant livestock type and density) and by cropping systems (i.e., classes of crop types) (Keller et al., 2001).
Figure 1
illustrates the main components of the model, which are programmed in MATLAB Version 5.3 (Mathworks, 1999). These modules are linked to the software UNCSAM (Janssen et al., 1994), a package for efficient Monte Carlo sampling in combination with regression and correlation analysis to perform sensitivity and uncertainty analyses.
 |
MATERIALS AND METHODS
|
|---|
Available Databases
PROTERRA-S is linked to agricultural statistics and soil information systems that are available in Switzerland. Since these agricultural databases are updated every five years, the model is designed to estimate average heavy-metal mass balances for five-year periods. The model input parameters are obtained from the following data sources: (i) agricultural statistics, (ii) soil maps and soil surveys, (iii) studies of experimental farms (e.g., heavy-metal concentration in manure), and (iv) guidelines and recommendations for agricultural practices (e.g., fertilization management, application of sewage sludge).
As a case study, we applied the model to the Sundgau, a rural region in the northwest of Switzerland. This region extends from the Swiss Jura in the southeast with hills to 800 m to the Rhine Valley in the northwest. Sundgau has an area of about 95 km2 and includes 201 farms that cultivate about 36 km2 of arable land. Pastures and meadows account for 49% of this area, agricultural crops 46% (dominated by bread and fodder cereals), and special crops (vineyards, orchards) 6%. Dairy farms and arable farms cover about 88% of the total agricultural area. The data used in this case study refer to the period from 1992 to 1997. The complete set of model input parameters used in this case study is available at http://www.ito.umnw.ethz.ch/SoilProt/staff/kellera or can be obtained from the corresponding author.
Derivation of Model Parameters
The Freundlich parameter, KF in Eq. [6], was determined from GIS mapped soil properties using the relationship given by Elzinga et al. (1999):
 | [7] |
where bkf, bpH, and bCEC are empirical regression coefficients and
is the regression error assumed to be normally distributed with mean zero and variance
2
(Draper and Smith, 1981). Table 1 gives the values of the regression coefficients calculated based on a wide range of soils for cadmium and zinc. Given the calculated KF values, we estimated the metal concentration in soil solution based on the Freundlich equation (see below for the estimation of the Freundlich coefficient n), and compared these estimated metal concentrations in soil solution with measured NaNO3extractable concentrations (Swiss Agency for the Environment, Forests and Landscape, 2001) from a regional soil survey by Bono (1999). For cadmium, the estimated concentrations were in the same range as the measured ones, while for zinc the regression function given by Elzinga et al. (1999) overestimated the respective measured values on average by a factor of five. Therefore, we adjusted the intercept, bkf, of the regression function to fit the average of these measurements.
The empirical parameter n in Eq. [6] was estimated according to Buchter et al. (1989). They studied the relationship between soil properties and the magnitude of the Freundlich parameters for 11 soils differing greatly in characteristics. For cadmium and zinc, Buchter et al. (1989) derived a regression function of the form:
 | [8] |
where bn and bnpH are empirical coefficients. Taking estimated regression coefficients from the literature (Table 1), here denoted as
i and
k, the total variance of the predicted Freundlich coefficients,
2Y, is given by:
 | [9] |
where Cik is the covariance matrix of the p regression variables and
2
is the regression error. Random realizations of
were drawn with Monte Carlo simulations to estimate
2
, taking into account the coefficient of determination of the regression functions, R2, according to:
 | [10] |
For each regression function the error term was assigned to the regression constants, that is, bkf and bn were treated as a normally distributed random variables with mean bkf and variance
2
and with mean bn and variance
2
, respectively.
Uncertainty Analysis
In this case study, we initially considered 42 parameters as uncertain, thus we treated them as random variables (RV). A preliminary sensitivity analysis using stepwise backward regression (Zar, 1984) showed that uncertainty of about half of these random variables had negligible effect on the uncertainty of the model results. Therefore, 20 of the most sensitive parameters were retained as uncertain random variables. These are listed in Tables 2 and 3. Table 2 contains the sources of data from which the probability distributions given in Table 3 were inferred. It should be noted that the classification of parameter uncertainty into the categories of informational and natural in Table 2 indicates the predominant uncertainty type of the random variables. A sharp distinction between such categories, however, is not possible as each variable has a measure of both uncertainties. The selection of parameters, considered here as random variables, depends on the agricultural characteristics (e.g., the predominant crop types and livestock production of the Sundgau region). In other regions, estimates of cadmium and zinc fluxes may be sensitive to a different set of parameters because other types of crop and manure may dominate the uncertainty of these estimates.
Probability distributions for spatially distributed parameters (pH, CEC, and ct) were derived from 1:5000 GIS-based soil maps and from regional soil surveys. The regional uncertainty of atmospheric deposition was inferred from bulk deposition measurements available from various locations of the Swiss Monitoring Network for Air Pollutants excluding measurements within urban areas.
Correlations between uncertain input parameters may significantly affect the uncertainty of the model results. However, the amount of data usually available for heavy-metal balance studies does not support such analyses. We estimated correlations between random variables for which we had obvious indication that they were correlated, as shown in the last column of Table 3. The Spearman correlation coefficients listed in Table 3 were obtained from soil surveys (correlation between ct and cCrop, pH, CEC) and from the literature (correlation between ct and metal concentration of fertilizers).
To propagate the uncertainties in the parameters through the model we used the Latin hypercube sampling technique (McKay et al., 1979). After some preliminary runs with different sample sizes of up to 600, we noticed that the results above 300 were very consistent. Thus, we made 400 model runs with PROTERRA-S. Based on the definitions of the probability functions, all possible combinations of soil characteristics and agricultural land-use types were taken into account by this approach. We performed the Latin hypercube sampling with and without taking parameter correlations into account, to determine the effect of such correlations on the model outputs. We used the procedure described by Janssen et al. (1992) for the Latin hypercube sampling with correlated parameters.
Determination of Uncertainty Contributions
To quantify how the uncertainties in model input parameters affect the uncertainty of the model results, various approaches have been proposed in the literature (e.g., Iman and Helton, 1992; Janssen, 1994). We employed a multiple regression model to evaluate the effect of each random variable Xi (listed in Table 3) on the regional net cadmium and net zinc balances (Y), of the form:
 | [11] |
where b0 and bi are the regression coefficients,
j are the residual errors, p is the number of random variables, and N is the number of data points. The residuals are assumed to be normally distributed with zero mean and variance
2
. To obtain a relative sensitivity measure that is independent of the scale of the random variables, standardized regression coefficients, SRC, were calculated as (Janssen et al., 1994):
 | [12] |
where
Y is the standard deviation of Y and is the standard deviation of the variable Xi. The SRC measures the relative change of
Y for each parameter with respect to its standard deviation, due to a relative change in
Xi with respect to its standard deviation, while keeping all other model parameters constant. To account for the influence of correlation between input parameters on the model output uncertainty, Janssen et al. (1994) introduced the so-called root of uncertainty, RTU:
 | [13] |
where rXiY is the correlation coefficient between parameter Xi and Y, and rXiXk is the correlation matrix of the model input parameters. If no correlation between the model input parameters exists, the RTUi equals approximately the SRCi. The uncertainty measure RTU is only meaningful if the regression function shows a good linear approximation (i.e., the adjusted coefficient of determination R2adj
1). Therefore, the assumptions of the regression analysis must be rigorously checked with respect to R2adj as well as the residual errors. While the R2adj gives information about the approximation of the linear regression to the model results, the residuals reveal if assumptions for the regression errors are valid. We checked the residuals by inspecting normal probability plots (quantiles of residuals versus quantiles of the standard normal distribution) and the plot of model results versus regression estimates.
Because the estimation of RTU by the least squares minimization method is very sensitive to outlier values, we applied the less-sensitive maximum likelihood (ML) estimation method in combination with a Huber-weighting function for the determination of the regression coefficients bi that minimize the following quantity (Huber, 1981):
 | [14] |
where
(
j) is an arbitrary function of the residuals
j, and s is a scale parameter. The residuals
j depend on the estimated regression coefficients bi, and are assumed to be independently distributed and follow the same distribution. The scale parameter s can be estimated by the (robust) median of absolute values, sMAV (Huber, 1981):
 | [15] |
The derivation of the function
(
j) leads to the ML estimator function
for bi of the form:
 | [16] |
where X is the matrix of the random variables
and is a function that assigns weights, wj, to the standardized residuals
j/sMAV. There are various
functions to obtain robust estimates for the bi (Hoaglin et al., 1983). We used the following
function proposed by Huber (1981):
 | [17] |
 | [18] |
where c is a tuning constant, which was set to a value of five in this study. The ML approach can be solved iteratively. In each iteration the ML estimator of the coefficients b given by
are calculated according to:
 | [19] |
where XT is the transpose of the matrix X and W is a diagonal weighting matrix with diagonal elements given by:
 | [20] |
The
j and sMAV are recalculated in each iteration.
Risk Analysis
For our purposes, we define risk in monetary terms as the expected cost associated with the probability of failure at a given year:
 | [21] |
where R(t) is the risk in year t, Cf(t) is the cost (Euro) that would arise as the consequence of failure in year t, and Pf(t) is the probability of failure in year t. In this study we define failure as the occurrence of a heavy-metal accumulation rate in excess of a critical rate. This critical rate was calculated as the tolerable metal accumulation during 200 years that would lead to a metal concentration in the plow layer equal to the Swiss guideline values (Swiss Agency for the Environment, Forests and Landscape, 2001).
The probability of failure was estimated by employing the same Monte Carlo simulation procedure as described above. The probability of failure can be reduced by decreasing the part of uncertainty that is associated with the informational components by improving data processing and additional surveys.
 |
RESULTS AND DISCUSSION
|
|---|
Heavy-Metal Fluxes
Model simulations showed an average enrichment of cadmium and zinc in the agricultural soils of the Sundgau region during the period considered (Table 4). The regional cadmium balance was mainly determined by atmospheric deposition (67% of total input) and by crop removal (43% of total input) (Fig. 2)
. Leaching of cadmium was of minor importance because of the large sorption capacity and high soil pH values in the region. Cadmium inputs from commercial fertilizers were found to be relatively small due to the predominance of manure as fertilizer. About 60% of the total P input (21.8 kg P ha-1 yr-1) was provided by manure and about 40% by commercial fertilizers. The amount of P supplied with sewage sludge was almost negligible due to the small amounts applied within the region. It should be emphasized that even for extensively used meadows and pastures, where atmospheric deposition is the only relevant cadmium source, the regional cadmium balance indicates a slight accumulation in soil. In contrast to cadmium, the regional zinc balance was mainly governed by animal manure (64% of total input), followed by atmospheric deposition (31%) and crop uptake (31%) (Fig. 3)
.

View larger version (35K):
[in this window]
[in a new window]
|
Fig. 2. Model results for the regional cadmium balance. (a) Input and output fluxes, where IMan = inputs by animal manure; ISe, ICF, and IPes = inputs with sewage sludge, commercial fertilizer, and pesticides, respectively; IAtm = atmospheric deposition; OCrop = output by crop removal; and OL = leaching. The boxes give the interquartile range (IQR). The horizontal bar in each box indicates the median. The length of the whiskers is 1.5 x IQR. Values outside these limits are plotted as points. (b) Net flux, where xmed = median flux. (c) Contribution of the uncertainty sources to total variance of the net flux, where RTU = root of uncertainty. For symbols of uncertainty sources, see Appendix 1.
|
|

View larger version (38K):
[in this window]
[in a new window]
|
Fig. 3. Model results for the regional zinc balance. (a) Input and output fluxes, where IMan = inputs by animal manure; ISe, ICF, and IPes = inputs with sewage sludge, commercial fertilizer, and pesticides, respectively; IAtm = atmospheric deposition; OCrop = output by crop removal; and OL = leaching. The boxes give the interquartile range (IQR). The horizontal bar in each box indicates the median. The length of the whiskers is 1.5 x IQR. Values outside these limits are plotted as points. (b) Net flux, where xmed = median flux. (c) Contribution of the uncertainty sources to total variance of the net flux, where RTU = root of uncertainty. For symbols of uncertainty sources, see Appendix 1.
|
|
The net metal fluxes in the region of our study were in good agreement with the regional and national balance studies in European agroecosystems. An overview of these studies has been given by Moolenaar (1998). According to these studies, net cadmium accumulation in European agricultural soils varied between 0.2 and 3.7 g ha-1 yr-1 and net zinc accumulation between -61 and 1083 g ha-1 yr-1. The median of the net metal fluxes in our study (Table 4) was about in the middle of these ranges.
Sources of Uncertainties
Propagation of uncertainties of the random variables through the model resulted in large variances of the regional aggregated net cadmium and net zinc fluxes. The coefficient of variation (CV) was about 120% for cadmium and about 50% for zinc.
The regression model used in this study to quantify the contribution of the uncertainty sources explained about 98% of the variance of the simulated net zinc flux (Fig. 4)
. Thus, the linear approximation was appropriate in this case. Since there were few outliers that were identified by the Huber function, the maximum likelihood estimator for the regression coefficients, bi, was quite similar for zinc to the corresponding least squares estimation. For cadmium, however, less of the variation was explained by the linear regression model (R2adj = 0.73), and the residuals showed a slightly nonlinear structure (see Fig. 4). In particular, the regression model overestimated the negative net cadmium fluxes. Down-weighting of the observed outliers in the robust regression led to significant changes in some regression coefficients and hence, to changes in the corresponding uncertainty measure RTU.

View larger version (32K):
[in this window]
[in a new window]
|
Fig. 4. Control plots of the residual analyses: normal probability plot (left graphs) and plot of the model estimate versus the robust regression estimates (right graphs) for the regional net flux of (a) cadmium and (b) zinc. The observations plotted in grey refer to weighted residuals (Huber tuning constant c = 5).
|
|
Table 5 shows that the uncertainty of the regional cadmium balance considering parameter correlations resulted primarily from the variation of cadmium concentration in crops (RTU = 18%), from uncertainty in Freundlich regression coefficients bkf and bn (16.7%), from atmospheric deposition (16.2%), and from uncertainty in leaching rates due to spatial variation of soil pH and CEC (15.2%). With respect to metal export by plant uptake, variation of cadmium concentration in grass on extensively used meadows (32% of regional area) contributed most to the total uncertainty of the net cadmium flux (RTU = 7.4%). In summary, informational uncertainty explained about 44% and natural uncertainty about 56% of the total variation of the cadmium net flux. This analysis indicates that substantial reduction in cadmium flux uncertainty can be achieved by additional surveys and investigations.
View this table:
[in this window]
[in a new window]
|
Table 5. Uncertainty measure of the regional net fluxes neglecting (RTU-0) and considering (RTU-1) parameter correlation. For an explanation of symbols see Appendix 1.
|
|
The sources of uncertainty in the cadmium balance were in general agreement with those identified by Tiktak et al. (1999) in a study on metal accumulation in Dutch soils. Considering variation of fundamental soil variables (pH, organic matter, clay, water content, water flux), crop variables (concentration and yield), and uncertainty resulting from the Freundlich regression functions, the uncertainty in the latter explained most of the total variance in soils with low pH. In contrast, in soils with high pH, organic matter content contributed the most to the total uncertainty.
Comparable results were obtained by Palm (1994), who coupled a hydrological model with a model for sorption and plant uptake to assess cadmium behavior for a soil profile. Cadmium uptake by plants was assumed to be proportional to root water uptake. He found that the partitioning of cadmium between sorption, leaching, and plant uptake was highly sensitive to sorption isotherms, whereas root distribution of plants and hydrological soil properties influenced the partitioning of these fractions only to a limited degree. However, the uncertainties in cadmium inputs were not taken into account in these two studies. As our case study shows, cadmium and zinc accumulation in agricultural soils may also be quite sensitive to the uncertainties in atmospheric deposition and fertilization management.
Our findings are also consistent with the findings of Boekhold and van der Zee (1991), who used a deterministic transport model to simulate cadmium accumulation in soil, leaching, and uptake by barley. They reported that leaching was very sensitive to the input rate, water flow, and sorption parameters. With time, however, this sensitivity decreased. Cadmium uptake by barley was very small, but in time became increasingly sensitive to flow and sorption parameters. From their analysis, Boekhold and van der Zee (1991) concluded that the input flux is the most sensitive parameter for short-term predictions of cadmium accumulation in soil, whereas for long-term predictions accumulation is also sensitive to variations of flow and sorption parameters. In our case, uncertainty of input fluxes explained about 42% of the total uncertainty of the regional net cadmium flux.
The uncertainty of the regional net zinc flux was primarily influenced by variations of zinc concentration in manure (RTU = 33.5%) and in crops (17.7%) as well as variations in atmospheric zinc deposition (15.8%), when parameter correlations were considered (Table 5). The large contribution from animal manure to the total uncertainty was mainly due to highly uncertain zinc concentrations in pig manure (23.1%) and in the manure from cattle and dairy cows (7.2%). Similar to uncertainty in cadmium uptake, uncertainty arising from metal uptake by grass grown on extensively used meadows yielded the largest RTU value (5.2%). In addition, a significant part of the uncertainty in the zinc balance resulted from spatial variation of soil pH and CEC (10.5%) and zinc concentration in soils (9.4%). The uncertainty resulting from regression of the Freundlich isotherms was less important for zinc than for cadmium. Both bkf and bn contributed only about 6.2% to the total variation of the regional net zinc flux. To sum up, about 58% of the total variation of the net zinc flux resulted from informational uncertainties, while 42% was caused by natural uncertainties. As in the case of cadmium, there is also a large potential for uncertainty reduction of the regional zinc balance by further detailed surveys and investigations.
Variation of heavy-metal concentrations in crops and in soils led to large uncertainties in the estimates of the crop uptake rate coefficients (Table 6). The largest rate coefficients were obtained for maize silage and intensively used meadows. For cadmium, the CVs of the rate coefficients ranged between 50 and 72%, while for zinc, values between 33 and 45% were found.
In response to the errors in the regression parameters and in the spatial uncertainty of soil pH and CEC, the CV of the Freundlich parameter KF was 58% for cadmium and 37% for zinc. For the CV of the Freundlich parameter n, a value of 10% was obtained for both cadmium and zinc (Table 6). Consequently, these uncertainties led to large variations in the leaching rate coefficient, kL, for cadmium and zinc. Furthermore, the distributions of the simulated kL values were strongly skewed. The standard deviation was about one order of magnitude larger than the median in the case of cadmium and about six times larger than the median in the case of zinc.
The Influence of Parameter Correlation
Correlations between the random variables had little effect on the regional net cadmium and net zinc fluxes. For individual land-use systems, however, their influence was quite significant. For example, for dairy farms (pastures and meadows, 0.21.0 cattle manure units [CMU] ha-1 cattle) the standard deviation of the net zinc flux decreased from 220 to 198 g ha-1 yr-1 (median 378 g ha-1 yr-1) when parameter correlations were incorporated in the simulations (Fig. 5)
.

View larger version (35K):
[in this window]
[in a new window]
|
Fig. 5. Frequency distributions of simulated net zinc fluxes for both with and without considering input parameter correlations. The left plot is is the regional balance, the right plot is for a specific land-use system (dairy farm types).
|
|
Quantification of the uncertainty contributions was more sensitive to parameter correlation than the net fluxes. Considering correlation between model input parameters, the RTU values with correlation (RTU-1) deviated considerably from those without considering correlations (RTU-0; see Table 5). Hence, even weak correlations can have a substantial influence on the contribution of input variables on the total uncertainty of the model output. With correlations, the contribution of crop concentrations to the total uncertainty of the net cadmium flux increased significantly, while the influence of the uncertainty in the Freundlich regression functions decreased. With respect to the net zinc balance, correlations between model input parameters were relevant for the uncertainty contributions of atmospheric deposition and of the zinc concentrations of crops and manure. These results suggest that the correlations between uncertain input variables should be analyzed more intensively. Presently, however, due to limited data sets, such correlations are rarely investigated and usually neglected. One possible approach to dealing with this problem is to treat the correlation coefficients between model input parameters also as random variables.
Model Uncertainty
The analysis of uncertainty propagation yields informational and natural uncertainties, which are conditional on the validity of the model apart from the calibration of its parameters. An additional source of uncertainty is the error in the model itself. Model uncertainty may arise by neglecting relevant processes such as linear nonlinear adsorption, choice of inadequate process descriptions such as steady-state or equilibrium, or specification of incorrect boundary conditions. In our case, metal fluxes with soil erosion or metal inputs by compost were not considered, although they may play an important role for the heavy-metal balance of certain agricultural soils. Moreover, simplification of the metal leaching processes (i.e., metal partitioning between the adsorbed and dissolved phases) and the choice of flux assessment methods for metal inputs through animal manure may lead to model uncertainties. Regarding the leaching process, the Freundlich regression functions used in this study tend to overestimate the metal concentrations in soil solution (see, for comparison, Elzinga et al., 1999). Therefore, we recommend calibrating these isotherms with local measurements if balance analyses indicate that model outputs are sensitive to these data.
In general, model uncertainty is assessed by comparison of different models or by comparison of model predictions with independent measurements (de Vries et al., 1998; Heuvelink, 1998). The latter validation approach was performed by Tiktak et al. (1998). They tested the validity of their model by comparing retrospective model predictions using historical cadmium inputs with present heavy-metal concentrations in soils of the Netherlands. They found that model predictions deviated, in general, by less than a factor of two from measured cadmium concentrations in soil.
Unfortunately, the data available are very limited for such a direct validation of regional scale flux models. For PROTERRA-S, this would require long-term monitoring, not only of soil conditions such as pH, organic matter content, and contaminant concentrations, but also of metal fluxes entering and leaving the soil. Consequently, model development and application must be linked to soil monitoring networks. Integration of modeling efforts into existing monitoring programs is also in the interest of such programs, as only models allow local extrapolation of metal inputs and outputs assessed at specific monitoring sites to other areas.
Reduction of Uncertainty in Heavy-Metal Balances
To reduce the uncertainty in metal balances, the two major types of uncertainty considered in this study are distinguished. While it is desirable to reduce the informational uncertainty, the natural uncertainty of parameters must be maintained for models to be realistic. To evaluate the potential gain that could result from the reduction of uncertainty in different parameters, we calculated the so-called uncertainty impact factor (UIF), defined as the ratio of the RTU to the metal fluxes (% of total input). For the regional cadmium balance the largest UIF values were found for the parameters pH, CEC (UIF = 1.3), and bkf and bn (UIF = 1.4), indicating that the contributions of these sources to the total uncertainty of the regional net cadmium flux were larger than the contribution of the leaching flux to the total input flux (Table 5). In contrast, the uncertainty of atmospheric deposition, which accounted for 67% of the total cadmium input, contributed only 16% to the total variation of the net cadmium flux (UIF = 0.3).
This finding illustrates that the UIF can be used as an indicator to express the contribution of the random variables to the total uncertainty, regardless of the magnitude of the affected metal flux. In our case, the uncertainty of the net cadmium flux could be reduced considerably if regional measurements were available to establish a local regression function for the prediction of the Freundlich parameters. In contrast, uncertainty resulting from the spatial variation of the soil pH and CEC can hardly be reduced. Since the agricultural data are not site specific, the model takes all possible combinations of soil characteristics and agricultural land use within the region into account. Site-specific agricultural data would provide heavy-metal balances for individual soil type/land-use combinations, and that would improve the accuracy of the simulated net fluxes substantially.
Risk Analysis
For a rough estimation of the risk involved in the current metal accumulation trends, we compared the simulated net metal flux of the regional balance and of specific land-use types with the critical accumulation rate. This estimation assumes that the rate of metal accumulation based on the five-year analysis holds for the next 200 years. This extrapolation only reveals the continuation of the status quo and does not include any future measures that would reduce or increase the rate of accumulation. In the case of zinc, a critical rate of 0.46 mg kg-1 yr-1 was obtained based on our calculations (Swiss guide value: 150 mg kg-1 HNO3extractable zinc). For a conversion from flux (g ha-1 yr-1) to critical value (mg kg-1 yr-1) we used a soil depth of 0.2 m and a soil bulk density of 1300 kg m-3. The distribution of the initial heavy-metal concentration (Table 3) for the region was based on measured values from data sources reported in Table 2. Figure 6
shows the predicted probability of failure based on different land-use systems. More details about the definition of the land-use classifications can be found in Keller et al. (2001). The probability of failure for the arable and dairy farms is zero, the probability of failure for the animal husbandry farms about 82%, and for the regional averaged balance about 3%. In an alternative presentation, metal accumulation in soil as estimated with the regional balance is compared with the soil treshold concentration (i.e., the Swiss guide value; Fig. 7) .

View larger version (23K):
[in this window]
[in a new window]
|
Fig. 6. Cumulative distribution function (cdf) of simulated zinc accumulation rates compared with the critical rate, where ccrit is the tolerable zinc accumulation during 200 years that would lead to concentrations equal to the Swiss guide value of 150 mg kg-1.
|
|

View larger version (15K):
[in this window]
[in a new window]
|
Fig. 7. Median zinc concentration in soil, ct, over time (solid line), based on the estimated regional accumulation rate compared with the soil threshold concentration (Swiss guide value). The 95% confidence interval is also given (small dashed line).
|
|
Estimating costs of failure in environmental studies is a very controversial issue as it may encompass various diverse factors such as land price, cleanup costs, and adverse health effects. However, presenting the risk in monetary terms makes estimates more tangible and is an important tool for communicating abstract ideas to other participants. For the purpose of this paper, however, we will limit ourselves to a very simple scenario. We assume that the current market value of agricultural land in the Sundgau region would decrease from 30 to 10 Euro per m2 if the zinc accumulation rate exceeded the critical value. The cost of failure would then equal to 20 Euro per m2 as a result of land depreciation. Multiplied by the probabilities given above, the expected risk R, will amount to about 22 million Euro for the entire Sundgau region (36 km2), and for a small 10-ha animal husbandry farm to about 1.6 million Euro. Since a large part of these expected risks originated from informational uncertainties, it would be worthwhile to invest in additional surveys and investigations to improve the informational databases.
 |
SUMMARY AND CONCLUSIONS
|
|---|
We applied the stochastic model PROTERRA-S to the Sundgau region and assessed the aggregated cadmium and zinc balances in agricultural soils and their uncertainty sources in the regional scale. The results indicated that cadmium and zinc inputs on average exceeded the outputs and, thus, led to net accumulation of these metals in soil during 19921997. The regional balance, however, provided only a general trend of metal fluxes. To take preventive measures, more detailed analyses relating the land-use systems and individual soil types within the region would be more useful. The model is designed to estimate fluxes for five-year periods; therefore, it estimates accumulation trends only for these periods. For temporal predictions beyond that, a dynamic model based on a more sophisticated description of the cause and effect relationships between metal cycling in agroecosystems and soil processes is required.
Taking into account informational and natural uncertainties of the input parameters resulted in large uncertainties of the heavy-metal balances. The standard deviations of the regional net cadmium and net zinc fluxes were of the same order of magnitude as their average values. Using a maximum likelihood regression approach for uncertainty analysis, we quantified the contribution of uncertainty in model input parameters to the total variation of the net heavy-metal fluxes. Uncertainty in the simulated net zinc flux originated mainly from informational uncertainty in zinc concentrations of manure and crops and from natural uncertainty in atmospheric deposition of zinc. In the case of cadmium, primarily the latter two sources of uncertainty as well as the natural uncertainty of the soil pH, CEC, and the uncertainty in the regression function to estimate the Freundlich parameters were the sources for the large variation of the regional net cadmium flux.
The modeling results gave valuable insights for (i) suitable preventive strategies for specific inputs to reduce heavy-metal enrichment in soils and (ii) reducing uncertainty in balances by further investigations to determine highly sensitive parameters more accurately. In order to identify the potential of parameters for efficient uncertainty reduction in the estimated metal balances, we defined the uncertainty impact factor (UIF). A large potential for uncertainty reduction in the net cadmium fluxes was found in the Freundlich regression functions. To reduce model uncertainty, we concluded that long-term monitoring data of soil conditions, metal inputs and outputs, and agricultural practices are required. Such data were not available for this study and therefore model uncertainty had not been accounted for. To improve the reliability and precision of metal balances, we suggest that model development and application should be integrated into soil monitoring programs in the future.
A simple risk analysis was conducted in order to exemplify the potential benefits of uncertainty analyses in the soil pollution monitoring. The analysis showed large risks to be associated with the current trends in accumulation of heavy metals in soils. However, since these risks resulted to a great degree from informational uncertainties, substantial risk reduction could be achieved by determining sensitive model parameters more accurately.
Notation of Symbols
bkf, bnFreundlich regression coefficients (unitless)cmetal content (mg kg-1)cttotal metal concentration in soil (mg kg-1)cssoluble metal concentration in soil (mg L-1)Ccovariance matrixCECcation exchange capacity (mmolc kg-1)Cfcost of failure (Euro)CVcoefficient of variation (%)
regression residualsImetal input flux (g ha-1 yr-1)IPphosphorus input flux (g ha-1 yr-1)kCcrop-uptake rate coefficient (yr-1)kLleaching-rate coefficient (yr-1)KF, nFreundlich parameters (L kg-1 if n = 1)mempirical parameter (unitless)Mmetal content in soil (g ha-1)Nnumber of data pointsOmetal output flux (g ha-1 yr-1)OPphosphorus ouput flux (g ha-1 yr-1)Pfprobability of failurepHsoil pH (unitless)qwwater flux (L m-2 yr-1)
volumetric water content (L m-3)
dry soil bulk density (kg m-3)rSpearman correlation coefficientR2coefficient of determinationRrisk (%)RTUroot of uncertainty (%)
2variance
standard deviationSRCstandardized regression coefficient (%)ttime (yr)UIFuncertainty impact factor (unitless)Xrandom variableYdry matter crop yield (kg ha-1 yr-1)zdepth of plow layer (m)
Process Indices
Agragricultural activitiesManmanureSesewage sludgeCFcommercial fertilizersPespesticidesAtmatmospheric depositionCropcropsLleachingIPtolphosphorus surplus
 |
ACKNOWLEDGMENTS
|
|---|
We are grateful to P.H.M. Janssen, National Institute for Public Health and Environmental Protection (RIVM), Bilthoven, the Netherlands, who provided the UNCSAM software package for efficient Monte Carlo simulations. We wish to thank the Swiss Federal Statistical Office, Bern, for providing the agricultural statistics, and the reviewers for their helpful comments.
 |
REFERENCES
|
|---|
- Abbaspour, K.C., R. Schulin, M.Th. van Genuchten, and E. Schläppi. 1998. Procedures for uncertainty analysis applied to a landfill leachate plume. Ground Water 36:874883.
- Andersson, A. 1992. Trace elements in agricultural soilsFluxes, balances and background values. Report no. 4077. Swedish Environ. Protection Agency, 75007 Uppsala, Sweden.
- Boekhold, A.E., and S.E.A.T.M. van der Zee. 1991. Long-term effects of soil heterogeneity on cadmium behavior in soil. J. Contam. Hydrol. 7:371390.
- Bono, R. 1999. Schwermetalle in Baselbieter Böden. Bau- und Umweltschutzdirektion Kanton Basel-Landschaft, Bodenschutzfachstelle, 4410 Liestal, Switzerland.
- Buchter, B., B. Davidoff, M.C. Amacher, C. Hinz, I.K. Iskandar, and H.M. Selim. 1989. Correlation of Freundlich Kd and n retention parameters with soils and elements. Soil Sci. 148:370379.
- De Vries, W., and D.J. Bakker. 1998. Manual for calculating critical loads of heavy metals for terrestrial ecosystemsGuidelines for critical limits, calculation methods and input data. Report 166. DLO Winand Staring Centre, Wageningen, the Netherlands.
- De Vries, W., J. Kros, C. van der Salm, J.E. Groenenberg, and G.J. Reinds. 1998. The use of upscaling procedures in the application of soil acidification models at different spatial scales. Nutr. Cycling Agroecosyst. 50:223236.
- Draper, N.R., and H. Smith. 1981. Applied regression analysis. 2nd ed. John Wiley & Sons, New York.
- Elzinga, E.J., J.J.M. van Grinsven, and F.A. Swartjes. 1999. General purpose Freundlich isotherms for cadmium, copper and zinc in soils. Eur. J. Soil Sci. 50:139149.
- Federal Office of Environment, Forests and Landscape. 1994. Guidelines for water protection in agricultureSubject: Farm manure. Federal Office for Agriculture (FOA), Federal Office of Environment, Forests and Landscape (FOEFL), 3003-Bern, Switzerland.
- Gerritse, R.G., W. van Driel, K.W. Smilde, and B. van Luit. 1983. Uptake of heavy metals by crops in relation to their concentration in the soil solution. Plant Soil 75:393404.
- Heuvelink, G.B.M. 1998. Uncertainty analysis in environmental modelling under a change of spatial scale. Nutr. Cycling Agroecosyst. 50:255264.
- Hoaglin, D.C., F. Mosteller, and J.W. Tukey. 1983. Understanding robust and exploratory data analysis. John Wiley & Sons, New York.
- Huber, J.P. 1981. Robust statistics. John Wiley & Sons, New York.
- Iman, R.L., and J.C. Helton. 1992. An investigation of uncertainty and sensitivity analysis techniques for computer models. Risk Analysis 8:7190.
- Janssen, P.H.M. 1994. Assessing sensitivities and uncertainties in models: A critical evaluation. p. 344361. In J. Grasman and G. van Straten (ed.) Predictability and nonlinear modelling in natural science and economics. Kluwer Academic Publ., Dordrecht, the Netherlands.
- Janssen, P.H.M., P. Heuberger, and O. Klepper. 1994. UNSCAM: A tool for automating sensitivity and uncertainty analysis. Environ. Software 9:111.
- Janssen, P.H.M., P. Heuberger, and R. Sanders. 1992. UNCSAM: A software package for sensitivity and uncertainty analysisManual. Report no. 959101004. Natl. Inst. of Public Health and Environ. Protection, Bilthoven, the Netherlands.
- Keller, A., B. von Steiger, S.E.A.T.M. van der Zee, and R. Schulin. 2001. A stochastic empirical model for regional heavy metal balances in agroecosystems. J. Environ. Qual. 30:19761989.[Abstract/Free Full Text]
- Kuboi, T., A. Noguchi, and J. Yazaki. 1986. Family-dependent cadmium accumulation characteristics in higher plants. Plant Soil 92: 405415.
- Mathworks. 1999. MATLAB user's manual. Version 5.3. The MathWorks, Inc., Natick, MA.
- McKay, D.M., R.J. Beckmann, and W.J. Conover. 1979. A comparison of three methods selecting values of input variables in the analysis of output from a computer code. Technometrics 21:239245.
- Moolenaar, S.W. 1998. Sustainable management of heavy metals in agro-ecosystems. Ph.D. thesis. Agric. Univ. of Wageningen, Wageningen, the Netherlands.
- Moolenaar, S.W., and Th.M. Lexmond. 1998. Heavy metal balances of agro-ecosystems in the Netherlands. Neth. J. Agric. Sci. 46:171192.
- Palm, V. 1994. A model for sorption, flux and plant uptake of cadmium in a soil profile: Model structure and sensitivity analysis. Water Air Soil Pollut. 77:169190.
- Reiner, I., C. Lampert, M. Piterkova, and P.H. Brunner. 1996. Stoffbilanzen landwirtschaftlicher Böden von ausgewählten Betriebstypen bei Verwendung von Klärschlamm und Kompost. BKK2-Endbericht. TU Wien. Institut für Wassergüte und Abfallwirtschaft (AWS)., Wien, Austria.
- Swiss Agency for the Environment, Forests and Landscape. 2001. Commentary on the Ordinance of 1 July 1998 relating to impacts on the soil. Swiss Agency for the Environment, Forests and Landscape (SAEFL), 3003 Bern, Switzerland.
- Schütze, G., and H.D. Nagel. 1998. Kriterien für die Erarbeitung von Immssionsminderungszielen zum Schutz der Böden und Abschätzung der langfristigen räumlichen Auswirkungen anthropogener Stoffeinträge. UBA Texte 19/98, Forschungsbericht 204 02 825. Umweltbundesamt Berlin, Germany.
- Tiktak, A., R. Alkemade, H. van Grinsven, and K. Makaske. 1998. Modelling cadmium accumulation on a regional scale in the Netherlands. Nutr. Cycling Agroecosyst. 50:209222.
- Tiktak, A., A. Leijnse, and H. Vissenberg. 1999. Uncertainty in a regional-scale assessment of cadmium accumulation in the Netherlands. J. Environ. Qual. 28:461470.[Abstract/Free Full Text]
- Van der Zee, S.E.A.T.M., H.N.M. Ferdinandus, and A.E. Boekhold. 1990. Long-term effects of fertilization and diffuse deposition of heavy metals on soil and crop quality. p. 323326. In M.L. van Beusichem (ed.) Plant nutritionPhysiology and applications. Proc. of the 11th Int. Plant Nutrition Colloquium, Wageningen, the Netherlands. August 1989. Kluwer Academic Publ., Dordrecht, the Netherlands.
- Van der Zee, S.E.A.T.M., and W.H. van Riemsdijk. 1987. Transport of reactive solute in spatially variable soil systems. Water Resour. Res. 23:20592069.
- Von Steiger, B., and P. Baccini. 1990. Regionale Stoffbilanzierung von landwirtschaftlichen Böden mit messbarem Ein- und Austrag. Nationales Forschungsprogramm "Boden". Rep. no. 38. 3097 Liebefeld-Bern, Switzerland.
- Von Steiger, B., A. Keller, and R. Schulin. 1998. Regional mass flux balancing for controlling gentle soil remediation operations. Nutr. Cycling Agroecosyst. 50:303306.
- Von Steiger, B., and J. Obrist. 1993. Available databases for regional mass balances in agricultural land. p. 3546. In R. Schulin et al. (ed.) Soil monitoringEarly detection and surveying of soil contamination and degradation. Birkhäuser Verlag, Basel, Switzerland.
- Zar, H.J. 1984. Biostatistical analysis. Prentice Hall, New Jersey.
This article has been cited by other articles:

|
 |

|
 |
 
Z. Miao, M. Trevisan, E. Capri, L. Padovani, and A. A. M. Del Re
Uncertainty Assessment of the Model RICEWQ in Northern Italy
J. Environ. Qual.,
November 1, 2004;
33(6):
2217 - 2228.
[Abstract]
[Full Text]
[PDF]
|
 |
|