2.1. Groundwater and Soil CO₂ Sampling

For demonstrating the performances of the developed techniques for ¹⁴C and δ¹³C analysis in DIC and soil CO₂, we carried out sampling and analysis of groundwater and soil CO₂ from the campus of Physical Research Laboratory (PRL) (Figure 1). Groundwater was collected from a borewell sucking water from a depth of ∼90 m from the surface. Water samples were collected in 60 mL polypropiline narrow mouth bottles (Tarson, Kolkata, India) directly from a pipe sucking the groundwater before falling on the reservoir from which water is distributed. For soil CO₂ collection, a homemade customized vacuum system with pumping arrangement was utilized as shown in Figure 2. A stainless-steel pipe (∼2 cm diameter) was introduced into the soil to the depth from which the soil air was to be collected. The soil material inside the pipe was removed and a new empty pipe was introduced into the hole. The empty spaces on the sides of the pipe were packed with soil materials to avoid contact with atmospheric air. The glass flask along with the pipeline was evacuated using a rotary vacuum pump. The flask was then filled to 1 atmosphere with the soil air from the desired depth, the pressure inside the flask was monitored using a pressure gauge. More details about the soil CO₂ extraction system are provided elsewhere (Mohanty et al. in preparation).

Figure 1 Location for the groundwater sampling site, the campus of the Physical Research Laboratory (PRL). Locations of the sampling sites whose groundwater samples were dated previously by Agarwal et al. (Reference Agarwal, Gupta, Deshpande and Yadava2006) are also shown (the site number were taken from their published article). Samples with radiocarbon ages neither modern nor beyond radiocarbon dating limit are listed here.

Figure 2 Soil air sampling setup. The soil air samples are collected from a specific soil depth introducing the pipe to that depth and filling the pre-evacuated high vacuum glass flask.

2.2. CO₂ Extraction from DIC of Groundwater

A schematic diagram for the extraction of DIC from water is shown in Figure 3. The extraction system was modified from some existing techniques (e.g., Carmi et al. Reference Carmi, Kronfeld, Yechieli, Yakir, Stiller and Boaretto2007; Macchia et al. Reference Macchia, D’Elia, Quarta, Gaballo, Braione and Maruccio2013; Ge et al. Reference Ge, Wang, Zhang, Luo and Xue2016; Takahashi et al. Reference Takahashi, Handa and Minami2021). About 10 mL commercial orthophosphoric acid (85%) was taken in a flask of ∼1 L capacity and was evacuated to ∼1 × 10^–2 mbar. Using a syringe, 20–200 mL water was introduced into the flask. The amount of water for reaction depends on the concentration of DIC. The flask was heated to ∼60ºC using a hot plate and a magnetic stirrer was used for fast reaction and release of CO₂. The reaction was allowed to proceed for 15–20 minutes. The liberated CO₂ was passed through a series of three U-tubes with first one immersed in alcohol-liquid nitrogen slush kept at –80ºC and the subsequent two in liquid nitrogen (LN₂) at –196ºC. The U-tube immersed in the slush, trapped most of the moisture while CO₂ passed through it and freezes in the two LN₂ traps. One of the LN₂ traps was a spiral tube to increase the CO₂ freezing efficiency (Figure 3). The extracted CO₂ was cleaned by freeze and thaw technique a few times while transferring from one trap to the other to remove remaining traces of water and other impurities. The amount of CO₂ liberated was measured by expanding it in a calibrated volume connected to a pressure gauge. Finally, pure CO₂ was collected in the pre-evacuated high vacuum glass ampule to load into the graphite preparation system. The complete extraction of CO₂ from a water sample takes almost an hour. The vacuum system takes an additional half an hour for cleaning by pumping between two samples.

Figure 3 Glass vacuum line used for extraction of carbon dioxide from dissolved inorganic carbon and soil CO₂. Water samples are injected into the sample flask, reacted with orthophosphoric acid under vacuum, passed through tubes immersed in slush (−80°C) and liquid nitrogen (−196°C) and finally collected in the sample bottle.

2.3. Estimation of CO₂ Yield from DIC

