Model-Based Integration and Analysis of Biogeochemical and Isotopic

Aug 10, 2013 - Faculty of Geosciences, Utrecht University, 3584 CB Utrecht, The ... School of Earth and Environment, University of Western Australia, ...
0 downloads 0 Views 2MB Size
Download PDF
Article pubs.acs.org/est

Model-Based Integration and Analysis of Biogeochemical and Isotopic Dynamics in a Nitrate-Polluted Pyritic Aquifer Yan-Chun Zhang,† Henning Prommer,*,‡,§,∥ Hans Peter Broers,⊥,#,▽ Caroline P. Slomp,† Janek Greskowiak,○ Bas van der Grift,⊥ and Philippe Van Cappellen†,◆ †

Faculty of Geosciences, Utrecht University, 3584 CB Utrecht, The Netherlands CSIRO Land and Water, Private Bag No. 5, Wembley WA6913 , Australia § School of Earth and Environment, University of Western Australia, Crawley WA6009, Western Australia ∥ National Centre for Groundwater Research and Training, Flinders University, Adelaide SA5001, South Australia, Australia ⊥ Subsurface and Groundwater, Deltares Research Institute, 3584 CB Utrecht, The Netherlands # TNOGeological Survey of The Netherlands, 3584 CB Utrecht, The Netherlands ▽ Critical Zone Hydrology Group, VU University, 1081 HV Amsterdam, The Netherlands ○ Hydrogeologie und Landschaftswasserhaushalt, University Oldenburg, 26129 Oldenburg, Germany ◆ Earth and Environmental Sciences, University of Waterloo, Waterloo, Ontario Canada N2L 3G1 ‡

S Supporting Information *

ABSTRACT: Leaching of nitrate from agricultural land to groundwater and the resulting nitrate pollution are a major environmental problem worldwide. Its impact is often mitigated in aquifers hosting sufficiently reactive reductants that can promote autotrophic denitrification. In the case of pyrite acting as reductant, however, denitrification is associated with the release of sulfate and often also with the mobilization of trace metals (e.g., arsenic). In this study, reactive transport modeling was used to reconstruct, quantify and analyze the dynamics of the dominant biogeochemical processes in a nitrate-polluted pyrite-containing aquifer and its evolution over the last 50 years in response to changing agricultural practices. Model simulations were constrained by measured concentration depth profiles. Measured 3H/3He profiles were used to support the calibration of flow and conservative transport processes, while the comparison of simulated and measured sulfur isotope signatures acted as additional calibration constraint for the reactive processes affecting sulfur cycling. The model illustrates that denitrification largely prevented an elevated discharge of nitrate to surface waters, while sulfate discharges were significantly increased, peaking around 15 years after the maximum nitrogen inputs.

■

carbon.9,13−18 Besides the competition with other reductants, important factors known to affect or control in situ rates of denitrification by pyrite include the pH of the groundwater18 and the reactivity and distribution of the pyrite in the aquifer.19,20 Linked to the release of sulfate, more hazardous contaminants, in particular trace metals (e.g., arsenic21 or nickel21,22) may be mobilized during pyrite oxidation. As a consequence, groundwater quality within capture zones of drinking water supplies and other water bodies may become and remain seriously deteriorated over extended periods. While in some highly industrialized countries changes in agricultural practices have resulted in declining nitrate loads to groundwater systems, a significant lag period typically exists until improve-

INTRODUCTION Leaching of nitrate from agricultural land to groundwater and the resulting nitrate pollution are a major environmental problem worldwide.1−5 In anoxic aquifers that contain sufficiently reactive electron donors, however, nitrate concentration may be significantly attenuated. For example, various field studies have demonstrated that microbially mediated6 (autotrophic) denitrification coupled to pyrite oxidation decreased nitrate concentrations while sulfate was released to the groundwater: 5FeS2 + 14NO3− + 4H+ → 5Fe2 + + 7N2 + 10SO4 2 − + 2H 2O

(1)

Where present, pyrite often is the dominant electron donor, that is, reductant7−10 even in cases where its potential reductive capacity11,12 is clearly exceeded by that of simultaneously present, but less reactive or inert sediment-bound organic © 2013 American Chemical Society

Received: Revised: Accepted: Published: 10415

May 29, 2013 August 9, 2013 August 10, 2013 August 10, 2013 dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422

Environmental Science & Technology

Article

