Published online 23 August 2007
Published in Vadose Zone J 6:668-678 (2007)
DOI: 10.2136/vzj2006.0148
© 2007 Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
ORIGINAL RESEARCH
Quantifying the Influence of Uncertainty and Variability on Groundwater Risk Assessment for Trace Elements
Sven Altfeldera,*,
Wilhelmus H. M. Duijnisvelda,
Thilo Streckb,
Gunnar Meyenburga and
Jens Utermanna
a Federal Institute for Geosciences and Natural Resources, Stilleweg 2, D-30655 Hannover, Germany
b Dep. of Soil Science and Land Evaluation, Univ. of Hohenheim, D-70593 Stuttgart, Germany
* Corresponding author (sven.altfelder{at}bgr.de).
All rights reserved. No part of this periodical may be reproduced or transmitted in any form or by any means, electronic or mechanical, including photocopying, recording, or any information storage and retrieval system, without permission in writing from the publisher.
Received 29 September 2006.
 |
ABSTRACT
|
|---|
The application of leaching models to predict field-scale heavy metal transport has been successful at several sites. Past studies involved site-specific sorption experiments to quantify the leaching process and to investigate its spatial variability. Uncertainty due to lack of knowledge was frequently ignored. Here we present a leaching model that is based on an extended Freundlich equation to describe sorption. The equation is derived on a nationwide scale and is applicable to arbitrary sites in Germany. Instead of relying on site-specific sorption experiments, it only requires information on the spatial variability of pH, organic carbon (Corg), and clay content. In addition to accounting for spatial variability, the model also considers the major sources of uncertainty influencing a leaching prediction. These sources involve uncertainty in the spatial distribution of pH, Corg, and clay content, as well as uncertainty of the extended Freundlich equation. The modeling strategy is based on a parallel soil column approach, in which a deterministic transport model is combined with a two-dimensional Monte Carlo method. The model's ability to predict downward movement is tested on two sites in Germany—a wastewater irrigation area and an area in the vicinity of a metal smelter—where cadmium leaching has been modeled previously using extended Freundlich equations derived from local sorption data. Without any parameter fitting involved, the predicted field-averaged cadmium profiles agree quite well with measured profiles in both cases.
Abbreviations: 2-D, two-dimensional • CDE, convection dispersion equation • CDF, cumulative density function • pdf, probability density functions • SCM, stochastic convective model.
|
INTRODUCTION
|
|---|
Environmental situations in which materials with low contamination are disposed of, or where soils are contaminated unintentionally, pose a risk of groundwater contamination. Trace elements that are frequently encountered in these materials are subject to leaching processes that should be assessed. We present a modeling strategy adapted to these elements that allows the prediction of their leaching potential in situations where information is scarce. In Germany, the need to provide such tools has been articulated as part of soil legislation passed in recent years (Federal Republic of Germany, 1998, 1999).
Most trace elements have a strong affinity to the solid phase of soils. Sorption is therefore a key process in modeling the transport of these elements. Using specific isotherm models such as Freundlich or Langmuir models to describe trace element sorption to soil is often not justified because detailed information on site-specific sorption behavior is either not available or its determination is too time consuming. To overcome this situation, a semi-empirical pedotransfer function may be applied to relate sorption or transport parameters to easily measurable soil properties. The pedotransfer function used here is an extended Freundlich isotherm relating trace element sorption to standard soil properties such as pH, clay content, or Corg (van der Zee and van Reimsdijk, 1987; Streck, 1993; Horn et al., 2004, Römkens et al., 2004). It belongs to a set of pedotransfer functions given by Utermann et al. (2005), who described the derivation of generally applicable extended Freundlich isotherms for Cd, Cr, Co, Cu, Mo, Ni, Pb, Sb, Tl, and Zn on a nationwide scale in Germany. The extended Freundlich isotherm has the advantage that by relating easily measurable soil properties to trace element sorption, the spatial variability of sorption can be incorporated easily in a leaching risk assessment. A disadvantage of the extended isotherm model is that it contributes considerably to uncertainty when used to describe sorption in a solute transport model, and it is less accurate than isotherm models such as the Freundlich or Langmuir models. Another nonnegligible source of uncertainty related to the use of the extended isotherm is the sampling uncertainty associated with the soil sampling strategy used to investigate the spatial variation of soil properties such as pH, clay content, and Corg at a given site.
Some solute-transport concepts that may be applied to field-scale transport problems with spatially variable soil properties are inappropriate with regard to the characteristics of field-scale solute flow, while others are too demanding with regard to data and computational requirements in a practical risk assessment. Applying the convection dispersion equation (CDE) at the field scale is only valid at travel times and distances that are sufficiently long (Jury et al., 1990), usually longer than encountered in reality. The stochastic convective model (SCM) is a better approximation of the physics in the initial phase of solute movement when transverse mixing is still minor. In theory, it converges toward the convective dispersive model asymptotically for linear sorption and longer transport distances (Jury and Roth, 1990, ch. 7). The SCM is often used for relatively mobile solutes, and although it is difficult to apply to layered soils, the SCM has been quite successful in the past (Simmons, 1982; Butters and Jury, 1989; Jury and Scotter, 1994).
For field- and regional-scale problems, van Dam et al. (2004) suggested that the SCM is in principle superior to the CDE but suffers from numerous practical limitations. Instead, a model of noninteracting parallel columns can be used (Dagan and Bresler, 1979) that may or may not involve the CDE at the column scale. In practice, the parallel columns model is similar to the stochastic-convective model but is more flexible (de Rooij, personal communication, 2005). Relatively recent developments such as the fractional advection dispersion equation or continuous time random walk methods have not yet been developed to a stage that allow a routine application to practical transport problems in the unsaturated zone (van Dam et al., 2004).
The parallel soil column model is often combined with Monte Carlo methods (van der Zee and van Reimsdijk, 1987; van der Zee and Boesten, 1991; Streck and Richter, 1997b). Basically, a certain number of simulations, in which each simulation is run deterministically with selected input variables randomly sampled from their respective distributions, yields probability density functions of soil concentrations as an output. Even though this approach is neither stochastic convective nor convective dispersive, it converges toward the SCM when the local column scale dispersion is small.
Our modeling approach is based on the parallel column model, with the investigated site divided into individual soil columns. Solute transport in each column is modeled with the CDE, taking into account sorption. Sorption coefficients are calculated with the extended Freundlich equation mentioned above. Uncertainty of the extended Freundlich equation, as well as spatial variability and sampling uncertainty, is explicitly considered using a two-dimensional (2-D) Monte Carlo technique (Cullen and Frey, 1999).
A guidance document on ecological risk assessment by the USEPA (2001) states that the different sources responsible for output variance should be distinguished within risk assessments; the 2-D Monte Carlo technique provides the capability of separating the effects of uncertainty and spatial variability.
Variability arises from true heterogeneity or variation in characteristics of the environment and/or receptors. Uncertainty, on the other hand, arises from incomplete knowledge of the world and represents lack of knowledge about certain factors. If necessary, variance due to uncertainty about the spatial distribution of soil properties may be decreased by a more detailed investigation of the site. Variance due to variability of soil characteristics, however, is an intrinsic property of a site and therefore is fixed.
Results of the 2-D Monte Carlo simulation are spatially averaged to yield a distribution of site- or field-averaged concentrations that reflects the uncertainty of the forecast methodology. To test the model's capability for risk assessment, we applied it to two case studies dealing with Cd transport at the field scale.
|
Theory
|
|---|
Deterministic Model Theory
Before describing the stochastic part of the chosen approach, we give a brief description of the deterministic transport model describing solute transport in an individual column of the parallel soil column model. The transport model was originally developed by Streck and Richter (1997b) and later modified by Ingwersen (2001). These authors successfully applied it in combination with a parallel soil column approach to predict downward Cd movement at the field and regional scale. The following text is a short introduction to the theory. More details can be found in Ingwersen (2001).
Water Transport
Since all investigated trace elements feature strong sorption and no degradation in soils, solute transport can be simulated assuming a stationary soil water regime (Swartjes, 1990; Stöfen, 2005). The one-dimensional local water balance is (Streck, 1993)
 | [1] |