The liberated CO₂ from the water samples was allowed to expand in a pre-calibrated volume attached to a pressure gauge (Setra DATUM 2000^TM) to measure the CO₂ content and hence DIC concentration in the sample. The gauge and attached volume were calibrated using pure Na₂CO₃ powder. Precisely weighed Na₂CO₃ was reacted with orthophosphoric acid in the flask at 60ºC for half an hour and the liberated CO₂ was extracted following the procedure discussed above. The amount of Na₂CO₃ reacted varied from 2 to 53 mg. The liberated CO₂ was expanded in the volume attached to the gauge and a calibration curve was prepared by plotting the amount of CO₂ in micromole against the gauge reading (Figure 4). The amount of liberated CO₂ in micromole was calculated from the actual weight of the Na₂CO₃ powder used in the reaction. The amount of CO₂ liberated from any sample with unknown concentration can be estimated directly from the gauge reading using the equation shown in Figure 4.

Figure 4 Calibration plot used for estimating the amount of CO₂ liberated from any sample.

2.4. Extraction of CO₂ from Soil Air

CO₂ from soil air for δ¹³C and ¹⁴C measurements was extracted using the same glass vacuum line used for extracting DIC (Figure 3). Instead of using water flask, the glass flask containing soil air was connected and CO₂ was quantitatively extracted from the soil air using cryogenic technique. From the total amount of CO₂ and the volume of the flasks, the concentration of CO₂ was measured with a precision of ∼1% (Mohanty et al. in preparation).

2.5. Preparation of Standards

In absence of water standards for radiocarbon dating, we used carbonate standards for estimating the accuracy of the radiocarbon measurements for the present technique. The International Atomic Energy Agency’s two carbonate standards viz. Carrara Marble (IAEA-C1) and Travertine (IAEA-C2) were used. The Carrara marble serves as background and Travertine has a radiocarbon content of 41.14 ± 0.03 pMC (Rozanski Reference Rozanski1991). The carbonate samples were finely powdered with agate-mortar and dried samples were loaded in the reaction flask followed by evacuation. About 20–30 mL of Milli Q water was introduced into the flask using a syringe through a rubber septum under vacuum condition. The carbonate powder was mixed with water to make the extraction procedure similar to that of water samples. About 5 mL orthophosphoric acid was then injected into the flask using a syringe. The CO₂ evolved was extracted following the same procedure as discussed for the water samples. This gives the estimate of the background and accuracy of the procedure followed for radiocarbon measurements in DIC samples.

2.6. Graphite Preparation and ¹⁴C Dating Using AMS

The purified CO₂ samples were converted into graphite using an in-house graphite preparation line (Figure 5). About 50 to 100 micromole CO₂ was transferred to a vacuum system in which pre-conditioned iron (∼5 mg) and zinc (∼25 mg) powders were loaded in the two arms of an L-shaped interconnected vacuum sealed quartz tubes (Jull et al. Reference Jull, Donahue, Hatheway, Linick and Toolin1986). Iron and zinc powders were heated to appropriate temperatures for CO₂ reduction to elemental carbon (graphite) which deposit on the surface of iron powder. Graphite precipitation using this technique takes 3 to 5 hours. The graphite-coated iron was then pressed into pellet and loaded into the source of the 1 MV accelerator mass spectrometer for ¹⁴C measurements at Physical Research Laboratory (PRL) Ahmedabad, India (Bhushan et al. Reference Bhushan, Yadava, Shah and Raj2019). Using ¹⁴C/¹²C or ¹⁴C/¹³C for samples, ages were estimated (Donahue et al., Reference Donahue, Linick and Jull1990). Along with the samples, measurements of international standards such as Oxalic Acid I, Oxalic Acid II and background materials (anthracite) were made to estimate the background subtracted activity of samples and measure their radiocarbon ages. The ages were corrected for the isotopic fractionation during sample preparation and within the AMS (Donahue et al. Reference Donahue, Linick and Jull1990).

Figure 5 Graphite preparation system. Fe: iron powder, Zn: zinc powder, PG: pressure gauge, VG: vacuum gauge, P: pump, SB: sample bottle to load CO₂ into the system.

2.7. Radiocarbon Dating in Groundwater

Conventionally radiocarbon age (using Libby half-life of 5568 yr) of any sample is calculated using the basic decay equation given by

(1) $$t = - {{5568} \over {\ln (2)}}\ln \left( {{A \over {{A_0}}}} \right)$$