that are (i) mostly located between 10 and 20 m b.g.l. and (ii) enriched in organic matter.33 This Pliocene aquifer is overlain by fluvial sands and gravels of Pleistocene age which were deposited by the Meuse River. Several faults cross the area which are revealed by horizontal jumps in hydraulic heads in the upper Pleistocene aquifer, one of them close to the research location.34 Distributed over the study site several multilevel sampling devices (MLS) (20 cm length, 3 cm diameter at 2 m depth intervals, up to 40 m depth b.g.l.) were installed in 1996. These were used to collect high-resolution geochemical data, including major ion and trace metals in 1996 and 2006. Sulfur isotope data of aqueous sulfate and solid phase sulfur were obtained in 2006 and 2007.33 At the same time 3H and 3He concentrations were measured to determine groundwater ages in samples collected from separate screens (2 in. diameter, 2 m length) available within the MLS setup. Our previous field studies showed a well-defined depth-dependent redox pattern, with (i) an oxic zone underlying the shallow unsaturated zone (ii) a relatively thin denitrifying zone at ∼10 m a.s.l. (above sea level) and a sulfate reducing zone in the lower part of this aquifer. Heavy metals (e.g., Ni, Zn, As) were detected in the drinking water wells in this area decades ago. Concentrations of Ni, As, and SO4 at the production well field increased between 1989 and 1996 and the well field was eventually closed. Earlier studies indicate that the presence of heavy metals is linked to the occurrence of pyrite oxidation by nitrate.16,35 The full details of the sampling and analytical techniques are reported elsewhere.16,33 In this paper, we mostly concentrate on employing the data collected at well 40, which is located on the selected flowpath. Modeling Study. A numerical reactive transport model was set up to integrate and interpret the comprehensive data set characterizing flow, transport, chemical reactions, and isotope patterns at the study site and to identify the time-scale of the response to changes in agricultural practices. Our approach was to (i) initially use the existing piezometric head and groundwater age data to constrain and quantify groundwater flow and conservative transport, before (ii) simulating the temporally and spatially changing biogeochemical transformations and the associated isotope signatures over a 100 year period from 1906 to 2006. Groundwater ages in most parts of the model domain, that is, except for the sections nearest to the base of the aquifer, were significantly younger than 100 years. Therefore the simulated spatial biogeochemical patterns, in particular those located in zones that affect solute discharge to surface waters, are not impacted by a potentially incomplete model spin-up (i.e., the model simulation period required until concentrations reach a quasi-steady-state). The key drivers for the temporal changes in biogeochemical transformations in the groundwater system are the variations in the groundwater recharge composition that are linked to changing agricultural practices and atmospheric deposition, particularly over the last 50 years. In contrast to the varying recharge composition, recharge fluxes, and groundwater flow conditions were assumed constant over the entire simulation period. By making the latter simplification we assume that, for example, the effects of localscale land-use changes such as the expansion of urban areas, variations in vegetation and the effects of climate variability are negligible. Even though we only have high-resolution data for a small portion of the field site, we have chosen to extend the model domain from the groundwater divide (upstream) all the way to the Meuse River (downstream). This approach allows for an assessment of the evolution of groundwater age patterns

ments are also seen in the surface waters that drain the polluted aquifers.21,23 To reconstruct and quantify the biogeochemical processes that have occurred in response to the historic changes in agricultural practices environmental tracers such as tritium (3H)/helium (3He), chlorofluorocarbons (e.g., CCl2F2) and stable isotope measurements have shown to provide useful additional constraints.24−26 For example, microbially mediated sulfate reduction leads to a measurable enrichment of the heavy isotope in the remaining sulfate because the lighter isotope (32S) is consumed preferentially. In contrast, pyrite oxidation may show no or only negligible sulfur isotope fractionation and can thus be used to elucidate sulfur sources through massbalance considerations.25,26 Beyond such mass-balance calculations, observed isotope patterns may also be analyzed and illustrated through reactive transport models.27−29 The numerical modeling approach allows for a fully integrated quantification of complex biogeochemical reactions mechanisms and the associated isotope signatures in aquifers characterized by both hydrogeological and geochemical heterogeneities. Previous modeling studies of denitrification coupled to pyrite oxidation12,20,30 have illustrated the role of pyrite as electron donor under steady state flow and constant31 or highly simplified12 nitrate input conditions. Also, Howden et al.32 recently presented a catchment-scale model to assess and predict the impacts of diffuse nitrate pollution, thereby assuming a fixed mean travel time. However, to our knowledge none of these studies have explored and quantified the longerterm spatial and temporal variations in turnover rates and isotopic signatures that result from changes in agricultural practices and corresponding nitrate loads to the groundwater. In this study we use a detailed data set from a shallow aquifer system in The Netherlands to constrain numerical simulations of groundwater flow, solute transport and biogeochemical dynamics associated with the transient, diffusive input of nitrate and its attenuation by pyrite oxidation. With the numerical model we integrate a multitude of data types and test whether our previously established conceptual models16,33 for the fate of nitrate and sulfate at the catchment scale hold in a fully quantitative context. As such it is the first modeling study that links a more complex sequence of hydrogeochemical and isotopic processes to changes in agricultural practices on a decadal time scale, and assesses the time-scales involved in the recovery of groundwater quality and the quality of groundwater discharge. The development and application of process-based quantification tools will contribute to a better understanding and eventually improved predictions of mass discharge characteristics to surface waters, while complementing the development of more simplistic models.

■

MATERIALS AND METHODS Study Site. The investigated field site is located in Oostrum, Limburg in the south of The Netherlands (Supporting Information (SI) Figure S1). Over the past decades the area was characterized by a high agricultural activity with resulting maximum nitrate concentrations in the recharge water of up to 4 mM in 1985.16 The shallow groundwater levels (3−5 m b.g.l., that is, below ground level) allow for a rapid transfer of pollutants from the ground surface through the vadose zone to the groundwater. The investigated aquifer has an approximate thickness of 50 m whereby impermeable clay layers at its bottom define the base of the aquifer. The main aquifer consists of coarse fluvial sands of Pliocene with unconnected clay lenses 10416

dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422

Environmental Science & Technology

Article

Biogeochemical Reaction Network. The reaction network was defined to consider all key water-sediment reactions affecting the redox and major ion chemistry over the simulation period. The reactions that were assumed to drive the most significant biogeochemical changes were the kinetically controlled reactions between the more oxidized recharge water and the aquifer’s major reductants pyrite, sedimentary organic matter (SOM) and siderite. Pyrite oxidation was assumed to follow a previously developed and applied rate expression7−9,42−44

