Modelling assessment of regional groundwater contamination due to historic smelter emissions of heavy metals

Available online at www.sciencedirect.com

Journal of Contaminant Hydrology 96 (2008) 48 – 68 www.elsevier.com/locate/jconhyd

Modelling assessment of regional groundwater contamination due to historic smelter emissions of heavy metals Bas van der Grift ⁎, Jasper Griffioen TNO Geological Survey of the Netherlands P.O. Box 80.015, 3508 TA Utrecht, The Netherlands Received 22 March 2005; received in revised form 5 October 2007; accepted 10 October 2007 Available online 18 October 2007

Abstract Historic emissions from ore smelters typically cause regional soil contamination. We developed a modelling approach to assess the impact of such contamination on groundwater and surface water load, coupling unsaturated zone leaching modelling with 3D groundwater transport modelling. Both historic and predictive modelling were performed, using a mass balance approach for three different catchments in the vicinity of three smelters. The catchments differ in their hydrology and geochemistry. The historic modelling results indicate that leaching to groundwater is spatially very heterogeneous due to variation in soil characteristics, in particular soil pH. In the saturated zone, cadmium is becoming strongly retarded due to strong sorption at neutral pH, even though the reactivity of the sandy sediments is low. A comparison between two datasets (from 1990 to 2002) on shallow groundwater and modelled concentrations provided a useful verification on the level of statistics of “homogeneous areas” (areas with comparable land use, soil type and geohydrological situation) instead of comparison at individual locations. While at individual locations observations and the model varies up to two orders of magnitude, for homogeneous areas, medians and ranges of measured concentrations and the model results are similar. A sensitivity analysis on metal input loads, groundwater composition and sediment geochemistry reveals that the best available information scenario based on the median value of input parameters for the model predicts the range in observed concentrations very well. However, the model results are sensitive to the sediment contents of the reactive components (organic matter, clay minerals and iron oxides). Uncertainty in metal input loads and groundwater chemistry are of lesser importance. Predictive modelling reveals a remarkable difference in geochemical and hydrological controls on subsurface metal transport at catchment-scale. Whether the surface water load will peak within a few decades or continue to increase until after 2050 depends on the dominant land use functions in the areas, their hydrology and geochemical build-up. © 2007 Elsevier B.V. All rights reserved. Keywords: Soil; Groundwater flow; Metal transport; Modelling; Catchment-scale

1. Introduction There are various reasons why the soil becomes contaminated with heavy metals and these leach to ⁎ Corresponding author. Tel.: +31 30 2564720. E-mail address: [email protected] (B. van der Grift). 0169-7722/$ - see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.jconhyd.2007.10.001

groundwater: impact of acid rain (Edmunds et al., 1992; Frei et al., 2000; Kjøller et al., 2004), long-term manureand fertilizer-borne input to agricultural soils (McBride et al., 1997; Keller et al., 2001; Miller et al., 2003; Xue et al., 2003), ore mining (Brown et al., 1998), dispersion and deposition from smelters (Wilkens and Loch, 1997; Fernandeze-Turiel et al., 2001; Seuntjens et al., 2002;

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Burt et al., 2003), forestation of cultivated soils (Römkens, 1998; Strobel et al., 2001), and diffuse atmospheric deposition (Tiktak, 1999). The fate of the contaminating heavy metals in the subsoil needs to be predicted, so that soil and water can be managed properly. For example, when the contaminant load is unacceptably high, it may be necessary to take measures to safeguard the abstraction of drinking water or the discharge to surface water. Reactive transport models are useful tools for the management of contaminated areas. Few studies have addressed the reactive transport of heavy metals in the field. For these studies, a clear distinction can be made between studies that address transport within the unsaturated zone, where crop uptake is relevant and transport is dominantly vertical, and studies of transport in the saturated zone where transport is three-dimensional. In the case of the unsaturated zone, Wilkens (1995) modelled transport of Cd and Zn in several field plots using a finite-difference approach with vertically variable distribution coefficients. Tiktak (1999) used a mass budget model to model accumulation of Cd in the top horizon in a regional context, as well as a onedimensional reactive transport model for the unsaturated zone of three field plots. Seuntjens et al. (2002) did a sensitivity analysis in order to rank the importance of spatially variable water flow and solute transport parameters affecting field-scale cadmium flux in a layered sandy soil. Keller et al. (2001) applied an empirical stochastic model to the plowed layer of agricultural soils, which considers heavy metal inputs from agricultural management and outputs from crop removal and leaching on a regional-scale. Brown et al. (1998) and Kjøller et al. (2004) modelled reactive transport of trace metals in the saturated zone one-dimensionally using multicomponent geochemical formulations for the hydrogeochemical reactions. Walter et al. (1994) used a two-dimensional vertical crosssection addressing both the major compounds as well as trace metals by taking into account multicomponent geochemical formulations for the hydrogeochemistry. Kent et al. (2000) and Curtis et al. (2006) applied twodimensional reactive groundwater transport models, Kent et al. (2000) incorporated semi-empirical surface complexation models to describe sorption, Curtis et al. (2006) used a semi-mechanistic surface complexation model. So far, regional-scale metal transport simulations have not been conducted. The objective of the study was to develop a regional approach, using reactive transport models, that addresses geochemical and hydrological

49

controls on subsurface transport of the trace metals Cd and Zn. The study was developed for three different catchments in the Kempen region in the Netherlands (Fig. 1), where the soil had been severely contaminated by atmospheric emissions from three Zn smelters. With this study we want to assess the history and to predict the future in terms of: (1) temporal and spatial variable leaching of metals from the unsaturated zone to groundwater; (2) developments in Cd and Zn concentrations in shallow and deeper groundwater, and (3) the metal loads of the surface water drainage network in relation to the hydrology and geology. 2. Site characterization 2.1. Study area In the study area (see Fig. 1), which is near the Dutch–Belgian border, there is heavy metal contamination from three zinc-ore smelters less than 10 km apart. The zinc ore most often used in these smelters is sphalerite (ZnS), which contains a wide range of other metals, the most abundant being manganese, cadmium, copper, arsenic, tin, gallium, antimony, and thallium (Levinson, 1977). From 1880 to 1974 the chimneys of these smelters emitted oxides of heavy metals that reached the soil either by dry deposition or with rainfall. The plants switched to an electrolytic process in 1974; since then, atmospheric emissions have diminished drastically. The historic input of heavy metals on the soil through atmospheric deposition has resulted in an excessive accumulation in the topsoil in the Kempen region, accompanied by increased leaching to the groundwater (Harmsen, 1977; Boekhold, 1992; Tiktak, 1999; Seuntjens et al., 2002, Sonke et al., 2002; CSO, 2001). The soils are mainly Spodosols and Earth soil which have developed in poor, fluvio-eolian, Pleistocene sands, and which are vulnerable to leaching because of the acidifying conditions and low contents of organic matter and clay (De Bakker and Schelling, 1989; Wilkens and Loch, 1997). In the brook valleys that intersect the sandy regions the soil varies from peaty via loamy to sandy. The Kempen region is flat lowland with surface levels decreasing from south to north, from about 40 m to 20 m above sea level. The difference in surface level between the brook valleys and the higher parts of the area is in the order of a few to 10 m. The climate of the area is humid and temperate with mean July and January temperatures of 19 and 4 °C, respectively, and annual precipitation averaging 700 mm. Geologically, the area can be divided into

50

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Fig. 1. Map showing the location of the Kempen area in the Netherlands, the Beekloop–Keersop, Buulder Aa and Tungelroijsche Beek catchments, the geological profile along the dotted line with the important Feldbiss Fault and Peel Boundary Fault, the locations of the drillings for geochemical analysis and the interpolated cadmium content in the topsoil for the year 1995 (CSO, 2001).

two: the southwestern part is the Kempen High bounded by the Feldbiss Fault with the Roer Valley Graben. This active tectonic subsidence area contains the thickest terrestrial fine-grained deposits with most complete record in the Netherlands (Schokker, 2003; De Mulder et al., 2003). The Late Pleistocene Boxtel Formation that

covers the Roer Valley Graben is a fine-grained fluvial and eolian sandy deposit intercalated with heterogeneous layers of loam and peat, and has a maximum thickness of 35 m. In the Kempen High, the top unit is the middle Pleistocene Sterksel Formation, which consists mainly of medium and coarse sand and some

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

