JEQ Grow Your Career With ASA
HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
QUICK SEARCH:   [advanced]


     

Published online 9 August 2006
Published in J Environ Qual 35:1702-1714 (2006)
DOI: 10.2134/jeq2005.0412
© 2006 American Society of Agronomy, Crop Science Society of America, and Soil Science Society of America
677 S. Segoe Rd., Madison, WI 53711 USA
This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (9)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Ingwersen, J.
Right arrow Articles by Streck, T.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Ingwersen, J.
Right arrow Articles by Streck, T.
Agricola
Right arrow Articles by Ingwersen, J.
Right arrow Articles by Streck, T.
Related Collections
Right arrow Municipal Wastes
Right arrow Spatial Variability
Right arrow Upscaling
Right arrow Solute Transport Models
Right arrow Heavy Metals

TECHNICAL REPORTS

Heavy Metals in the Environment

Modeling the Environmental Fate of Cadmium in a Large Wastewater Irrigation Area

Joachim Ingwersen* and Thilo Streck

Univ. of Hohenheim, Institute of Soil Science and Land Evaluation, Biogeophysics Section, D-70593 Stuttgart, Germany

* Corresponding author (jingwer{at}uni-hohenheim.de)

Received for publication October 31, 2005.

   ABSTRACT
TOP
 ABSTRACT
INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
The fate of cadmium in soils is governed by spatially heterogeneous processes that proceed from decades to centuries. This study aimed at modeling the fate of Cd within the wastewater irrigation area (WIA) of Braunschweig (Germany). The sandy soils (mainly Dystric Cambisol or Typic Haplumbrept) at this site (28 km2) have received considerable loads of heavy metals by irrigation of municipal wastewater for up to 40 yr. The soils of the WIA are in agricultural use. The main crops are sugar beet (Beta vulgaris L.), potato (Solanum tuberosum L.), and wheat (Triticum aestivum L.). As a result of asparagus (Asparagus officinalis L.) cropping, about 15% of the soils have been converted to Rigosols. In 1996, we measured the vertical distribution (0 to 1.2 m) of soil pH, organic carbon content, and the EDTA–extractable content and the solution phase concentration of Cd at 153 sites. At sites not used for asparagus cultivation, Cd has migrated on average to a depth of about 0.5 m. Due to deep plowing, which accelerates migration, Cd has been displaced on average to about 0.7 m at the Rigosol sites. To model the fate of Cd at the scale of the WIA, we used different parallel soil column approaches. In each column the local model SEFAH was used to simulate both displacement and plant uptake of Cd. The model was fed with measured or randomly generated soil data. The results of retrospective simulations from 1957 to 1996 agreed well with observed Cd profiles. The better the spatial variability of sorption was described, the better the performance. Our simulation results show that Cd pollution of soil at first affects the soil-plant pathway. The breakthrough of Cd to the groundwater is dampened and is delayed for many decades.

Abbreviations: BMWA, Braunschweig Municipal Wastewater Association • EDTA, ethylene-diamine-tetra-acetic acid • GM, grid model • MC, Monte-Carlo • OC, organic carbon • PSC, parallel soil column • PTF, pedotransfer function • WIA, wastewater irrigation area


   INTRODUCTION
TOP
 ABSTRACT
 INTRODUCTION
MATERIALS AND METHODS
 MODELING
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
HEAVY METALS, such as Cd, enter the environment via atmospheric deposition and fertilization as well as by the application of sewage sludge, wastewater, or compost. The relevant spatial scale of this environmental issue is often the field or the regional scale. From field-scale studies we know that the transport of sorbing solutes is dominantly governed by retardation, which strongly depends on the influence of sorption properties of the soil such as pH and organic carbon (OC) content. Van der Zee and van Riemsdijk (1987) elucidated this point in a theoretical work. Streck and Richter (1997a, 1997b) investigated a field (120 by 72 m) that had been irrigated with municipal wastewater for 29 yr. Streck and Richter (1997a, 1997b) found that displacement was spatially highly variable and that field-scale dispersion of Cd and Zn could be attributed to a large extent to the spatial variability of sorption properties of soil.

The direct use of local process equations and parameters at larger spatial scales is in general not feasible (Refsgaard and Butts, 1999). Appropriate upscaling procedures are needed to transfer local models to the field or regional scale. A widely used upscaling procedure is the parallel soil column (PSC) approach, which has also been referred to as the "stream tube" model (Jury and Roth, 1990). In this approach, the region is represented as an ensemble of vertical, not connected soil columns, where transversal mixing is excluded. The solute transport in each column is described either as convective or as convective-dispersive. To account for spatial variability of sorption both Monte-Carlo (MC) methods and the grid model (GM) have been used. In the MC method (Ammoozegar-Fard et al., 1982; Seuntjens et al., 2002; Streck and Richter, 1997b), each soil column is parameterized with random input variables. The random variables can be correlated or uncorrelated. The field- or regional-scale results are obtained by averaging the local model outputs of the PSC ensemble. In the GM approach, simulations are performed at points in the region where the parameters are known from measurements (Anlauf, 1988; Streck and Richter, 1997b).

Besides solute transport, plant uptake is another important process in the environmental fate of heavy metals. To simulate plant uptake, several models have been developed during the last few decades. These models differ in their conceptual approach, degree of complexity, and consequently, in data input needed. An example of a model that incorporates a very detailed description of the most fundamental mechanisms is the Barber-Cushman model (Barber and Cushman, 1981). The model describes the radial flow of nutrients by diffusion and convection through soil to plant roots. In combination with a Michaelis-Menten kinetically controlled active solute uptake, nutrient depletion zones around the root can be modeled. Models that treat root solute uptake in a simplified manner have been published by Behrendt et al. (1995), Christensen and Tjell (1984), Schoups and Hopmans (2002), and Ingwersen and Streck (2005). These models do not explicitly consider the underlying mechanisms governing root solute uptake. In these models, uptake equals the product of water transpired, solute concentration in soil solution, and a plant-specific empirical parameter. This simplification is accepted for the benefit of less model data input.