upstream of and within the study site and eliminates the need to define the transient and depth-varying water compositions at the model’s upstream boundary. Groundwater flow was simulated with the USGS code MODFLOW,36 whereas the subsequent reactive transport simulations were carried out with PHT3D.37 Model Geometry and Setup. The model was constructed as a 2D vertical transect along a streamline passing through the local study site. The model domain has a lateral extension of 10 000 m and extends vertically between −20 and 40 m a.s.l. over a thickness of 60 m. The conceptual hydrogeological model was derived on the basis of earlier field investigations.16,33 The key hydrogeological features (see also Figure 2 and SI Table S2) affecting the groundwater flow within the model domain are (i) a layer of heterogeneously distributed clayey material (clay lenses) separating the sandy aquifer into an upper Pleistocene and a lower Pliocene section (ii) a vertical fault that is hypothesized to act as a flow barrier in the upper aquifer, leading to younger, more polluted groundwater reaching deeper aquifer sections and (iii) two surface water courses that allow the groundwater to discharge. The exact distribution of the clay lenses is not well characterized. Therefore we investigated a range of possible distributions and options to represent the clay layers in the numerical model. The three main variants tested were (i) a clay layer of constant thickness with holes at regular intervals, (ii) a continuous clay layer with an increased hydraulic conductivity relative to the first option, and (iii) a random conductivity field with specified horizontal and vertical correlation length. The latter option provided the most plausible results and was therefore employed as a basis for the subsequent model development. Two water courses modeled as drains to which groundwater can discharge, are located at a distance of 3000 and 8000 m from the upstream groundwater divide (SI Figure S1), in accordance with a previously established and calibrated larger-scale (regional) groundwater flow model for this area.38 The drain at x = 8000 m represents a meander of the Oostrum brook which has cut deeply in the fluvial terraces of the Meuse river, the drain at 3000 m represents a smaller headwater which discharges in the Oostrum brook. Groundwater Age. The reconstructed transient recharge water composition included the varying 3H (tritium) input (SI Figure S2), and its distinct peak in the 1960s.39 In the model simulations the advective-dispersive transport of 3H and its transformation to 3He (helium) was considered, allowing fieldmeasured depth-profiles to be used as constraints for the flow model calibration. A half-life of 12.32 yrs was used to represent the decay of 3H and production of 3He: dC 3 H = −k 3HC 3H dt

0.67 ⎤ ⎡ ⎛ A pyr ⎞⎛ C ⎞ ⎥ rpyr = ⎢(CO2 0.5 + f2 C NO3−0.5)C H+−0.11⎜10−10.19 ⎜ ⎟ ⎟ ⎢ V ⎠⎝ C0 ⎠ ⎥⎦ ⎝ ⎣ pyr

(4)

where r is the reaction rate, c refers to the concentration of either O2, NO3− or H+, (C/C0)pyr reflects the change in the surface area of pyrite due to their dissolution, Apyr/V represents a ratio of pyrite surface area to solution volume, and f 2 is a constant. The kinetically controlled oxidation of SOM using either oxygen, nitrate or sulfate as electron acceptor8,9 was assumed to follow ⎡ CO2 rSOM = ⎢k O2 + kNO3− ⎢⎣ 2.94 × 10−4 + CO2 C NO3− + k SO4 2− 1.55 × 10−4 + C NO3− ⎤ ⎥ + CSO4 2− ⎥⎦

CSO4 2− 1.0 × 10−4

(5)

where kO2, kNO3, and kSO4 are reaction rate constants. The dissolution of iron carbonates (FeCO3) was modeled as rsiderite = ksiderite(1 − SR siderite)

(6)

where SRsiderite is the saturation ratio of siderite and ksiderite is a reaction rate constant. Other potentially present mineral phases such as quartz and Fe-silicates were assumed to have a relatively low reactivity compared to the considered solid phase reactants and were therefore not included in the simulations. As in some previously described, comparable modeling studies, iron was decoupled from the global redox equilibrium and oxidation of dissolved ferrous to ferric iron was modeled as8 −

dFe = k[OH−]2 (PO2 + C NO−3 )C Fe2+ dt

(7)

where k is the rate constant, [OH−] is the hydroxyl ion activity. PO2, CFe2+, and CNO3− represent oxygen partial pressure and the concentrations of Fe2+ and nitrate, respectively. The conceptual model and its numerical implementation were successively adapted and improved during the model development and calibration procedure. All kinetic reactions were implemented into the PHREEQC-2 standard database utilized by PHT3D. The employed reaction rate parameters are listed in SI Table S3. Sulfur Isotope Signatures. The evolution of the sulfur isotope signature was simulated by separately tracking the fate of the light (32S) and heavy (34S) sulfur isotopes.27,45 This allowed the consideration of the different isotope ratios of the sulfate contained within the recharge water, estimated to be ∼2.6‰33, and that of the sulfate released during pyrite

(3)

where C3H is the 3H concentrations and k3H is the first-order rate constant. If degassing of 3He occurred,40 this could have affected the 3He concentrations. However, no correction for degassing has been applied in the present study. The groundwater age was also determined through a direct age simulation following the approach of Goode et al.41 to create a full spatial illustration of the groundwater age patterns that prevail within the model domain. In this approach the groundwater age was simulated by employing a 0th-order reaction (production) in the advection-dispersion-reaction equation. 10417

dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422

Environmental Science & Technology

Article

oxidation (∼−4.2‰33). The oxidation rate of the lighter pool of pyrite is computed from r 32 pyr = rpyr

32