Where A is the activity measured in a sample and A ₀ is the initial activity. These activities are corrected for isotopic fractionation also called normalized activities (Stuiver and Polach Reference Stuiver and Polach1977) “Closed system” is a basic requirement for radiocarbon age determination i.e., the system should remain closed to subsequent gains or losses of the carbon except through radioactive decay. However, most of the natural samples particularly subsurface water samples experience many physical and chemical processes and mixing with many different water channels during transport and within the aquifers altering the radiocarbon contents (Plummer and Glynn, Reference Plummer and Glynn2013 and references therein). Decomposition of old organic matter or dissolution of carbonates may cause incorporation of DIC into groundwater in addition soil CO₂. These two components might have different ¹⁴C content compared to the original groundwater. In such cases actual ¹⁴C content present in a groundwater sample can be different from that attained just by decay of the original ¹⁴C content of the sample. Though, there is no standard way to correct for the incorporation of carbon from the surroundings, some schemes for radiocarbon age corrections in groundwater samples are discussed below.

If DIC in groundwater is solely derived from atmospheric CO₂ and/or the modern soil CO₂ and no mixing with other flow paths or incorporation of carbon from other sources take place, the determined age will be close to its exact value. It is to be noted that the determined age should be corrected for age of the soil CO₂ particularly when the ages of CO₂ in the unsaturated zone are significantly old or if a significant fraction of carbon is originated from the post 1950 CE atmosphere (bomb carbon). For example, ¹⁴C content of the soil CO₂ at Ti Tree (Australia) at a depth of 10 m is ∼50 pMC though it is close to 100 pMC near the surface (Wood et al. Reference Wood, Cook and Harrington2015). Low ¹⁴C contents in soil CO₂ in deeper soil layers were also observed at Saskatoon (Canada), North Dakota and Navada (USA) (Bacon and Keller Reference Bacon and Keller1998; Thorstenson et al. Reference Thorstenson, Weeks, Haas and Fisher2016; Haas et al. Reference Haas, Fisher, Thorstenson and Weeks2016). On the other hand, presence of bomb carbon in unsaturated zone CO₂ or in organic matter (e.g., Blume et al. Reference Blume, Guitard, Crann, Orekhov, Amos, Clark, Lapen, Blowes, Ptacek, Craiovan and Sunohara2022; Laskar et al. Reference Laskar, Yadava and Ramesh2012) may lead to underestimation of groundwater ages if not corrected. In most of the studies related to the groundwater dating, ¹⁴C content in CO₂ in the unsaturated zone is assumed to be modern which may lead to significant discrepancies between the actual and measured ages. The ¹⁴C contents are also expressed as percent modern carbon pMC = (A/A₀)100% or in the form of delta notation such as D¹⁴C = (A/A₀–1)1000‰ and or fraction modern F¹⁴C = A/A₀. A sample with pMC greater than 100% or positive D¹⁴C or F¹⁴C greater than 1 indicate presence of bomb carbon.

Incorporation of carbon with different ¹⁴C signature in groundwater DIC pool in some cases can be corrected using stable carbon isotope ratio (δ¹³C) in DIC as follows. Average δ¹³C values of C3 and C4 vegetation are ∼ –27‰ and –13‰, respectively (Deines Reference Deines1980; Khon Reference Kohn2010). The organic carbon present in soil under a particular vegetation type has similar or slightly enriched δ¹³C values (Laskar et al. Reference Laskar, Sharma, Ramesh, Jani and Yadava2010, Reference Laskar, Yadava, Sharma and Ramesh2013, Reference Laskar, Yadava and Ramesh2016). The enrichment is due to incorporation of ¹³C depleted CO₂ in the atmosphere during the industrial era (Suess Effect) (McCarrol and Loader Reference McCarroll and Loader2004; Paul et al. Reference Paul, Balesdent and Hatté2020) and preferential removal of lighter carbon during decomposition by microbes (Laskar et al. Reference Laskar, Yadava and Ramesh2016). Soil CO₂ has a δ¹³C value controlled by the decomposition of organic matter, root respiration and diffusive fractionation. As soil CO₂ mainly controls the DIC contents in the groundwater, its δ¹³C value can be used to estimate its fraction in DIC taking into account the appropriate fractionation factor (Szaran Reference Szaran1997). Therefore, by measuring δ¹³C values in soil CO₂ in recharge region and in DIC, contribution from other sources (bedrock) can be estimated and corrected for the ¹⁴C activity in the sample. If ¹⁴C dilution is caused by dissolution of carbonate of marine origin, the correction can be relatively easily made using δ¹³C values as these carbonates are free from ¹⁴C and their δ¹³C values are 0 ± 2‰ (Ripperdan Reference Ripperdan2001) which is very different from soil CO₂ (∼ –20‰). On the other hand, if dissolution takes place from soil carbonates formed by evaporation or lacustrine carbonates it may be difficult to estimate the contributions simply using δ¹³C of DIC due to weak constraints and possible overlapping in the range of their δ¹³C values with that dissolved from soil CO₂. In addition, groundwater does not flow like the flow in a tube (piston flow), but encounter dispersion, advection and mixing making it a more complex system. Therefore, appropriate models are required to handle those mixing issues as discussed in Section 2.8.

