Groundwater residence time estimates obscured by anthropogenic carbonate

See allHide authors and affiliations

Science Advances  21 Apr 2021:
Vol. 7, no. 17, eabf3503
DOI: 10.1126/sciadv.abf3503


Groundwater is an important source of drinking and irrigation water. Dating groundwater informs its vulnerability to contamination and aids in calibrating flow models. Here, we report measurements of multiple age tracers (14C, 3H, 39Ar, and 85Kr) and parameters relevant to dissolved inorganic carbon (DIC) from 17 wells in California’s San Joaquin Valley (SJV), an agricultural region that is heavily reliant on groundwater. We find evidence for a major mid-20th century shift in groundwater DIC input from mostly closed- to mostly open-system carbonate dissolution, which we suggest is driven by input of anthropogenic carbonate soil amendments. Crucially, enhanced open-system dissolution, in which DIC equilibrates with soil CO2, fundamentally affects the initial 14C activity of recently recharged groundwater. Conventional 14C dating of deeper SJV groundwater, assuming an open system, substantially overestimates residence time and thereby underestimates susceptibility to modern contamination. Because carbonate soil amendments are ubiquitous, other groundwater-reliant agricultural regions may be similarly affected.


Billions of people are reliant upon groundwater as a source of drinking water (1, 2) and for irrigating nearly half of all global agriculture (3). Although the vast majority of Earth’s groundwater is ancient, having been recharged centuries to millennia ago (4, 5), it remains vulnerable to anthropogenic contamination (6) and depletion (710) due to increasing groundwater extraction. Accurately determining groundwater residence times is important for understanding groundwater contaminant sources and transport pathways, calibrating models of groundwater flow that aid in sustainable management, and contextualizing human impacts on groundwater systems (11).

Here, we investigate groundwater residence time distributions in the central-eastern San Joaquin Valley (SJV) using a combination of geochemical modeling and analytical constraints informed by multiple independent age tracers. We measured radiocarbon (14C), tritium (3H), and a comprehensive suite of groundwater geochemical parameters from a network of 17 wells (36.65°N to 36.75°N and 117.55°W to 117.85°W) south of Fresno, California, which are or have been surrounded by irrigated farmland over the past 50 years (fig. S1). We also measured noble gas age tracers 39Ar (n = 3) and 85Kr (n = 4), which have historically been underused for groundwater dating (12, 13), in a subset of wells to independently explore potential chemical controls on 14C. The SJV is California’s largest and most productive agricultural region, and it is particularly vulnerable to contamination and long-term declines in groundwater reserves due to drought, increased rates of extraction, and seepage from irrigation (10, 1416). The natural groundwater system in the SJV has been fundamentally altered by agricultural activities over the past several decades including both enhanced discharge and recharge due to groundwater pumping and irrigation seepage, respectively (10, 1619). Given the influence of modern agriculture on groundwater in the SJV, this region serves as an ideal field-scale laboratory to assess the impact of anthropogenic perturbations to groundwater chemistry on residence time estimates.

The primary aim of this study is to understand natural and anthropogenic controls on the radiocarbon content of dissolved inorganic carbon (DIC), which is the most widely used groundwater age tracer (20). Measured trends in inorganic carbon chemistry in the SJV are modeled and interpreted in combination with multiple independent age tracers, from a subset of wells to precisely isolate geochemical controls on 14C activity. The use of multiple age tracers that are sensitive across a range of time scales, including inert dating tools like 39Ar, provides key constraints for disentangling the impacts of groundwater mixing and DIC chemistry on 14C. Using these data, we model the extent to which agricultural perturbations have affected 14C of DIC and, in turn, affected 14C-based groundwater dating. Our results have broad implications for the susceptibility of deep groundwater to modern contamination.

Radiocarbon dating of groundwater

Radiocarbon activity of DIC is a ubiquitous groundwater dating tool because of its relative ease of sampling, inexpensive analytical cost, and applicability across a wide temporal range. 14C has not only a 5730-year half-life (T1/2) but also a modern atmospheric excess from 20th century bomb testing (21), rendering it dually sensitive on time scales of both thousands of years and recent decades. However, the direct translation of measured 14C activity to groundwater residence time is frequently complicated by 14C-free inputs of DIC (22), isotopic exchange of DIC with soil CO2 (23, 24), and by the mixing of groundwater water parcels of different age (2527). Using 14C to estimate groundwater residence times therefore requires conceptual models both for DIC chemistry and for groundwater flow and mixing.