Groundwater Recharge and Recharge Water Compositions. The estimated groundwater recharge rates varied spatially within the model domain, depending on predominant land use but were assumed to remain constant over time (SI Table S2). The input of nitrate to the groundwater and other constituents of the recharge water composition (spatially constant) has varied significantly over the past decades in this region. Nitrate input at farmlands began to decrease after ∼1985 when a national manure law was introduced for groundwater and drinking water quality protection.16 The modeled temporal changes in the input of major compounds NO3, SO4, K, Mg, Ca, and Na over the last 50 years were initially adopted from detailed investigations for a nearby area which was subject to similar agricultural activities and changes thereof (SI Table S1).50 For the period from 1906 to 1956, that is, for the model spin-up period, the recharge water composition was assumed to be similar to that reconstructed for 1956, except for nitrate. Nitrate concentration before 1920 were assumed to be constant at 0.24 × 10−3 mol L−1 before linearly increasing to 0.86 × 10−3 mol L−1 in 1956. In order to better represent the local conditions for the Oostrum site, the nitrate input function between 1956 and 2006 (SI Figure S2) was scaled such that the maximum nitrate concentration in the groundwater recharge was 4 × 10−3 mol L−1 (256 mg L−1) which is within the typical range observed for similar types of land use in the southern part of The Netherlands.

S

C 32 pyr C32 pyr + C 34pyr

(8)

where C32pyr and C32pyr are the concentrations of the light and heavy isotopes, respectively. The rate of the heavy isotope is correspondingly computed from r 34pyr = rpyr

C 34pyr C 32 pyr + C 34pyr

(9)

In addition, the kinetic isotope effect on sulfur isotope ratios was modeled. Similar to Prommer et al.27 the isotope fractionation linked to the oxidation of SOM under sulfatereducing conditions was included. The rate of sulfate removal during SOM oxidation for the light isotope is computed from r 32SO24− = rSO24−

C 32SO24− C 32SO24− + C 34SO24−

(10)

where C32SO42− and C34SO42− are the dissolved sulfate concentrations corresponding to the light and heavy isotopes, respectively and rSO2− is the sulfate removal rate which can be 4 computed from the sulfate-related term in eq 5, scaled by a stoichiometric coefficient. The sulfate removal rate for the heavy isotope is computed from r 34SO24− = rSO24−

C 34SO24− C 32SO24− − + C 34SO24−

■

RESULTS AND DISCUSSIONS Simulated Groundwater Flow and Age Patterns. Using head measurements and measured 3H/3He concentration profiles as joint constraints, the model calibration resulted in plausible groundwater flow patterns with mostly small deviations between simulations and measurements. The model calibration indicated that the simulated age patterns were strongly affected by the two surface water features (drains), which is in agreement with previous studies in lowland areas.51,52 The observed vertical head gradient could only be reproduced when vertical conductivities of the clayey aquifer horizon were assumed sufficiently low, confirming the results of earlier studies.38,52 Figure 1 shows the modeled 3He and 3H depth profiles for well 40 (MLS 40) for the years 1986, 1996,

(ε + 1) (11)

where ε is the enrichment factor, which was adopted from previously reported values27,45−47 and set to 13‰. In the absence of any specific, laboratory-derived information about a possible reaction-rate dependency of the enrichment factor for sulfate reduction, as discussed by Sikora et al.48 and dePaolo49 for calcite precipitation, we omitted to include such a dependency. Modeled Ambient Hydrogeochemistry. Initial mineral and SOM concentrations were adopted from Zhang et al.16 The reductants pyrite, SOM, and siderite were assumed to be abundant throughout the model domain below ∼10 m a.s.l., thus reflecting the measured total S and total C contents of the sediments within this depth range.35 The initial concentrations of dissolved compounds employed for the model simulations were constructed by (i) setting the water composition in the oxidized part of the aquifer (above the pyrite and SOM containing zone) similar to the recharge water composition that was employed for the simulation period between 1906 and 1956 and (ii) by assuming a more reducing water composition for the deeper parts of the aquifer. The latter was generated on the basis of the recharge water composition, whereby oxygen and nitrate concentrations were decreased to zero before charge-balancing and equilibrating the remaining solution with respect to the prevailing mineral composition. All initial concentrations and associated initial isotope signatures employed in the model are listed in SI Table S1. The defined initial water composition was successively changing during the first years of the model simulations before an almost stable concentration pattern had evolved. Therefore, the first 50 years of simulation time acted as model spin-up period.

Figure 1. Simulated 3H, 3HeTrit and apparent age depth profiles for the years 1986 (dotted), 1996 (dashed) and 2006 (solid lines) in comparison with the measured data (squares) for well 40 from 2006. 10418

dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422

Environmental Science & Technology

Article

patterns developing from purely physical transport in the absence of reactive processes. Under these conditions the simulated temporal and spatial variations of, for example, nitrate and sulfate concentrations mainly reflect the trend of historical variations in the groundwater recharge composition. For nitrate this leads to a peak in the concentration profile at ∼5−10 m a.s.l., where it reached ∼3 × 10−3 mol L−1 (Figure 3). This maximum corresponds to the estimated historical peak of the agricultural activity at the Oostrum site, showing the potential vertical travel distance of nitrate in the absence of reactions and the downward concentration decrease that resulted from mixing/dispersion alone. The lower nitrate concentrations in the younger water above this peak resulted from the introduction of a national manure law that limited the application of nitrate.50 However, nitrate inputs in 2006 were still significantly higher than those under pristine conditions, that is, prior to the start of intensive agricultural activities. In this conservative transport scenario the simulated sulfate concentrations do barely show a pronounced concentration peak (Cmax ∼0.7 × 10−3 mol L−1), that is, the input peak of ∼1.8 × 10−3 mol L−1 was mostly attenuated within the aquifer. Sulfur Dynamics. Compared to the conservative transport scenario the reactive transport simulations clearly illustrate the effect of pyrite oxidation on oxygen and nitrate concentrations, which decrease sharply at a depth of ∼10 m a.s.l., whereas sulfate and dissolved ferrous iron concentrations increase correspondingly (Figure 3). At greater depth sulfate concentrations decrease, partially due to the lower historic sulfate input in the recharge water associated with a decrease in atmospheric deposition after 1965, but more importantly due to (i) the lower historic nitrate input and the accordingly lower sulfate production by nitrate reduction and (ii) sulfate reduction by sedimentary organic carbon in the deeper parts of the aquifer.33 This becomes evident when comparing the reactive simulation profile with the depth profiles for the nonreactive simulation and a modified reactive case in which only SOM degradation was switched off (Figure 3). Simulated sulfate concentrations from the fully reactive model over-