where q is the water flux density (m d–1), z is depth (m), and W(z) is the root water uptake. Assuming that root water uptake is proportional to the relative root density (RRD = 1 – e–
z) as a function of z and that evapotranspiration (ET in [mm yr–1]) is completely taken up by the roots, one may write after integration of Eq. [1]:
 | [2] |
where N (mm yr–1) is precipitation, ET (mm yr–1) is the actual evapotranspiration, and is a parameter describing root density distribution with depth. The factor on the right side of Eq. [2] converts (mm yr–1) to (m d–1). If the source of contamination lies below the root zone, water uptake by roots becomes unimportant and Eq. [2] simplifies to
 | [3] |
Solute Transport
One-dimensional transport of sorbing and nondegradable contaminants in the model is described with the convection-dispersion equation:
 | [4] |
where C is the dissolved concentration (µg L–1), Ds = 
q is the apparent dispersion coefficient (m2 d–1) with proportionality constant (m) representing the dispersion length, S is the sorbed concentration (µg kg–1), 
is the volumetric water content, 
is the soil bulk density (kg L–1), and t is time (d). The sorbed amount, S, is assumed to be a function of C and soil properties like pH, Corg, and clay content. The derivation of this relationship for several trace elements from nationwide data in Germany is described in Utermann et al. (2005). More details, especially for Cd, are outlined below.
The initial condition is
 | [5] |
where S0 is equal to the natural background sorbed concentration for an uncontaminated soil. The boundary condition at the upper boundary is a flux-averaged concentration of the infiltrating source, Cinf (µg l–1) so that
 | [6] |
At the lower boundary, we assume
 | [7] |