The application of 14C of DIC as a groundwater dating tool requires precise knowledge of the initial 14C activity of DIC (C0 14) at the time of recharge. The measured 14C activity (Cm 14) in groundwater (in the absence of mixing between different-aged waters) is linked to residence time (τ) by the decay equationC14m=C140eλτ=ΓC14atmeλτ(1)where λ is the decay constant [0.121 ka−1 (thousand years)], Catm 14 is the atmospheric 14C activity at the time of recharge, and Γ is a scalar that relates C0 14 to Catm 14 such that C0 14=Γ Catm 14. While the history of Catm 14 over the relevant time range for groundwater dating is well constrained (21, 28), input of 14C-free DIC sources to groundwater can dilute the overall 14C concentration in DIC (by lowering the 14C/12C ratio) to values considerably lower than Catm 14. In this sense, Γ represents the effective lowering of groundwater C0 14 from the contemporary atmospheric value [~100 pmc (percent modern carbon) before the 1950s and >100 pmc since then], which is essential for separating signals controlled by DIC chemistry from decay-related 14C signals and thereby accurately quantifying residence times.

Groundwater DIC is intrinsically linked to exchange between dissolved and gaseous soil CO2. The respiration of soil organic carbon generates steady-state CO2 partial pressures (pCO2) above the water table that often exceed pCO2 in the atmosphere by an order of magnitude or more (22). This large reservoir of soil CO2 contributes to groundwater DIC both directly, by dissolving in water, and indirectly, through the generation of acidity that leads to weathering of minerals present in soil and aquifer rocks. Differences in isotopic composition between these sources of DIC partially govern both the stable (13C and 12C) and radioactive (14C) carbon isotopic compositions of DIC. However, (i) gas exchange between DIC and CO2 and (ii) equilibrium fractionation among DIC species further modify the isotopic composition of DIC. Where carbonate minerals are present in soils and rocks, it is commonly assumed that DIC in groundwater is derived from two primary sources: 14C-rich soil CO2 and 14C-free carbonate minerals. Soil CO2 is characteristically depleted in 13C (δ13C between −15 and −30‰VPDB (22)) and is typically assumed to have a roughly atmospheric 14C activity arising from the respiration of recently formed organic matter. The 14C activity of CO2 and DIC in the vadose zone can be somewhat lower due to the production of CO2 from deeper and older organic material (29) and from drought-related impacts on water-sediment interaction (30). Conversely, carbonate minerals are generally 14C-free, with δ13C near 0‰VPDB (20, 22). For a closed system in which carbonate dissolution occurs beneath the water table in isolation from soil CO2, the dissolution of carbonates tends to increase DIC concentrations and δ13C while decreasing 14C activity by adding 14C-free DIC with δ13C of ~0‰ to the DIC pool (31, 32). Thus, relatively young groundwater that has experienced closed-system carbonate dissolution may have lower 14C activity and higher δ13C than DIC acquired purely from dissolution of soil CO2. This assumption forms the basis of conventional closed-system dissolution models for groundwater 14C dating that are constrained by observed values of δ13C (31, 32) or DIC (33). For an open system in which carbonate dissolution occurs in the presence of soil CO2, isotopic equilibration between soil CO2 and DIC leads to a near atmospheric 14C activity of recently recharged groundwater (22). More complex models that involve complete or partial equilibrium isotope exchange and fractionation (23, 24, 34) and account for multiple reactions and processes (35) that alter 14C have also been developed. In each conceptual model, the nature of interaction between soil CO2, groundwater DIC, and mineral weathering collectively determine the isotopic composition of DIC, pH and total DIC content.


Systematic controls on radiocarbon activity

Across the 17 SJV wells sampled in this study, we observe unexpected, yet coherent, relationships between groundwater 14C, DIC, and pH (Fig. 1). Specifically, 14C activity increases with DIC, from ~40 pmc at ~2 mM DIC to ~120 pmc at ~10 mM DIC, in direct contrast to the common inverse relationship expected for closed-system dissolution of 14C-free carbonates (33), in which the progressive dissolution of carbonates in isolation from soil CO2 yields an anticorrelation between 14C activity and DIC, pH, and δ13C (36). In the SJV samples, 14C activity consistently decreases with increasing pH, and DIC is positively correlated with Mg, Ca, and Sr concentrations (Fig. 1 and fig. S8), consistent with the general expectation for carbonate dissolution. Mean 3H activities (T1/2 = 12.5 years) indicate the presence of predominantly decade-old water (3H activity >1 TU; where 1 TU = 1 tritium unit = 0.118 becquerel per liter) in wells where high DIC and low pH (<7.5) are observed. Despite the close correspondence between 14C, DIC, and pH, we observe no coherent relationship between any of these parameters and the stable carbon isotopic composition of DIC (δ13C). Together, closed-system dilution of DIC by 14C-free carbonate dissolution—the common expectation for natural systematic 14C variability in many global aquifers (37)—cannot explain these observed trends in the SJV.

Fig. 1 Relationships between carbon-14 of DIC and key geochemical parameters.