gravel deposits from the rivers Rhine and Meuse. Deeper geological formations in the area with increasing age are: the lower and middle Pleistocene Stramproy Formation of eolian and local terrestrial origin, the lower Pleistocene Waarle Formation, and the Late Miocene and Pliocene Kiezeloöliet Formation. The latter two are deposits from predecessors of the river Rhine. The Roer Valley Graben is bounded on the northeast by the Peel Boundary Fault which marks the transition to the Peel Horst (which lies outside the study area). This paper focuses on three catchments that differ in their hydrogeological and geochemical characteristics: the Beekloop–Keersop on the Kempen High and Buulder Aa and Tungelroijsche Beek catchments in the Roer Valley Graben. The upper streams of these catchments are located in immediate vicinity of the smelters. The three catchments encompass an area that is over 330 km2 in size. 2.2. Groundwater quality In the period 2003–2004 groundwater in the Kempen area was sampled and analyzed. At 62 sampling fields the uppermost meter of groundwater was sampled using temporary monitoring wells. At each of the sampling fields, four wells were drilled in a square plot of about 100 ⁎ 100 m to a depth of 1 m below the actual groundwater level. These groundwater levels vary between 0.2 and 4.5 m depth. Groundwater samples of the four wells at one location were collected, mixed and analyzed for macro-chemistry (pH, electrical conductivity, alkalinity, Cl, NO3, SO4, PO4, NH4, Na, K, Mg, Ca, Al, Fe), heavy metals (As, Cd, Ni, Cr, Cu, Pb, Zn) and DOC using routine techniques (ICP-MS, Ion chromatography, Element Analyser, TOC Analyzer). The monitoring wells were removed after sampling. Additionally, 51 samples were taken from existing permanent monitoring wells; these ranged in depth from 1.5 to 30 m. These data were augmented with data from a monitoring program performed in 1990 in the Kempen area, in which the uppermost groundwater had been sampled and analyzed for Cd, Zn and pH (Tauw, 1991). We used the major groundwater chemistry data as input variables for the Cd and Zn sorption isotherms in groundwater environment, as explained later. The groundwater data on Cd and Zn were not input into the model but were instead used for model evaluation. Here, the groundwater chemistry data are grouped according to the concept of “areas of homogeneous groundwater composition”, i.e., areas within which the chemical composition of groundwater varies less than it does between such areas (Broers, 2002). Therefore, dominant

51

controls on the chemical groundwater composition had to be assigned. In the Netherlands, these dominant controls are input loads, which are directly related to land use, and hydrologic regimes (Broers, 2002). We established the homogeneous areas for two land uses (agriculture and nature areas) and three hydrologic regimes (infiltration, intermediate and discharge areas). Nature areas are areas without an intensive agricultural or residential function and exist mainly of pine forest and moor land. About 25% of the land use in the study area is nature, rising to 29% within 5 km of the smelters. Infiltration areas are higher-lying parts of the landscape between the brooks. The water table is deep (on average, 2–6 m) and there are no drainage ditches. Discharge areas correspond with the brook valleys in the lower parts of the landscape; here, groundwater exfiltrates and the water table is close to the surface. In between are the intermediate areas, which are characterized by the presence of locally recharged flow systems that discharge into the local ditches and small brooks, and in which the water table is at about 1–3 m depth (Broers and Van der Grift, 2004; Broers, 2004). 2.3. Sediment geochemistry Data about the aquifer geochemistry was obtained from 11 deep drillings to about 100 m depth and 6 shallow drillings in the Boxtel Formation in the Roer Valley Graben (see Fig. 1). From the deep drillings, a number of samples were selected for analysis on the basis of the lithostratigraphy. Thus, within every geological formation we selected a series of samples at 10 m intervals, and also at the transition between two formations. From the Boxtel Formation drillings, samples were taken every meter (Schokker, 2003). In total, 298 sediment samples were analyzed for bulk geochemistry (X-ray Fluorescence Spectrometry and Thermogravimetric Analyzer) to determine the contents of reactive components: organic matter, total Al2O3 (as proxy for clay fraction), carbonates, Stotal (as proxy for pyrite), Fereactive (as proxy for iron oxides and hydroxides and siderite). The results for the reactive components were averaged per geologic unit and used as input variables for sorption affinity of Cd and Zn in the saturated zone model (Table 1). 3. Modelling approach For this study, we developed a modelling approach coupling an unsaturated zone leaching model with a groundwater flow and transport model (Fig. 2). With the unsaturated zone model we calculate the spatially and

52

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Table 1 Median values of reactive component in the geological formations derived from analyzed samples from 17 drillings in the Kempen area

Boxtel Formation Beegden Formation Sterksel Formation Stramproy Formation Waalre Formation Breda Formation

Origin

N

Al2O3 (%)

Stotal (%)

Fereactive (%)

OM (%)

Carbonates (%)

Local and eolian Fluvial Fluvial Local and eolian Fluvial Marine

204 5 45 31 1 12

4.24 1.96 3.33 5.15 6.73 4.15

0.11 0.03 0.10 0.18 0.09 0.80

0.03 0.00 0.07 0.41 0.00 1.63

0.94 0.25 0.34 2.55 1.33 1.55

0.30 0.24 1.34 2.29 0.75 1.53

temporally variable leaching of heavy metals to groundwater. We subsequently applied a fully distributed threedimensional reactive transport model to the saturated zone to assess subsurface transport and surface water load. 3.1. Unsaturated zone model HYDRUS-1D was used to model the leaching of cadmium and zinc to the water table for two periods: 1880 – present and present – 2050. The first period served for historic modelling and verification purposes, and the second for prediction and assessment of vulnerability. Breakthrough curves of Zn and Cd were calculated for unique combinations, considering 76 soil types, 2 land use functions, 10 classes of depths of water tables and 11 classes of metal input loads. In all there were 16,720 unique combinations of these

variables. The modelling was performed for yearly averaged conditions of the precipitation excess and groundwater level, and free drainage as bottom boundary condition. Soil chemical and physical parameters are needed to model metal leaching from topsoil to groundwater. Note that the sandy soils in the Kempen area are mainly aerobic and mildly acid (pH 4–6.5) and therefore precipitation of Cd or Zn as carbonates or oxides plays no role in this region. We, therefore, described leaching by sorption together with advective/dispersive transport. In this area, sorption on organic matter and clay minerals must be key processes controlling reactive transport in the unsaturated zone, because oxide contents are relatively low (Wilkens and Loch, 1997). We used an empirical-partitioning model from Römkens et al. (2004) to obtain adsorption isotherms for the individual

Fig. 2. Structure of the coupled model.

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

soils horizons in the area. The individual soil types were derived from Steur and Heijink (1991). The Freundlichtype partition equations have been derived from multiple-linear regression analysis of 1450 soil–soil solution data of 60 very diverse Dutch soils, including the soils in the study area: h i 0;54 log QCd =CCd ¼ 5:01 þ 0:27 log ½k clay þ 0:65 log ½k SOM þ 0:29 pH R2 ¼ 0:77

ð1Þ

seðY Þ ¼ 0:37

h i 0;78 log QZn =CZn ¼ 4:96 þ 0:36 log ½k clay þ 0:51 log ½k SOM þ0:52 pH R2 ¼ 0:85

ð2Þ

seðY Þ ¼ 0:41

where QMe is the metal content in soil (mol/kg) using 0.43 N HNO3 destruction, CMe is the total metal concentration in the soil solution (mol/m3), SOM is sedimentary organic matter and pH is the acidity of the soil (0.01 M CaCl2 extraction). Input parameters for the empirical-partitioning equations are pH, organic matter and clay content. We used national data on the physical– chemical characterization of the soil types in the Netherlands (De Vries, 1999) as input. This characterization provides the characteristics of the most common soil horizons to a depth of 120 cm, including the soils in the study area. Groundwater levels were derived from the regional groundwater flow model (Buma et al., 2002) (see Section 2.3). The Van Genuchten parameters for modelling the unsaturated water flow in the 76 soil profiles were calculated with pedo-transfer functions from the soil texture (clay, loam, organic matter content and the median of the sand fraction) and bulk density (Wösten et al., 2001). These parameters were also available from the physical–chemical characterization of Dutch soil types (De Vries, 1999). A precipitation excess of 259 mm/year for nature areas and 277 mm/year for agricultural land was used for the unsaturated zone model, which takes annual precipitation and land use dependent evapotranspiration into account. 3.2. Soil surface load The following three sources of Cd and Zn were considered: (1) atmospheric deposition, (2) fertilizers

53