where G is the depth of the groundwater table (m). The CDE (Eq. [4]) is approximated numerically using the Crank–Nicholson implicit–explicit finite difference scheme (Smith, 1985) with a grid resolution of 1 cm and a time step of 10 d (Streck, 1993). The nonlinear system of equations is solved using the Newton–Raphson method (Press et al., 1988). More details regarding the numerical procedures may be found in Streck (1993).
Stochastic Model Theory
Incorporating Variability and Uncertainty
Uncertainty and variability are considered in the framework of a parallel soil column model. The use of this model to describe trace element transport at the field scale is possible based on the assumption that transverse solute mixing can be neglected. If the transport distance is short compared to the spatial variation of the variables governing sorption and transport, this assumption is likely to be fulfilled (Dagan and Bresler, 1979). Unlike its conventional form, the parallel soil column model used here does not consider the stochastic distribution of local pore water velocities. Instead, only the sorption properties are assumed to vary in space. The choice is based on the following reasoning. First, according to Streck and Piehler (1998), the mean travel time of strongly sorbing solutes is insensitive to variations in the water content . Hence, instead of assessing the probability density functions (pdf) of v (e.g., by means of a tracer experiment), modeling of the transport of strongly sorbing solutes can be based on the pdf of the water flux density q and assuming constant water content. Second, the variation of the water flux density pdf is usually small compared to the variation of soil sorption properties. This was illustrated by Streck and Richter (1997b), who showed that ignoring the spatial variation of q had only a minor effect on the prediction of Cd transport at the field scale. As will be shown later, if uncertainty and variability are considered, the Freundlich sorption coefficient may vary by several orders of magnitude.
Preferential flow that might cause the rapid transport of small solute fractions is not relevant in the present case. Soil properties relevant for solute movement are varied over the ensemble of parallel soil columns to account for spatial variability and sampling uncertainty of these properties. An additional important source of uncertainty, described in more detail in the Appendix, is the limited precision of the extended Freundlich equation that relates soil parameters like pH, Corg, and clay content to sorption.
Variability and uncertainty in the model input are explicitly considered by assigning distributions to the respective model parameters to account for the two sources of input variance. These parameter distributions are then propagated through the deterministic model within a Monte Carlo approach, leading to an output that is itself a distribution that may be evaluated with regard to threshold values.
To discriminate what effect input variability and uncertainty have on the model's output distribution, the Monte Carlo method can be extended to two statistical dimensions (2-D Monte Carlo). The two dimensionality of this so called "double looping" method (Cullen and Frey, 1999) is achieved by
- running the model at k discrete locations of the investigated site (variation of input parameters in the statistical dimension "variability"),
- running the model l times for each of the k discrete locations (variation of input parameters in the statistical dimension "uncertainty").
Step 1 leads to a set of k local results reflecting one realization of site variability, while step 2 leads to l sets with k local results where the ensemble of sets reflects uncertainty in the knowledge of variability. The k local results at a certain depth obtained from step 1 may be displayed as a cumulative density function (CDF) visualizing variance due to variability. Step 2 allows the calculation of an ensemble of l CDFs displaying uncertainty in the knowledge on spatial variation. Another advantage of this method is that it allows the calculation of field- or site-averaged results by averaging only in dimension "variability" without losing information on uncertainty. This is useful as it is expected that decision making with regard to threshold values will be based on field-averaged concentrations instead of local ones. These concentrations may as well be displayed in form of a CDF, visualizing uncertainty of the field-averaged result.
Column Generation
The risk assessment takes place at the field scale. Before a transport prediction is performed within a risk-assessment framework, it is assumed that the unsaturated zone of an investigated site has been characterized with regard to the soil properties governing transport. As outlined above, only those properties that control sorption are considered for highly sorbing trace elements. The result of this characterization will be depth-dependent distributions of each property based on a random sample derived from soil profile sampling. Our modeling approach considers normal or lognormal distributions with given mean and standard deviation. Further, they may be truncated to minimum and maximum values. In addition to its depth-dependent distribution, a property may be correlated with itself over depth or with other properties. Correlation is characterized by the correlation matrix. Distribution parameters and correlations are uncertain because they are based on sampling. While uncertainty in the correlation matrix is neglected, bootstrap sampling based on the property-specific parametric distribution and sample size (Cullen and Frey, 1999) is used to estimate uncertainty in the mean and standard deviation of the distributions.
Random columns are generated with soil properties distributed according to the type of distribution, the minimum and maximum, the correlation matrix, and uncertainty in the mean and standard deviation. Random number generation is done using a JAVA routine based on the method of mixed linear congruential random number generation (Knuth, 1981). Values exceeding minimum or maximum are replaced by resampling. The Cholesky transform of the correlation matrix is used to correlate the multivariate normal random numbers. After correlating the random numbers, remaining values exceeding minimum or maximum are corrected again, this time by setting them to the respective minimum or maximum value to keep shifts in correlations at a minimum. To retain correlation and distribution characteristics of the multivariate random sample, the limiting minimum–maximum values must not be too close to the mean.
The layerwise generation of the soil properties in the random columns is achieved in the following manner. First, based on the statistics of the layer specific sampling distribution of each property (i.e., mean, standard deviation, sample size), a mean and standard deviation pair is generated via bootstrap sampling for each soil layer. Second, the set of layer-specific mean and standard deviation pairs are combined with the corresponding minimum and maximum of the layer-specific sampling distribution and the correlation matrix to generate k random soil columns representing one realization of the soil variability at the investigated site (first dimension of the 2-D Monte Carlo Method). This procedure is repeated l times, each time with a new set of mean and standard deviation pairs obtained through bootstrap sampling, to quantify uncertainty (second dimension of the 2-D Monte Carlo Method).
Figure 1 shows CDFs of two soil properties, pH and Corg, in a single soil layer of a set of random columns generated in this manner. The distribution data used to generate these two CDFs are from a study on Cd transport at a wastewater irrigation site (Streck, 1993). This site is one of two sites where the model approach was tested. Model testing is described in detail in a later section. It can be seen that in this specific case, the variability of pH and Corg (a single CDF) has a much larger influence on the overall variance than uncertainty (confidence intervals of the ensemble of CDFs).
View larger version (32K):
[in this window]
[in a new window]
|
FIG. 1. Generated cumulative density functions (CDF) of pH and Corg with uncertainties in a single soil layer of a set of random columns.
|
|
Generation of Freundlich Coefficients
Antilog and rearrangement of the extended Freundlich equation (Eq. [A2]) derived in the Appendix leads to
 | [8] |