SJV groundwater 14C activity (A) increases as a function of DIC and Mg2+ + Ca2+ concentrations and (B) decreases with pH such that tritiated samples (3H activity > 1 TU) tend to have pH lower than 7.5 and 14C greater than 100 pmc, while (C) shows no coherent relationship with δ13C. Note the linear and log scales of color bars in (A) and (B), respectively. (note that DIC data are unavailable for well 180-4). The gray arrows indicate the conventional expected direction of change for each variable shown in these panels due to closed-system carbonate dissolution, which is clearly inconsistent with the data in (A) and (C).

The highest observed 14C activities (~120 pmc) in tritiated modern samples clearly indicate that DIC of recently recharged groundwaters experienced near-complete isotopic equilibration with soil CO2, as the mean northern hemisphere atmospheric 14C activity between 1955 and 2010 was ~125 pmc (21). The δ13C values associated with high-14C (>100 pmc) samples fall within a −16.5 to −13.5‰ range, also consistent with open-system equilibrium exchange between DIC and biogenic soil CO2 after accounting for equilibrium isotopic fractionation between DIC and CO2 (appendix S2). While δ13C values in lower-14C samples are slightly below those of the higher-14C samples, we note that the predicted sensitivity of δ13C to open- and closed-system dissolution is small relative to potential end-member variability (appendix S2 and fig. S5). Given the high potential for major changes in carbonate and soil CO2 end-member δ13C changes driven by fundamental reconfiguration of the SJV groundwater system due to irrigation-derived recharge over the past century, we suggest that shifts in end-member δ13C have likely overprinted meaningful signals related to open- or closed-system dissolution. However, the conspicuous relationships between 14C and DIC, and 14C and pH, with marked concave-down curvatures, strongly indicate that a systematic process besides progressive closed-system carbonate dissolution—and not radioactive decay—controls much of the observed 14C variability in this system. Here, we investigate whether DIC chemistry and physical transport and mixing processes can account for the observed relationships between 14C, DIC, and pH and evaluate the implications for the reliability of 14C-based groundwater residence times.

Typically, initial 14C activity (C0 14 in Eq. 1) is inferred from various combinations of measured DIC, δ13C, and major and trace ion concentrations by using models of varying complexity (24, 31, 32, 34, 38). Among existing single-sample models used to estimate C0 14 (38), carbonate dissolution is assumed to occur under either entirely open- or closed-system conditions, whereby a portion of DIC isotopically exchanges with soil CO2 (23, 24, 34) or carbonate minerals (23, 24), while the remainder of DIC does not. Conventionally, an appropriate model is chosen and validated based on observed trends in 14C, δ13C, and DIC in modern groundwater with high 3H activity (39). However, if fundamental changes to the DIC system have occurred over recent decades due to human activity, observed trends in modern groundwater (e.g., increases in 14C activity with DIC) may not be representative of the preceding natural processes in older groundwater. In the specific context of SJV measurements (Fig. 1), it is therefore critical to understand whether the inference of near-atmospheric 14C0 for modern groundwater (Γ ~ 1 based on the >100 pmc 14C activities of tritiated samples) is representative of natural, predevelopment 14C0. In other words, if it is assumed that the 14C0 of deeper, 3H-free groundwater is similarly near-atmospheric (Γ ~ 1), then the lowest observed 14C activities imply mean residence times of thousands of years; however, if this assumption is invalid (e.g., Γ < 1), then 14C measurements would imply considerably shorter residence times.

Insights from multiple age tracers

Radioactive noble gas isotopes such as 39Ar (T1/2 = 269 years) and 85Kr (T1/2 = 10.8 years) are inert and thus insensitive to any chemical process linking 14C, DIC, and pH, making them ideal for understanding groundwater flow and mixing processes (13). While 85Kr is sensitive on the same decadal time scale as tritium, 39Ar is a unique intermediate age tracer, filling a gap between tritium and 14C (40). Because of their conceptual simplicity as monatomic, unreactive dissolved gas tracers with well-constrained input functions (40, 41), 85Kr and 39Ar offer valuable physical constraints but have remained underused due to analytical challenges (12, 13), with only a handful of studies having used this set of age tracers in the same groundwater samples [e.g., (42, 43)]. In our SJV study area, these noble gas age tracers constitute a powerful tool to independently evaluate whether the observed low 14C activities in low-DIC, high-pH, 3H-free groundwaters (Fig. 1) arise from systematic chemical processes that loweredC0 14. In this study, we sampled a subset of six wells for analysis of inert age tracers in addition to 14C and 3H to test for consistency among independent age tracers under physically realistic mixing scenarios.