Regional-scale studies on the fate of heavy metals in soil are scarce. Studies of Keller et al. (2001) and Tiktak et al. (1998) focused on the prediction of heavy metal accumulation in topsoil in a regional context and are based on a mass-balance calculation. In both studies, leaching is calculated only as convective transport from topsoil to subsoil and metal transport in the subsoil is not considered. Plant uptake is either calculated with approaches similar to the simplified ones mentioned above (Keller et al., 2001) or even ignored because of the limited effect of harvesting on the heavy metal balance of topsoil (Tiktak et al., 1998). The core of both models is a pedotransfer function (PTF) that is used to estimate the Freundlich coefficient of the sorption isotherm based on soil properties such as OC content, clay content, and pH. Soil properties are assessed by analyzing soil maps using geographic information systems. Tiktak et al. (1998) evaluated the model validity by testing whether it was possible to reconstruct the present Cd contents of the Netherlands' topsoils (33 800 km2) using estimates of historical emission rates. Despite the many simplifications and the large uncertainties of model input data such as land use, soil properties, or pH dynamics, current Cd contents were only slightly underestimated by the model. Keller et al. (2001) compared estimated Cd and Zn balances of different land use systems with reported metal balance studies on experimental farms within a rural region (95 km2) in northwestern Switzerland and found that they were in good agreement.

The study area of the present work was the wastewater irrigation area (WIA) of Braunschweig, which is the largest of its kind in Europe. The sandy soils within this region (29 km2) have received considerable loads of heavy metals by irrigation using municipal wastewater for up to 40 yr. A detailed description of the WIA and its operation history is given in Ingwersen (2001). In a previous paper (Ingwersen and Streck, 2005), we introduced, calibrated, and tested the plant uptake model used in the present study. The plant model assumes that uptake of Cd is proportional to mass flow, i.e., the product of water transpired, Cd concentration in soil solution, and a plant-specific empirical parameter. The model was specifically calibrated for the conditions of the WIA using extensive field data. Within the WIA, we sampled soil and plant material from 40 potato, 40 sugar beet, and 32 winter wheat fields. After site-specific calibration, simulations agreed well with the observed Cd contents in plant material. The model explained between 66% and 87% of the observed variance.

The overall goal of this study was to simulate the environmental fate of Cd in the soils of the WIA employing different PSC approaches. Based on an extensive field data set, model approaches were parameterized and tested. For model parameterization it was necessary (i) to quantify the spatial variability of soil properties governing sorption, (ii) to set up a site-specific extended sorption isotherm, and (iii) to reconstruct the Cd input history. Finally, we used the model to carry out scenario simulations with regard to leaching and plant uptake.


   MATERIALS AND METHODS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
MODELING
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
The Wastewater Irrigation Area of Braunschweig
The WIA (43 km2) is located in the northwest of the city of Braunschweig (250 000 residents), Germany. It is part of the city's sewage management system. Today, this system consists of three components: the sewage plant "Steinhof," a sewage farm, and the WIA (Fig. 1 ). After a mechanical and a primary and secondary biological treatment, one-third of the effluent of the sewage plant is pumped into the sewage farm and the remaining two-thirds are pumped to the WIA, where they are used for sprinkling irrigation. In 1996, about 15 million m3 were applied to an area of 28.7 km2, which corresponds to an annual irrigation of about 523 mm. The irrigated area is in agricultural use. The main crops are sugar beet, potato, winter cereal, and summer cereal. The remaining 14.3 km2 of the WIA are mainly forests and urban areas. The irrigation area is divided into four pumping station districts. In 1957, wastewater irrigation started in district 1. In the other three districts irrigation started within the following 9 yr. In the early 1980s, the annual Cd loads to the WIA reached maximum values of about 200 g ha–1. Today annual Cd loads are below 5 g ha–1. Cadmium contents (determined by aqua regia extraction) in topsoil range from 0.10 to 1.49 mg kg–1 and are on average 0.36 mg kg–1 (Ingwersen, 2001).


Figure 1
View larger version (38K):
[in this window]
[in a new window]
  		Fig. 1. Sketch of the wastewater irrigation area (WIA) of Braunschweig, most important installations, and sampling sites.

 
Until the opening of the sewage plant in 1979, raw sewage was used for irrigation. The raw sewage contained high amounts of particles. To avoid clogging of the nozzles of the sprinklers, most of the sludge was removed by sedimentation in the settling tanks of each pumping station. The collected sludge was pumped into nearby sludge basins, where it was composted for 2 to 3 yr. After compostation, the sludge was used to fertilize fields within the WIA. In 1975, new irrigation devices with larger nozzles were installed, so that higher amounts of sludge could be used for irrigation. Subsequently, the sludge basins of districts 1 and 2 were shut down. Until the opening of the sewage plant in 1979, the sludge of these two districts was pumped to districts 3 and 4, where the additional sludge was composted and applied to fields (Abwasserverband Braunschweig, 1979).

The dominant soil type in the WIA is a Dystric Cambisol (Food and Agriculture Organization, 1990) or Typic Haplumbrept (Soil Survey Staff, 1994). In areas with a high groundwater table, Cambic Gleysols (US Soil Taxonomy: Typic Haplaquept) are also present. Toward the north, the portion of Humic Cambisols (US Soil Taxonomy: Typic Hapludoll) increases.