where K* is the intrinsic Freundlich coefficient, Clay is clay content, a, b, and c are regression coefficients, n is the Freundlich exponent accounting for sorption nonlinearity, and 
is the residual error remaining after regression. Clay and Corg are in weight percentage. Knowing the empirical regression coefficients, the spatial variability and uncertainty of pH, Corg, and clay, and the residual error distributions of the on-site error 1 and the between-sites error 2 (see Appendix for details), Freundlich coefficients K may be calculated for each column and layer by using Eq. [8]. The empirical regression coefficients of the extended Cd–Freundlich equation, as well as the parameters of the residual error distributions, are given in Table 1.
View this table:
[in this window]
[in a new window]
|
TABLE 1. Empirical regression coefficients and residual error distribution parameters of the extended Cd–Freundlich equation.
|
|
An example distribution of Freundlich coefficients K calculated from Eq. [8] for a single soil layer (the same layer used to draw Fig. 1) is given in Fig. 2. The distribution of K combines the uncertainty of all input parameters. Comparing Fig. A5 in the Appendix and Fig. 2 reveals that the residual error (= 1 + 2) dominates the overall variance and uncertainty of K. While has a span of approximately 2.5 orders of magnitude on the log scale, K has a span of three orders of magnitude. Only about one-sixth of the overall variance of K is related to variability and uncertainty of soil properties (limited here to pH and Corg; see Case Study 1), while the rest results from the uncertainty of the extended Freundlich equation. If the dominant influence of the extended Freundlich equation on the variance and uncertainty of K can be confirmed for additional sites, a decision to ignore spatial variation and sampling uncertainty of soil properties governing sorption may be justified. This would be helpful in situations where depth-specific soil profile samples are not available, which is not uncommon for many contaminated sites. In this case, field averages of the soil properties that govern sorption may suffice for transport modeling when a non-site-specific pedotransfer functions is used to predict sorption. This would simplify the risk assessment of trace elements with regard to the necessary site-specific measurements. However, the 2-D Monte Carlo approach would still be needed to allow the separate treatment of the on-site error 1 and the between-sites error 2.
View larger version (16K):
[in this window]
[in a new window]
|
FIG. 2. Cumulative density function (CDF) of Freundlich coefficient K (log scale, left panel; carteisian, right panel) in one soil layer of the random soil columns with uncertainty (same layer as in Fig. 1).
|
|
Modeling Solute Transport
Deterministic transport model simulations are run for each soil column of the parallel soil column model, using the Freundlich coefficients calculated for the column layers and appropriate boundary conditions. The Freundlich coefficients integrate variability and uncertainty of pH, Corg, and clay, as well as the residual error of the extended Freundlich equation according to Eq. [8]. Local column scale dispersion, if not available for the investigated site, may be derived from literature (for a review of dispersivity length, see Vanderborght and Vereecken, 2007). Column scale dispersion has no significant effect on trace element transport at the field scale because at this scale, solute spreading is caused to a high degree by the spatial variability of sorption (Streck and Piehler, 1998). The water content of the horizons is set to field capacity.
Figure 3 shows the procedure for aggregating the transport simulation results from the local column scale to the field scale. At first, concentration depth profiles are simulated in k parallel soil columns. The k simulations are repeated l times to account for uncertainty in the spatial variation represented by the confidence band of the CDFs in Fig. 2.
View larger version (47K):
[in this window]
[in a new window]
|
FIG. 3. Aggregation scheme for quantifying the distribution of field-averaged concentration depth profiles from local profiles.
|
|
The concentration depth profiles belonging to each of the l sets can be averaged for each set to obtain l field-averaged depth profiles of the contaminant. Variation of these profiles reflects the remaining uncertainty with which contaminant movement at the field scale can be predicted (see upper four panels in the right column of Fig. 3). Distribution parameters like median or certain percentiles can also be calculated for the ensemble of field-averaged profiles to quantify this uncertainty (bottom right panel in Fig. 3).
|
Case Studies
|
|---|
The displacement of trace elements in soil is a very slow process (Dowdy and Volk, 1983). For Pb, Maskall et al. (1995, 1996) reported migration rates at historical smelter sites of 0.72 to 0.75 cm yr–1 in sandy materials and 0.07 to 0.31 cm yr–1 in clays. Because measurable transport distances develop only within decades (Swartjes, 1990), it is difficult to measure transport in field experiments. Sites that have a long history of contamination are an ideal test ground for the modeling approach outlined above, presuming that the contamination history can be reconstructed (Hawkins et al., 1995; Streck and Richter, 1997b; Ingwersen and Streck, 2006). At two Cd-contaminated sites where these conditions are met, the present depth distribution of Cd in soil is predicted from historical data to test the model approach. A good agreement between the predicted and measured present profiles may be regarded as an indication that the model captures the dominant transport processes and may be used to predict future developments.
The first site is a field in a wastewater irrigation area near Braunschweig, Germany. When data for the field were collected by Streck (1993), it already had a 29-yr history of Cd pollution due to irrigated wastewater from the municipal wastewater treatment plant nearby. The downward movement of Cd at the site has been successfully modeled before by Streck and Richter (1997b), using an extended Freundlich equation derived from local Cd-sorption data. Briefly, Streck (1993) took a total of 480 soil samples at 48 locations on a regular 11 by 7 grid (120 m long and 72 m wide) with 10 subsamples at each location (first sample 0–30 cm, followed by 9 samples in 10-cm increments) down to a depth of 120 cm. On these samples, Streck (1993) measured the Cd concentration in the solution and sorbed phases as well as the soil properties pH and Corg. From these data, he derived an extended Freundlich equation relating measured sorption equilibria to pH and Corg. More details may be found in Streck and Richter (1997a).
The second site is located close to a metal smelter in the city of Nordenham, Germany. The site has an 85-yr history of Cd emissions from the smelter. The Cd transport at this site has been modeled before by Beyer (2002), using an extended Freundlich equation relating measured sorption data to pH and Corg. Data on Cd content, pH, and Corg on a 100 by 100 m plot, with 25 regular grid points and 5 depth increments down to a depth of 80 cm were used to derive a local transfer function for transport prediction (Beyer, 2002).
Streck (1993) and Beyer (2002) derived their site-specific extended Freundlich equations by relating Cd sorption to pH and Corg, but not to clay. This approach seems suitable on an intermediate scale (e.g., the field scale) for a relatively homogeneous soil texture because a pronounced variation of clay content is unlikely. The influence of clay on sorption is more or less constant and independent of location at a given site and can be treated as part of the intrinsic sorption coefficient. As a consequence, the measurement of spatial clay content distribution was omitted in both investigations.
The lack of clay content information poses a problem when applying the universal extended Freundlich equation to the two sites. The function was derived from data gathered on the national scale, and on this scale the variation of clay content cannot be neglected (Utermann et al., 2005). We therefore face the task of incorporating clay content at the two sites in our model predictions even though its spatial distribution was not investigated in either study.
Case Study 1: Wastewater Irrigation
Figure 4 shows the boundary conditions for the wastewater irrigation site near Braunschweig. The precipitation rate is the 30-yr average from a nearby weather station. The evapotranspiration rate is based on the site-specific crop rotation. The irrigation rate was retrieved from irrigation records, while the Cd concentration in the irrigation water was reconstructed from various data sources (Streck, 1993). Layerwise distribution parameters of pH and Corg reported by Streck (1993) are given in Table 2. The minimum and maximum values were set to the sample limits of the respective layer.
View this table:
[in this window]
[in a new window]
|
TABLE 2. Depth-specific distribution parameters of soil properties pH and Corg (clay content is assumed constant) at the wastewater site near Braunschweig, Germany.
|
|
As mentioned above, the spatial variation of clay content was not investigated at this site. However, a preliminary screening study by Streck (1993) suggests that clay content at the site has a nearly constant and rather low concentration of only about 1% of mass. Because the contribution of clay content to overall sorption is small when clay content is low, a clay content of 1% was assumed in each soil layer (see Table 2).
For the 2-D Monte Carlo simulation, a total of 40,000 columns were generated—200 columns for variability times 200 realizations of uncertainty. Column generation was performed using the parameters listed in Table 2 and covariance matrices (not shown) describing the correlations between Corg and pH in each layer. The generated pH and Corg distributions for the 50- to 60-cm layer are presented in Fig. 1. The generated pH and Corg distributions and the constant clay content were combined with the parameters and residual errors listed in Table 1 to calculate layerwise Freundlich coefficients for the 40,000 columns. The generated K distribution shown in Fig. 2 is also for the 50- to 60-cm layer.
The boundary conditions shown in Fig. 4 were used in combination with the Freundlich coefficients of a single column to calculate Cd displacement down to a depth of 120 cm in each of the columns. Evaporation from the soil surface was neglected, and infiltration at the surface is equal to precipitation plus irrigation. The rate of water taken up by roots is approximated by the potential evapotranspiration due to the high irrigation rates (Streck and Richter, 1997b). The empirical parameter = 5.87 m–1 describing the relative root density and the dispersion length = 0.01 m were taken from Streck and Richter (1997b). Water content was set to the steady state depth distribution of , corresponding to the flux rates q instead of field capacity, to be in line with Streck and Richter (1997b). The effect of plowing was simulated by perfect mixing of the Cd content in the plow layer (0.3 m) each year, using the same approach as Streck and Richter (1997b). The results were aggregated as illustrated in Fig. 3 and compared to measured data. In Fig. 5 the measured field-averaged dissolved and sorbed Cd concentrations are plotted versus depth. The measured concentrations are layer-wise arithmetic means of the 48 grid samples. In addition the figure shows the depth distribution predicted with the local extended Freundlich equation (Streck, 1993) as well as the depth distribution predicted in this study. The median of the predicted concentrations matches the measured field averaged depth profiles of the dissolved and sorbed concentration remarkably well. It is also close to the concentrations predicted with the local extended Freundlich equation.
View larger version (20K):
[in this window]
[in a new window]
|
FIG. 5. The sorbed and dissolved Cd-concentration profiles measured at the site in Braunschweig, Germany, sorbed and dissolved concentrations predicted with the local extended Freundlich equation (Streck and Richter, 1997b) as well as the depth distribution (including the 95% confidence interval) predicted with the extended Cd–Freundlich equation derived by Utermann et al. (2005).
|
|
Case Study 2: Metal Smelter
Figure 6 shows the boundary conditions for the metal smelter site in Nordenham. The average yearly precipitation is assumed to be 740 mm, based on data from the Deutscher Wetterdienst (1996). An average yearly evaporation of 610 mm is estimated for the site using the method of Renger and Wessolek (1990) for pastures as described by Hennings (2000). The Cd-emission concentration was reconstructed by a proportional distribution of the total mass of Cd recovered on the site based on the production numbers of the smelter.
Because no site-specific data on clay content is available, we searched the database of the soil information system of Lower Saxony (NIBIS) (Heineke and Eckelmann, 1998) for nearby clay content measurements. We retrieved four soil profiles with a close distance to the center of the site (ranging from 280 to 580 m) where the clay content was measured down to depth of 80 cm. All profiles are in the same mapping unit as the modeled site (Calcaric fluvisol, according to the World Reference Base [WRB, 1998]). The clay content in the different layers of these profiles is shown in Table 3. Down to a depth of 40 cm, the clay content stays relatively constant. Below 40 cm, there seems to be a shift to a material containing less clay. The omission of clay in the Beyer (2002) derivation of the site-specific extended Freundlich equation, based on the assumption that variations in clay content are small, is apparently incorrect at this site.
We derived distribution parameters for each layer by assuming that these profiles are representative for the site-specific clay distribution. Because these four soil profiles represent an area that is larger than the investigated site, it may be that the layerwise variation of clay is more pronounced than on the site itself. Correlations between clay contents in different layers and correlations of clay content with the two other properties (pH, Corg) were not considered. Assuming a lognormal distribution for clay, minimum and maximum clay content were estimated by taking the log transformed clay content and adding or subtracting twice the standard deviation to the mean.
In addition to the distribution characteristics of clay, Table 4 lists the layerwise distribution parameters and distribution types of pH and Corg reported by Beyer (2002). The minimum and maximum values for pH and Corg were set to the sample limits for each layer. The empirical parameter = 7.34 m–1 was estimated on a measured root density under pasture (Aboling, 1997). Dispersion length was set to 0.01 m based on the study of Cernik et al. (1994). The water content was set to field capacity, estimated at 0.41 (dimensionless) for this site based on a method reported in Hennings (2000). As in Case Study 1, evapotranspiration is assumed to be completely taken up by roots. The transport modeling procedure with the columns is the same as described for Case Study 1.
View this table:
[in this window]
[in a new window]
|
TABLE 4. Depth specific distribution parameters of soil properties pH, Corg, and clay at the site near Nordenham, Germany.
|
|
Figure 7 shows the field-averaged depth profile of the sorbed-Cd concentration for the site. The dissolved concentration was not measured. It also shows the depth distribution predicted with the local extended Freundlich equation (Beyer 2002) as well as the results from this study. The median of the predicted concentrations on the left matches the field-averaged depth profile of the sorbed concentration. It is also close to the concentrations predicted with the local extended Freundlich equation.
View larger version (25K):
[in this window]
[in a new window]
|
FIG. 7. The sorbed Cd-concentration profile measured at Nordenham, Germany, sorbed concentrations predicted with the local extended Freundlich equation (Beyer, 2002) as well as the depth distribution (including the 95% confidence interval) predicted with the extended Cd–Freundlich equation derived by Utermann et al. (2005).
|
|
|
Conclusions
|
|---|
We have successfully demonstrated that a relatively simple modeling approach, applicable when data are scarce, allowed the prediction of trace element transport in soil at two sites with a known history of Cd contamination. The modeling approach relies on a combination of (i) an extended Freundlich equation applicable nationwide to describe trace element sorption from easily measurable soil properties, (ii) a parallel soil column model, and (iii) a 2-D Monte Carlo technique. Instead of calculating a single deterministic result, the approach estimates a distribution of concentrations in the solid and liquid soil phase by explicitly accounting for the spatial variability and uncertainty in the input data as well as model uncertainty of the extended Freundlich equation. By taking into account input variances, the model is a valuable tool for groundwater pollution risk assessment of trace elements at sites where available data are sparse. The predicted probability distributions of sorbed and dissolved concentrations are as precise as possible under the given conditions and allow a direct evaluation of the certainty of a result. This is especially useful in situations where more simple methods do not allow a clear yes-or-no decision with regard to the transgression of legislative threshold values.
|
Appendix
|
|---|
Uncertainty of the Extended Freundlich Equation
The set of extended Freundlich equations derived by Utermann et al. (2005) for different trace elements is based on experiments with soil samples gathered at 133 sites in Germany. From this set, the extended Freundlich equation for Cd relating pH, Corg and clay content (Clay) to Cd-sorption has the following form:
 | [A1] |