and (3) animal manure. Rozemeijer (2002) made a reconstruction of the historic atmospheric deposition data for Cd and Zn originating from the smelters across the Kempen region. Rozemeijer (2002) reconstructed the time-dependent historic atmospheric deposition rate of Cd and Zn in the period 1880–1975 from present measured contents of these metals in forest soils samples at 19 locations in the study area with varying distance from the smelters, taking developments in historic zinc production into account. As it can be assumed that in the 20th century the yearly atmospheric emission and deposition of Cd and Zn increased concomitantly with the zinc production rate of the smelters, the temporal trend in annual deposition rate in the period 1880–1975, has been linearly correlated to the zinc production of the Budel smelter (Makaske et al., 1995). Rozemeijer (2002) used HYDRUS-1D (Šimůnek et al., 1998) to calibrate the level of the 1880–1975 atmospheric deposition to the analyzed cadmium and zinc contents in the soil. The metal content ranged from b0.2 mg/kg Cd to b 20 mg/kg Zn at locations over 30 km from the smelters to 5.3 mg/kg Cd and 733 mg/kg Zn at a location within two km from the Budel smelter. Non-linear Freundlich adsorption isotherm coefficients (KF) for the HYDRUS model were calculated from the pH, organic matter and clay content of the analyzed soil samples, as described above. The results of Rozemeijer (2002) indicated that the historic atmospheric deposition of the metals due to emission from the smelters decreased strongly with distance from the smelters. For cadmium, the modelled average atmospheric deposition at the 19 locations in the period 1880–1975 varied between 15 and 642 g/(ha yr) depending on the distance from the smelters. For zinc, the levels varied between 730 and 25000 g/(ha yr). Using the distance from the smelters and the predominant wind direction as input variables, we spatially interpolated the deposition at the 19 locations as calculated by Rozemeijer (2002) to a map of atmospheric deposition rates (Fig. 3). For the period 1975– 2005 regional values for atmospheric deposition of metals were used (RIVM/CBS, 2002). The load from agriculture sources was determined from data on local and national use of manure and fertilizer salts in the period 1950–2000 (Broers and Van der Grift, 2004). These data were augmented with data on the metal composition of manure and fertilizers (Hotsma 1997; Raven and Loeppert 1997; Moolenaar and Lexmond, 1998), in order to calculate the historic cadmium and zinc loads. Far less heavy metals originate from agriculture than from the historic atmospheric deposition in the vicinity of the smelters. The average

54

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Fig. 3. Average cadmium deposition in period 1880–1975 as used for the unsaturated zone model (in g/ha/year).

agricultural cadmium and zinc loads in the 1880–2000 period were 1.5 g/(ha yr) and 606 g/(ha yr), respectively. Peak levels around 1985 were 4.3 and 1830 g/(ha yr), respectively. The soil surface load of Cd and Zn expected in 2005– 2050 was calculated from data on the current values of atmospheric deposition of metals (RIVM/CBS, 2002) and the statutory norms for the use of manure and fertilizers in the Netherlands (NMP-RIVM, 2002). 3.3. Saturated zone groundwater flow model The saturated zone transport model was based on a MODFLOW finite-difference groundwater flow model (McDonald and Harbaugh, 1988) and an MT3DMS transport model (Zheng and Wang, 1999). By use of telescopic mesh refinement, detailed models for the three individual catchments were clipped out from an existing supraregional groundwater flow model (Buma et al., 2002) that covers over 13,000 km2 and thus extends over the entire province of Noord–Brabant (Fig. 1). This model has natural boundaries located in surrounding provinces and in Belgium. The supraregional model has 9 geohydrological layers covering the geological formations and contains over 9 ⁎ 106 grid cells. It has been calibrated on transmissivity and hydraulic layer resistance using a representer-based inverse method for groundwater flow (Valstar et al.,

2004). All groundwater abstractions above 1000 m3/yr were included in the model, plus their exact depths of abstraction. Hydraulic heads from the supraregional model were used as boundary conditions for the three catchment flow models. These models were adjusted to perform solute transport simulations for each individual catchment and have there boundaries outside the natural boundaries of the catchment. The Tungelroijsche Beek catchment does not cross the Belgian boarder and only small parts of the Beekloop–Keersop and Buulder Aa catchments are situated in Belgium. The catchment models contain only the first geohydrological layers of the original supraregional model to a depth of about 70 to 100 m, because Zn and Cd are not transported beyond these depths. In order to minimize numerical dispersion, the catchment models have 13 layers with increasing thickness from surface: 3 m for the upper 5 layers, 5 m for the next 4 layers to about 15 m for the deepest layers. Transmissivity of each layer of the new layers was assigned on basis of there depth towards the original layers and the conductivity and resistance of these layers. The horizontal grid discretization of the models is 100 ⁎ 100 m. The calculation time steps were in the order of 13 days, fulfilling the Courant condition. For example, the Buulder Aa catchment groundwater model contains 546,000 grid cells encompassing an area of 24 ⁎ 17.5 km.

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68 Table 2 Median values of groundwater composition in homogeneous areas at three depth intervals derived from the groundwater quality monitoring network in the area Clay minerals Cation-exchange capacity 10 meq/kg Amorphous Fe-hydroxide (Dzombak and Morel, 1990) Capacity, strongly selective 200 mmol/mol Fe Capacity, weakly selective 5.0 mmol/mol Fe Surface area 600 m2/g FeOOH Humic acid Sorption capacity, carboxyl 2.54 mol/kg humic acid Sorption capacity, phenol 3.75 mol/kg humic acid Heterogeneity 0.54

3.4. Saturated zone groundwater transport model We applied the transport model on grid cells inside the natural boundaries of catchment. Grid cells outside this natural boundary and grid cells inside the catchment boundary in Belgium were set to inactive. For retarding solutes like Cd and Zn neglecting the lateral crossboarder transport is not a serious limitation. Due to long groundwater travel times, metals that leached to the groundwater in Belgium will not reach the Netherlands groundwater and surface water in relevant amounts in the timeframe our model. We used results of the unsaturated zone HYDRUS model runs at the groundwater table as recharge concentrations for the reactive groundwater transport modelling of Cd and Zn. This output was averaged per 5 years to match input for the groundwater transport model. The transport of cadmium and zinc in groundwater was calculated for the period from 1950 to 2050; note that prior to 1950 no breakthrough at the water table was simulated with the unsaturated zone model. Initial groundwater concentrations for Cd and Zn in 1950 were set at zero. Sorption of Zn or Cd was described for each individual grid cell by general Freundlich isotherms type curves. However, the empirical Freundlich isotherms used for the unsaturated zone could not be used for the groundwater transport model. This unsaturated zone approach was developed to understand the relationship between the soil and the soil solution. This empirical model could not be extrapolated to sandy aquifer material, which contains significantly lower contents of sorbents. Therefore, we used Freundlich isotherm functions derived by a meta-model from Griffioen et al. (1998). This approach was specially developed for aquifer environments. These functions calculate the fraction of the sorbent occupied by the trace metal of interest as a function of pH. First, the

55

activity of the free trace metal (e.g. Cd2+) was calculated from a groundwater composition, taking into account the major inorganic aqueous complexes and complexation with dissolved organic acids (Griffioen et al., 1998). The median concentrations of macro-chemical components and the pH of the groundwater were input in the aqueous complexation and sorption models. These medians were derived per homogeneous area, as explained before. We distinguished four depth intervals when interpreting groundwater quality data: uppermost groundwater (0–6 m), moderately deep groundwater (6–15 m), deep groundwater (15–35 m) and very deep groundwater (N 35 m) and they were assumed to be constant over the period modelled. For the pH, which can be considered as the dominant factor controlling aqueous complexation and sorption of metals, this assumption is justified. As soil pH changes are buffered by several chemical weathering reactions, the excess atmospheric acid input of the last decades has not lead to lower groundwater pH values but to higher concentrations of elements like aluminium, calcium and magnesium (Appelo and Postma, 1993; Kros and Mol, 2001; Broers, 2002). Second, sorption to heterogeneous sediment was considered using a meta-model approach. We considered three types of sorbents and assumed that the sorption to these three individual sorbents was additive (Fest et al., 2005). The three sorbents were clay minerals, iron hydroxides and organic matter. The sorption capacities for the three different sorbents are presented in Table 2 together with other relevant data. The Gaines–Thomas selectivity coefficients for cationexchange to clay minerals refer to values for Camp Berteau montmorillonite (Bruggenwert and Kamphorst, 1982). The intrinsic binding constants for amorphous Fe-hydroxide were obtained from Dzombak and Morel (1990). The intrinsic binding constants of H+, Ca2+ and Cd2+ for sorption to humic acids as organic matter were obtained from Kinniburgh et al. (1996); non-published data were used for Zn2+, Al3+, AlOH2+ and Al(OH)22+ (Table 3). Cation-exchange to clay minerals and surface

Table 3 Characteristics of sorbents considered Cation

Carboxyl group

Phenol group

log KMe

Non-ideality

log KMe

Non-ideality

Zn2+ Al3+ AlOH2+ Al(OH)+2

1.55 0.51 0.95 3.11

0.48 0.36 0.32 0.32

4.40

0.64

56

Cl (mmol/l)

NO3 (mmol/l)

SO4 (mmol/l)