and 2006 in comparison with the measured data from 2006, demonstrating a good agreement between simulated and observed 3H and 3He concentrations at this location and thus the apparent groundwater age. Figure 2 shows the correspond-

Figure 2. Results of the direct age simulation in the model domain and computed steady state streamlines for the end of the simulation period (2006), whereby “age” refers to the time in years since the start of the model simulation in 1906. The location of MLS 40 and two drains are indicated by the yellow bars.

ing age contours for the end of the simulation period (2006), as derived from the direct groundwater age simulation. While generally simulated groundwater ages increase with depth, this pattern is interrupted near the points where groundwater from the shallow aquifer discharges to surface drainage points (Figure 2). In the deeper parts of the aquifer flow velocities decrease locally in the proximity of stagnation points, resulting in increased simulated water ages. The measured 3H and 3He concentrations were used to constrain the transverse vertical dispersivity (best fit for αT = 5 mm), which controls vertical mixing by hydrodynamic dispersion. The longitudinal dispersivity employed for the simulations was 5 m. However, results are relatively insensitive to the exact value. Nonreactive Simulations. Multispecies nonreactive transport simulations were carried out to illustrate the concentration

Figure 3. Depth profiles of selected simulated concentrations (in mol L−1) at MLS 40 and isotope ratios (in ‰) after 100 years of simulation time in comparison with values observed in 2006 (squares). Results of reactive simulations are shown by solid gray lines and the corresponding nonreactive results by dashed gray lines. Dotted red lines represent the results of a comparative simulation in which SOM mineralization under sulfate reducing condition was not considered. 10419

dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422

Environmental Science & Technology

Article

Figure 4. Simulated evolution of nitrate (top panel) and sulfate (middle panel) concentrations (in mol L−1) as well as the sulfur isotope ratio (bottom panel, in ‰) after 50, 75, 90, and 100 years simulation time.

literature-based enrichment factor of 13‰ associated with sulfate reduction by SOM oxidation. A comparative model run in which SOM mineralization and thus sulfate reduction was deactivated shows a depth profile in which the simulated isotope ratios remained constant below the clay layer (Figure 3). An increase in sulfate reduction rates, which would have minimized the discrepancy between simulated and observed sulfate concentrations, would not have only resulted in an overestimation of alkalinity and pH (as discussed above), but also an overestimation of the 34S enrichment (data not shown). Figure 4 illustrates how the increase of sulfate released by pyrite oxidation resulted in a slight decrease of δ34S within and below the pyrite oxidation front between 1956 and 2006. These patterns are interrupted at locations where upward flow toward drainage points occurs. Below the clay layer, where oxygen and nitrate were generally fully depleted, sulfate reduction coupled to the oxidation of SOM resulted in the model in a steady increase of the δ34S with depth (see Figure 3). The measured data show a greater variability, which can most likely be attributed to either (i) smaller-scale geochemical heterogeneities that are not explicitly resolved in the model or (ii) to the fact that a spatially constant recharge water composition was assumed, which neglected the influence of the now ∼1 km long section near Drain 1 that is covered by an urban area and a section downstream of the study site, which is covered by a forested area. Linkages between Recharge, Discharge and Reactive Transformation Rates. The calibrated model was used to analyze the temporal evolution of mass budgets and mass transfer rates and how they responded to the dynamic changes in the groundwater recharge quality. Figure 5 illustrates the major sources and sinks that affect the nitrate and sulfate mass within the aquifer and the corresponding simulated discharge to drains and across the model’s downstream boundary. Simulated sulfate and nitrate input into the aquifer via groundwater recharge show maximum values in 1965 and 1985, which coincide with peaks in atmospheric sulfur deposition and in agricultural nitrate inputs, respectively. Interestingly the peak in nitrate discharge occurred already approximately in 1990, that is, a relatively short time after the input maximum. The reason for this short lag-time is that only the small fraction of the

estimate somewhat those measured concentrations that correspond to the period in which nitrate input had peaked (at approximately −2.5 m a.s.l.). The discrepancy could either stem from a slightly overestimated nitrate input or to some extent also from the incomplete transformation of sulfides to solid phase elemental sulfur rather than dissolved sulfate, a reaction pathway that was suggested by Zhang et al.16 but was not further explored in the present model. An increase in the modeled sulfate reduction rates could have also reduced the simulated sulfate concentrations. However, in additional model scenarios that assumed higher sulfate reduction rates the associated additional mineralization of organic carbon (i.e., SOM) resulted in a clear overestimation of both alkalinity and pH, compared to the measured values (simulation results not shown). Figure 4 illustrates the historical dynamic changes of the sulfate plumes that were mostly created by pyrite oxidation. It can be seen that in 1956 (i.e., after 50 years simulation time) sulfate concentrations were not significantly raised as a result of pyrite oxidation within the clayey aquifer horizon. Even in 1981, around the time when nitrate input peaked, sulfate concentrations were only modestly raised. The highest sulfate concentrations occurred in 2006, despite the fact that nitrate input had decreased. The main reason for this behavior is the approximately 15 years lag period between the maximum nitrate input reaching the groundwater table and its arrival at the deeper, pyrite-containing zones, where sulfate release coupled to nitrate-driven pyrite oxidation occurs. Simulated Sulfur Isotope Signature. As illustrated in Figure 4, the simulated sulfur isotope signature varied strongly in space, but also, though much less pronounced, with time. While in the shallow aquifer parts the simulated sulfur isotope ratio remained at ∼3‰33, corresponding to the ratio that persisted in the sulfate supplied by groundwater recharge (∼2.6‰33), mixing with the sulfate released during pyrite oxidation resulted in a decrease of δ34S to ∼ −4‰ by 2006. This value is close to the mean value of δ34S (−4.2‰) assumed for the solid phase, that is, the δ34S of the pyrite prevailing in the clay layer through which the groundwater penetrates. The observed isotopic enrichment trend in the deeper part of the aquifer was correctly reproduced by using in the simulation the 10420

dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422

Environmental Science & Technology

■

Article

AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS Financial support was provided by The Netherlands Organisation for Scientific Research (NWO-water and Vidi grant 86405.004 to CPS). We thank A. Visser, J. Rozemeijer, A. Dale, D. van de Meent, H. de Waard, G. Nobbe, G. Laruelle, M. Verheul, and G. Klaver for help with the field work. We thank TNO for making available data from an earlier project.

■

Figure 5. Simulated, spatially integrated recharge, discharge and transformation rates of sulfate and nitrate between 1940 and 2040. Marked lines (red ∇ and gray ∇) represent simulated sulfate and nitrate input rates into the aquifer via groundwater recharge. Solid lines represent simulated sulfate and nitrate discharge rates to the two drains and across the model’s downstream boundary. Dashed lines represent simulated sulfate and nitrate discharge rates for the nonreactive model. The cross-dashed line shows the sulfate release rates by pyrite oxidation coupled to nitrate reduction while the open circle line shows the corresponding nitrate consumption rates. The arrows mark the time of the maximum values in sulfate and nitrate recharge and discharge rates.

(1) Strebel, O.; Duynisveld, W. H. M.; Böttcher, J. Nitrate pollution of groundwater in western Europe. Agr. Ecosyst. Environ. 1989, 26 (3− 4), 189−214. (2) Böhlke, J.-K. Groundwater recharge and agricultural contamination. Hydrogeol. J. 2002, 10 (1), 153−179. (3) Zhang, W. L.; Tian, Z. X.; Zhang, N.; Li, X. Q. Nitrate pollution of groundwater in northern China. Agr. Ecosy. Environ. 1996, 59 (3), 223−231. (4) Puckett, L. J.; Tesoriero, A. J.; Dubrovsky, N. M. Nitrogen Contamination of Surficial Aquifers–A Growing Legacy. Environ. Sci. Technol. 2011, 45 (3), 839−844. (5) Burow, K. R.; Nolan, B. T.; Rupert, M. G.; Dubrovsky, N. M. Nitrate in groundwater of the United States, 1991−2003. Environ. Sci. Technol. 2010, 44 (13), 4988−4997. (6) Jørgensen, C. J.; Jacobsen, O. S.; Elberling, B.; Aamand, J. Microbial oxidation of pyrite coupled to nitrate reduction in anoxic groundwater sediment. Environ. Sci. Technol. 2009, 43 (13), 4851− 4857. (7) Descourvieres, C.; Prommer, H.; Oldham, C.; Greskowiak, J.; Hartog, N. Kinetic Reaction Modeling Framework for Identifying and Quantifying Reductant Reactivity in Heterogeneous Aquifer Sediments. Environ. Sci. Technol. 2010, 44 (17), 6698−6705. (8) Wallis, I.; Prommer, H.; Simmons, C. T.; Post, V.; Stuyfzand, P. J. Evaluation of Conceptual and Numerical Models for Arsenic Mobilization and Attenuation during Managed Aquifer Recharge. Environ. Sci. Technol. 2010, 44 (13), 5035−5041. (9) Prommer, H.; Stuyfzand, P. J. Identification of TemperatureDependent Water Quality Changes during a Deep Well Injection Experiment in a Pyritic Aquifer. Environ. Sci. Technol. 2005, 39 (7), 2200−2209. (10) Tesoriero, A. J.; Puckett, L. J. O2 reduction and denitrification rates in shallow aquifers. Water Resour. Res. 2011, 47 (12), W12522. (11) Hartog, N.; Van Bergen, P. F.; De Leeuw, J. W.; Griffioen, J. Reactivity of organic matter in aquifer sediments: Geological and geochemical controls. Geochim. Cosmochim. Acta 2004, 68 (6), 1281− 1292. (12) Descourvières, C.; Hartog, N.; Patterson, B. M.; Oldham, C.; Prommer, H. Geochemical controls on sediment reactivity and buffering processes in a heterogeneous aquifer. Appl. Geochem. 2010, 25 (2), 261−275. (13) Postma, D.; Boesen, C.; Kristiansen, H.; Larsen, F. Nitrate reduction in an unconfined sandy aquifer - water chemistry, reduction processes, and geochemical modeling. Water Resour. Res. 1991, 27 (8), 2027−2045. (14) Pauwels, H.; Kloppmann, W.; Foucher, J.-C.; Martelat, A.; Fritsche, V. Field tracer test for denitrification in a pyrite-bearing schist aquifer. Appl. Geochem. 1998, 13 (6), 767−778. (15) Tesoriero, A. J.; Liebscher, H.; Cox, S. E. Mechanism and rate of denitrification in an agricultural watershed: Electron and mass balance along groundwater flow paths. Water Resour. Res. 2000, 36 (6), 1545− 1559.