where K* is the intrinsic Freundlich coefficient, Clay is clay content, a, b, and c are regression coefficients, n is the Freundlich exponent accounting for sorption nonlinearity, and is the residual error remaining after regression. Clay and Corg are in weight percentage. Utermann et al. (2005) demonstrated that knowledge of pH, Corg, and clay is the minimum information required for a successful application of these functions. Figure A1 shows measured and predicted (using Eq. [A1]) sorbed-Cd concentrations from the study of Utermann et al. (2005). The residual error between the two quantities in Fig. A1 is up to the order of one magnitude and must be considered in a risk assessment.
View larger version (19K):
[in this window]
[in a new window]
|
FIG. A1. Measured and predicted Cd concentrations in the sorbed phase for the 133 nationwide monitoring sites in Germany investigated by Utermann et al. (2005). Sorbed concentration (S) predicted was obtained by multiple linear regression between the dependant variable S measured and the independent variables pH, Corg, and clay using Eq. [A1] (Utermann et al., 2005).
|
|
In risk assessment, site- or field-averaged concentrations are compared to threshold values. Given that the extended Freundlich equation is based on isotherm measurements gathered at 133 soil monitoring sites of a nationwide network, the residual error is likely to combine an error that is site-specific (on-site error) and an error occurring between different sites. When field-average concentrations are calculated, on-site error will influence the absolute value of the field average but not its uncertainty. Uncertainty will only be influenced by the error between sites. To separate the on-site error from the between-site error, an ANOVA of the error residuals was performed. A prerequisite for an ANOVA is a constant sample error variance. Figure A2 shows a box plot of the residual error in the predicted sorbed concentration S based on Corg, pH, and clay. The residual error is plotted against the site (sample) ID. To increase confidence in the results of the ANOVA, the analysis was limited to residuals from sites where isotherms were measured in at least four horizons, leaving 36 sites for analysis.
View larger version (27K):
[in this window]
[in a new window]
|
FIG. A2. Box plot of the on-site errors in the predicted sorbed concentration S (extended Cd–Freundlich equation based on Corg, pH, and clay).
|
|
Figure A2 shows that the prerequisite of a nonvarying on-site error variance is approximately fulfilled. Two sites appear to have a relatively large error variance (BW71 and NS57). Figure A3 shows four additional diagnostic plots that characterize the distribution of the on-site residuals (or error variance) in more detail. The two plots on the left show the absolute residuals and the square root of the standardized residuals versus the fitted ANOVA-mean for each site. Both plots allow a visual test of constant variance and influence of outliers. Outliers are more prominent in the upper plot, while in the lower plot resolution of the ordinate is increased because of standardization. Both plots include a moving average (averaged over 0.2 units on the abscissa) to check for constant variance. The higher resolution in the lower plot leads to more fluctuations but neither trend nor prominent peaks are visible in the moving average. The upper-right plot, where the standardized residuals are compared to quantiles of a normal distribution, reveals a good agreement between the two. In the lower-right plot, the Cook's distance (Cook and Weisberg, 1982), a measure for the influence of outliers on an ANOVA-fit, is plotted for every single observation. No prominent outlier is visible in the plot, which means that no single observation has a prominent influence on the ANOVA-fit.
Results of the ANOVA are presented in Table A1. The F value suggests that there is a difference between the site-specific means at the 1% significance level. Comparing the sums of squares reveals that on-site variance is approximately 70% of the total variance and between-site variance is approximately 30% of the total variance. These variance fractions are approximate because the between-sites sum of squares is a biased estimator for the error variance. However, the on-site sum of squares is an unbiased estimator and can readily be used to estimate the standard deviation of the normally distributed on-site random error with a mean of zero.
A complete statistical description of the residual error is achieved by providing the standard deviation of the on-site error and a statistical distribution of the site-specific error means that characterize the between-sites error. The left panel of Fig. A4 displays a histogram and summary statistics of site-specific group means obtained from the ANOVA procedure, and the right panel shows the corresponding quantile–quantile plot. Figure A4 suggests that the site-specific means are approximately described by a normal distribution but may deviate moderately in the tails. We approximate the distribution of site-specific means in the following as a normal distribution with mean and standard deviation as given in the left panel of Fig. A4.
The residual error term in the extended Freundlich equation (Eq. [A1]) is broken into two normally distributed error terms, representing contributions from on-site error, 1, and between-sites error, 2, yielding
 | [A2] |
