Kinetic model parameters used in the SEAM-3D model.
Nitrate contamination of groundwater systems used for human water supplies is a major environmental problem in many parts of the world. Fertilizers containing a variety of reduced nitrogen compounds are commonly added to soils to increase agricultural yields. But the amount of nitrogen added during fertilization typically exceeds the amount of nitrogen taken up by crops. Oxidation of reduced nitrogen compounds present in residual fertilizers can produce substantial amounts of nitrate which can be transported to the underlying water table. Because nitrate concentrations exceeding 10 mg/L in drinking water can have a variety of deleterious effects for humans, agriculturally derived nitrate contamination of groundwater can be a serious public health issue.
The Central Valley aquifer of California accounts for 13 percent of all the groundwater withdrawals in the United States . The Central Valley, which includes the San Joaquin Valley, is one of the most productive agricultural areas in the world and much of this groundwater is used for crop irrigation. However, rapid urbanization has led to increasing groundwater withdrawals for municipal public water supplies. That, in turn, has led to concern about how contaminants associated with agricultural practices will affect the chemical quality of groundwater in the San Joaquin Valley . Crop fertilization with various forms of nitrogen-containing compounds can greatly increase agricultural yields. However, leaching of nitrate from soils due to irrigation has led to substantial nitrate contamination of shallow groundwater . That shallow nitrate-contaminated groundwater has been moving deeper into the Central Valley aquifer since the 1960s . Denitrification can be an important process limiting the mobility of nitrate in groundwater systems . However, substantial denitrification requires adequate sources of electron donors in order to drive the process. In many cases, dissolved organic carbon (DOC) and particulate organic carbon (POC) are the primary electron donors driving active denitrification in groundwater. The purpose of this chapter is to use a numerical mass balance modeling approach to quantitatively compare sources of electron donors (DOC, POC) and electron acceptors (dissolved oxygen, nitrate, and ferric iron) in order to assess the potential for denitrification to attenuate nitrate migration in the Central Valley aquifer.
1.1. Interactions of dissolved organic carbon, oxygen, and nitrate
There are at least three distinct compartments present in groundwater systems that store natural organic carbon capable of reacting with dissolved oxygen, nitrate, and other electron acceptors. Dissolved organic carbon (DOC) is present at varying concentrations in all groundwater systems [5,6,7], and this dissolved compartment can store substantial amounts of organic carbon. In addition to DOC, many groundwater systems contain particulate organic carbon (POC) in various stages of diagenesis [5,6,8,9,10]. Microbial degradation of POC can be an additional source of DOC to soil interstitial water  and groundwater . Finally, silicate, iron oxyhydroxide, and other minerals present in aquifer solids have a capacity to adsorb DOC , removing it from the dissolved phase [13,14,9,15]. These adsorption processes are partially reversible, so that desorption of organic carbon from aquifer materials is also a potential source of DOC . A mass balance model of organic carbon dynamics in groundwater systems, therefore, will need to account for each of these carbon-storing compartments and their interactions.
In contrast to the complexity inherent in the multiple sources, sinks, and composition of DOC, atmospheric oxygen carried through the unsaturated zone by infiltrating precipitation is the sole source of dissolved oxygen (DO) to groundwater systems which lack active photosynthesis. In addition, DO’s relatively limited solubility in fresh water (10.1 mg/L at 15 °C) provides a convenient upper limit to concentrations of DO that can be delivered to the water table. These characteristics will be useful in constraining a quantitative mass balance between DOC and DO.
The interaction of DO with the three compartments of organic carbon present in groundwater systems determines the transformation or lack of transformation of nitrate. The usual ecological succession of electron acceptor utilization in groundwater systems (oxygen>nitrate>Fe(III)> sulfate> carbon dioxide  implies that once concentrations of DO drop below approximately 0.5 mg/L, reduction of nitrate will occur and may coincide with any of the succeeding predominant terminal electron-accepting processes. Constructing a mass balance between the sources and bioavailability of DOC, DO, and nitrate, therefore, is central to assessing the fate and transport of nitrate.
2.1. Study area
The San Joaquin Valley occupies the southern two-thirds of the Central Valley of California (Figure 1), a large northwest-trending, asymmetric structural trough filled with marine and continental sediments up to 10 km thick . East of the valley the Sierra Nevada rise to an elevation of more than 4,200 m. West of the valley are the Coast Ranges, a series of parallel ridges of moderate elevation. Streams in the northern part of the San Joaquin Valley drain through the San Joaquin River northward to the San Francisco Bay. During predevelopment, groundwater generally moved toward the center of the valley where it discharged to the San Joaquin River. However, extensive development of groundwater for agriculture and public water supply has substantially altered the flow system.
The hydrologic system in the Modesto area (Figure 1) is complex, in part because of the heterogeneous nature of the hydrogeological setting. The primary aquifers in the study area are a complex sequence of overlapping alluvial fan deposits that have been eroded from the Sierra Nevada (Figure 2). These alluvial fan deposits consist of coarse-grained sands and gravels with discontinuous clayey silts and clays . A relevant characteristic of these sediments for the present study is that they contain very low amounts of organic carbon, typically in the range of 0.01 to 0.1 weight percent organic carbon (Figure 2). The low sediment organic carbon content reflects the generally arid climate in the recent geologic past, and the overland transport of alluvial fan deposits prior to final deposition.
The relatively low amounts of available organic carbon have resulted in groundwater, that is predominantly oxic  (Figure 1). Most of the individual analyses (Figure 1) of groundwater are either oxic or mixed oxic and anoxic. The specific criteria used to distinguish redox conditions are described in . In general, concentrations of dissolved oxygen tend to decrease and concentrations of dissolved iron tend to increase as groundwater approaches the discharge area near the San Joaquin River.
2.2. MODFLOW model of the San Joaquin Valley aquifer
The numerical model in this study is based on a regional model of groundwater conditions in the Central-Eastern San Joaquin Valley [19, 21]. The U.S. Geological Survey (USGS) three-dimensional finite-difference code MODFLOW-2000  was used to simulate groundwater flow and water-level distributions across the study area (Fig. 1). This regional model was constructed using a three-dimensional grid consisting of 153 rows and 137 columns and 16 layers. The 400 m by 400 m cells were uniform in dimension, and model layers varied from 0.5 m to 16 m above layer 8 and from 20 to 74 m below layer 8. The model area extends from the Coast Ranges to the Sierra Nevada foothills, although the area west of the San Joaquin River was not simulated. The external boundaries of the regional model were a no-flow boundary on the northeastern boundary and general head boundaries on the other three sides. Hydraulic conductivities were estimated based on sediment texture documented in drilling logs. Complete details of the model, including input parameters, calibration, flow budget and travel times are described in .
2.3. SEAM-3D model of the San Joaquin Valley aquifer
The Sequential Electron Acceptor Model in three dimensions (SEAM-3D) was used to construct a quantitative mass balance between electron donors and acceptors in this study . Because of the computational complexity of this mass balance, it was not feasible to model the entire area of the aquifer (Fig. 1). Rather, a cross-sectional model approximately 25 km in length (Fig. 1) was constructed from the regional model of  from the Tuolumne River to the San Joaquin River (Fig. 1). Only the top 11 of the 16 layers of  were included in the SEAM cross-sectional simulations; the bottom 5 layers were below a regional confining unit (Corcoran clay) and were likely to contain water that is too old for evaluating changes in redox conditions as a result of anthropogenic processes. Specified heads cells were used at the eastern and western boundaries of the model to approximately reproduce the configuration of the flow system prior to development. The USGS program MODPATH was used to simulate groundwater travel times. MODPATH is a particle tracking post-processing model that computes three dimensional flow paths using output from groundwater flow simulations based on MODFLOW . The simulated head distribution and approximate times-of-travel for the cross-sectional model are shown in Figure 3. In general, the simulated flow system reflects recharge near the Sierra Nevada in the vicinity of the Tuolumne River and discharge at the San Joaquin River. The undulations of the flowpaths shown in Figure 3 reflect lithological heterogeneities that cause variations in the distribution of hydraulic conductivity. The triangles shown on each of the flowpaths delineated in Figure 3 show the lateral distance traveled by groundwater in 10 years. The shallowest flowpaths have travel times on the order of 50 years and the deepest on the order of hundreds of years.
Equations for the mass balance of bioavailable organic carbon and electron acceptors (EAs) assume that DOC serves as the primary electron donor (carbon/energy source) for a heterotrophic microbial population in the aquifer system. Physical and biogeochemical processes incorporated in equations of transport include advection, dispersion, microbially-mediated biotransformation, rate-limited sorption and desorption. For example, the mass balance equation for bioavailable DOC is given as
where is the concentration of bioavailable DOC in the aqueous phase [M L-3]; is the concentration of bioavailable DOC in the solid phase [M M-3]; is the bulk density of the subsurface material [M L-3];is the average pore water velocity [L T-1]; Dij is the tensor for the hydrodynamic dispersion coefficient [L2 T-1]; is a biodegradation sink term dependent on the mode of respiration [M L-3 T-1]; is the DOC concentration of the source or sink flux [M L-3]; θ is aquifer porosity [-];is distance along the respective Cartesian coordinate axis [L]; t is time [T]; and qs is the volume flow rate per unit volume of aquifer representing fluid sources (positive) and sinks (negative) [T-1].
Mass balance equations of the aqueous phase EAs (DO and nitrate, respectively) are
where and are the aqueous phase concentrations [M L-3] of DO and nitrate, respectively;and are the DO and nitrate concentrations [M L-3] of source or sink fluxes, respectively; and and are the EA biodegradation sink terms [M L-3 T-1], respectively. This treatment assumes the effects of sorption on nitrate transport are small. The consumption of the bioavailable Fe(III) concentration [M M-3] in the solid phase, , is expressed as
Biodegradation of DOC is a function of EA availability and is described using modified Monod kinetics . In summary, the overall approach is to write an equation of mass balance for each individual solute and solid-phase constituent considered and then to solve these equations simultaneously in order to compute a true mass balance as a function of time and space.
2.4. Parameter estimation
The SEAM-3D model was initially parameterized to reproduce the approximate steady-state distribution of dissolved oxygen prior to large-scale agricultural development. A relatively high (6 mg/L) concentration of dissolved oxygen was assumed to enter the aquifer at the eastern boundary, simulating recharge from the Sierra Nevada foothills (Figure 4). In addition, concentrations of sediment organic carbon were fixed at 0.02 weight percent (wt%) throughout the model domain, concentrations of DOC entering the model with recharge were assumed to be 1 mg/L, the half-saturation constant (Ks) and maximum oxygen utilization rate (Vmax) were initially set at 5.0 g m-3 and 0.002 d-1, respectively. The steady-state distribution of dissolved oxygen given these assumptions is shown in Figure 4. Simulated dissolved oxygen concentrations ranged from 6.0 to about 1.0 mg/L and generally decreased with depth.
The next step in constraining the model was to use parameter estimation (PEST) to further refine the monod kinetic parameters . First, water-quality data estimated from predevelopment observations of dissolved oxygen and nitrate were used to refine the Vmax values (Table 1). In a second step, water-quality observations measured along the flowpath of the cross sectional model  (Figure 1) were used to constrain the Vmax values (Table 1). In general, the PEST approach tended to lower Vmax values for the O2-DOC, NO3-DOC, and Fe-DOC reactions. The final Vmax values as constrained by the flowpath water-quality data are listed in Table 1.
|Parameter||Initial Estimate||Predevelopment observations (PEST)||Flowpath observations (PEST)|
|O2 half saturation constant (g m-3)||1.0||1.85||3.0|
|NO3 half saturation constant (g m-3)||1||1||1|
|O2-DOC Vmax (d-1)||0.0002||0.00068||0.0032|
|Fe-DOC Vmax (d-1)||0.00002||0.00002||0.000013|
|Sediment organic carbon (wt%)||0.02||0.02||0.02|
|Organic carbon dissolution rate (d-1)||0.01||0.01||0.0041|
|DOC solubility (mg/L)||1||1||1|
2.5. Sediment organic carbon measurements
The sediments used in this study to quantify concentrations of particulate and adsorbed organic carbon were collected by the USGS using a wire-line core barrel driven into the bottom of a borehole drilled with mud rotary coring methods. Cores were collected from two bore holes named MRWA and MREA located west and east of the city of Modesto respectively (Figure 2). The cores were collected on site, logged, and stored in core boxes prior to transporting them to a storage facility. The cores were subsampled for analysis of organic carbon approximately 4 years after collection. In the storage facility, the cores were broken in half and the center of the core subsampled for analysis in order to minimize the effects of residual drilling mud. Concentrations of sediment organic carbon (Fig. 2) were analyzed with an organic carbon analyzer (Costech Analytical Technologies, Inc., Valencia, California) with a detection limit 0.001 w% organic carbon.
2.6. Modeling organic carbon dynamics in groundwater systems
The soil science literature has extensively considered the dynamics of DOC formation and transport in soils , and this literature is a useful starting point for considering DOC dynamics in groundwater systems. The standard conceptual model of DOC dynamics in soils  considers that the total amount of carbon present at any given time and place reflects DOC and POC delivery from surface sources, production of DOC from bioavailable POC, biodegradation of DOC and POC, and the adsorption/desorption of DOC on soil particles. This conceptual model, in turn, can be used to build a quantitative mass-balance model of organic carbon dynamics in groundwater systems.
DOC mobilized from surface soils is one source of bioavailable DOC to groundwater, but DOC can also be generated from POC present in aquifer material as well . Of all the particulate organic carbon (POCtot) present in an aquifer, only a fraction is available to support microbial metabolism:
where f is the fraction of bioavailable POC. Microbial metabolism of POCbio can generate additional DOC, and this process can be conceived of as a mass-transfer from the particulate to the dissolved phase according to the equation:
Where is the bulk density of the subsurface material [M L-3], is the rate constant for DOC production [T-1]; is the aqueous concentration of DOC in equilibrium with POC at any point in the system [M L-3]. The sorption-desorption of DOC onto aquifer materials such as silicate minerals, ferric oxyhydroxides, and POC, can be approximated using a simple linear isotherm:
where is the concentration of DOC adsorbed to aquifer material; and is the distribution coefficient [L3 M-1]. Once the fraction of bioavailable POCtot is specified, the numerical model used in this paper uses equations 2-7 to iteratively calculate concentrations of bioavailable DOC and POC in the system as a function of time and space. This bioavailable DOC then drives the sequential utilization of electron acceptors (EAs) such as dissolved oxygen (DO), nitrate, and Fe(III). Note that this approach yields a closed mass balance, as envisioned by the model of , for the total amount of organic carbon stored in all three compartments as a function of time. This approach builds on previous numerical methods used to simulate of redox processes in groundwater systems .
3. Results and discussion
The sediments analyzed for particulate organic carbon in this study (Fig. 2) were predominantly coarse to fine-grained poorly sorted sands with inter-bedded lenses of silt and clay. Many of the silt and clay sediments showed visible remains of roots indicating that they represented fossilized soils, which developed in between sedimentation events. The modern surface soils contained between 0.2 and 0.5 wt% organic carbon. The alluvial fan sediments, however, showed very low concentrations of organic carbon, typically between 0.01 and 0.05 wt%. Modern fluvial sediments typically contain well in excess of 0.1 wt% organic carbon. Thus, the San Joaquin alluvial fan deposits contain unusually low amounts of organic carbon. That characteristic, in turn, will affect the fate and transport of nitrate.
The simulated transport of nitrate through the Central Valley aquifer, using the flowpath PEST-estimated kinetic parameters (Table 1) is shown in Figure 5. This simulation assumes that recharge coming in from the agricultural areas east of the San Joaquin River contains 20 mg/L of nitrate. This assumption is consistent with ambient conditions that existed approximately in the year 2000 . The starting point for the simulations that follow, therefore, may be thought of as beginning in 2000.
After ten years of simulation, nitrate at concentrations of about 10 mg/L have moved into the shallow portions of the aquifer, which is consistent with observed nitrate contamination in this system . After 50 years of simulation, the shallowest portion of the aquifer shows nitrate concentrations ranging from 15 to 20 mg/L, and nitrate has reached the discharge area near the San Joaquin River. After 100 years of simulation, nitrate concentrations in portions of the aquifer near the discharge area have risen above 10 mg/L. The results of these simulations predict that, given the kinetic parameters estimated by PEST, nitrate will penetrate deeper into the aquifer for the foreseeable future. The results also suggest that concentrations of nitrate in the deeper portions of the Central Valley aquifer may increase above the 10 mg/L maximum concentration limit (MCL) established by the U.S. Environmental Protection Agency.
Because the kinetic parameters used in the model are subject to uncertainty, the next step in the analysis assessed the sensitivity of the results to changes in those parameters. For this step, the focus was on the cells at the discharge area of the San Joaquin river at the western terminus of the cross-sectional model in layer 6. The first of these simulations assumed that the Vmax for nitrate reduction coupled to DOC oxidation was increased by a factor of 20. This scenario reflects the possibility that the PEST-derived Vmax values could be underestimated. The results are summarized in Figure 6. The decrease in DO concentrations reflects the increase in the the O2-DOC Vmax from 0.0002 to 0.0032 d-1 indicated by the flowpath PEST. As was the case in the simulation shown in Figure 5, nitrate is predicted to begin arriving at the discharge area after about 40 years of simulation. The maximum nitrate concentrations (~10 mg/L) predicted in the simulation of Figure 6, however, were about 20% lower than those shown in the simulation of Figure 5. So, an increase in the NO3-DOC Vmax by a factor of 20 does increase the attenuation of nitrate. However, nitrate concentrations are still predicted to increase substantially at the discharge area over time.
Note the simulated behavior of nitrate and dissolved iron indicated in Figure 6. As dissolved oxygen concentrations decrease (due to the PEST-indicated increase in O2-DOC Vmax), iron concentrations at the discharge area begin to increase. This is consistent with observed detections of dissolved iron near the discharge area (Figure 1). However, as nitrate encroaches on the discharge area, nitrate metabolism begins to replace iron metabolism and iron concentrations begin to decrease. This reflects the design of SEAM-3D which uses the most efficient available electron acceptor, in this case nitrate, preferentially to Fe(III).
In addition to increasing nitrate concentrations in recharge water, agricultural practices such as fertilization, tilling, and irrigation have the capacity to increase concentrations of DOC as well. The next simulation, therefore, explored the sensitivity of the model to increasing concentrations of DOC in recharge water from 1 to 10 mg/L. The higher NO3-DOC Vmax used in Figure 6 was also used. The results are shown in Figure 7 and indicate additional attenuation of nitrate, with maximum concentrations decreasing from 10 mg/L (Figure 6) to about 6 mg/L. DOC concentrations at the discharge area also increase from 1 to 4 mg/L. This, in turn, indicates that increasing concentrations of DOC entering the aquifer with recharge may indeed increase nitrate attenuation. However, the model results still indicate substantial nitrate concentrations reaching the discharge area.
4. Summary and conclusions
Here, a numerical mass-balance modeling approach was used to simulate the long-term fate and transport of agriculturally-derived nitrate in the aquifer system underlying the San Joaquin Valley in California. The SEAM-3D code (Sequential Electron Acceptor Model in three dimensions) used in this study couples the oxidation of dissolved organic carbon (DOC) to the reduction of dissolved oxygen (DO), nitrate (NO3), and ferric iron (Fe(III)) using monod kinetics and including inhibition functions to force the utilization of electron acceptors in the order DO> NO3>Fe(III). A cross-sectional model 25 kilometers in length was constructed by taking the hydraulic conductivity distribution from a calibrated regional model and providing boundary conditions that approximate the historical steady-state distribution of hydraulic head. The model was initially parameterized by matching model-derived concentrations of DOC and DO to historically measured concentrations. The parameterization was then refined using parameter estimation (PEST) on measured point concentrations of DOC, DO, NO3, and dissolved iron (Fe(II)). The parameterized model was then used as a hypothesis-testing tool to evaluate different possible future scenarios of nitrate transport in the Central Valley aquifer. Model simulations using the PEST-derived model parameters indicate that the amount of dissolved and particulate organic carbon available in the aquifer is inadequate to consume the amount of DO that typically recharges the aquifer, and that NO3 derived from agricultural activities will be drawn deeper into the aquifer in the foreseeable future. Model simulations that increase the assumed rate of nitrate reaction with DOC by a factor of 20 decrease simulated concentrations of nitrate near the discharge area of the aquifer by about 20 percent, but simulated concentrations are still substantial.
Nitrate concentrations have increased substantially in shallow groundwater of the San Joaquin Valley in recent years. This phenomenon has led many public supply well operators to tap progressively deeper groundwater in order to avoid nitrate contamination. That practice, in turn, has led to concern that elevated nitrate concentrations will continue to be drawn deeper into this groundwater system over time. The current study used a mass balance modeling approach to assess the possible transport and attenuation of nitrate in this system. The results indicate that the amount of available organic carbon in this system, either DOC or particulate organic carbon in aquifer solids, is inadequate to fully attenuate the nitrate that is now entering the shallow portion of the aquifer. This finding, in turn, suggests that nitrate contamination in the Central Valley aquifer will continue to move deeper into the system, and may eventually reach the discharge area near the San Joaquin River.
This research was supported by the National Water-Quality Assessment (NAWQA) program of the U.S. Geological Survey. Use of trade names is for identification purposes only and does not constitute endorsement by the U.S. Government.