Prior work in the eastern SJV has demonstrated that the unconfined hydrogeologic setting is well characterized by two distinct regimes: (i) a deeper (3H free) portion of the aquifer recharged naturally by precipitation and river runoff before major agricultural and urban development and (ii) a shallower, mostly tritiated (i.e., nonzero 3H activity), portion recharged over the past 50 to 100 years at substantially higher rates by local irrigation seepage (16). We refer to these regimes as predevelopment and postdevelopment, respectively. Notably, the modern rates of irrigation-derived recharge in the SJV are known to exceed natural recharge by roughly sevenfold (10, 18). This substantial influx of agriculturally affected water also flows laterally (19), therefore influencing groundwater chemistry in more urban regions near irrigated farmland (17). On the basis of 3H, the interface between these two regimes in the study area appears to fall between 50 and 100 m below the land surface (fig. S2b), similar to previous findings to the north (16). Here, we first focus on three wells for which 39Ar was measured, which are screened entirely in the predevelopment regime and thus provide the simplest opportunity to compare inert tracers with 14C without needing to account for mixing between pre- and postdevelopment groundwater.

Direct comparison of 14C-derived and 39Ar-derived age estimates requires physical assumptions about mixing of different aged groundwaters in each well (appendix S1). All three deep wells are 3H-free, with 39Ar in the range ~23 to 50 pmAr (percent of modern atmospheric 39Ar activity) and 14C in the range ~ 43 to 69 pmc (table S1). It is readily apparent that residence time estimates based on 39Ar and 14C that neglect mixing (i.e., Eq. 1) are incompatible with a 14C scaling factor Γ of 1, as both 39Ar and 14C observations across each well fall in the range of ~50% modern activity, although the half-lives of the two tracers differ by an order of magnitude. However, even when groundwater mixing is rigorously accounted for by using a physical model, Γ must be well below unity in each well to reconcile 14C with 39Ar (Fig. 2). During pumping, water is drawn into each well through a screened interval, inducing mixing that leads to a distribution of groundwater ages in the sample. In the deeper predevelopment portion of the aquifer, an exponential age-depth profile is expected based on the unconfined nature of the aquifer (16, 44, 45). Our parameterization of mixing assumes a partial-exponential distribution of groundwater ages (e.g., in Fig. 2B) that represents the fraction of each water mass sampled by a given well. We then (i) convolve this distribution with the 39Ar and Catm 14(21, 28) atmospheric histories (e.g., in Fig. 2A) using different scaling factors for C0 14 and (ii) account for radioactive decay to predict present-day measured 39Ar and 14C activities in these wells. We constrain the age distributions based on known well geometry and for maximum agreement between simulated and measured 39Ar such that an optimal 14C scaling factor Γ can then be independently determined by comparing simulated 14C with observations for each well (appendix S1). Figure 2 shows the results of this analysis for well 180-1, as an example. For all deep wells (figs. S7, S10, and S11), we find that all optimal scaling factors Γ are substantially below unity (0.44 to 0.73). Critically, we find that the assumption of no 14C dilution (Γ = 1) yields mean residence times between 3200 and 8300 years for these wells, whereas 39Ar-optimized Γ values (Γ < 1) yield smaller mean residence times of 240 to 840 years (table S2).

Fig. 2 Multiple tracer age comparison for Well 180-1.

(A) Initial groundwater 14C activity histories for differing 14C dilution factors (Γ; Eq. 1) values are convolved with (B) a modeled distribution of groundwater ages, constrained by 39Ar, to (C) simulate tracer activities for direct comparison with measurements. Tracer activities are shown in units of pmc, pmAr, dpm per cm3, and TU for 14C, 39Ar, 85Kr, and 3H, respectively. Open and closed circles show simulated 14C activities associated with atmospheric (Γ = 1) and subatmospheric (Γ < 1) 14C0, respectively.

Two shallower wells with screened intervals from 49 to 95 m and 66 to 69 m below the surface were found to have detectable 85Kr activity and 3H above 1 TU, indicating substantial contributions from modern recharge (figs. S8 and S9). Therefore, these wells likely represent more complicated bimodal distributions of waters from both predevelopment and postdevelopment regimes (16), and a single-distribution modeling approach is thus overly simplified. However, comparison of 85Kr with 14C provides qualitative insights into potential temporal changes in Γ factors. Following the same partial exponential mixing (PEM) modeling approach for the deeper wells but using 85Kr instead of 39Ar to tune age distributions, we find considerably higher Γ factors for these wells (0.85 and 0.91). Therefore, our inert tracer-14C age comparison case study points toward a shift from low C0 14 to higher C0 14 in both shallow and deep wells. This is interpreted to be associated with the transition from the natural, predevelopment to modern recharge systems. Because a failure to account for this modern shift (i.e., by assuming Γ = 1) would yield 14C-derived residence times in the deeper wells that are thousands of years older than those independently implied by 39Ar, it is essential to understand the underlying mechanisms.