Note that by deriving a local extended Freundlich equation from site-specific sorption data the between-sites error 2 could be neglected—such a function is unlikely to exhibit a site-specific bias. In contrast, the on-site error 1 may only be reduced by choosing a more accurate process-oriented sorption model. A consequence is that in the context of our model approach 2 is treated like uncertainty (reducible), while 1 is related to variability (irreducible). In addition, since the random error 1 varies on the site scale (the field scale), it contributes to the absolute value of the field-average concentration but not to its uncertainty (which is only influenced by 2). This is important considering that risk assessments will be based on field-averaged concentrations.
In the context of the 2-D Monte Carlo approach, the statistical distribution of 2 is used to generate l mean residual error values, comparable to different realizations of a site-specific bias. The mean residual error 2 is kept constant in the layerwise calculation of Freundlich coefficients in each of the k columns of a single set describing on-site variability, but it is varied over the l sets describing uncertainty. In contrast, the statistical distribution of 1 is used to generate an on-site random error that varies for each soil layer in the k columns that describe on-site variability. The values of 1 and 2 are assumed statistically independent, in contrast to the soil properties.
An actual simulation run of one set of k columns is therefore based on layerwise Freundlich coefficients calculated with a varying error 1 but a constant error 2; 2 is only varied between the l different sets. Figure A5 shows the distribution of the CDFs for the combined residual error 1 + 2. Because the residual error distributions are independent of soil layer and depth, one CDF belonging to this distribution applies to each soil layer within one set of k columns.
|
ACKNOWLEDGMENTS
|
|---|
This study was partly financed by the German Federal Ministry for Education and Science (BMBF) within the Federal Research Program Sickerwasserprognose.
|
REFERENCES
|
|---|
- Aboling, S. 1997. Untersuchungen zu Vegetation, Wurzellängendichte und Futterqualität intensiv und extensiv bewirtschafteter Rinderweiden mit besonderer Berücksichtigung der Randbereiche. Ph.D. diss. Univ. Hannover, Germany.
- Beyer, C. 2002. Sickerwasserprognose für einen schwermetallimmissionsbelasteten urbanen Raum unter besonderer Berücksichtigung von Unsicherheit und Variabilität. Diploma thesis. Tech. Univ. Carolo-Wilhelmina, Braunschweig, Germany.
- Butters, G.L., and W.A. Jury. 1989. Field scale transport of bromide in an unsaturated soil: 2. Dispersion modeling. Water Resour. Res. 25:1583–1589.
- Cernik, M., P. Federer, M. Borkovec, and H. Sticher. 1994. Modeling of heavy metal transport in a contaminated soil. J. Environ. Qual. 23:1239–1248.[Web of Science]
- Cook, R.D., and S. Weisberg. 1982. Residuals and influence in regression. Chapman and Hall, London.
- Cullen, A.C., and H.C. Frey. 1999. Probabilistic techniques in exposure assessment: A handbook for dealing with variability and uncertainty in models and inputs. Plenum Press, New York.
- Dagan, G., and E. Bresler. 1979. Solute dispersion in unsaturated heterogeneous soil at field scale: I. Theory. Soil Sci. Soc. Am. J. 43:461–467.
- Deutscher Wetterdienst. 1996. Klimadaten von Deutschland, Zeitraum 1961–1990. Selbstverlag des Deutschen Wetterdienstes, Offenbach am Main, Germany.
- Dowdy, R.H., and V.V. Volk. 1983. Movement of heavy metals in soils. p. 229–240. In D.W. Nelson, D.E. Elrick, and K.K. Tanji (ed.) Chemical mobility and reactivity in soil systems. SSSA Spec. Pub. 11. ASA, SSSA, Madison, WI.
- Federal Republic of Germany. 1998. Federal soil protection act of 17 March 1998 (Federal law gazette I.: P. 502): Act on protection against harmful changes to soil and on rehabilitation of contaminated sites (Federal soil protection act—BBodSchG). Available at http://www.umweltbundesamt.de/altlast/web1/berichte/pdf/bbodschg.pdf (verified 29 May 2007). Federal Republic of Germany.
- Federal Republic of Germany. 1999. Federal soil protection and contaminated site ordinance (BbodSchV). Available at http://www.umweltbundesamt.de/altlast/web1/berichte/pdf/bbodschv-engl.pdf (verified 29 May 2007). Federal Republic of Germany.
- Hawkins, J.L., M.I. Sheppard, and S.S. Jorgensen. 1995. Predicting soil lead migration: How can ancient church roofs help? Sci. Total Environ. 166:43–53.[CrossRef]
- Heineke, H.J., and W. Eckelmann. 1998. Development of soil information systems in the Federal Republic of Germany: Status report. p. 125–132. In H.J. Heineke, W. Eckelmann, A.J. Thomasson, R.J.A. Jones, L. Montanarella, and B. Buckley (ed.) Land information systems: Developments for planning the sustainable use of land resources. European Soil Bureau Research report no. 4. Office for Official Publications of the European Communities, Luxembourg.
- Hennings, V. 2000. Methodendokumentation Bodenkunde. Auswertungsmethoden zur Beurteilung der Empfindlichkeit und Belastbarkeit von Böden, Sonderhefte Reihe G- Geologisches Jahrbuch. Vol. 1. Bundesanstalt für Geowissenschaften und Rohstoffe, Hannover, Germany.
- Horn, A.L., R.-A. Düring, and S. Gäth. 2004. Sorption of Cd in soils: Pedotransfer functions for the parameters of the Freundlich sorption isotherm. Water Air Soil Pollut. Focus 4:61–71.
- Ingwersen, J. 2001. The environmental fate of cadmium in the soils of the waste water irrigation area of Braunschweig: Measurement, modelling, and assessment. Ph.D. diss. Tech. Univ. Carolo-Wilhelmina, Braunschweig, Germany. Available at http://www.digibib.tu-bs.de/?docid=00001266 (verified 26 Mar. 2007).
- Ingwersen, J., and T. Streck. 2006. Modeling the environmental fate of cadmium in a large wastewater irrigation area. J. Environ. Qual. 35:1702–1714.[Abstract/Free Full Text]
- Jury, W.A., and K. Roth. 1990. Transfer functions and solute movement through soils. Theory and applications. Birkhäuser, Basel, Switzerland.
- Jury, W.A., D. Russo, G. Streile, and H. Elabd. 1990. Evaluation of volatile organic chemicals residing below the soil surface. Water Resour. Res. 26:13–20.[CrossRef]
- Jury, W.A., and D.R. Scotter. 1994. A unified approach to stochastic-convective transport problems. Soil Sci. Soc. Am. J. 58:1327–1336.[Web of Science]
- Knuth, D.E. 1981. Seminumerical algorithms. The art of computer programming. Vol. 2. 2nd ed. Addison-Wesley, Reading, MA.
- Maskall, J., K. Whitehead, C. Gee, and I. Thornton. 1996. Long-term migration of metals at historical smelting sites. Appl. Geochem. 11:43–51.[CrossRef]
- Maskall, J., K. Whitehead, and I. Thornton. 1995. Heavy metal migration in soils and rocks at historical smelting sites. Environ. Geochem. Health 17:127–138.[CrossRef][Web of Science]
- Press, W.H., B.P. Flannery, S.A. Teukolsky, and W.T. Vetterling. 1988. Numerical recipes in C: The art of scientific computing. Cambridge Univ. Press, Cambridge, UK.
- Renger, M., and G. Wessolek. 1990. Auswirkungen von Grundwasserabsenkung und Nutzungsänderung auf die Grundwasserneubildung. p. 295–307. In Folgen anthropogener Einflüsse auf den Wasserhaushalt und die Wasserbewirtschaftung. Vol. 386. Mitt. Inst. für Wasserwesen, Univ. der Bundeswehr München, Germany.
- Römkens, P.F.A.M., J.E. Groenenberg, L.T.C. Bonten, W. de Vries, and J. Bril. 2004. Derivation of partition relationships to calculate Cd, Cu, Ni, Pb, and Zn solubility and activity in soil solutions. Report no. 305. Alterra, Wageningen, the Netherlands.
- Simmons, C.S. 1982. A stochastic-convective transport representation of dispersion in one-dimensional porous media systems. Water Resour. Res. 18:1193–1214.
- Smith, G.D. 1985. Numerical solution of partial differential equations: Finite difference methods. 3rd ed. Clarendon Press, Oxford, UK.
- Stöfen, H. 2005. Entwicklung eines Verfahrens für Sickerwasserprognosen im Sinne der Bundes-Bodenschutz- und Altlastenverordnung. Ph.D. diss. TU Hamburg-Harburg, Germany.
- Streck, T. 1993. Schwermetallverlagerung in einem Sandboden im Feldmaßstab: Messung und Modellierung. Ph.D. diss. Tech. Univ. Braunschweig, Germany.
- Streck, T., and H. Piehler. 1998. On field-scale dispersion of strongly sorbing solutes in soils. Water Resour. Res. 34:2769–2773.[CrossRef]
- Streck, T., and J. Richter. 1997a. Heavy metal displacement in a sandy soil at the field scale: II. Measurements and parameterization of sorption. J. Environ. Qual. 26:49–56.[Web of Science]
- Streck, T., and J. Richter. 1997b. Heavy metal displacement in a sandy soil at the field scale: II. Modeling. J. Environ. Qual. 26:56–62.
- Swartjes, F.W. 1990. Numerische Simulation der eindimensionalen Schwermetallverlagerung im homogenen, gesättigten/ungesättigten Boden. Ph.D. diss. Technische Universität Berlin, Germany.
- USEPA. 2001. Risk assessment guidance for superfund. Vol. 3, Part A. Process for conducting probabilistic risk assessment. EPA/540-R-02-002. USEPA, Washington, DC.
- Utermann, J., G. Meyenburg, S. Altfelder, H.E. Gäbler, W. Duijnisveld, A. Bahr, and T. Streck. 2005. Entwicklung eines Verfahrens zur Quantifizierung von Stoffkonzentrationen im Sickerwasser auf der Grundlage chemischer und physikalischer Pedotransferfunktionen. Final Report, BMBF-Forschungsvorhaben 02WP0206.
- van Dam, J.C., G.H. de Rooij, M. Heinen, and F. Stagnitti. 2004. Concepts and dimensionality in modeling unsaturated water flow and solute transport. p. 1–36. In R.A. Feddes, G.H. de Rooij, and J.C. van Dam (ed.) Hydrology of the unsaturated zone: Progress, challenges, and applications. Proc. of the Int. Wageningen UR Frontis Symp. Wageningen, the Netherlands. 2–5 October 2004. Kluwer Academic, Dordrecht, the Netherlands.
- Vanderborght, J., and H. Vereecken. 2007. Review of dispersivities for transport modeling in soils. Vadose Zone J. 6:29–52.[Abstract/Free Full Text]
- van der Zee, S.E.A.T.M., and J.J.T.I. Boesten. 1991. Effects of soil heterogeneity on pesticide leaching to groundwater. Water Resour. Res. 27:3051–3063.[CrossRef]
- van der Zee, S.E.A.T.M., and W.H. van Reimsdijk. 1987. Transport of reactive solutes in spatially variable soil systems. Water Resour. Res. 23:2059–2069.
- WRB. 1998. World reference base for soil resources. World Soil Resources report no. 84. ISSS/ISRIC/FAO, Rome.
TOP
ABSTRACT
INTRODUCTION