Na (mmol/l)

K Mg (mmol/l) (mmol/l)

Ca (mmol/l)

Al (mmol/l)

Fe (mmol/l)

DOC (mmol C/l)

0.00

0.25

0.24

0.29

0.26

0.02

0.04

0.05

0.21

0.00

0.43

16 4.51 180

0.00

0.40

0.17

0.60

0.37

0.07

0.17

0.23

0.04

0.00

0.21

20 5.71 208 18 4.25 288

0.10 0.00

0.38 0.51

0.01 0.27

0.47 0.69

0.38 0.42

0.04 0.05

0.09 0.06

0.25 0.17

0.00 0.31

0.07 0.01

0.13 0.73

5 5.26 300

0.06

0.60

0.01

0.77

0.92

0.07

0.21

0.63

0.00

0.21

0.38

1 5.81 97 15 4.98 446

0.10 0.03

0.27 0.47

0.01 1.87

0.19 0.57

0.27 0.35

0.01 0.49

0.03 0.42

0.17 1.00

0.00 0.04

0.02 0.01

0.07 1.16

17 5.31 596

0.05

0.96

0.27

0.90

0.68

0.32

0.40

1.14

0.01

0.00

0.32

8 5.55 646 58 5.67 554

0.10 0.24

1.08 0.64

0.02 1.33

1.67 0.75

0.89 0.62

0.09 0.70

0.82 0.59

1.48 1.11

0.00 0.02

0.15 0.01

0.16 2.11

59 5.76 554

0.24

0.95

0.06

0.89

0.88

0.10

0.42

1.12

0.00

0.09

0.46

21 5.98 421 18 6.05 477

0.41 0.79

1.20 0.60

0.01 0.52

0.80 0.63

0.68 0.65

0.06 0.50

0.24 0.46

0.90 1.06

0.00 0.02

0.18 0.02

0.27 2.38

6 6.75 646

3.31

0.91

0.01

0.56

0.51

0.04

0.24

1.85

0.00

0.12

0.38

3 7.00 469

4.00

0.41

0.17

0.27

0.47

0.07

0.26

1.88

0.00

0.09

0.32

Type of homogeneous Depth interval area

N

Nature recharge

40 4.31 167

Nature intermediate

Agriculture recharge

Agriculture intermediate

Discharge

Uppermost (0–1 m) Moderately deep (7–10 m) Deep (20–25 m) Uppermost (0–1 m) Moderately deep (7–10 m) Deep (20–25 m) Uppermost (0–1 m) Moderately deep (7–10 m) Deep (20–25 m) Uppermost (0–1 m) Moderately deep (7–10 m) Deep (20–25 m) Uppermost (0–1 m) Moderately deep (7–10 m) Deep (20–25 m)

pH

Ec Alk (μS/cm) (mmol/l)

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Table 4 Process parameters used in NICA-DONNAN model for sorption of Zn and Al to humic acid (where KMe refers to ion-specific, average intrinsic binding constant; cf. Milne et al., 2003)

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

complexation to amorphous Fe-hydroxide were modelled using PHREEQC (Parkhurst and Appelo, 1999). Surface complexation to humic acid was modelled using ECOSAT (Keijzer and Van Riemsdijk, 1998). General single-solute isotherm type curves were established for these three sorbents from a representative range of fresh groundwater systems. For example, sorption to iron oxides and organic matter was described by a Freundlich isotherm:  n bMe ¼ K Me2þ ð3Þ with log K ¼ a þ b pH þ c pH2 þ d pH3 n ¼ a Vþ bV pH þ c VpH2 þ d V pH3

57

from an empirical fit of data points that relate the total sorbed amount, SMe, to the aqueous concentration of the trace metal. Four trace metal concentrations were considered: the actual concentration, the drinking water limit, and one tenth and ten times the value of the latter. Now, total concentrations instead of the associated aqueous activity of Zn or Cd are used in order to obtain single-solute sorption isotherms that can be used in MT3DMS (Zheng and Wang, 1999). A typical range of Freundlich K values for Zn in the uppermost groundwater is 2.36E-04 to 6.13E-03 (based on concentration in water in μg/l and adsorbed on soil in mg/kg), n values varies between 0.82 and 0.99. 4. Results and discussion 4.1. Regional groundwater contamination

where βMe is the fraction of the sorbent occupied by the metal, [Me2+] is the free activity of Me2+, K is the sortion affinity and n is an exponential term. For each of the three sorbents, the fraction occupied by a trace metal was calculated for a series of fresh water compositions: (1) from pH 4 to 7, (2) from unpolluted (soft, slightly mineralized water) to hard, mineralized but still fresh groundwater, (3) one tenth of the maximum tolerable concentration for drinking water for the trace metal of interest, to ten times that value. This yielded a series of 63 groundwater compositions that is representative for groundwater composition in Pleistocene, sandy areas in the Netherlands (Van den Brink et al., 2007). For each pH value in the range of 4 to 7 with an interval of a half pH value, we fitted a Freundlich isotherm (log k and n) through the set of trace metal free activities and calculated fractions occupied. This resulted in 7 pH-dependent isotherms. For each sorbent, the values for the third-order polynomial pH-function were subsequently obtained from a fit through the series of 7 log k and n values using GENSTAT 5 (1987). The third step was to calculate the total amount sorbed as the sum of the products of the fraction sorbed and the amounts of individual sorbents present for a given activity of a trace metal and a given sediment composition: SMe ¼ Sclay bMe;clay þ Sox bMe;ox þ SSOM bMe;SOM

ð4Þ

where SMe refers to the amount of sorbent present in a geological unit and βMe is the fraction of the sorbent occupied by Me. As amount of sorbent present we used the medians per geological unit of measured sorbent contents (Table 1). The Freundlich isotherm used in the groundwater transport model was calculated per grid cell

Our field assessment shows that in the nature areas there is regional contamination of groundwater at the water table. In the 2003 monitoring program, 31 of the 60 samples of uppermost groundwater exceeded the statutory Dutch intervention limit for cadmium in groundwater (6 μg/l) and 33 exceeded the intervention limit for zinc (800 μg/l). A clear pH-dependency of groundwater contamination with Cd and Zn is observed, which is here not presented. The monitoring programs in the Kempen revealed that the median concentration of Cd is statistically significantly higher than those calculated for surrounding sandy soils areas in the province Noord– Brabant: 6.64 μg Cd/l versus 1.3 μg Cd/l. Note that these values are also anthropogenically influenced. Table 4 presents the median values of the macrochemical components and pH for the homogeneous areas. In all areas, the pH rises with increasing depth. Agricultural areas show higher concentrations of major cations and anions than nature areas. The higher pH in agricultural areas is attributable to the use of manure and lime fertilizers, which boosts the concentrations of the major cations and anions compared with those observed in nature areas. 4.2. Transport through the unsaturated zone and leaching to groundwater Fig. 4 shows the modelled cadmium concentration in the uppermost groundwater (0–3 m depth) in 2005. The dominant pattern is that after starting with very high values close to the smelters, cadmium concentration decreases with increasing distance. The concentrations in the uppermost groundwater are more variable and more scattered than the metal contents in the topsoil

58

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Fig. 4. Simulation results of cadmium concentration (μg/l) in uppermost groundwater (0–3 m depth) in 2005.

(Fig. 1) and atmospheric deposition pattern that is used as upper boundary for the model (Fig. 3). The factors mainly responsible for the short scale spatial variation in concentrations in cadmium content are differences in retardation for various soil types, and also depth of water table. The spatial variability in loading at the water table was found to be high. Fig. 5 shows the

median, 25 and 75 percentile breakthrough curves of cadmium at groundwater level in nature infiltration areas and agricultural infiltration areas, i.e., areas homogeneous in terms of their land use groundwater regime. The median breakthrough curve for the nature infiltration areas peaks about 30 years earlier than for the agricultural infiltration areas, but there is large variation:

Fig. 5. Median, 25 and 75 percentile breakthrough curves of cadmium (μg/l) at the phreatic surface in nature infiltration and agriculture infiltration areas derived from all modelled breakthrough curves in those areas.

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

the 25 percentile breakthrough curve of agricultural infiltration areas occurs before the 75 percentile curve of the nature infiltration areas. The variation in breakthrough curves within homogenous areas is mainly

59

attributed to variation of soil pH and groundwater levels within homogeneous areas. Other factors controlling the variation are differences in soil organic matter and content and differences in the metal load.

Fig. 6. Measured and modelled cadmium concentrations (μg/l) in the upper meter of groundwater in nature infiltration, nature intermediate and exfiltration areas as boxplots with median values, P25–P75 and non-outlier range (N monitoring: 1990 nature infiltration = 13, 2003 nature infiltration = 33, 1990 nature intermediate = 7, 2003 nature intermediate = 46, 1990 discharge = 6, 2003 discharge = 15).