Midcentury shift to open system dissolution

To explain the apparent shift in C0 14 in groundwater over recent decades, we propose a simple conceptual DIC model that largely reproduces observed trends in 14C, DIC, and pH, as well as the findings of low 14C scaling factors Γ in deeper wells. The model assumes that groundwater DIC is acquired through varying proportions of open-system (i.e., in the vadose zone, exposed to soil CO2) and closed-system (i.e., below the water table) carbonate mineral dissolution (appendix S2). The model also assumes that when carbonate minerals dissolve in the vadose zone, the isotopic composition of resultant carbonate-derived DIC is nearly atmospheric (Γ ~ 1) due to rapid isotopic equilibration with soil CO2 via gas exchange [(22) and appendix S3]. Conversely, when carbonate dissolution occurs in isolation from soil CO2 under closed-system conditions (i.e., beneath the water table), the resultant 14C-free carbonate-derived DIC dilutes the 14C activity of the DIC pool, thereby driving down C0 14.

In the DIC model, the extent of “system openness” is parameterized on the basis of the total DIC content of a groundwater parcel at the water table (DICiso), acquired under open-system conditions, before it is isolated from soil CO2 by displacement beneath the water table due to subsequent recharge. Before any carbonate dissolution, the equilibrium dissolution of soil CO2 in water provides a minimum DIC (DIC0) in solubility equilibrium with soil CO2. If carbonate dissolution occurs rapidly enough to achieve equilibrium at or above the water table, before isolation, then groundwater will obtain a maximum open-system DIC (DICeq,open) in isotopic equilibrium with soil CO2. However, if groundwater is isolated from soil CO2 before reaching DICeq,open, then additional closed-system carbonate dissolution below the water table will consume acidity and dilute the DIC pool with 14C-free input without replenishment by soil CO2, thereby limiting the equilibrium DIC (DIC < DICeq,open) and raising the equilibrium pH. Here, we define a system openness parameter (f) as the fraction of DIC that is acquired under open-system conditions before isolation from overlying soil CO2fDICisoDIC0DICeq,openDIC0(2)

The openness parameter (f) ranges from 0 (fully closed system) to 1 (fully open system). Figure 3 shows measured groundwater data alongside the modeled relationship between 14C activity and pH for varying system openness (f). At greater f, 14C is higher, a larger fraction of DIC exists as CO2(aq), and, correspondingly, the pH is lower. In Fig. 3, we demonstrate that this direct relationship between f and pH—rather than radioactive decay—can account for the inverse relationship between 14C and pH observed in Fig. 1. Notably, our idealized carbon model predicts a positive, concave-down curved relationship between 14C and DIC, similar to the data plotted in Fig. 1 (fig. S4) with a small range of variability in δ13C relative to source endmember uncertainties (fig. S5). f is quasi-conservative with respect to mixing such that the DIC and 14C activity of an admixture of waters with different f values will closely approximate the model-predicted DIC and 14C activity for their average f (fig. S7). In Fig. 3, comparison of measured data with the model curves suggests a shift at f ~ 0.5 such that groundwaters with substantial modern recharge fractions (3H > 1 TU) are best explained by mostly open-system dissolution (f > 0.5) with near-modern atmospheric 14C of soil CO2 (120 pmc). In contrast, older samples (3H < 1 TU) are best explained by mostly closed-system dissolution (f < 0.5) with a predevelopment 14C of soil CO2 (100 pmc).

Fig. 3 Measured groundwater 14C activity and modeled 14C activity curves versus pH (or CO2(aq)/DIC) for a soil pCO2 of 21,000 μatm.

Most tritiated samples (<50 years mean recharge age) fall on the modeled curve for high values of f (>0.5) and an initial 14C of mean modern atmospheric 14C (120 pmc). Most 3H-free samples fall onto or below the modeled curve, indicating greater degrees of closed-system equilibration and contributions from radioactive decay. Tritiated samples with high DIC and low pH are interpreted to result from modern open-system equilibrium at higher soil pCO2.

We hypothesize that the central-eastern SJV experienced a shift from mostly closed-system to mostly open-system carbonate dissolution regimes after the mid-1950s (Fig. 4). This conceptual framework suggests a major shift in C0 14 from as low as 50% atmospheric (Γ ~ 0.5) before the 1950s to nearly 100% atmospheric (Γ ~ 1) since the 1950s. Thus, for observed predevelopment (deeper) groundwater 14C activities between ~40 and 70 pmc, the contribution of radioactive decay to lowering 14C is estimated to only be ~1 to 10 pmc, resulting in 14C-derived age estimates that are closer to those implied by 39Ar (table S2).

Fig. 4 Overview of system openness hypothesis.