A simple approach for correction of ¹⁴C age in groundwater due to contamination from carbonate dissolution and decomposition of old organic matter is as follows. Suppose a groundwater sample is diluted due to incorporation of carbon from bedrocks and q is the dilution factor (see Equation 3 below), then the ¹⁴C activity in the recharge area would be qA_o, and the corresponding radiocarbon age of such a groundwater, as discussed elsewhere (e.g., Clark and Fritz Reference Clark and Fritz1997; Aggarwal et al. Reference Agarwal, Gupta, Deshpande and Yadava2006; Cartwright Reference Cartwright, Currell, Cendón and Meredith2020) would be

(2) $$t = - {{5568} \over {\ln (2)}}\ln \left( {{A \over {q{A_0}}}} \right)$$

Unfortunately, determination of the dilution factor q is challenging in most of the cases. In case of dilution by marine carbonate under closed condition, the dilution factor q can be calculated using the following two-component isotopic mass balance equation:

(3) $${\rm{{\delta ^{13}}{C_{DIC}} = {\rm{ }}{\delta ^{13}}{C_{carb}} \times q + {\delta ^{13}}{C_{recharge}} \times \left( {1-q} \right)}}$$

where δ¹³C_DIC, δ¹³C_recharge and δ¹³C_carb are the δ¹³C values observed in the DIC in the groundwater, DIC of recharging water and carbonate rocks in the aquifers, respectively. δ¹³C_recharge can be obtained from the δ¹³C value of soil CO₂ as δ¹³C in bicarbonate (most of the DIC at neutral pH) is ∼7.4‰ higher than δ¹³C value of the soil CO₂ at 30ºC (Vogel et al. Reference Vogel, Grootes and Mook1970; Mook et al. Reference Mook, Bommerson and Staverman1974; Deines et al. Reference Deines, Langmuir and Harmon1974; Szaran Reference Szaran1997). In the above isotopic mixing calculation (Equation 3), contribution from atmospheric CO₂ is neglected because of the fact that most of the dissolution of CO₂ in groundwater takes place in the unsaturated zone due to high concentration of CO₂ in the latter. Soil CO₂ from decomposition/oxidation of old organic matter if significant can be corrected using the ¹⁴C in soil CO₂ in the recharge areas. The values of q can also be calculated from major ion geochemistry, graphical techniques and isotopic mass balance as discussed by many previous researchers (e.g., Vogel and Ehhalt Reference Vogel and Ehhalt1963; Vogel et al. Reference Vogel, Grootes and Mook1970; Pearson and Hanshaw Reference Pearson and Hanshaw1970; Mook Reference Mook1972; Tamers Reference Tamers1975; Fontes and Garnier Reference Fontes and Garnier1979; Clark and Fritz Reference Clark and Fritz1997; Geyh Reference Geyh2000; Agarwal et al. Reference Agarwal, Gupta, Deshpande and Yadava2006; Coetsiers and Walraevens Reference Coetsiers and Walraevens2009; Blaser et al. Reference Blaser, Coetsiers, Aeschbach-Hertig, Kipfer, Camp, Loosli and Walraevens2010; Clark Reference Clark2015; Han and Plummer Reference Han and Plummer2016; Cartwright et al. Reference Cartwright, Currell, Cendón and Meredith2020).