60

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

4.3. Verification of historic modelling The model results were evaluated by comparing the measured and modelled Cd and Zn concentrations for groundwater at the water table. The two datasets of shallow groundwater samples from 1990 to 2003/2004 were used. Input variables for the transport model were derived from regional maps and surveys and hydrologic models. Therefore, since the input parameters were incorporated per area that is homogeneous in terms of land use, soil type and hydrogeology, it is most useful to compare the modelled and measured concentrations per homogeneous area statistically. Fig. 6 shows the median and range of modelled and measured concentrations for 1990 and 2003/2004 as boxplots for the nature infiltration areas, nature intermediate areas and discharge areas. Here, the statistics for modelled concentrations were obtained from all layer-1 grid cells within a homogeneous area. In both the measured and modelled concentrations there is a clear increase from 1990 to 2003. It seems that the medians, the range of concentrations and its temporal changes have been well simulated. Here, it must be noted that the reactive transport model was developed without any calibration. The comparison also reveals that the highest median values and the largest range in Cd and Zn concentrations occur in the uppermost groundwater in nature infiltration areas. In 2003 the modelled and measured median Zn concentrations were 1799 and 1713 μg/l in these areas, compared with 1249 and 900 μg/l in “intermediate” nature areas and 155 and 61 in discharge areas. The high concentration in the nature infiltration areas is attributable to the unsaturated zone being prone to leaching because of its low soil pH. Furthermore, in these areas, the shallow groundwater is not influenced by presence of old pristine exfiltrating groundwater. Large differences of up to two orders of magnitude for individual measurement locations should also be noted, which are attributed to local-scale variation in soil pH, etc. When we compare our modelling fit with other studies as cited below, we note that a match within 50% of observed values for most observation points has not been reached in any 3D field modelling study, even those where process parameters were calibrated. The difference between measured and modelled values is often several factors when individual points are compared, a regional approach cannot cope with individual locations. The discrepancy is attributed to e.g. errors in mass estimate (Schirmer et al., 2000) or problems with specifying the contaminant boundary condition in the source area (Brun et al., 2002). Ideally, physically-based distributed models are applied and all parameters are independently obtained

in field application of transport models. This can be reached for systems that are rather homogenous in their spatial characteristics, have a temporally constant input and where a limited amount of processes is operational. However, when models are used at a regional or larger scale this can seldom be reached, especially in cases of a 3D modelling approach as in this study. Two kinds of approaches are now recognized. First, values for process parameters are obtained by means of calibration. Schäfer and Therrien (1995), Keating and Bahr (1998), Brun et al. (2002) and Prommer and Stuyfzand (2005) used this approach in 3D, local field studies where groundwater transport happened within a single geological stratum and a limited amount of process parameters needed to be calibrated. Mostly, one or a few biodegradation constants were calibrated. A problem with calibration is that non-unique calibration solutions may exist even for these relatively simple field systems (Schäfer and Therrien, 1995; Saaltink et al., 2003). Another problem is that objective calibration criteria have often not been used and the reproduction of patterns is generally felt satisfactory. A minimum requisite for calibration would be use of optimization software and even this software needs subjective weighing criteria for the calibration parameters involved (Van Breukelen et al., 2004). The second kind of approach is based upon use of best available information (Brusseau, 1991). Historical modelling with a predictive approach is performed using independently estimated or measured parameters to test the performance of the model (Kent et al., 2000; Schirmer et al., 2000; Schreiber et al., 2004; Curtis et al., 2006; Thorsen et al., 2001). Sensitivity analysis is performed: 1. to investigate model sensitivity to parameters and 2. to explore the influence of measurements error or system spatial and temporal variability. After this, the model is suitable for future scenario predictions. This approach is more useful for underconstrained problems, where non-unique calibration solutions will certainly exist and unrealistic values for parameters may become obtained. We performed a sensitivity analysis on historic surface loads, groundwater composition and sediment geochemistry. These are the most dominant input parameters controlling the subsurface metal concentration. We modelled a scenario with half of the calculated surface input load and a double surface input load. The statistics of modelled zinc concentrations in nature infiltration areas are plotted in Fig. 7. The difference between modelled median, 25 percentile and 75 percentile of the basis scenario (respectively, 1799, 859 and 2867 μg/l) and those measured (1713, 878 and

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

61

Fig. 7. Boxplots with median values, P25–P75 and non-outliner range of modelled zinc concentrations in uppermost groundwater in nature infiltration areas as a result of a sensitivity analysis on Zn surface load.

3784 μg/l) is smaller than that of the double surface load scenario (4003, 2493 and 6049 μg/l) and half surface load scenario (765, 293 and 1328 μg/l). The second and third sensitivity scenarios were done with the 25 percentile and 75 percentiles concentrations of the individual macro-chemical groundwater components and the sorbent contents of the geological units, respectively. As explained before these input values were derived by taking statistics per homogeneous area of groundwater quality and for sediment geochemistry per geologic unit. The best available information scenario was done with median values from Tables 1 and 4. Because the

Fig. 8. Modelled zinc retardation factor in uppermost groundwater in agricultural area as result of a sensitivity analysis on groundwater chemistry (hydrochemistry) and sediment composition (geochemistry). The retardation was modelled using median, 25 percentile and 75 percentile values of individual macro-chemical groundwater element (including pH) and sediment sorbents (clay minerals, organic matter and iron oxides).

inverse relation between ionic strength of a solution and the relative sorption of an individual trace metal to soil particles, the 25 percentile scenario was done with a solution that has a lower pH than the basis scenario and higher concentrations of macro-chemical elements. Fig. 8 gives the result for these sensitivity scenarios as calculated retardation of zinc in shallow groundwater in agricultural areas. Not surprisingly, the mobility of metals is higher in the 25 percentile scenario and lower in the 75 percentile scenario. However, there is an obvious difference in model sensitivity towards groundwater composition and the sediment geochemistry. As a result of variation in groundwater chemistry the retardation varies between 23 and 34. This is a minor change in mobility comparing to sediment geochemistry scenario. The retardation in this scenario varies between 1 and 131. There are two reasons for this difference in sensitivity. First, the large uncertainty in sediment composition within a geological formation compared to the uncertainty in groundwater concentrations within a homogeneous area: the 25 percentile of the clay content in the Boxtel Formation is lower than the detection limit of 0.1% while the 75 percentile is 12.6%, the 25 percentile of the pH is the uppermost groundwater in agricultural areas is 4.84 and the 75 percentile is 5.45. Second, Cd and Zn are not strong complex-forming metals in fresh groundwater systems. This results in a behavior in which the free metal activity is not strongly dependent on the aqueous composition. Curtis et al. (2006) came to the opposite conclusion for uranium(VI) sorption due to formation of U(VI)-carbonate complexes in groundwater with a higher alkalinity and pH.

62

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Fig. 9. Depth profiles of average values of modelled cadmium and zinc concentrations in homogeneous areas and measured concentration for all areas.

After leaching from the unsaturated zone, the fate of metals depends on the flow of the groundwater and sorption capacity within the groundwater-saturated zone. In infiltration areas metal-contaminated groundwater was found at a depth of 18 m but in intermediate or exfiltration areas clean pristine groundwater was also found at depths near surface. Fig. 9 presents the depth profiles of average Cd and Zn concentrations from all groundwater samples in the area. Unfortunately, there were insufficient samples to allow us to statistically analyze the groundwater samples grouped according to the different homogeneous areas. Concentrations drop strongly from 3 to 13 m below surface. The average groundwater quality is influenced by contaminants leaching from the surface down to a depth of 13.5 m. Increased concentrations were occasionally found below this depth in the monitoring survey as well as in the model. The model predicted no increased metal concentration deeper that 18 m for the present situation. This is in agreement with the results of the monitoring survey. There is a clear difference in the depths of the Cd and Zn front between the homogeneous areas. As expected, the metal front is much deeper in the nature areas than in the agricultural areas, due to difference in pH and associated retardation factor in the shallow subsurface. The average values in nature infiltration areas between 0–3 and 3–6 m depth are almost the same (2210 and 1965 μg/l for Zn, and 29 and 24 μg/l for Cd), indicating that leaching to shallow groundwater is at moment around its maximum in these areas.