Illustration of hypothesized shift from natural carbonate dissolution (left) within the saturated zone (predominantly closed-system DIC exchange) to modern conditions where agricultural-lime influenced dissolution leads to predominantly open-system DIC exchange (right) through the vadose zone. In both panels, red arrows indicate the flux of soil CO2 to groundwater DIC, gray arrows indicate the flux of carbonate dissolution (natural or anthropogenic) to the DIC pool, and the dashed arrows indicate the direction(s) of isotope exchange. Under closed-system conditions, there is no isotopic exchange between DIC and soil CO2 during carbonate dissolution; hence, this arrow is unidirectional on the left but bidirectional on the right. Dots on the right panel represent finely ground agricultural carbonate particles that dissolve throughout the unsaturated zone and at the water table in the presence of CO2.

We caution, however, that while this idealized carbon model predicts coherent trends that are consistent with system-wide observations, it neglects processes that may be important for individual wells. For example, (i) silicate weathering likely accounts for some portion of DIC, (ii) observations of low modern pH imply heterogeneity in modern soil pCO2 such that some locations exceed the uniform value (2.1%) prescribed in the model, (iii) enhanced open-system dissolution of natural carbonates driven by acidic agricultural infiltration could lead to excess DIC with modern atmospheric 14C activity, and (iv) soil 14C activities may have been somewhat lower than the atmosphere in the predevelopment system because modern irrigation seepage likely more efficiently delivers modern organic carbon to the vadose zone. For example, a greater contribution of soil CO2 produced from respiration of older organic matter deeper in the vadose zone may have further lowered predevelopment 14C0 (29).

Our model is similar to several existing models for partially open-system dissolution of carbonate minerals (23, 24, 34) but differs in a subtle but important way. These existing models assume that DIC is composed of (i) dissolved CO2 in isotopic equilibrium with soil CO2 and bicarbonate derived from carbonate mineral dissolution that is either (ii) in isotopic equilibrium with CO2 or (ii) isolated from CO2 (thereby retaining the carbon isotopic composition of the carbonate minerals). These models parameterize an absolute concentration (q) of bicarbonate in equilibrium with soil CO2. Unlike these models, our approach defines a fraction (i.e., f) of the maximum equilibrium concentration of DIC that can be acquired via open-system carbonate dissolution (DICeq,open) relative to the minimum concentration of DIC in which only soil CO2, but no carbonate minerals, dissolves in the unsaturated zone. Our approach hence enables prediction of the evolution of carbon isotopic composition and total DIC as a function of f (i.e., system openness), whereas previous models do not directly parameterize total DIC as a function of q. The curved relationship between 14C and DIC that our model predicts is an important clue in interpreting these groundwater data. However, it is crucial to note that, like previous models, our idealized approach accounts for just one of many potentially important processes. We suggest that idealized models and simple graphical approaches [e.g., (39)] are helpful in identifying dominant processes and complement more complex geochemical models [e.g., (35)] that include multiple mechanisms but are often underconstrained.

The indication of a general modern shift toward greater system openness points to an anthropogenic influence on groundwater chemistry, which we suspect is linked to irrigation seepage into shallow groundwater. While such a shift from closed- to open-system conditions directly affects the DIC of wells that draw mixtures of pre- and postdevelopment groundwater (i.e., wells screened throughout the ~50- to100-m interface between these two regimes; fig. S2), it also substantially affects the interpretation of 14C in deep wells that are unaffected by direct mixing with modern groundwater. That is, conventionally inferring that a near-atmospheric 14C0 (Γ ~ 1), consistent with observations of high 14C in shallow modern wells, applies to deeper wells would yield erroneously long residence times. In this sense, this modern shift to system openness masks the closed-system nature of deeper, predevelopment groundwater and thereby obfuscates 14C dating of these wells without the additional insight from inert age tracers and consideration of broader trends in parameters related to DIC chemistry.