There can be more potential sources of DIC with different ¹⁴C content compared to the soil CO₂ such as soil carbonate, carbonate nodules and addition of bomb carbon in groundwater which can alter the q value. Addition of such carbon in DIC would result in over or underestimation of groundwater residence time if not corrected. Estimation of the input of ¹⁴C-free DIC or any input into the DIC pool of an aquifer with different initial ¹⁴C content than the aquifer water is a major challenge in interpreting radiocarbon data in groundwater (Vogel and Ehhalt Reference Vogel and Ehhalt1963; Fontes Reference Fontes1992; Aravena et al. Reference Aravena, Wassenaar and Plummer1995; Kalin Reference Kalin, Cook and Herczeg2000; Coetsiers and Walraevens Reference Coetsiers and Walraevens2009; Han et al. Reference Han, Plummer and Aggarwal2012; Plummer and Glynn Reference Plummer and Glynn2013; Han and Plummer Reference Han and Plummer2016; Nydahl et al. Reference Nydahl, Wallin, Laudon and Weyhenmeyer2020; Cartwright et al. Reference Cartwright, Currell, Cendón and Meredith2020). As shown in Equation 3, for a two-component mixing, values of q may be calculated from the δ¹³C values of soil CO₂, DIC in groundwater and carbonate rocks (Pearson and Hanshaw Reference Pearson and Hanshaw1970; Tamers Reference Tamers1975; Clark and Fritz Reference Clark and Fritz1997; Geyh Reference Geyh2000; Coetsiers and Walraevens Reference Coetsiers and Walraevens2009; Clark Reference Clark2015; Han and Plummer Reference Han and Plummer2016).

The ages estimated using decay equations (1 and 2) need to be converted into calendar ages using radiocarbon calibrations (e.g., Reimer et al. Reference Reimer, Austin, Bard, Bayliss, Blackwell and Bronk Ramsey2020; Hogg et al. Reference Hogg, Heaton, Hua, Palmer and Southon2020). Calibration is necessary to compare the obtained ages with other paleoclimate datasets. Other factors (e.g., hard-water effect) may cause erroneous estimation of radiocarbon ages and need to be carefully accounted before calibrations. The calibrated ages of groundwater can be considered as the residence times in some limited cases if all the water in an aquifer is recharged at the same time without much dispersion and advection during its transport. Hardly such simple systems exist as groundwater flows along paths of varying lengths and undergoes hydrodynamic dispersion and advection causing ¹⁴C content modification. In other words, a groundwater sample may contain many fractions of water from different origins and have a range of residence times rather than being of a single age (Maloszewski and Zuber Reference Maloszewski and Zuber1982; Maloszewski Reference Maloszewski2000; Cox Reference Cox2003; Suckow Reference Suckow2014).

2.8. Estimation of Residence Time of Groundwater Using Radiocarbon

There is no standard way to account all the complexities discussed above. Researchers address these complexities using different models (Seltzer et al. Reference Seltzer2021; Starn et al. Reference Starn, Kauffman, Carlson, Reddy and Fienen2021; Fareze et al. Reference Fareze, Rajaobelison, Ramaroson, Razafitsalama and Rakotomalala2022). One approach is to use Lumped Parameter Models (LPM) (Maloszewski and Zuber Reference Maloszewski and Zuber1982; Maloszewski and Zuber Reference Maloszewski and Zuber1992; Maloszewski Reference Maloszewski2000; McGuire and McDonnell Reference McGuire and McDonnell2006; Jurgens et al. Reference Jurgens, Böhlke and Eberts2012; Cartwright et al. Reference Cartwright, Currell, Cendón and Meredith2020) to estimate the groundwater residence times. LPM are mathematical models that deal with the transport of water from recharge zones to discharge regions based on aquifer geometry and flow configurations. The model’s approach involves convolution of tracer input and decay functions with a weight factor and the obtained output is matched to a known time series (Maloszewski and Zuber Reference Maloszewski and Zuber1982; Maloszewski et al. Reference Maloszewski, Stichler and Zuber2004). The models allow the variability of the ¹⁴C input, define a recharge zone, time taken by the water to penetrate the soil zone and the ratios of advection to dispersion. In other words, the model assumes more realistic and flexible flow path geometries, and the effects of dispersion may be taken into account to a large extent.