As in other sandy areas of the so-called Geest in Lower Saxony (northwest Germany), asparagus is cultivated in the WIA. Before asparagus is planted, the humus and nutrient-rich topsoil is transferred to a depth of about 0.3 to 0.6 m by deep plowing. This anthrophic soil developed from deep plowing for asparagus cultivation is denoted according to the German soil taxonomy as a Rigosol (FAO: Hortic Anthrosol; US Soil Taxonomy: Typic Plagganthrept). The total area converted to Rigosols within the WIA is about 15% (Ingwersen, 2001).

Sampling
The investigation area was divided into 500 by 500 m blocks. In each block, a field was randomly selected where we took at one location, one soil sample in the following way. First, we excavated a pit in the Ap horizon (about 0.5 by 0.5 by 0.3 m). From the excavated material we took a subsample of about 3 kg of field-wet soil. At the bottom of the pit, we took a soil monolith of 0.07-m diameter and 0.9-m depth using a mechanical soil auger (Cobra, Atlas Copco Inc., Essen, Germany). Each monolith was divided into 9 soil layers (0.3–0.4, 0.4–0.5, ..., 1.1–1.2 m). In each layer a sample was taken from the longitudinal center. The large diameter of the probe guaranteed a spreading-free sampling. Both top and subsoil samples were stored in plastic cups, dried at 40°C, and passed through a 2-mm sieve. Altogether 153 soil monoliths were taken, of which 14 were Rigosols (Fig. 1).

Measurements
All soil samples were analyzed for pH, OC content, EDTA–extractable Cd content, and Cd concentration in 2.5 mM CaCl2. We assumed that the EDTA–extractable content represents the fraction of Cd that participates in sorption–desorption reactions (Filius et al., 1991; Gäbler et al., 1999). Furthermore, we assumed that a 24 h equilibration of soil with 2.5 mM CaCl2 yields the concentration in soil solution to a good approximation (Streck and Richter, 1997a).


   MODELING
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 
The Local Model SEFAH
To simulate solute transport and plant uptake at the local scale we used the 1D SEFAH (Simulating the Environmental Fate of Heavy Metals) model (Streck and Richter, 1997b).

Cadmium fate was described with the convection-dispersion equation

Formula 1[1]
where C (kg m–3) denotes the concentration in soil solution, S (kg kg–1) the sorbed phase concentration of the solute, DS (m2 s–1) the apparent dispersion coefficient, q the water flux density (m s–1), {theta} (m3 m–3) the volumetric water content, and {rho} (kg m–3) the bulk density. The sink term {phi} (kg m–3 s–1) accounts for the uptake of solutes by plant roots. Symbols t (s) and z (m) are time and space coordinates, respectively. Nonlinear sorption was described by a Freundlich sorption isotherm

Formula 2[2]
where k (kg–m m3m) is the Freundlich coefficient and m is the Freundlich exponent.

Assuming that plant uptake is mass–flow-governed and that root water uptake is proportional to root length distribution, where the latter is described by an exponential function, the sink term {phi} may be written as a function of depth

Formula 3[3]
Here, {omega} (m–1) is an empirical parameter describing the root distribution with depth, {eta} is a crop- and cultivar-specific empirical parameter, and TP (m s–1) stands for the total transpiration. The transpiration was computed with the approach of Bierhuizen and Slatyer (1965)

Formula 4[4]
where {Delta}e (Pa) is the average saturation deficit of the air during the main vegetation period, kp (Pa) is a crop-specific constant, {rho}w (kg m–3) is the density of water, and Y is the dry matter production (kg m–2 s–1). Dividing a soil in n layers of thickness {Delta}z, integrating Eq. [3], and taking into account the partitioning of absorbed Cd between processed (e.g., grain) and unprocessed plant parts (e.g., straw), the Cd content of processed plant part may be written as:

Formula 5[5]
Here, QHM denotes the Cd content ratio and QY the dry matter yield ratio between unprocessed and processed plant parts. This approach has been successfully tested for the conditions of the WIA (Ingwersen and Streck, 2005).

In case of strongly sorbing solutes, the water regime can be modeled in a simplified way because displacement of Cd in soil is insensitive to volumetric water content (Streck and Piehler, 1998; Seuntjens et al., 2002). Therefore, the model runs with constant volumetric water content with depth. Moreover, short-term variations of water flux have a negligible effect on heavy metal movement and the water regime can be simulated under steady-state conditions (Swartjes, 1990). The water flux density can then be written as

Formula 6[6]
where I (m s–1) and P (m s–1) denote the irrigation and precipitation rates (rainfall and snow), respectively, and E (m s–1) is the evaporation.

The initial distribution of the soil solution concentration was chosen as

Formula 7[7]
whereby C0 denotes the background concentration. The Cd concentration at the upper boundary Cinf (kg m–3) was calculated by

Formula 8[8]
where Cw (kg m–3) denotes the Cd concentration in wastewater and Cp (kg m–3) the concentration in precipitation.

The model SEFAH contains a routine that simulates the metal displacement by annual plowing and deep plowing. Because heavy metal leaching is very slow, annual plowing can be regarded as a soil homogenizing process (Schimmack and Bunzl, 1986). We assumed that after plowing, all compartments of the Ap horizon have the same total concentration. The new equilibrium after plowing between solution and sorbed phase concentration was calculated by numerically solving the mass-balance equation using the Newton–Raphson algorithm.

Besides annual plowing of the Ap horizon, deep plowing was only used on the newly installed asparagus plantations. Deep plowing transfers soil material from the Ap horizon (0 to 0.3 m) to a depth of about 0.3 to 0.6 m and vice versa. Depending on the depth of deep plowing, m layers below the Ap horizon are affected. After deep plowing the total metal concentration in the Ap horizon CT,0new may be written as