Like many other major agricultural regions worldwide (46, 47), human activities have substantially influenced groundwater in the SJV since the mid-20th century. Groundwater overdraft has induced meters of land subsidence (15), and land-use changes have both increased groundwater salinity (14) and led to contamination by agriculture (48, 49). Here, we propose that irrigation seepage has directly influenced groundwater DIC in the SJV through the delivery of carbonate soil amendments, known as agricultural lime, to shallow groundwater. Agricultural lime (finely ground calcite, dolomite, and/or magnesite; hereafter referred as AL) is ubiquitously used as a buffer for acidic soils on farmlands in the SJV (50) and globally. Inefficient tillage and limited depth mobility of AL in crop soils (51) can lead to surface excesses of these particulate carbonate minerals, which can be advected by precipitation or irrigation seepage into the subsurface. Crucially, the availability of easily dissolvable carbonate in the vadose zone, where pCO2 is far higher than the atmosphere even at shallow depths, is key to raising the equilibrium groundwater DIC content via open-system dissolution. The prevalence of farmland throughout the study area (fig. S1) raises the possibility that the spatial distribution of recharge may have changed over recent decades, providing widespread pathways for alteration of groundwater chemistry and delivery of AL to the vadose zone. Previous work has identified prominent geochemical signals associated with AL dissolution in groundwater recharged by irrigation seepage (52). Because AL is designed to dissolve rapidly, due to its high surface area, transport of this particulate carbonate through the vadose zone would result in enhanced rates of open-system carbonate dissolution relative to the natural predevelopment system, in which carbonate dissolution in subsurface rocks and soils was slower and thus occurred predominately beneath the water table. Although the 14C activity of AL (sourced from old, naturally occurring, 14C-free carbonate minerals) is negligible, rapid isotopic equilibration of DIC with the large reservoir of soil CO2 via gas exchange would efficiently reset the 14C activity of DIC to a near-atmospheric value over a time scale of weeks to months (appendix S3). Figure 4 presents a schematic of pre-1950s closed-system and post-1950s open-system carbonate dissolution, where the mid-20th century shift is driven by anthropogenic carbonate input to groundwater.

The hypothesized agricultural influence on modern groundwater chemistry is further supported by measurements of elevated Sr2+, NO3, and SO42− accompanying high DIC, Mg2+, and Ca2+ in low-pH, tritiated groundwaters (fig. S8). While NO3 and SO42− are common indicators of agricultural seepage (49, 53), Sr2+ is associated with AL dissolution due to its compatible valence and thus frequent incorporation into carbonate minerals (52). Irrigation seepage carrying particulate agricultural carbonates to the vadose zone thus provides a compelling mechanism for a modern shift toward greater open-system carbonate dissolution, leading to near atmospheric 14C activities and higher DIC in modern tritiated groundwaters (Fig. 1). Critically, if this widely used soil amendment (i.e., AL) has led to masking of natural closed-system dilution of 14C throughout the SJV, corresponding to thousands-of-year overestimations of groundwater age, it likely affects the use of 14C as an age tracer in other global aquifers that are relied upon for major agricultural production worldwide.

Implications for groundwater dating

In this study, we observe two independent lines of evidence indicating a major mid-20th century shift in groundwater DIC in the central-eastern SJV: a residence-time modeling case study constrained by inert age tracers and a wider consideration of geochemical trends corroborated by an idealized DIC model. These two approaches point to a transition from mostly closed-system to mostly open-system dissolution of carbonates, which we hypothesize is due to anthropogenic input of AL. This modern shift in DIC chemistry has masked underlying natural closed-system carbonate dissolution, which led to low C0 14 under predevelopment conditions. A key consequence of this finding is that without multiple age tracers and a mechanistic consideration of regional controls on inorganic carbon chemistry, apparent 14C-based groundwater residence times in the deeper, predevelopment portion of the aquifer (under the assumption that Γ = 1, as suggested by modern tritiated water), would be overestimated by hundreds to thousands of years. Specifically, we find that inert 39Ar-derived mean age estimates in deeper wells, which range from ~240 to 840 years, are an order of magnitude younger than those implied by 14C (if Γ = 1), which range from ~3200 to 8300 years.

Erroneously long estimates of recharge time scales would invariably lead to underestimation of the susceptibility of deeper groundwater to modern contamination and can affect calibration of groundwater flow models that are important for sustainable management. A new understanding that the deeper wells considered in this study have mean residence times of centuries, rather than millennia, raises the concern that continued groundwater withdrawal in the SJV and depression of the regional water table may soon allow modern irrigation-derived recharge to reach the depths of these wells. These findings have broad relevance to other global aquifers that have experienced fundamental changes from predominantly natural to anthropogenic recharge over the past century, driven by large-scale agricultural development (46, 47). We therefore suggest that a combined inert-gas and inorganic carbon system modeling approach represents the best means of improving estimates of groundwater residence times in groundwater systems susceptible to DIC perturbations by recent agricultural activity, without the added uncertainties introduced by more complex models. In addition, we emphasize that the agricultural regions most prone to masking of natural dilution of C0 14 by irrigation seepage are often those most heavily reliant on groundwater and, therefore, most vulnerable to declines in groundwater quality and quantity.


Groundwater measurements