4.4. Predictive transport modelling: hydrogeological and geochemical controls As described before, sorption of metals on clay minerals, iron oxides and organic matter is the dominant process controlling reactive transport at shallow depths. The pH of the infiltrating groundwater increases to near neutral with increasing depth (see Table 1), resulting in an increase of sorption. Therefore, cadmium and zinc become strongly retarded in the saturated zone, despite the low reactivity of the sandy sediments. Acidification of groundwater in the past decades has also been limited to no more than several tenths of pH units at shallow depths, and thus no strong mobilization has occurred at regional-scale. The relative importance of sorptive control versus groundwater travel time is thus of prime interest to relate future Cd and Zn surface water load to leaching and transport in groundwater zone. Table 5 provides the mass balances of cadmium and zinc for the three catchments. The mass balances were calculated with the following equations: ðStoredunsat þ Storedsat Þt ¼ Int  Outt

ð5Þ

Rcht ¼ Int  Storedunsat þ Storedsat þ Outt

ð6Þ

Storedunsat ¼ Storedsoil moisture þ Storedsorbed

ð7Þ

Storedsat ¼ Storedgroundwater þ Storedsorbed

ð8Þ

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

63

Table 5 Calculated mass budgets for cadmium and zinc of the Beekloop–Keersop, Buulder Aa and Tungelroijsche Beek catchment (soil surface load in 1000 kg and other budgets in percentage of the total soil surface load) Cadmium

Zinc

Soil surface load

Stored in soil

1000 kg

%

Cumulative recharge flux

Stored in saturated zone Solute

Cumulative discharge flux

Sorbed

Beekloop–Keersop 1965 33.9 98.8 1990 51.7 95.7 2005 52.1 86.2 2050 52.5 52.6

0.9 4.2 13.8 47.3

0.1 0.8 2.3 4.9

0.7 2.7 9.3 32.3

0.1 0.8 2.3 11.0

Buulder Aa 1965 44.7 1990 68.2 2005 68.6 2050 69.1

98.7 90.6 80.2 53.1

1.3 9.4 19.9 46.8

0.3 2.6 4.1 4.6

1.0 5.9 13.7 33.7

0.1 0.9 2.1 8.6

Tungelroijsche Beek 1965 33.9 98.5 1990 51.9 92.7 2005 52.3 86.2 2050 52.8 65.2

1.4 7.3 13.9 34.9

0.4 1.9 2.3 3.3

0.8 3.6 7.7 17.0

0.3 1.8 3.8 14.6

where Int is the soil surface load, Outt is the cumulative discharge flux, Rcht is the cumulative recharge flux, Storedunsat is the metal amount in the unsaturated zone and Storedsat is the amount in the saturated zone. All three catchments show an increase of the ratio sorbed to dissolved in the saturated zone over time, which is a result of the downward movement of the contamination front with associated near neutral pH conditions in the deeper subsurface. The model predicted that in 2050, for all the three catchments, 87% of the cadmium in the saturated zone will be sorbed on the soil particles (compared with 70% in 1990). However, there is a substantial increase in the total amount of cadmium and zinc present in the saturated zone in the period 1990–2050. The model predictions for the three catchments need further attention to illustrate the differences on hydrogeological and geochemical controls. In Beekloop– Keersop intensive leaching of metals to the groundwater starts a few years later than in Buulder Aa and Tungelroijsche Beek. In 1990 4% of the total cadmium load was leached to the groundwater in Beekloop– Keersop; in Buulder Aa and Tungelroijsche Beek this was 9 and 7%. From 1990 to 2005 the leaching in Beekloop–Keersop increases by more than a factor 3.3, compared with comparable factors of 2.1 and 1.9 for

Soil surface load

Stored in soil

1000 kg

%

Cumulative recharge flux

Stored in saturated zone Solute

Sorbed

Cumulative discharge flux

Beekloop–Keersop 1932 93.6 3079 84.5 3219 71.6 3342 29.6

6.4 15.5 28.4 70.4

0.9 2.0 3.7 7.2

4.6 10.0 18.1 44.0

1.0 3.6 6.7 19.3

Buulder Aa 2949 94.6 4614 83.2 4750 69.9 4896 30.0

5.4 16.8 30.1 70.0

1.4 4.0 5.8 6.8

3.4 10.2 19.6 48.7

0.6 2.6 4.7 14.6

Tungelroijsche Beek 2010 94.2 5.9 3218 83.7 16.4 3378 73.8 26.2 3576 39.0 61.0

1.1 2.8 3.1 2.9

3.4 8.7 14.5 29.6

1.3 4.9 8.6 28.6

Buulder Aa and Tungelroijsche Beek. In the end, Tungelroijsche Beek shows the slowest leaching (65% of cadmium still present in the unsaturated zone in 2050, this is both 53% for Beekloop–Keersop and Buulder Aa) and least metals present in the saturated zone (only 58% of the leached cadmium is present in the saturated zone in 2050, this is 79 and 82% for Beekloop–Keersop and Buulder Aa). The soils in the Beekloop–Keersop catchment are most vulnerable to leaching: the soil parent material in the Kempen High area consists of medium and coarse sand and some gravel deposits characterized by lower pH and smaller contents of organic matter and clay fractions than in the Roer Valley Graben. Less intensive leaching in the period before 1990 is explained by the fact that the zinc smelters are not located within the catchment. This catchment does not contain the areas with the highest level of contamination around the smelters. Because non-linear sorption behavior, the retardation factor is inversely related to increasing concentrations, leaching is less intensive. Despite the higher reactivity of the soil parent material, the leaching in the Buulder Aa catchment is as intensive as in the Beekloop–Keersop catchment. The Buulder Aa catchment contains the largest proportion of nature areas. As mentioned before, retardation within the unsaturated zone is greater in agricultural soils than

64

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Fig. 10. Simulated cadmium load by seepage of groundwater to surface waters in the Beekloop–Keersop and Tungelroijsche Beek catchments, flux in kg/year and concentration in μg/l.

in nature soils because the pH is higher due to liming. The Tungelroijsche Beek catchment has the highest retardation of cadmium in the unsaturated zone. Due to the parent material and agriculture as dominant land use, this catchment is less vulnerable for metal leaching to groundwater. Fig. 10 provides the simulated cadmium load of surface water by groundwater exfiltration and drainage in the Beekloop–Keersop and Tungelroijsche Beek catchments as flux in kg/year and average concentration in μg/l. Ditches and drains – particularly the latter – are widespread in the catchments studied and part of the surplus precipitation drains away rapidly via them. The surface water load was calculated by multiplying the concentration in a grid cell by the river and/or drainage flux for that specific grid cell and summation of all grid cell in a catchment. Qj;t ¼

n X

Ci;j;t d Frivi;j;t

ð9Þ

i¼1

Where Qt,j is the surface water load in catchment j and at time t (kg yr− 1), Ci,j,t is the metal concentration in grid cell i (kg m3) and Friv, i,j,t is the river or drainage flux in grid cell i (m3 yr− 1). The concentrations were derived by dividing the mass flux (kg yr− 1) by the water flux (m3 yr− 1) of

exfiltrating groundwater. Here, the concentrations were used without considering geochemical transformation processes in the stream sediment. The modelled concentrations of exfiltrating groundwater are in the same range of several μg/l in both catchments, and will increase in the coming decades. The simulations show that for the Beekloop–Keersop catchment a peak will be reached at 2030, whereas the flux and concentration will gradually increase until 2050 for Tungelroijsche Beek. To establish the hydrogeological control, Fig. 11 presents the relative concentration of a conservative solute in the surface water system as a result of groundwater exfiltration. In this scenario the recharge groundwater contains a contamination block front with a concentration of 1 during 10 years. In the remaining 90 years the recharge concentration was set at 0. The maximum concentration in exfiltrating groundwater is reached in both catchments after 10 years, which is at the end of the block front when the concentration of the recharge water switched from 1 to 0. This means that there is an instantaneous response of the surface water load to changes in groundwater recharge concentration. Note that transport in the unsaturated zone was not considered in this conservative simulation, but for a conservative solute this takes about 1 to 2 years. The two main differences between the Beekloop– Keersop and Tungelroijsche Beek catchments are the

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

65

Fig. 11. Simulated surface water concentration as result of transport of a non-reactive contamination block front input (recharge concentration of 1 during 10 years followed by 0 during 90 years).

maximum relative concentration at 10 years after the start of the simulation and the concentrations from 30 years on. The maximum concentration for the Tungelroijsche Beek catchment is 0.56 and that for Beekloop–Keersop catchment is 0.45, which is a difference of 23%. After 30 years the concentration of the Tungelroijsche Beek has decreased to almost zero. For the Beekloop–Keersop, the relative concentration increases slightly to a second peak of 0.068 at 50 years. Short streamlines feeding the drainage system are more important for the Tungelroijsche Beek than for the Beekloop–Keersop. As well as being fed by the drainage system, the Beekloop–Keersop is also fed by baseflow, which impacts on the surface water system for more than 100 years. Regarding exfiltration of solutes from groundwater to surface water, the results of the reactive transport simulation are opposite to the results of conservative transport (compare Figs. 10 and 11). For a non-reactive tracer, the Tungelroijsche Beek is the quickest responding system, whereas the Beekloop–Keersop is the quickest for reactive transport modelling. This reversal is caused by differences in soil type, dominant land use and geology between the areas: due to dominantly nature land use combined with the low reactivity of the soil parent material, the Beekloop–Keersop catchment is more vulnerable to leaching of metals to the surface water than the Tungelroijsche Beek catchment.