Formula 9[9]
where h0 (m) denotes the thickness of the Ap horizon, {rho}0 (kg m–3) the bulk density of the Ap horizon, and {alpha}i (kg m–2) describes the soil exchange between the Ap horizon and the i-th layer. If the sum of {alpha}i equals h0{rho}0 the topsoil is completely exchanged by subsoil. The total concentration in the i-th layer may be calculated by

Formula 10[10]
Due to deep plowing, not only Cd, but also OC and protons are exchanged between the Ap horizon and subsoil. Changes in the OC content and soil pH caused by deep plowing were also calculated using Eq. [9] and Eq. [10]. The new Freundlich coefficient of each layer was estimated by PTF 1 (see Model Parameterization section). Finally, the new equilibrium between solution and sorbed phase concentration was computed.

Deep plowing forces the OC content in Ap horizon into a new, nonequilibrium state. For simplicity, we assume that (i) after plowing the OC content linearly approaches the former OC content, (ii) it takes 25 yr to reach the former equilibrium state (Gäth et al., 1997), and (iii) the OC content of layers below the Ap horizon is constant during simulation.

Regional Modeling
Grid Model
In the GM approach the local model SEFAH was applied at points in the region where measurements had been taken, and the Freundlich coefficient k could be estimated from measured soil pH and OC content by the PTF:

Formula 11[11]
We used two versions of this PTF. In PTF 1, the intrinsic Freundlich coefficient k* and the empirical parameters a and b were assumed to be constant. In PTF 2, k* is assumed to vary horizontally, that means, for each soil monolith a site-specific k* was used and the PTF becomes

Formula 12[12]
where ki* is the intrinsic Freundlich coefficient of the i-th monolith. These two approaches are denoted as GM 1 and GM 2.

The Cd load ICd (mg m–2) of each sampled soil monolith was estimated by

Formula 13[13]
where CEDTA (mg kg–1) and BGEDTA (mg kg–1) denote the measured EDTA–extractable Cd content in 1996 and the estimated EDTA–extractable Cd background content before irrigation started, respectively. The value of BGEDTA was estimated by applying the PTF 1 with the background concentration C0.

Monte-Carlo Simulations
MC 1
The MC 1 modeling approach uses random input variables. Based on our field observations, we assumed that OC contents and Cd loads are lognormally distributed and pH values are normally distributed. In MC methods the Freundlich coefficient k was estimated by PTF 1. Preliminary experiments had shown that at least 400 simulations were necessary to derive first and second moments of the model output. A set of 500 random soil columns with z-correlated random numbers from a multivariate normal distribution was generated. The vector of means and the covariance matrix were estimated from the data set of non-Rigosols (N = 139). The random numbers were generated using the statistic program Statgraphics for DOS Version 7.0 (Manugistics, 1999). We confined the analysis of the simulation result to the first and second moments of the model output.

MC 2
The MC 1 approach ignores the spatial structures of soil pH, OC content, and Cd load. The MC 2 modeling approach incorporates these two characteristics by combining the GM and the geostatistical method of "Conditional Simulation" (Journel and Huijbregts, 1991). A conditional simulation is one "possible realization" of a random field, which reproduces the first two spatial moments of the measured data (mean and variogram, as well as the histogram). Moreover, the random field is conditioned on the observations, i.e., at the sampling points, simulated and observed values are the same. The input variables considered random were generated as follows:

  1. A pool of random soil columns (N = 4000) with multivariate normally distributed pH and multivariate lognormally distributed OC content was generated using the vector of means and covariance matrix as estimated from the measurements. As in MC 1, each soil column consists of 10 layers (0–0.3 m, 0.3–0.4 m, ..., 1.1–1.2 m).
  2. Variogram models were fitted to experimental variograms of pH and OC content at layers with the largest spatial correlation length (pH, layer 8 [0.9 to 1.0 m]; OC content, layer 1 [0 to 0.3 m]) and of the Cd load (Fig. 2 ).
  3. A total of 500 (x,y) coordinates were randomly drawn from the arable area of the WIA. At each point the depth of groundwater table was taken from a map of interpolated groundwater tables.
  4. The Cd load, pH value of layer 8, and the OC content of layer 1 were computed at each of the 500 (x,y) coordinates by the method of conditional simulation using the geostatistical software package GStat Version 2.0 g (Pebesma and Wesseling, 1998). Three two-dimensional random fields of horizontally correlated Cd loads, pH, and OC contents were generated.
  5. Finally, soil columns from the ensemble generated in step 1 were randomly assigned to the two-dimensional fields. A matching soil column had to show, within defined tolerances, a similar pH and a similar OC content as the corresponding two-dimensional random field (pH ± 0.05 units; OC content ± 5% deviation from generated value).
Because the 1D random soil columns and the 2D random fields are generated separately, the ensemble of soil columns produced is not truly but quasi 3D correlated.


Figure 2
View larger version (12K):
[in this window]
[in a new window]
  		Fig. 2. Experimental and fitted semivariograms of (A) soil pH in 0.9 to 1.0 m layers (exponential model: sill = 0.365, range = 937 m), (B) organic carbon content in Ap horizons (spherical model: nugget = 0.053, partial sill = 0.154, range = 2192 m), and (C) Cd loads (exponential model: nugget = 0.145 mg2 m–2, partial sill = 0.219 mg2 m–2, range = 429 m). Organic carbon contents (% by mass) and Cd loads were log-transformed before analysis.
 
Model Testing, Model Assessment Criteria, and Sensitivity Analysis
The validity of the model was tested by hindcast simulations. We performed retrospective simulations from 1957 to 1996 and compared the final simulated Cd profiles with those observed in 1996. To compare the modeling results we applied two criteria: modeling efficiency EF (Loague and Green, 1991) and displacement index DI. The EF is defined as the proportion of the total variance of observed data explained by the model:

Formula 14[14]
Here, Pi and Oi denote predicted and observed values, respectively, n is the number of samples, and Formula 14 is the mean of the observed data. If measured and predicted data are all equal, EF is unity. If EF is less than zero, the predicted values are less accurate than simply using the observed mean.

The DI is defined as the ratio between predicted and observed EDTA–extractable Cd in subsoil (0.3 to 1.2 m)

Formula 15[15]
where CPi (kg kg–1) denotes the predicted EDTA–extractable Cd content and COi (kg kg–1) the observed EDTA–extractable Cd content of the i-th subsoil layer. Hence, the DI is a measure of the extent to which the model overestimates or underestimates the Cd displacement. If DI is less than unity, for example, the model underestimates Cd displacement.

In the sensitivity analysis we calculated 1% scaled sensitivities (S0.01) using a forward–difference approximation (Hill, 1998). A 1% scaled sensitivity indicates the change in the prediction produced by a 1% change in the parameter value. Each S0.01 was calculated as

Formula 16[16]
where c is the unperturbed parameter value, {Delta}c is the perturbation (1% of c), and y(c) and y(c + {Delta}c) are the simulated values. The sensitivities were calculated using the simulated Cd profiles and plant Cd contents of the year 1996.

Model Parameterization
Table 1 gives an overview of the vertical distribution of soil pH and OC content within the WIA. At the non-Rigosol sites, the variability of pH increases with increasing depth. The OC content decreases continuously with depth. The Rigosols are characterized by lower pH values. Particularly between 0.3 and 0.6 m, the average soil pH is up to 0.4 units lower than at the non-Rigosol sites. In these layers the OC contents are also very different from those of the non-Rigosols. The average OC content in layers 0.4 to 0.5 m and 0.5 to 0.6 m is about 2 times higher than at the non-Rigosol sites.


View this table:
[in this window]
[in a new window]
  		Table 1. Vertical distribution of pH and organic carbon (OC) content in the soils of the wastewater irrigation area of Braunschweig.

 
Distribution of soil pH and OC contents were checked for normality and lognormality with the Kolmogorov–Smirnov (KS) test modified by Lilliefors (1967). In each layer down to a depth of 1 m, OC contents were lognormally distributed at {alpha} = 0.05. Soil pH values of the Ap horizons were normally distributed ({alpha} = 0.05). Below the Ap horizon, soil pH values showed a slightly left-skewed distribution (Fig. 3 ). The coefficient of skewness decreases from 0.03 in the 0.3 to 0.4 m layer to –1.22 in the 1.1 to 1.2 m layer. Because of the left-skewness the KS test ({alpha} = 0.05) rejected the null hypothesis of both normal and lognormal distribution. The model, however, assumes normally distributed pH values in all soil layers. The consequences for the validity of simulation results are discussed below.


Figure 3
View larger version (22K):
[in this window]
[in a new window]
  		Fig. 3. Histograms of soil pH in (A) Ap horizons and (B) 0.9–1.0 m layers. The curves are fitted normal distributions.

 
Table 2 compiles the model input parameters and variables with regard to weather, solute transport, and plant uptake. For non-Rigosols, simulations were run with a crop rotation of sugar beet, wheat, wheat, potato, wheat, and wheat. In hindcast simulations, Cd taken up by plants was completely returned to the Ap horizon at the end of each year. In forecast simulations, only Cd of the unprocessed plant parts remained at the site and were returned to the Ap horizon. Rigosol sites were assumed to be cropped with asparagus for 10 yr (Wonneberger et al., 2004) followed by 20 yr of recovery (Krug et al., 2002) during which it is cropped with the rotation mentioned above. After the recovery phase the crop rotation starts again with a 10-yr period of asparagus cropping. In 1996, about 6% of the area of the WIA was cropped with asparagus. According to interviews with employees of the BMWA, the average irrigation at asparagus sites is 100 mm yr–1. Because of the lack of experimental data, the uptake of Cd by asparagus was ignored in the model.


View this table:
[in this window]
[in a new window]
  		Table 2. Model input parameters and variables for the SEFAH model.

 
The parameters of the PTF 1 in Eq. [11] were derived by estimating the coefficients of log-transformed extended Freundlich isotherms:

Formula 17[17]
The coefficients were estimated by multiple linear regression using 616 complete data sets of S, C, OC content, and soil pH. In PTF 2 (Eq. [12]), we estimated for each soil monolith a site-specific k* value. In this case, multiple linear regression was performed using the "blocking" procedure (Draper and Smith, 1966). Equation [17] was expanded by N dummy variables Zi and site-specific coefficients {xi}i.

Formula 18[18]
If the observation is from the i-th soil monolith Zi = 1 otherwise Zi = 0. After multiple regression, the intrinsic Freundlich coefficient of the i-th soil monolith may be calculated as:

Formula 19[19]
The estimated parameters of the two PTFs are given in Table 3. The high modeling efficiencies show that both PTFs allow an acceptable estimation of Cd sorption.


View this table:
[in this window]
[in a new window]
  		Table 3. Parameters of pedotransfer functions (PTF) used in this study. The parameters were estimated by multiple linear regression (N = 616).

 
For 1980 to 1996, data of Cd loads to the WIA were available from routine measurement of the BMWA. Before 1980, however, Cd loads to the WIA were not measured. Thus, for 1957 to 1979 the Cd input history was reconstructed from our measurements (Fig. 4 ). We used a simple mass-balance approach taking into account measured Cd loads (see Eq. [13]) and different operation periods of districts as well as the sludge transfer between districts from 1975 to 1979. Because of the negligible effect on the mass-balance, the Cd flux into groundwater and Cd output via harvest were ignored.