recharged nitrate which is transported along shallow flow paths reaches the drains. In other words, only nitrate that was recharged within close proximity to the drains was captured by the drains, whereas denitrification completely removed nitrate during its passage across the reactive zone. Figure 2 illustrates this point, showing that most of the simulated path lines cross the clay horizon. This occurs because hydraulic conductivities in the deeper aquifer parts are higher than in the shallower parts. All nitrate that discharges directly from the shallow, oxidized part of the aquifer does so relatively shortly after infiltration. This is in accordance with the results of Velde et al.,53 who reported a fast nitrate discharge response in lowland areas with an intensive drainage system. In contrast to nitrate, a much larger lag-time exists for the discharge peak of sulfate. The simulated sulfate discharge peak in 2000 occurs ∼15 years after nitrate concentrations peaked in the groundwater recharge. Although this peak is mainly linked to sulfate production by pyrite oxidation, the simulated peak value is clearly lower than the simulated sulfate release rate by pyrite oxidation. The reasons for this are (i) that a large fraction of the generated sulfate is still stored in the aquifer (also illustrated in the sulfate subplot for 2006 in Figure 3) and (ii) that sulfate reduction occurred at a relatively constant rate. The model simulations suggest that, while at this particular site nitrate and sulfate discharge rates will successively further decline over the next 30 years, the largest portion of the response to the changes in agricultural practises has already occurred by now.

■

REFERENCES

ASSOCIATED CONTENT

S Supporting Information *

Additional Tables and Figures are provided in the Supporting Information. This material is available free of charge via the Internet at http://pubs.acs.org. 10421

dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422

Environmental Science & Technology

Article

denitrification in a sandy aquifer. Chem. Geol. 2012, 300−301, 123− 132. (34) Ernst, L. F.; de Ridder, N. A. High resistance to horizontal groundwater flow in coarse sediments due to faulting. Geologie en Mijnbouw 1960, 39, 66−85. (35) Broers, H. P.; Buijs, E. A. The origin of trace metals and arsenic at the drinking water well field. Report TNO-NITG97−198-A (in Dutch). 1997. (36) McDonald, J. M.; Harbaugh, A. W. MODFLOW. A Modular 3D Finite Difference Ground Water Flow Model; U.S.G.S.: Reston,VA, 1988. (37) Prommer, H.; Barry, D. A.; Zheng, C. MODFLOW/MT3DMSbased reactive multicomponent transport modeling. Ground Water 2003, 41 (2), 247−257. (38) Van Rooij, R. Regional modelling of groundwater flow around the Oostrum production well field; TNO: 1995. (39) Schlosser, P.; Stute, M.; Dörr, H.; Sonntag, C.; Münnich, K. O. Tritium/3He dating of shallow groundwater. Earth Planet. Sc. Lett. 1988, 89 (3−4), 353−362. (40) Visser, A.; Schaap, J. D.; Broers, H. P.; Bierkens, M. F. P. Degassing of 3H/3He, CFCs and SF6 by denitrification: Measurements and two-phase transport simulations. J. Contam. Hydrol. 2009, 103 (3−4), 206−218. (41) Goode, D. J. Direct simulation of groundwater age. Water Resour. Res. 1996, 32 (3), 289−296. (42) Williamson, M. A.; Rimstidt, J. D. The kinetics and electrochemical rate-determining step of aqueous pyrite oxidation. Geochim. Cosmochim. Ac. 1994, 58 (24), 5443−5454. (43) Eckert, P.; Appelo, C. A. J. Hydrogeochemical modeling of enhanced benzene, toluene, ethylbenzene, xylene (BTEX) remediation with nitrate. Water Resour. Res. 2002, 38 (8), 1130. (44) Wallis, I.; Prommer, H.; Simmons, C. T.; Post, V. E. A.; Stuyfzand, P. J. Evaluation of Conceptual and Numerical Models for Arsenic Mobilization and Attenuation during Managed Aquifer Recharge. Environ. Sci. Technol. 2010, 44 (13), 5035−5041. (45) Gibson, B. D.; Amos, R. T.; Blowes, D. W. 34S/32S Fractionation during Sulfate Reduction in Groundwater Treatment Systems: Reactive Transport Modeling. Environ. Sci. Technol. 2011, 45 (7), 2863−2870. (46) Bottrell, S. H.; Hayes, P. J.; Bannon, M.; Williams, G. M. Bacterial sulfate reduction and pyrite formation in a polluted sand aquifer. Geomicrobiol. J. 1995, 13 (2), 75−90. (47) Spence, M. J.; Bottrell, S. H.; Thornton, S. F.; Richnow, H. H.; Spence, K. H. Hydrochemical and isotopic effects associated with petroleum fuel biodegradation pathways in a chalk aquifer. J. Contam. Hydrol. 2005, 79 (1−2), 67−88. (48) Visser, A.; Broers, H. P.; van der Grift, B.; Bierkens, M. F. P. Demonstrating trend reversal of groundwater quality in relation to time of recharge determined by 3H/3He. Environ. Pollut. 2007, 148 (3), 797−807. (49) Sikora, E. R.; Johnson, T. M.; Bullen, T. D. Microbial massdependent fractionation of chromium isotopes. Geochim. Cosmochim. Acta 2008, 72, 3631−3641. (50) DePaolo, D. J. Surface kinetic model for isotopic and trace element fractionation during precipitation of calcite from aqueous solution. Geochim. Cosmochim. Acta 2011, 75, 1039−1056. (51) Broers, H. P. The spatial distribution of groundwater age for different geohydrological situations in the Netherlands: implications for groundwater quality monitoring at the regional scale. J. Hydrol. 2004, 299 (1−2), 84−106. (52) Modica, E.; Reilly, T. E.; Pollock, D. W. Patterns and Age Distribution of Ground-Water Flow to Streams. Ground Water 1997, 35 (3), 523−537. (53) van der Velde, Y.; de Rooij, G. H.; Rozemeijer, J. C.; van Geer, F. C.; Broers, H. P. The nitrate response of a lowland catchment: on the relation between stream concentration and travel time distribution dynamics. Water. Resour. Res. 2010, 46 (11), W11534.