5. Conclusions Using coupled unsaturated and saturated zone flow and reactive transport models, a remarkable difference was found in the behavior of both cadmium and zinc within three catchments studied in the Cd- and Zn-contaminated Kempen area in the south of the Netherlands. This resulted in a difference in surface water load due to exfiltration of contaminated groundwater between the catchments. It is attributable to differences in the geohydrologic and geologic structure, soil type and dominant land use. The intensively drained Tungelroijsche Beek catchment, which in terms of its hydrogeology should be the most vulnerable for solute discharge from groundwater to the surface water, turns out to be the least vulnerable for metals. This is due to its soil parent material that is strongly correlated to geological build-up of the Roer Valley Graben and dominant agricultural land use that retards metal leaching to groundwater for a long time. The Beekloop–Keersop catchment which in terms of its hydrogeology should be the least vulnerable turns out to be the most vulnerable for the metals. This is mainly due to the soil parent material in this catchment on the Kempen High. The subsurface is as a result of its low pH and organic matter content extremely vulnerable for leaching. Overall, it can be concluded that geochemical controls are dominating the subsurface metal transport.

66

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Comparison of statistical results of the historical modelling with those of actual data at the level of homogenous areas is a proper way for evaluation of a regional subsurface transport model: the modelling approach predicts for homogeneous areas cadmium and zinc concentrations in shallow groundwater that are in the same order of magnitude and have the same range in concentrations as measured values. Depth profiles of concentrations are also comparable between model results and field measurement. As concentration depth profiles typically include all the important features of subsurface solute transport that are soil input load, groundwater flow and geochemical processes, these are very useful for verification of regional subsurface transport models. A sensitivity analysis on metal surface input loads, groundwater composition and sediment geochemistry reveals that the best available information scenario predicts the range in observed concentrations very well. Due to the wide range in measured organic matter, clay and iron oxide content of the sediment the model is sensitive to the paramerisation of the sediment geochemistry. Acknowledgements The authors are grateful to ‘Actief Bodembeheer de Kempen’, for funding a major part of the research. Ebel Smidt is gratefully acknowledged for stimulating discussions. The research was co-financed by TNO research program for the Ministry of Environment. Peter Venema is thanked for providing non-published data for binding constants to humic acids. Joy Burrough advised on the English. The authors thank the editor-in-chief, Dr. E.O. Frind, and three anonymous reviewers for their useful critical comments on the draft versions of this paper. References Appelo, C.J.A., Postma, D., 1993. Geochemistry, groundwater and pollution. A.A. Balkema. Boekhold, A.E., 1992. Field scale behavior of cadmium in soil. Ph.D thesis, Agricultural University Wageningen, the Netherlands. Broers, H.P., 2002. Strategies for regional groundwater quality monitoring, Ph.D thesis, Utrecht University, the Netherlands. Broers, H.P., 2004. The spatial distribution of groundwater age for different geohydrological situations in the Netherlands: implications for groundwater quality monitoring at the regional-scale. J. Hydrol. 299, 84–106. Broers, H.P., Van der Grift, B., 2004. Regional monitoring of temporal changes in groundwater quality. J. Hydrol. 296, 192–220. Brown, J.G., Bassett, R.L., Glynn, P.D., 1998. Analysis and simulation of reactive transport of metal contaminants in groundwater in Pinal Creek Basin, Arizona. J. Hydrol. 209, 225–250. Bruggenwert, M.G.M., Kamphorst, A., 1982. Survey of experimental information on cation exchange in soil systems. In: Bolt, G.H.

(Ed.), Soil Chemistry B. Physico-chemical models. Elsevier, Amsterdam, pp. 141–204. Brun, A., Engesgaard, P., Christensen, Th.H., Rosbjerg, D., 2002. Modelling of transport and biogeochemical processes in pollution plumes: Vejen landfill, Denmark. J. Hydrol. 256, 228–247. Brusseau, M.L., 1991. Application of a multiprocess nonequilibrium sorption model to solute transport in a stratified porous medium. Water Resour. Res. 27, 589–595. Buma, J.T., Kremers, A.H.M., Van der Meij, J.L., Te Stroet, C.B.M., Vernes, R.W., 2002. Watergoals — desired ground- and surface water regime. Development model instrument and exploring calculations of bottlenecks, measures and overall solutions (in Dutch). TNO report NITG. Burt, R., Wilson, M.A., Keck, T.J., Dougherty, B.D., Strom, D.E., Lindahl, J.A., 2003. Trace element speciation in selected smeltercontaminated soils in Anaconda and Deer Lodge Valley, Montana, USA. Adv. Environ. Res. 8, 51–67. CSO, 2001. Actief Bodembeheer de Kempen - Gevalsafbakening (in Dutch), CSO rapport nr. 00.298, Bunnik, the Netherlands. Curtis, G.P., Davis, J.A., Naftz, D.L., 2006. Simulation of reactive transport of uranium (VI) in groundwater with variable chemical conditions. Water Resour. Res. 42 (art. no. W04404). De Bakker, H., Schelling, J., 1989. A system of soil classification for the Netherlands. The higher levels. Winand Staring Centre, Wageningen, the Netherlands. De Mulder, E.F.J., Geluk, M.C., Ritsema, I., Westerhoff, W.E., Wong, T.E., 2003. The subsurface of the Netherlands (in Dutch). Netherlands Institute of Applied Geoscience TNO, Utrecht, the Netherlands. De Vries, F., 1999. Characterization of Dutch soil on physical– chemical characteristic (in Dutch). DLO-Staring Centrum, report no. 654. Wageningen, the Netherlands. Dzombak, D.A., Morel, F.M.M., 1990. Surface complexation modeling. Hydrous Ferric Oxide. John Wiley & Sons, New York. Edmunds, W.M., Kinniburgh, D.G., Moss, P.D., 1992. Environ. Pollut. 77, 129–141. Fernandez-Turiel, J.L., Acenolaza, P., Medina, M.E., Llorens, J.F., Sardi, F., 2001. Assessment of a smelter impact area using surface soils and plants. Environ. Geochem. Health 23, 65–78. Fest, E.P.M., Temminghoff, E.J.M., Griffioen, J., Van Riemsdijk, W.H., 2005. Proton buffering and metal leaching in sandy soils. Environ. Sci. Technol. 39, 7901–7908. Frei, M., Bielert, U., Heinrichs, H., 2000. Effects of pH, alkalinity and bedrock chemistry on metal concentrations of springs in an acidified catchment (Ecker Dam, Harz Mountains, FRG). Chem. Geol. 170, 221–242. Genstat 5, 1987. Genstat 5 Reference Manual. Clarendon Press, Oxford, U.K., 749 pp. Griffioen, J., Lourens, A.L., Venema, P., te Stroet, C.B.M., Minnema, B., Laeven, M.P., Stuyfzand, P.J., van Beek, C.G.E.M., Beekman, W., 1998. An integrated model for predicting and assessing the development of groundwater quality. In: Peters, J., et al. (Ed.), Artificial Recharge of Groundwater; Proc. Intern. Symp., TISAR 98. Amsterdam, the Netherlands, pp. 419–421. Harmsen, K., 1977. Behaviour of heavy metals in soils. Ph.D thesis, Agricultural University Wageningen, the Netherlands. Hotsma, P.H., 1997. Heavy metal contents in manure (in Dutch). Information and Knowledge Centre Agriculture. Ede, the Netherlands. Keating, E.H., Bahr, J.M., 1998. Using reactive solutes to constrain groundwater flow models at a site in northern Wisconsin. Water Resour. Res. 34, 3561–3571.

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68 Keijzer, M.G., Van Riemsdijk, W.H., 1998. Ecosat user manual. Agricultural University, Wageningen. Keller, A., Von Steiger, B., Van der Zee, S.E.A.T.M., Schulin, R., 2001. A stochastic empirical model for regional heavy-metal balances in agroecosystems. J. Environ. Qual. 30, 1976–1989. Kent, D.B., Abrams, R.H., Davis, J.A., Coston, J.A., LeBlanc, D.R., 2000. Modeling the influence of variable pH on the transport of zinc in a contaminated aquifer using semiempirical surface complexation models. Water Resour. Res. 36, 3411–3425. Kinniburgh, D.G., Milne, C.J., Benedetti, M., Pinheiro, J.P., Filius, J., Koopal, L.K., Van Riemsdijk, W.H., 1996. Metal ion binding by humic acid: application of the NICA-Donnan Model. Environ. Sci. Technol. 30 (5), 1687–1698. Kros, J., Mol, J.P., 2001. Historic pH and nitrogen availability in forest and nature areas. Alterra, report no 217, Wageningen, the Netherlands. Kjøller, C., Postma, D., Larsen, F., 2004. Groundwater acidification and the mobilization of trace metals in a sandy aquifer. Environ. Sci. Technol. 38, 2829–2835. Levinson, A.A., 1977. Introduction to exploration geochemistry. Applied Publishing Ltd. Makaske, G.B., Vissenberg, H.A., van Grinsven, J.J.M., Tiktak, A., Sauter, F.J., 1995. Metras: a one-dimensional model for assessment of leaching of trace metals from soil. RIVM, report no. 715501005. McBride, M.B., Richards, B.K., Steenhuis, T., Russo, J.J., Sauve, S., 1997. Mobility and solubility of toxic metals and nutrients in soil fifteen years after sludge application. Soil Sci. 162 (7), 487–500. McDonald, M.G., Harbaugh, A.W., 1988. A modular three-dimensional finite-difference groundwater flow model. U.S.G.S. Tech. Water-Res. Inv., Book 6, Chapter Al. Miller, C.W., Foster, G.D., Majedi, B.F., 2003. Baseflow and stormflow metal fluxes from two small agricultural catchments in the Coastal Plain of the Chesapeake Bay Basin, United States. Appl. Geochem. 18, 483–501. Milne, C.J., Kinniburgh, D.G., Van Riemsdijk, W.H., Tipping, E., 2003. Generic NICA-Donnan model parameters for metal-ion binding by humic substances. Environ. Sci. Technol. 37, 958–971. Moolenaar, S.W., Lexmond, Th.M., 1998. Heavy metal balances of agroecosystems in the Netherlands. Neth. J. Agric. Sci. 46, 171–192. Parkhurst, D.L., Appelo, C.A.J., 1999. User's guide to PHREEQC (version 2) — A computer program for speciation, batch-reaction, one-dimensional transport, and inverse geochemical calculations/ U.S. Geological Survey Water-Resources Investigations Report 99-4259. 312 pp. Prommer, H., Stuyfzand, P.J., 2005. Identification of temperaturedependent water quality changes during deep well injection experiment in a pyritic aquifer. Environ. Sci. Technol. 39, 2200–2209. Raven, K.P., Loeppert, R.H., 1997. Trace element composition of fertilizers and soil amendments J. Environ. Qual. 26, 551–557. RIVM (Ed.), 2002. Minerals Accounting System and the environment. Balance and outlook. National Institute for Public Health and the Environment, Bilthoven, the Netherlands. RIVM/CBS, 2002. Milieucompendium, The environment in numbers. Internet (www.milieucompendium.nl). National Institute for Public Health and the Environment, Statistics Netherlands, Bilthoven, the Netherlands. Römkens, P.F.A.M., 1998. Effects of land use changes on organic matter dynamics and trace metal solubility in soils. PhD thesis, Groningen University, the Netherlands, 156 pp. Römkens, P.F.A.M., Groenenberg, J.E., Bonten, L.T.C., de Vries, W., Bril, J., 2004. Derivation of partition relationships to calculate Cd, Cu, Ni, Pb and Zn solubility and activity in soil solutions. Alterrareport no 305, Wageningen, the Netherlands.