Figure 4
View larger version (18K):
[in this window]
[in a new window]
  		Fig. 4. Reconstructed cadmium input to the districts of the wastewater irrigation area of Braunschweig.

 
The lower boundary of each soil column was set to the depth of the groundwater table. The regional distribution of groundwater table depth was derived by overlaying a digital terrain model with a groundwater elevation map (Ingwersen, 2001). Groundwater elevation had been interpolated from levels of 26 groundwater-monitoring wells supplied by the BMWA. The groundwater elevation map was produced by Universal Kriging with a linear trend removal. The average depth of interpolated water table was 1.8 m with a coefficient of variation of 40%. When the depth of the groundwater table was below the sampling depth, the sorption characteristics of the unsampled soil zone were approximated by that of the deepest sampled layer. In GMs as well as in the MC 2 approach, soil columns were georeferenced and the depth of groundwater table were taken directly from the map. In the MC 1 approach soil columns were not georeferenced. Here, for each soil column the depth of groundwater table was randomly drawn from the empirical distribution obtained from the map of groundwater table.

From the late 1960s to the early 1980s, in Lower Saxony the plowing depths increased on average by 0.1 m (Nieder et al., 1993). Since then, they have been kept more or less constant. Based on this information we assumed that in the WIA the plowing depth increased linearly from 0.2 to 0.3 m from 1967 to 1984.

The amounts of soil exchanged by deep plowing were estimated from the measured OC profile. First, the OC profile before deep plowing was approximately reconstructed by assuming that the highest OC content observed below the Ap horizon is the former OC content of the Ap horizon, OC0,old. The former OC content of layers below the Ap horizon OCi,old was set to the OC content measured in the 0.9 to 1.0 m layer. Under these assumptions the soil exchange {alpha}i (kg m–2) can be calculated by

Formula 20[20]
At each Rigosol site, the year of deep plowing was obtained from the farmers.


   RESULTS AND DISCUSSION
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS AND DISCUSSION
CONCLUSIONS
 REFERENCES
 
Hindcast Simulation for Non-Rigosols
For the non-Rigosols, Fig. 5 shows observed and simulated Cd profiles in 1996, after up to 40 yr of wastewater irrigation. About 16% of the EDTA–extractable Cd was detected in the subsoil. Five percent of the EDTA–extractable Cd was found below 0.5 m indicating that, on average, displacement has not progressed much deeper. Streck and Richter (1997a) also used this 5% threshold value in their study to characterize Cd displacement. The authors investigated field-scale displacement of Cd through a soil within district 2 of the WIA. They found 5% of EDTA–extractable Cd below 0.7 m indicating that at the study site of Streck and Richter (1997a) Cd transport was distinctly faster than on the regional average. The faster transport is probably due to a lower pH and OC content of the investigated soil (Ap horizon: pH = 5.25, OC content = 0.72% by weight).


Figure 5
View larger version (26K):
[in this window]
[in a new window]
  		Fig. 5. Observed and simulated average Cd concentration profiles at non-Rigosol sites. (A) solution phase concentrations and (B) EDTA–extractable content.

 
Despite uncertainties of model input data and model simplifications such as assuming stationary soil pH, the agreement between observed and simulated data is good (Fig. 5). The differences between modeling approaches are moderate (Table 4). All simulations overestimate the overall Cd displacement by about 20%. This is mainly caused by the disagreement between model and observation in the 0.3 to 0.4 m layer. With the PTF with horizontally variable k*, the GM performs slightly better than with constant k*. In the MC simulations, taking into account the horizontal correlation between soil pH, OC content, and Cd loads improved the performance. The EF of the MC 2 simulation with regard to the solution phase concentration is 0.12 points higher than that of the MC 1 approach (Table 4). The output of MC 2 simulations is georeferenced. Hence, different from the MC 1 approach, simulation results may be used, for example, to delineate areas that are particularly prone to leaching or soil-plant transfer of Cd.


View this table:
[in this window]
[in a new window]
  		Table 4. Modeling efficiencies (EF) and displacement indices (DI) for simulations performed either with grid models (GM) or Monte-Carlo (MC) methods. Modeling efficiency and DI are defined in Eq. [14] and Eq. [15], respectively. The symbol C denotes the solution phase concentration and CT the EDTA–extractable content. The symbol N stands for the number of local simulations.

 
Cd displacement is highly variable within the WIA (Fig. 6 ). At some non-Rigosol sites, even after 40 yr of wastewater irrigation, increased Cd contents are limited to the Ap horizon and the adjacent subsoil layer. At other non-Rigosol sites, Cd has been displaced to below 1.2 m.


Figure 6
View larger version (15K):
[in this window]
[in a new window]
  		Fig. 6. Cadmium concentrations at non-Rigosol sites after up to 40 yr of wastewater irrigation. (A) All data and (B, C) two selected profiles.

 
The model presumes normally distributed pH values in all layers, although in subsoil measured pH values were slightly left-skewed. For the simulations presented, the simplifying assumption of normally distributed pH values in subsoil appears to have a negligible effect on the model outcome. In subsoil, simulated solution phase concentrations of GM and MC method are very similar. However, applying the model for a long-term forecast of the Cd seepage concentration, the MC method would probably tend to underestimate the transport velocity of Cd.

Hindcast Simulation for Rigosols
The displacement pattern of the Rigosols is completely different from that of the non-Rigosols (Fig. 7 ). Whereas at the non-Rigosol sites Cd concentrations continuously decrease with depth, the Rigosol profiles show a concentration peak between the 0.3 and 0.6 m depth. These are the layers affected by deep plowing. About 65% of the EDTA–extractable Cd was recovered in the subsoil, and 5% of the EDTA–extractable Cd was found below 0.7 m. Measured concentration profiles are matched fairly well by the GM simulations, which suggests that the calculation of soil exchanged by deep plowing from the measured OC profile is suited to model the transfer of Cd between top and subsoil.