In the LPM, at steady state, the ¹⁴C content of groundwater ¹⁴C (t) at a time t in any outlet position in the aquifer can be calculated from the corrected ¹⁴C content of DIC in the recharging water (¹⁴C_i) using exit age distribution function and the radioactive decay via the following convolution integral (Maloszewski and Zuber, Reference Maloszewski and Zuber1982):

(4) $${}^{14}{C_t} = \int\limits_0^\infty {{C_i}(t - \tau )} g(\tau )\exp ( - \lambda \tau )d\tau $$

Where τ is the mean residence time, t- τ is the time when the water recharged, λ is the decay constant and g(τ) is the response function describing the distribution of flow paths and residence times in the aquifer. Depending on the aquifer geometries, different LPMs can be applied (e.g., Maloszewski and Zuber, Reference Maloszewski and Zuber1982; Maloszewski and Zuber Reference Maloszewski and Zuber1992; Amin and Campana, Reference Amin and Campana1996; Cook and Bohlke Reference Cook, Bohlke, Cook and Herczeg2000; Maloszewski Reference Maloszewski2000; Jurgens et al. Reference Jurgens, Böhlke and Eberts2012).

Here we used a tracer LPM program developed by Jurgens et al. (Reference Jurgens, Böhlke and Eberts2012) to estimate the residence times of groundwater extracted from the PRL Campus. We also applied the LPM program on previously measured ¹⁴C ages in groundwater from different parts of Gujarat, Western India. We applied three different LPMs viz. piston flow model (PFM), dispersion model (DM) and binary mixing model (BMM) to show the variations in residence times under different dispersion and mixing parameters. Schematics of the PFM and DM are shown in Figure 6.

Figure 6 Schematics of (A) piston flow model (after Bethke and Johnson, Reference Bethke and Johnson2008) and (B) dispersion model (after Maloszewski and Zuber Reference Maloszewski and Zuber1982).

The PFM assumes that a tracer (¹⁴C) travels from the recharge area to the outlet position (e.g., a well) without mixing or hydrodynamic dispersion. For PFM, radiocarbon output for a constant input (¹⁴C_i), the solution of Equation (4) can be expressed as (Zuber and Maloszewski Reference Maloszewski2000):

(5) $${C_{\left( t \right)}} = {C_i}\left( {t - \tau } \right){\rm{exp}}\left( { - {\rm{\lambda \tau }}} \right)$$

The PFM is applicable to hydrogeologic settings with low dispersion, high linear velocity or short flow path from recharge to discharge region. ¹⁴C tracer can be assumed to follow piston-flow behaviour in shallow, short screened monitoring wells in unconfined aquifers or in confined aquifers with a relatively small recharge area.

The DM is based on a solution of a one dimensional advection dispersion equation for a semi-infinite medium with an instantaneous injection and detection of the tracer in the fluid flux. A solution for the Equation (4) for a constant input (¹⁴C_i) with advection and dispersion can be expressed as (Zuber and Maloszewski Reference Maloszewski2000):

(6) $${\rm{^{14}{C_t}{ = ^{14}}{C_i}\ exp\left[ {{{\left( {2{P_D}} \right)}^{-1}} \times \{ 1 - {{\left( {1 + 4\lambda {P_D}\tau } \right)}^{1/2}}\} } \right]}}$$

Where P_D is the dispersion parameter which is the reciprocal of Peclet number describing the relative importance of dispersion and advection. The Peclet number depends on the velocity of flow field and characteristics length of the system (Rapp Reference Rapp2017). The DM can give an approximate description of age distributions in samples from multiple aquifer configurations.

Binary mixing model (BMM) describes two-component mixtures in which each component can be obtained by a model. For example, a “BMM-PFM-DM” model describes a binary mixture in which one component of the mixture is modelled by using PFM (¹⁴C₁) and the other component by DM (¹⁴C₂) and the corresponding fractions are f1 and f2(=1-f1). The output in a BMM is given by:

(7) $${\rm{^{14}{C_{out}} = {f_1}^{14}{C_1} + {\left( {1 - {f_1}} \right)^{14}}{C_2}}}$$

Binary mixing models can be appropriate for wells screened across multiple aquifer units and aquifers with short-circuit pathways that result in age mixtures of significantly different mean ages.