67

Rozemeijer, J.C., 2002. Heavy metals in the unsaturated zone (in Dutch). TNO report NITG 02-081-A. Saaltink, M.W., Ayora, C., Stuyfzand, P.J., Timmer, H., 2003. Analysis of a deep well recharge experiment by calibrating a reactive transport model with field data. J. Contam. Hydrol. 65, 1–18. Schäfer, W., Therrien, R., 1995. Simulating transport and removal of xylene during remediation of a sandy aquifer. J. Contam. Hydrol. 19, 205–236. Schirmer, M., Molson, J.W., Frind, E.O., Barker, J.F., 2000. Biodegradation modelling of a dissolved gasoline plume applying independent laboratory and field parameters. J. Contam. Hydrol. 46, 339–374. Schokker, J., 2003. Patterns and processes in a Pleistocene fluvioaeolian environment, Roer Valley Graben, south-eastern Netherlands. Ph.D thesis, Utrecht University. Schreiber, M.E., Carey, G.R., Feinstein, D.T., Bahr, J.M., 2004. Mechanisms of electron acceptor utilization implications for simulating anaerobic biodegradation. J. Contam. Hydrol. 73, 99–127. Seuntjens, P., Mallants, D., Šimůnek, J., Patyn, J., Jacques, D., 2002. Sensitivity analysis of physical and chemical properties affecting field-scale cadmium transport in a heterogeneous soil profile. J. Hydrol. 264, 185–200. Šimůnek, J., Šejna, M., van Genuchten, M.Th., 1998. The Hydrus-1D Software Package for simulating the one-dimensional movement of water, heat and multiple solutes in variable-saturated media, Version 2, U.S. Salinity Laboratory. Sonke, J.E., Hoogewerff, J.A., van der Laan, S.R., Vangonsveld, A., 2002. Chemical and mineralogical reconstruction of Zn-smelter emissions in the Kempen region (Belgium), based on organic pool sediment cores. Sci. Total Environ. 292, 101–119. Steur, G.G.L., Heijink, W., 1991. Soil map of the Netherlands, scale 1:50.000, DLO-Staring Centrum, Wageningen, the Netherlands. Strobel, B.W., Hansen, H.C.B., Borggaard, O.K., Andersen, M.K., Raulund-Rasmussen, K., 2001. Cadmium and copper release kinetics in relation to afforstation of cultivated soil. Geochim. Cosmochim. Acta 65, 1233–1242. Tauw, 1991. Additional investigation of diffuse groundwater contamination in the Kempen. (in Dutch) Tauw, report no. 3106691, Deventer, the Netherlands. Tiktak, A., 1999. Modeling non-point source pollutants in soils. Applications to the leaching and accumulation of pesticides and cadmium. PhD thesis, Agricultural University Wageningen, the Netherlands, 230 pp. Thorsen, M., Refsgaard, J.C., Hansen, S., Pebesma, E., Jensen, J.B., Kleeschulte, S., 2001. Assessment of uncertainty in simulation of nitrate leaching to aquifers at catchment scale. J. Hydrol. 242, 210–227. Valstar, J.R., McLaughlin, D.B., te Stroet, C.B.M., van Geer, F.C., 2004. A representer-based inverse method for groundwater flow and transport. Water Resour. Res. 40, 1–12. Van Breukelen, B.M., Griffioen, J., Röling, W.F.M., Van Verdeveld, H.W., 2004. Reactive transport modelling of biogeochemical processes and carbon isotope geochemistry inside a landfill leachate plume. J. Contam. Hydrol. 70, 249–269. Van den Brink, C., Frapporti, G., Griffioen, J., Zaadnoordijk, W.J., 2007. Statistical analysis of anthropogenic versus geochemicalcontrolled differences in groundwater composition in the Netherlands. J. Hydrol. 336, 470–480. Walter, A.L., Frind, E.O., Blowes, D.W., Ptacek, C.J., Molson, J.W., 1994. Modelling multicomponent reactive transport in groundwater. 2. Metal mobility in aquifers impacted by acidic mine tailings discharge. Water Resour. Res. 30, 3137–3148.

68

B. van der Grift, J. Griffioen / Journal of Contaminant Hydrology 96 (2008) 48–68

Wilkens, B.J., 1995. Evidence for groundwater contamination by heavy metals through soil passage under acidifying conditions. Ph. D thesis, Utrecht University, the Netherlands, 146 pp. Wilkens, B.J., Loch, J.P.G., 1997. Accumulation of cadmium and zinc from diffuse emission on acid sandy soils, as a function of soil composition. Water Air Soil Pollut. 96, 1–16. Wösten, J.H.M., Veerman, G.J., de Groot, W.J.M., Stolte, J., 2001. Water retention and conductivity characteristics of upper en deeper soil layers in the Netherlands (in Dutch). Alterra, report 153, Wageningen the Netherlands.

Xue, H., Nhat, P.H., Gächter, R., Hooda, P.S., 2003. The transport of Cu and Zn from agricultural soils to surface water in a small catchment. Adv. Environ. Res. 8, 69–76. Zheng, C., Wang, P.P., 1999. MT3DMS; A modular three-dimensional multispecies transport model for simulation of advection, dispersion, and chemical reactions of contaminants in groundwater systems; Documentation and user's guide. Department of Geological Sciences, University of Alabama.