The wells selected for this study were sampled either as a part of this work or in previous U.S. Geological Survey (USGS) groundwater measurement programs (54). Groundwater samples were collected following standard USGS protocols and analyzed at USGS and other laboratories using established analytical methods (55) as part of the USGS National Water-Quality Assessment Enhanced Trends Network program (54). The carbon isotopic composition of DIC was analyzed at the Woods Hole Oceanographic Institution National Ocean Sciences Accelerator Mass Spectrometry laboratory (56); major and trace elements were measured via inductively coupled plasma mass spectrometry/atomic emission spectroscopy, atomic absorption spectroscopy, and ion chromatography at the USGS National Water Quality Laboratory (57, 58); water isotopes at the USGS Reston Stable Isotope Laboratory (59); tritium at the USGS Menlo Park Tritium Laboratory (60); and noble gasses at USGS Denver Noble Gas Laboratory (61). Measurements of 39Ar in three deep wells were made at Pacific Northwest National Laboratory via ultra-low background proportional counting (62). Radioactive Kr isotopes (85Kr and 81Kr) were measured in four wells spanning multiple depth intervals at the same site (wells 180-1 to 180-4) at Argonne National Laboratory via Atom Trap Trace Analysis (13). DIC concentrations and speciation were determined from measured pH and alkalinity using inorganic carbon system equilibrium constants presented in appendix S2.

Lumped parameter modeling

Residence time distributions and associated tracer simulations were carried out in MATLAB using multiple age distribution models presented in TracerLPM (45). Atmospheric histories of 14C activity (21, 28), 3H activity (45), and 85Kr activity (41) were prescribed, and atmospheric 39Ar activity was assumed to be constant over the 3-ka time period considered in this study. For each of the five wells analyzed for inert age tracers (39Ar and 85Kr) in the case study (appendix S1), the PEM model was used to estimate residence time distributions. Perforated depth intervals (table S1) were used to constrain the PEM model, reducing the residence time distribution to a function only of mean residence time (τ). For each well, τ was first estimated by maximizing agreement between simulated and measured inert age tracer activities (39Ar or 85K). Next, using intertracer constrained PEM model was used to simulate 14C activity and Γ was optimized to minimize disagreement between measured and simulated 14C activities. Because all deep wells were 3H free (indicating predominance of pre-1950s water), it was assumed that DIC was invariable across the age mixture for each well. Last, Γ values were compared with τ and pH, indicating a consistent shift toward lower Γ with increasing τ and pH. Appendix S1 provides additional details on lumped parameter modeling.

Carbonate dissolution model

The idealized inorganic carbon model used in this study solves for 14C activity and δ13C of DIC and for concentrations of CO2(aq), CO3, and HCO3, pH, and Me2+ (where Me2+ = Mg2+ + Ca2+) as a function of system openness, f, based on prescribed values of soil pCO2, soil and carbonate mineral 14C and δ13C, and temperature. CO2 solubility, reaction constants, and carbon equilibrium isotopic fractionation factors were prescribed based on literature values (22). The model solves for system parameters at three conceptual stages: (stage 1) dissolution of CO2 before carbonate dissolution occurs, (stage 2) progressive carbonate dissolution under open-system conditions until the point of isolation of groundwater from overlying soil CO2, and (stage 3) subsequent carbonate dissolution in the saturated zone under closed-system conditions. The model assumes equilibrium between DIC species at each stage and isotopic equilibrium between soil CO2 and DIC for stages 1 and 2. Appendix S2 provides additional details on the idealized carbonate dissolution model, including example simulations under varying conditions.


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We thank F. Pavia for helpful discussions about carbonate dissolution and acknowledge K. Belitz for designing the national-scale approach for the ETN program. We also thank B. Little and the City of Fresno for allowing the instrumentation and collection of data at City wells. We are grateful to three anonymous reviewers for helpful feedback. Funding: This work was conducted as a part of the USGS National Water Quality Assessment Program (NAWQA) Enhanced Trends Project ( Measurements at Argonne National Laboratory were supported by Department of Energy, Office of Science under contract DE-AC02-06CH11357. Measurements at Pacific Northwest National Laboratory were part of the Ultra-Sensitive Nuclear Measurements Initiative conducted under the Laboratory Directed Research and Development Program. PNNL is operated by Battelle for the U.S. Department of Energy under Contract DE-AC05-76RL01830. This work was also partially supported by NSF award OCE-1923915 (to A.M.S. and P.H.B. at WHOI). Author contributions: J.T.K. and A.M.S. designed the study and collected groundwater samples. A.M.S., J.T.K., D.V.B., P.H.B., K.E.D., and B.J. interpreted the data. A.M.S. carried out data analysis and modeling with input from J.T.K. and B.J. E.K.M. and C.E.A. measured 39Ar. J.C.Z. and P.M. measured 85Kr and 81Kr. A.M.S. wrote the manuscript with edits from all authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. All measurements and model output in this study are included as a Supplementary Dataset and can be publicly accessed online through the USGS National Water Information System (NWIS; using the well station identification numbers provided in Supplementary Dataset. Additional data related to this paper may be requested from the authors.

Stay Connected to Science Advances

Navigate This Article