Figure 7
View larger version (29K):
[in this window]
[in a new window]
  		Fig. 7. Observed and simulated average Cd concentration profiles of Rigosols. (A) solution phase concentrations and (B) EDTA–extractable content.

 
Similar to the non-Rigosol sites, the use of the PTF with horizontally variable k* produces better results than the simulation with a constant k*. The EF of the solution phase concentration increases from 0.47 to 0.68 (Table 4). This finding indicates that additional factors, which are not explicitly considered in the PTF, such as contents of iron oxide or clay, might have a significant impact on heavy metal sorption. The MC simulations reproduce the EDTA–extractable Cd content but systematically underestimate the solution phase concentration. The reason for the underestimation is that the random numbers were obtained from the vector of means and covariance matrix of the non-Rigosol sites. Subsoil pH values of the Rigosols, however, are on average, about 0.4 units lower and more variable than those of non-Rigosols (Table 1). Consequently, although EDTA–extractable Cd is fairly well modeled, the solution phase concentration is underestimated. This explanation is supported by the results of the sensitivity analysis (Table 5). With regard to leaching, the most sensitive variable/parameter are soil pH and the PTF coefficient a. Because both parameters are perfectly linearly correlated, they have the same sensitivity (see Eq. [11] and Eq. [12]). The sensitivity of the intrinsic Freundlich coefficient k* is in the same order of magnitude. The plant uptake parameters (kP, {eta}, and RAp), the average saturation deficit of the air, and the volumetric water content, have the lowest sensitivity with regard to leaching. Concerning plant uptake, the most sensitive model inputs are soil pH and the empirical PTF parameter a. This finding agrees with results of many previous studies that have shown that soil pH plays a key role in the plant uptake of Cd (Grant et al., 1998; Welp and Brümmer, 1999).


View this table:
[in this window]
[in a new window]
  		Table 5. One-percent scaled sensitivities of selected input parameters of the SEFAH model.

 
The estimated Freundlich coefficients used in the GM 1 and GM 2 simulations ranged from 28 to 1020 µg1–m Lm kg–1. That is a typical range for sandy soils of northwest Germany. Springob and Böttcher (1998) studied Cd sorption isotherms of sandy soils of the Fuhrberg catchment, which is located about 90 km northwest of the WIA. The investigated topsoils had a comparable soil pH range but the OC content was on average 2.4% by mass, which is moderately higher than within the WIA. Springob and Böttcher (1998) measured Freundlich coefficients in the range between 36 and 1275 µg1–m Lm kg–1.

In all simulations we assumed that the Freundlich exponent m of the sorption isotherm is constant. This assumption has been used in many previous studies (van der Zee and van Riemsdijk, 1987; Boekhold and van der Zee, 1992; Streck and Richter, 1997a, 1997b; Elzinga et al., 1999). Chardon (1984) reported that m was relatively constant over a wide range of experimental conditions. In contrast, Buchter et al. (1989) found a close correlation between soil pH and m values. Springob and Böttcher (1998) also observed a significant linear correlation between soil pH and m, although the coefficient of correlation was low (r2 = 0.3).

Scenario Simulations

Scenario 1: How will Cd concentration in seepage develop if the Cd loads are kept at the level of the year 1996?
Simulations were run from 1996 to 2246. We assumed that during this period of time, Cd loads by wastewater irrigation were at the level of year 1996 (1.1 g ha yr–1). The total atmospheric deposition was set to 0.9 g ha–1 yr–1 (Dämmgen et al., 1995) and the simulations were run with the MC 2 approach. Climatic input variables were set to the averages of the period from 1957 to 1996. Crop yields were kept constant at the levels of year 1996.

Figure 8 shows the average Cd profiles in soil for years 1996 and 2046. The solution phase concentration in the Ap horizon decreases from 1.7 µg L–1 in 1996 to about 0.9 µg L–1 in 2046, and Cd is displaced to a depth below 1.2 m. The Cd concentration in the seepage water gradually increases from about the background level in 1996 to a maximum concentration of 0.65 µg L–1 around year 2100 (Fig. 9 ). This maximum concentration is distinctly lower than the WHO guideline value for drinking water (3 µg L–1; WHO, 1993) or the limit of the German drinking water ordinance (5 µg L–1; TrinkwV, 2003). Moreover, the forecasted maximum concentration in seepage is about three times lower than the present average solution phase concentration in Ap horizons.


Figure 8
View larger version (17K):
[in this window]
[in a new window]
  		Fig. 8. Simulated regionally-averaged solution phase concentrations and EDTA–extractable Cd contents in 1996 and 2046. From 1996 to 2046, Cd load was kept at the level of 1996. Simulations were performed with the MC 2 approach.

 

Figure 9
View larger version (14K):
[in this window]
[in a new window]
  		Fig. 9. Predicted average Cd concentrations in seepage water at the depth of water table. From 1996 to 2246, Cd load was kept at the level of 1996. Simulations were performed with the MC 2 approach.

 

Scenario 2: What would have happened if Cd loads had not been reduced in the mid eighties?
In the mid 1980s the monitoring of illicit emissions to the municipal sewage canal system had been enforced and the sewage plant was equipped with a sewage press. Due to these two measures, the Cd loads into the WIA have been markedly reduced. Scenario 2 assumes that these countermeasures were not taken in the mid 1980s. It is assumed that from 1985 to 1996 the Cd load of each district remained on the average level between 1980 and 1984. Scenario 2 was computed using the GM 1 approach.