(16) Zhang, Y.-C.; Slomp, C. P.; Broers, H. P.; Passier, H. F.; Van Cappellen, P. Denitrification coupled to pyrite oxidation and changes in groundwater quality in a shallow sandy aquifer. Geochim. Cosmochim. Acta 2009, 73 (22), 6716−6726. (17) Baker, K. M.; Bottrell, S. H.; Hatfield, D.; Mortimer, R. J. G.; Newton, R. J.; Odling, N.; Raiswell, R.; Tellam, J. H. Reactivity of pyrite and organic carbon as electron donors for biogeochemical processes in the fractured Jurassic Lincolnshire limestone aquifer, UK. Chem. Geol. 2010, 332, 26−31. (18) Torrentó, C.; Cama, J.; Urmeneta, J.; Otero, N.; Soler, A. Denitrification of groundwater with pyrite and Thiobacillus denitrificans. Chem. Geol. 2010, 278 (1−2), 80−91. (19) Spiteri, C.; Slomp, C. P.; Tuncay, K.; Meile, C. Modeling biogeochemical processes in subterranean estuaries: Effect of flow dynamics and redox conditions on submarine groundwater discharge of nutrients. Water Resour. Res. 2008, 44 (2), W02430. (20) Miotliński, K. Coupled Reactive Transport Modeling of Redox Processes in a Nitrate-Polluted Sandy Aquifer. Aquat. Geochem. 2008, 14 (2), 117−131. (21) Broers, H. P. Release of Ni, Co and As from aquifer sediments at the Oostrum well field. In Groundwater Quality, Remediation and Protection, Herbert, M.; Kovar, K., Eds. Tübingen, Germany, 1998, 36, 106−108. (22) Larsen, F.; Postma, D. Nickel mobilization in a groundwater well field: Release by pyrite oxidation and desorption from manganese oxides. Environ. Sci. Technol. 1997, 31 (9), 2589−2595. (23) Hansen, B.; Thorling, L.; Dalgaard, T.; Erlandsen, M. Trend reversal of nitrate in Danish groundwater − A reflection of agricultural practices and nitrogen surpluses since 1950. Environ. Sci. Technol. 2011, 45 (1), 228−234. (24) Böhlke, J. K.; Denver, J. M. Combined use of groundwater dating, chemical, and isotopic analyses to resolve the history and fate of nitrate contamination in two agricultural watersheds, Atlantic coastal plain, Maryland. Water Resour. Res. 1995, 31 (x), 2319−2339. (25) Böhlke, J. K.; Wanty, R.; Tuttle, M.; Delin, G.; Landon, M. Denitrification in the recharge area and discharge area of a transient agricultural nitrate plume in a glacial outwash sand aquifer, Minnesota. Water Resour. Res. 2002, 38 (10), 1−26. (26) Moncaster, S. J.; Bottrell, S. H.; Tellam, J. H.; Lloyd, J. W.; Konhauser, K. O. Migration and attenuation of agrochemical pollutants: insights from isotopic analysis of groundwater sulphate. J. Contam. Hydrol. 2000, 43 (2), 147−163. (27) Prommer, H.; Anneser, B.; Rolle, M.; Einsiedl, F.; Griebler, C. Biogeochemical and Isotopic Gradients in a BTEX/PAH Contaminant Plume: Model-Based Interpretation of a High-Resolution Field Data Set. Environ. Sci. Technol. 2009, 43 (21), 8206−8212. (28) D’Affonseca, F. M.; Prommer, H.; Finkel, M.; Blum, P.; Grathwohl, P. Modeling of the long-term and transient evolution of biogeochemical and isotopic signatures in coal tar contaminated aquifers. Water Resour. Res. 2011, 47, W05518. (29) Green, C. T.; Böhlke, J. K.; Bekins, B. A.; Phillips, S. P. Mixing effects on apparent reaction rates and isotope fractionation during denitrification in a heterogeneous aquifer. Water Resour. Res. 2010, DOI: 10.1029/2009WR008903. (30) Engesgaard, P.; Kipp, K. L. A Geochemical Transport Model for Redox-Controlled Movement of Mineral Fronts in Groundwater Flow Systems: A Case of Nitrate Removal by Oxidation of Pyrite. Water Resour. Res. 1992, 28 (10), 2829−2843. (31) Frind, E. O.; Duyisveld, W. H. M.; Strebel, O.; Boettcher, J. Modeling of multi-component transport with microbial transformation in groundwater: The Fuhrberg case. Water Resour. Res. 1990, 26 (8), 1707−1719. (32) Howden, N. J. K.; Burt, T. P.; Mathias, S. A.; Worrall, F.; Whelan, M. J. Modelling long-term diffuse nitrate pollution at the catchment-scale: Data, parameter and epistemic uncertainty. J. Hydrol. 2011, 403 (3−4), 337−351. (33) Zhang, Y.-C.; Slomp, C. P.; Broers, H. P.; Bostick, B.; Passier, H. F.; Böttcher, M. E.; Omoregie, E. O.; Lloyd, J. R.; Polya, D. A.; Van Cappellen, P. Isotopic and microbiological signatures of pyrite-driven 10422

dx.doi.org/10.1021/es4023909 | Environ. Sci. Technol. 2013, 47, 10415−10422