The results of the simulations are compared with those of a hindcast simulation in Fig. 10 . The fluctuation of simulated Cd grain content is caused by different annual saturation deficits of the air. In a mass–flow-governed plant uptake model the saturation deficit of the air affects the transpiration and thus the uptake of Cd. Due to the reduced Cd loads to the WIA in the hindcast simulation, the regionally averaged Cd grain content levels off in the 1990s and approaches the regulatory limit of 0.24 mg kg–1 but does not go beyond this content. However, locally the regulatory limit is exceeded. In the period from 1990 to 1996 at 36 of 153 sites, the simulated Cd grain content was, on average, higher than 0.24 mg kg–1. Simulated Cd grain contents of Scenario 2 exceed the regulatory limit after 1989. The outcome of the scenario simulation demonstrates that the decision of the BMWA in the early 1980s to reduce Cd inputs was useful.


Figure 10
View larger version (21K):
[in this window]
[in a new window]
  		Fig. 10. Cadmium content in wheat grain for Scenario 2 and a hindcast simulation. Scenario 2 assumes that Cd loads were not reduced in the mid 1980s. In the hindcast simulation the reconstructed ("real") Cd input history was used.

 
The simulated Cd contents of sugar beet (hypocotyl) and potato (tuber) showed similar results. Because of the lack of experimental data no statement about the Cd content of asparagus, the only vegetable that has been cropped in the WIA, is possible.


   CONCLUSIONS
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS AND DISCUSSION
 CONCLUSIONS
REFERENCES
 
In this study we presented data on regional-scale displacement of Cd in agricultural soils after up to 40 yr of wastewater irrigation. Despite favorable leaching conditions (high water flux densities, sandy soils, and low organic carbon content) Cd migration has progressed on average not much deeper than about the 0.5-m depth at the non-Rigosol sites. The situation is very different at sites that had been cropped with asparagus. Here, deep plowing is an important displacement process, which instantaneously exchanges heavy metals between topsoil and subsoil and thus significantly accelerates metal migration. Under the conditions of the WIA both GM and MC methods were found to be suited for upscaling Cd migration in soils from the local to the regional scale. The adequate description of the spatial variability of both Cd input and soil properties that govern Cd sorption is a prerequisite to obtain sound simulation results. Therefore, the GM 2 (grid model with site-specific intrinsic Freundlich coefficients) and, if sufficient soil data are available, the MC 2 (quasi-3D Monte-Carlo simulations) methods are favorable for regional Cd transport modeling. Our results show that the irrigation of soils with municipal wastewater over a period of up to 40 yr has caused Cd accumulation in soil. Currently, the soil-plant pathway is most strongly affected. The breakthrough of Cd to the groundwater is damped and is retarded for many decades. Our study demonstrates that after a site-specific parameterization and validation, modeling can be a useful tool to assess and forecast the effect of land management practices on the environmental fate of Cd.


   ACKNOWLEDGMENTS
 
We thank T. Eggers, F. Seeßelberg, and H. Blickwede of the BMWA for their cooperation. We thank B. Heine, C. Marheineke, D. Müntner, and A. Küsters for laboratory assistance. The authors thank two anonymous reviewers for critical comments and suggestions, which helped to improve the manuscript. This research was funded by the Deutsche Forschungsgemeinschaft (German Research Foundation, Grant no. RI 269/36-1 and RI 269/36-2).


   REFERENCES
TOP
 ABSTRACT
 INTRODUCTION
 MATERIALS AND METHODS
 MODELING
 RESULTS AND DISCUSSION
 CONCLUSIONS
 REFERENCES
 




This article has been cited by other articles:


Home page
 J. Environ. Qual.Home page
T. W. Biggs and B. Jiang
Soil Salinity and Exchangeable Cations in a Wastewater Irrigated Area, India
J. Environ. Qual., March 25, 2009; 38(3): 887 - 896.
[Abstract] [Full Text] [PDF]

Home page
 Vadose Zone JHome page
S. Altfelder, W. H. M. Duijnisveld, T. Streck, G. Meyenburg, and J. Utermann
Quantifying the Influence of Uncertainty and Variability on Groundwater Risk Assessment for Trace Elements
Vadose Zone J., August 23, 2007; 6(3): 668 - 678.
[Abstract] [Full Text] [PDF]

This Article
Right arrow Abstract Freely available
Right arrow Figures Only
Right arrow Full Text (PDF) Free
Right arrow Alert me when this article is cited
Right arrow Alert me if a correction is posted
Services
Right arrow Similar articles in this journal
Right arrow Similar articles in Web of Science
Right arrow Similar articles in PubMed
Right arrow Alert me to new issues of the journal
Right arrow Download to citation manager
Citing Articles
Right arrow Citing Articles via HighWire
Right arrow Citing Articles via Web of Science (9)
Right arrow Citing Articles via Google Scholar
Google Scholar
Right arrow Articles by Ingwersen, J.
Right arrow Articles by Streck, T.
Right arrow Search for Related Content
PubMed
Right arrow PubMed Citation
Right arrow Articles by Ingwersen, J.
Right arrow Articles by Streck, T.
Agricola
Right arrow Articles by Ingwersen, J.
Right arrow Articles by Streck, T.
Related Collections
Right arrow Municipal Wastes
Right arrow Spatial Variability
Right arrow Upscaling
Right arrow Solute Transport Models
Right arrow Heavy Metals

HOME HELP FEEDBACK SUBSCRIPTIONS ARCHIVE SEARCH TABLE OF CONTENTS
The SCI Journals Agronomy Journal Crop Science
Journal of Natural Resources
and Life Sciences Education		 Vadose Zone Journal
Soil Science Society of America Journal Journal of Plant Registrations The Plant Genome