Skip to main content


Front. Water, 28 September 2022
Sec. Water and Hydrocomplexity
This article is part of the Research Topic Advances and Emerging Methods in Tracer Hydrogeology View all 5 articles

Explicit simulation of environmental gas tracers with integrated surface and subsurface hydrological models

  • 1Centre for Hydrogeology and Geothermics (CHYN), Faculty of Science, University of Neuchâtel, Neuchâtel, Switzerland
  • 2Climate and Environmental Physics and Oeschger Center for Climate Change Research, Faculty of Science, University of Bern, Bern, Switzerland
  • 3Hydrogeology, Department of Environmental Sciences, University of Basel, Basel, Switzerland
  • 4Eawag, Swiss Federal Institute of Aquatic Science and Technology, Dübendorf, Switzerland

Environmental gas tracers allow inferring groundwater travel times and mixing ratios. Their concentrations are commonly interpreted with simplified and indirect approaches that are conceptually at odds with the high degree of complexity found in natural systems. However, the information content of the tracers can potentially be fully explored through the explicit simulation of an advection-dispersion transport equation, for example using integrated surface-subsurface hydrological models (ISSHMs). These integrated models can be used to explicitly simulate environmental tracers in complex environments. ISSHMs are usually variably saturated flow models. However, these models do not explicitly simulate gas partitioning with the aqueous phase, restricting explicit simulation of gas tracers to fully saturated conditions or to tracers with very low solubilities. We propose a mathematical formulation for the production of environmental gas tracers that are emanated in the subsurface. The production is scaled according to gas/water partitioning and water saturation, which is already computed by the model. Therefore, ISSHMs can now be used to their full potential to explicitly simulate tracer concentrations under variably saturated and dynamic conditions. The new formulation has been successfully verified against reference simulations provided with a multi-phase flow and transport model. In addition, explicit simulation of 222Rn and 37Ar groundwater concentrations in a synthetic alluvial river-groundwater system was demonstrated, for the first time, with an ISSHM.


Hydrologic environmental tracers are widely distributed in the near-surface environment of the Earth, such that their abundance facilitates their consideration for studying a range of subsurface flow processes (Cook and Herczeg, 2012; Partington et al., 2020; Dwivedi et al., 2022). Many environmental tracers (e.g., CFCs, N2, O2, Ar, Ne, Kr, He, Xe, CH4, CO2, 222Rn, and SF6) are gases that dissolve into water. These environmental gas tracers can be either natural (e.g., 39Ar, 85Kr) or anthropogenic (e.g., CFCs, SF6) and their concentrations can either remain stable, increase (production) and/or decrease (decay) in the subsurface as a function of time.

The residence time of water in the subsurface is the sum of the transit times through the unsaturated and saturated zones. Residence times can be obtained from the decay of radioactive tracers and from tracers that are produced (accumulating tracers) in the subsurface. Therefore, their measured concentrations can be used to determine pathways and timescales (i.e., ages) of groundwater (Kazemi et al., 2006).

Tracer concentrations are usually considered for estimating apparent ages that are then used in the calibration of hydrological models (Van Huijgevoort et al., 2016; Schilling et al., 2017a, 2019, 2022). However, these simplified and indirect approaches are conceptually at odds with the high degree of complexity found in natural systems (McCallum et al., 2015, 2017; Peel et al., 2022). The explicit simulation of tracer mass-transport through the solution of an advection-dispersion equation is a promising pathway that can potentially increase the information contained in the tracers (Goode, 1996; Turnadge and Smerdon, 2014). A zero-order source term (or boundary condition) can be used to consider environmental tracers that are produced in the subsurface.

Gases that are emanated in the unsaturated media will partition between the aqueous and gaseous phases as a function of their solubility in water and water saturation. Therefore, a rigorous way to model accumulating tracers in a numerical flow model requires the consideration of both gas and liquid phases to allow exchanges between the two. For gases with high solubilities (e.g., 222Rn), deep unsaturated zones, and for tracers where the concentration is a function of depth (e.g., Ar isotopes), the consideration of a gas phase in the model is a prerequisite for robust tracer simulation.

The consideration of both gas and liquid phases requires the use of multiphase fluid flow and transport models that solve the advection-dispersion equation. HYDRUS, NUFT, and Getflows (Nitao, 1998; Tosaka et al., 2000; Simunek et al., 2016) are examples of multiphase flow and transport models. However, these mechanistic models are computationally demanding and are therefore usually preferred for controlled experiments and to model mass-transport in unsaturated soils (Guillon et al., 2016).

In groundwater hydrology, the simulation of the unsaturated zone along with the groundwater flow domain and the surface water flow domain can be done using integrated surface and subsurface hydrologic models (ISSHM) (Sebben et al., 2013; Paniconi and Putti, 2015). These integrated models are particularly well-suited to study complex hydrodynamic behaviors such as river-groundwater interactions that characterize shallow subsurface environments (Brunner et al., 2017). HydroGeoSphere (Brunner and Simmons, 2012; Aquanty Inc., 2022), ParFlow (Kollet and Maxwell, 2008), or OpenGeoSys (Kolditz et al., 2012) are examples of such integrated models. Some of these models allow the explicit simulation of tracers in the saturated zone of the aquifer. For instance, explicit modeling or a radiogenic tracer such as 222Rn has been recently proposed with HydroGeoSphere (Gilfedder et al., 2019). However, most of the integrated models are single-phase flow, preventing the consideration of a gas phase in the unsaturated zone. Thus, the zero-order source, which allows a tracer to be produced in the subsurface, does not consider gas/water partitioning, potentially leading to an overly large concentration in the unsaturated zone, as the produced tracer accumulates in a smaller volume of water compared to what would be the case under full saturation.

To overcome this issue we propose a simplified mathematical formulation for the production of environmental gas tracers that are emanated in the subsurface, which allows the solute mass production rate for variably saturated conditions to be scaled according to the degree of saturation and the gas/water partitioning coefficient. The proposed formulation is designed for single-phase surface and subsurface flow and transport in integrated hydrological models that are capable of simulating variably saturated subsurface flow. In the present study this formulation was implemented in the ISSHM HydroGeoSphere (HGS) and is called via the “Zero-order source with partitioning” command (Aquanty Inc., 2022). The new formulation in HGS was validated against an explicit multiphase flow and transport reference case implemented in the software HYDRUS-1D (Simunek et al., 2016). The capability of the new formulation is illustrated for a river-groundwater flow system where 222Rn and 37Ar are explicitly simulated.

Mathematical formulation

An instant-equilibration mass balance

The present section introduces the derivation of the proposed instant-equilibration mass-balance to scale the production of dissolved gas tracers according to the degree of water saturation of the subsurface. An underlying assumption is that any produced amount of tracer will partition immediately into both the gas and aqueous phases. In a small test volume Vtest, we can write:

δQtracer,totprodδt=PsatVpore     (1)

where, δQtracer,tot prod [M] is the amount of tracer produced in a small time increment δt [T], Psat [M·L−3·T−1] is the total tracer production rate (assumed independent of the degree of water saturation), and Vpore [L3] is the pore volume of the matrix (porosity × Vtest).

If we assume that the amount of tracer δQtracerprod is partitioned into both gas and aqueous phases, we can write the following mass balance relationship:

δQtracerprod=δQtracer,waterprod+δQtracer, airprod=δCtracer,waterprodVwater                  + δCtracer,airprodVair    (2)

where, δQtracer,waterprod and δQtracer,airprod are the amounts of produced tracer [M] which go into the liquid and gas phases, respectively. Similarly, δCtracer,waterprod and δCtracer,airprod are changes in tracer concentrations from the zero-order source [M·L−3] in the given volumes of water Vwater [L3] and air Vair [L3].

Assuming a solubility law of the form Hcc=CwaterCair with instant equilibration, we can write:

δQtracerprod=Hcc δCtracer,airprodVwater+ δCtracer,airprodVair                  =δCtracer,airprod(HccVwater+Vair)    (3)

If we divide both sides by Vpore, we get:

1VporeδQtracerprod=δCtracer,airprodHccVwater+VairVpore                                    =δCtracer,airprod(HccSw+Sa)    (4)

where, Sw and Sa = 1−Sw are respectively, the degree of water and air saturations [–]. Therefore, by combining Equation (4) with Equation (1) we can write:

Psat= 1VporeδQtracerprodδt=δCtracer,airprodδt(HccSw+[1Sw])    (5)

Which can also be written as follows:

δCtracer,airprodδt= PsatHccSw+[1-Sw]    (6)

If we define the effective production rate in water Peff [M·L−3·T−1] as the amount of tracer input only into the water phase per unit time, we can write:

Peff=δCtracer,waterprodδt=Hcc δCtracer,airprodδt=HccPsatHccSw+[1-Sw]    (7)

This formulation (Equation 7) is equivalent to the original zero-order source as implemented in HGS (production per unit aquifer volume) in the limit Hcc → ∞ (i.e., no partitioning) and for Sw = 1 regardless of Hcc. This means that all produced species remain in the aqueous phase. In the case Hcc → 0, the tracer is immediately volatilized. Mean groundwater age, defined with the age-mass concept (Goode, 1996), can be estimated directly by setting Hcc = 1 with a solute mass production rate equal to the aquifer porosity and with a decay constant of zero (no decay of mean age).

Equation 7 assumes that every produced atom in the unsaturated zone partitions immediately into both gas and aqueous phases according to a fixed partitioning coefficient. Thus, the proposed formulation is equivalent to an equilibrium solute transport model where the mass of volatilized tracer, which in theory would interact with groundwater in case of variations in water saturation, is not considered. The phase-equilibrium partitioning coefficient Hcc, is considered constant, synonymous with isothermal conditions. Moreover, Hcc is a well-constrained parameter as it strictly depends on water temperature and the considered species. Interested readers should refer to Cook and Herczeg (2012) and Schubert et al. (2012) for default values of Hcc for 222Rn and 37Ar, respectively.

The hypothesis of saturation-independent emanation is supported by the fact that 222Rn (and likely 37Ar) emanation has been shown to reach its maximum value at moderate degrees of water saturation, typically between 0.2 and 0.3 (e.g., Sun and Furbish, 1995; Zhuo et al., 2006). Indeed, the presence of a thin moisture film on the surface of mineral grains will prevent produced atoms recoiled into the pore space from embedding into adjacent grains, as recoil ranges of atoms in water are orders of magnitude lower than those in air. Only in very dry conditions and low residual water contents in the porous medium can this hypothesis no longer be justified.


To verify the proposed instant-equilibration mass balance as newly implemented in HGS, simulations of 222Rn in an unsaturated medium using HGS are compared with explicit multi-phase simulations using HYDRUS-1D. HYDRUS-1D is based on a finite element model that solves Richard's equation for variably-saturated subsurface flow and Fickian-based advection-dispersion equations for heat and mass transport.

In the proposed experiment, HYDRUS-1D was used to solve the advection-dispersion solute transport equations with a linear equilibrium between the aqueous and the gaseous phases. Diffusion in both aqueous and gas phases were considered to be instantaneous (i.e., equilibrium transport). Moreover, we did not consider the effects of adsorption-desorption and we have assumed no immobile fraction of water. Thus, the advection-dispersion equation can be written in the following form:

θct+agt= x(θDw × cx)+ x(aDg × gx)                          - qcx-(μwθc+ μgag)+ (γwθ+ γga)    (8)

where, c and g are solute concentrations in the aqueous [M·L−3] and gaseous [M·L−3] phases, respectively; θ [L3·L−3] is the water content, ρ [M·L−3] is the soil bulk density, a [L3·L−3] is the air content, Dw [L2·T−1] is the dispersion coefficient for the aqueous phase, Dg [L2·T−1] is the diffusion coefficient for the gas phase, q [L·T−1] is the volumetric fluid flux density, μw and μg [T−1] are first-order rate (decay) constants for solutes in the aqueous and gas phases, respectively; γw and γg [M·L−3·T−1] are zero-order rate constants (or source) for the aqueous and gas phases, respectively.

Furthermore, equilibrium exchange between water and gas concentrations of the solute in the soil system was assumed. The adsorption isotherm was considered linear. The concentrations in the aqueous and in the gas phases (c and g) are related by a linear expression of the form:

g=kg × c (9)    (9)

where, kg [-] is an empirical constant equal to (KH·Ru·TA)−1, in which KH [M·T2·M−1·L−2] is Henry's law coefficient, Ru [M·L2·T−2·K−1·M−1] is the universal gas constant and TA [K] is the absolute temperature. The advection-dispersion equations were solved in the gas and in the aqueous phases, thus permitting the simulation of solute transport simultaneously in both the aqueous and gas phases.

222Rn production

222Rn is produced as an intermediary product in the 238U decay chain, and is released into groundwater as a result of the alpha decay of matrix-bound 226Ra (Cecil and Green, 2000). This radioactive noble gas with a half-life of 3.82 days is commonly used in groundwater hydrology as a tracer of surface water—groundwater interactions, as 222Rn activities in groundwater are often orders of magnitude higher than those in surface water (Cecil and Green, 2000). Upon infiltration of surface water with low 222Rn activities into the subsurface, 222Rn concentrations will gradually increase until they reach secular equilibrium after ~15–20 days (Hoehn and Von Gunten, 1989). Therefore, it has been extensively used to study very young groundwater dynamics and mixing [i.e., timescales of several hours to a few weeks (Cook et al., 2008; Bourke et al., 2014; Cartwright and Gilfedder, 2015; Peel et al., 2022)]. Assuming spatially homogenous 222Rn production, the relationship between production rate per unit pore volume γRn [M·L−3·T−1] and equilibrium 222Rn activity concentration in groundwater CRn,eqSat [M·L−3] is given by the following equation:

CRn,eqSat= γRnλRn (10)    (10)

where, λRn [T−1] is the 222Rn decay constant. The value for CRn,eqSat, which is site specific since it is closely related to aquifer geochemistry and mineral texture, has been considered homogenous for the present study, and set to 15 Bq·l−1 (e.g., Popp et al., 2021).

Model set-up

As a verification experiment, a one-dimensional soil column 5-m thick uniformly discretized with 0.05 m vertical intervals is considered for both HYDRUS-1D and HGS.

In both models, the soil column was assumed homogeneous and variably saturated flow was solved using the van Genuchten—Mualem model with the fitting parameters α [L−1] and β [–] set to 3.48 and 1.75, respectively, simulating a typical alluvial gravel (Dann et al., 2009; Schilling et al., 2020). The residual soil water content was set to 0.095 and the porosity to 0.41. A value of 1 m·d−1 was used for the saturated hydraulic conductivity. The tortuosity parameter in the conductivity function was set to 0.5. The longitudinal horizontal dispersivity was set to 0.1 m, while the transverse dispersivities were set to 0.01 m. A pore pressure of 0 was applied to nodes at the bottom (z = −5m) of the column, making the bottom of the soil model a constant water table elevation.

For HYDRUS-1D, the bulk density was set to 1.593 g·cm−3. In Equation 9, which controls the interaction between the liquid and the gas concentrations, the empirical constant kg was set to 2.155 (-), which is the inverse of the partitioning coefficient Hcc as defined in Equation 3. The first-order rate (decay) constant for 222Rn was set to 0.1814 d−1 for all phases. The γw, and γg zero-order rate constants were set to 2721 Bq·m−3·d−1, which has been derived from Equation 10, assuming a saturated 222Rn activity CRn,eqSat of 15 Bq·l−1.

Flow and transport through the soil column were simulated according to two scenarios: (i) drainage conditions and (ii) recharge conditions. The soil column was initially considered fully saturated with water. For scenario (i), the model was run without any infiltration fluxes at the top so that the soil column is progressively drained until a pseudo steady-state condition was reached after 365 days (Figure 1B). We refer to pseudo steady-state as transient simulations that reaches a dynamic equilibrium state, that is, a point at which there is minimal drift in the model states (similar to actual steady-state conditions). Starting from this new pseudo steady-state condition, for scenario (ii), a constant input flux of 0.1 m·d−1 was then considered at the top of the column, simulating an extremely high infiltration flux. The model was then run for another 365 days until a new pseudo steady-state condition was reached (Figure 1D). To evaluate the performance of the HGS simulations in transient conditions, snapshots of dissolved 222Rn activities were also taken five days after the beginning of each simulation scenario (Figures 1A,C). For a majority of locations, the simulated saturation profiles of HGS and HYDRUS-1D were nearly identical (Supplementary Figure A1), therefore allowing direct comparisons of the tracer concentrations simulated with the two different software.


Figure 1. Concentrations of 222Rn in water for a synthetic one-dimensional soil column model under (A) drainage driven transient (t = 5 days) flow, (B) drainage driven steady-state flow, (C) recharge driven transient (t = 5 days) flow, and (D) recharge driven steady-state flow.


Under drainage driven flow and transport (scenario (i); Figures 1A,B), a perfect match of the simulation of 222Rn after 365 days was reached with HYDRUS-1D (under consideration of explicit multi-phase flow) and HGS (under consideration of the new partitioned zero-order source production), demonstrating the consistency of the simplified approach under pseudo steady-state conditions. After five days of drainage (i.e., highly transient conditions), the HGS-based simulations with partitioned zero-order source production slightly overestimated the 222Rn concentrations compared to the multi-phase flow simulations with HYDRUS-1D. This reflects the ability of HYDRUS-1D to explicitly account for 222Rn in both the gas and liquid phases. When the saturation drops, dissolved 222Rn concentrations re-equilibrate with those in the gas phase; in other words, a drop in saturation will lead to the volatilization of some of the dissolved 222Rn to the gas phase and therefore a reduction of the concentration in the aqueous phase. In HGS, the gas phase is not simulated explicitly and the 222Rn concentration in water is therefore not directly affected by changes in saturation; only production is scaled, which leads to an overestimation of 222Rn activities in water after five days of transient simulation of drainage.

Under recharge driven flow (scenario (ii); Figures 1C,D), the opposite mechanism occurs, and 222Rn concentrations in water are slightly underestimated with HGS compared to HYDRUS-1D in highly transient conditions. As in the drainage scenario, the difference is more pronounced under transient conditions compared to the pseudo steady-state, when the simulated concentration profiles of both codes almost match. The reason for the mismatch lies in the fact that once the saturation increases within the soil column because of the constant input flux from the top, dissolved 222Rn concentrations start to increase as a result of soil matrix production, but also because of the dissolution of 222Rn from the gas phase into the water. Again, as in HGS, the gas phase is not explicitly simulated, dissolved 222Rn concentrations in the soil column tend to be underestimated compared to the explicit multi-phase simulations with HYDRUS-1D.

Under both drainage and recharge conditions, more realistic simulations were obtained with HGS when the zero-order source was partitioned with the water saturation (red lines in Figure 1). When the production of 222Rn was not scaled according to water saturation (i.e., unpartitioned zero-order source), the 222Rn concentrations in the unsaturated soil column were significantly overestimated by HGS because the same amount of total produced tracer mass accumulated in a comparatively smaller volume of water (orange lines in Figure 1). Although some discrepancies between HYDRUS-1D and HGS with the proposed partitioned zero-order source along a one-dimensional unsaturated soil column could be identified, the improvement achieved with the partitioned zero-order source implementation is significant, resulting in nearly matching pseudo steady-state conditions under both drainage and recharge conditions. The proposed partitioned zero-order source approach is thus a highly efficient and simplified, yet realistic, way to include simulations of environmental gas tracers in the unsaturated zone for applications involving both surface and the subsurface flow processes, especially where the application of computationally demanding and complex two-phase flow models is not desirable or practical.

Illustrative example

The simulation of environmental gas tracer concentrations with a single-phase ISSHM is illustrated below using a more realistic and complex model configuration. For this purpose, a synthetic alluvial plain model was used for the simulation of 222Rn and 37Ar mass-transport. The consideration of both 222Rn and 37Ar was motivated by their short half-lives, making these tracers particularly well-suited to track recently infiltrated surface water in shallow alluvial aquifers (Brunner et al., 2017; Schilling et al., 2017a).

37Ar production

37Ar is a rare radioactive isotope of argon with a half-life of 34.95 ± 0.08 days (Renne and Norman, 2001). A originally 37Ar free water parcel that enters the subsurface reaches 88% of the secular production-decay equilibrium after three half-lives of 37Ar, i.e., after approximately 110 days. In the shallow underground, 37Ar atoms are primarily produced by cosmogenic neutron activation of calcium (40Ca(n)37Ar), and to some extent spallation of potassium (39K(n,p2n)37Ar) (Fabryka-Martin, 1988; Loosli and Purtschert, 2005; Riedmann and Purtschert, 2011). The flux of fast cosmogenic neutrons is exponentially attenuated with depth on a length scale, which is proportional to the density of the medium without of its elemental composition (Gosse and Phillips, 2001). Other production reactions, e.g., induced by muons or by neutrons from (α, n) reactions can be neglected in the depth range considered here (Spannagel and Fireman, 1972; Heisinger et al., 2002; Sramek et al., 2017). The depth dependent production rate 37P(z)r in (atoms·yr−1·cm−3) has been parameterized through a production function equation (Guillon et al., 2016) and can be described as follows:

 37P(z)r=( [Ca]+0.38 [K])× Xsc×P0                                  ×ezΛ [(1n)ρg+ n ρw Sw]    (11)

where, Ca and K are calcite and potassium contents in (mg·kg−1), XSC (-) is a scaling factor for latitude and altitude, Po (atoms·cm−3·yr−1·mgCa−1·kgrock) is the surface production normalized for sea level and high latitude (SLHL) (Fabryka-Martin, 1988) and for Ca and K content. Λ (−2) is the attenuation length for fast neutron in soils/rocks. z (g·cm−2) is the product of depth below ground surface (cm) and water density (−3), n (-) is the porosity. ρg (g·cm−3) is the grain density and Sw (-) is the water saturation.

From the above production function, the daily production of 37Ar in water 37P(z)w (atoms·m−3·d−1) can be computed as follows:

 37P(z)w=  37P(z)r ×[ϵn]× 106365    (12)

where, ϵ (-) is the emanation factor, which was determined by irradiation experiments (Johnson et al., 2021; Musy et al., 2022). The attenuation length of cosmogenic fast neutrons depends on the density of the material crossed. Along a given column of water, the saturation and porosity is 1 in Equation 11, resulting in a lower density and higher attenuation length compared to a same column of rocks. Thus, 37Ar production below a stream (37Pstream) is reduced compared to the production in the aquifer at the same depth. The reduction factor depends on the height of the water column h (cm) and the water density ρw (g·cm−3) such that:

 37Pstream=  37P(z)w h ×ρwΛ    (13)

Most of the parameters related to the production of 37Ar should be measured in the field so that realistic concentrations of 37Ar can be generated. Therefore, measurements from a pre-alpine alluvial aquifer were here considered to define the values for the parameters introduced in this section (Table 1). As the attenuation length parameter itself cannot be measured directly in the field, for the illustrative model it has been manually estimated such that the secular equilibrium profile in a virtual background piezometer (Figure 2) reproduces a 37Ar secular equilibrium profile determined from measurements of 37Ar activity concentrations in a sandy-gravel aquifer in the Swiss pre-Alps (Schilling et al., 2017a).


Table 1. Parameter values involved in the production of 37Ar from fast cosmic ray neutrons.


Figure 2. Conceptual model for the synthetic alluvial plain model. The bright blue area represents the stream, which is in losing condition. Five wells separated by a distance of 100 m are located at 50 m from the stream. The virtual background piezometer has been used to estimate the attenuation length for 37Ar production function. The observation point has been used to generate one-dimensional concentration profiles as illustrated in Figure 4. Solid black lines are groundwater head equipotential, illustrating losing river conditions and the influence of pumping wells.

Model set-up

The HGS model for the synthetic alluvial river-groundwater system has a spatial extent of 300 × 500 × 30 m along the x-, y-, and z-directions, respectively (Figure 2). The alluvial plain was gently inclined toward the stream outlet with a slope of 0.003 m/m in the y-direction and discretized with a 2D unstructured triangular mesh with lateral internodal spacing ranging from 4 m along the stream to 7 m on the alluvial plain. The 2D mesh was generated with AlgoMesh (HydroAlgorithmics Pty Ltd, 2016). The model consists of 41 layers of elements used to discretize the alluvial plain in the z-direction, with an exponential range of layer thicknesses varying from 0.05 m at the surface without exceeding more than 1 m at depth (Figure 2). This relatively fine vertical discretization has been chosen to appropriately simulate the vertical and exponential decrease of 37Ar emanation rates. The resulting 3D mesh has a total of 2,87,615 nodes for 5,44,000 elements.

The stream was conceptualized as a 20 m wide and 2 m deep channel (Figure 2). As in all ISSHM, streamflow was simulated explicitly, allowing infiltration of stream water and exfiltration of groundwater to arise naturally, that is, without the need to specify additional internal boundary conditions. Stream water inflow at the upstream end of the model was conceptualized as a second-type (specified flux) boundary condition (BC) with a value of 2 m3·s−1. Surface water outflow was implemented as a critical depth BC to nodes in the stream. Due to very rough surface usually found in alluvial streams, a high value for Manning's n (1.7 × 10−6 d·m−1/3) was used. To ensure pressure continuity between the surface and the subsurface, the coupling length (i.e., the first-order exchange coefficient), was set to a low value of 0.001 m (Liggett et al., 2012).

The aquifer was conceptualized as a homogeneous sandy gravel porous medium with a van Genuchten α of 3.40 m−1, a van Genuchten β of 1.71 (-), and a residual water saturation of 0.05 (-). The hydraulic conductivity (K) and porosity (n) of the homogeneous aquifer were set to 500 m·d−1 and 0.15, respectively, mimicking typical values for sandy-gravel aquifers (Schilling et al., 2017a, 2020). A head-dependent flux (Cauchy-type) BC was used for both the upstream and the downstream groundwater BC. Constant hydraulic heads equal to 99.5 and 93.20 m were considered at the upstream (y = 0 m) and downstream boundary (y = 500 m) of the model, respectively, with a conductance of 5.8 m2·s−1. These boundary conditions were manually adjusted to favor natural losing river conditions. Additionally, five pumping wells were implemented as one-dimensional vertical line elements (along the x-axis = 235 m) screened between 80 and 90 m. From each well, 45 l·s−1 where abstracted via a nodal flux BC set at an elevation of 82 m (Figure 2). Focusing on losing stream conditions was also preferred to better illustrate the use of these natural gas tracers to track the recent infiltration of surface water into a groundwater body. Note that pumping in the wells was not considered when the attenuation length for the 37Ar production function was manually adjusted to reproduce the secular equilibrium profile of 37Ar. This choice was made to ensure that the 37Ar background piezometer (Figure 2) was in secular equilibrium. The lateral model boundaries (x = 0 m and x = 300 m), and the bottom boundary, were assumed to be impermeable. No riverbed was considered, and it was assumed that the stream is fully connected to the aquifer [see Brunner et al. (2009) and Schilling et al. (2017b) for implications]. A meteoric groundwater recharge BC was applied (300 mm·y−1) at the top nodes of the subsurface flow domain to favor vertical flow in the unsaturated zone. The result of this flow model configuration is a shallow water table aquifer (≃ 2.5 m thick unsaturated and 27.5 m thick saturated zone) with a losing river all along the model and localized drawdown at the pumping wells (Figure 2). The proposed model design was inspired by the common configuration of bank filtration wells in mountainous river corridors.

The lateral transport BC, at x = 0 m and x = 300 m, and bottom boundaries, were set as zero mass gradients. The upstream boundary of the model, at y = 500 m, was set as a third-type BC with fixed concentrations at secular equilibrium. Surface nodes act as discharge and/or recharge points regarding the mass gradient between the surface and the subsurface. Due to the negligible atmospheric concentrations as well as the volatility of both 222Rn and 37Ar, a zero concentration BC was specified for the surface nodes. Based on the scaling relationships of Schulze-Makuch (2005), the longitudinal dispersivity of the current transport model was set to 5 m, while the transverse dispersivities were set to 0.5 m. Within the matrix, radiogenic production of 222Rn and 37Ar was set according to the new partitioning zero-order source, while the disintegration was simulated with a first-order decay. The 222Rn aquifer production rate per unit pore volume was computed from Equation 10. The production of 37Ar was computed from Equations 11, 12. For the elements below the stream, the production of 37Ar is slightly attenuated according to Equation 13, assuming a constant water column with a thickness of 1 m. The phase-equilibrium partitioning coefficient Hcc was assigned a value of 0.35 for 222Rn and 0.04182 for 37Ar in mole fraction ratio as defined in Equation 3 for 10°C water temperature.


The model was forced with constant boundary conditions and run long enough to reach secular equilibrium for both 222Rn and 37Ar. A comparison of the simulated concentrations of 37Ar at the virtual background piezometer shows that the simulations of 37Ar are in agreement with the secular equilibrium profile estimated from the real-world sandy gravel alluvial aquifer that was used for calibration of the attenuation length parameter (Supplementary Figure A2). The fast neutron attenuation length was estimated at 98 g·cm−2, which is of the same order as empirically established values (Phillips et al., 2001). Note that the apparent attenuation length observed at the background piezometer is the combined result of depth-dependent production and dispersive smoothing.

The simulated concentrations of 222Rn and 37Ar with the synthetic alluvial plain model under pumping conditions are illustrated in Figure 3. For both gas tracers, for a given depth, groundwater concentrations increase with distance from the stream. This is an expected pattern as the freshly infiltrated water from the stream lowers the groundwater concentrations of 222Rn and 37Ar. Minimum concentrations are located directly underneath the losing stream, corresponding to the location where surface water infiltrates into subsurface. The effect of depth-dependent 37Ar production is clearly visible in the groundwater concentrations.


Figure 3. Simulated steady-state 37Ar and 222Rn concentrations (Bq·l−1) for the synthetic alluvial plain model with the proposed zero-order source feature. Cross sections below are taken at y = 20 m.

The pumping at the wells induces faster flow paths from the edge of the riverbank toward the wells (Figure 3). Thus, high groundwater flow velocities are located between the riverbank and the well corridor, where drawdown is induced by the pumping wells. Within these high flow paths or capture zones, recently infiltrated water spends less time in upper layers and is thus less concentrated in 37Ar compared to other stretches along the stream. The low production rates located in the deeper layers are not sufficient to increase the concentrations of 37Ar and the low concentration signals allow tracking the recently infiltrated surface water (Figure 3). In contrast, the production of 222Rn was assumed independent of depth and therefore, concentrations are only related to by residence times in the alluvial aquifer (Figure 3).

Additionally, at the location of the observation point inside the model domain (see Figure 2), vertical one-dimensional profiles of 37Ar and 222Rn tracer concentrations can be drawn across the 27.5 m saturated portion of the aquifer, thus focusing on groundwater tracer concentrations (Figure 4). The main interest is to evaluate the impact of simulating tracer concentrations in the unsaturated zone on tracer concentrations in groundwater, which correspond to measurements usually gathered in the field. In this section, in addition to the standard (non-partitioning) zero-order source and the proposed new zero-order source with partitioning, a (non-partitioning) zero-order source with saturation threshold (which neglects the production of tracers below a given saturation threshold, assumed here at 99%) is also considered. A saturation threshold fixed at 99% restricts the production of tracers to the saturated zone, essentially neglecting the production of tracers in the unsaturated zone. Groundwater concentration profiles are thus computed for 37Ar and 222Rn with these three approaches.


Figure 4. Simulated pseudo steady-state 37Ar and 222Rn groundwater concentrations (Bq·l−1) profiles for the synthetic alluvial plain model with the (non-partitioning) zero-order source (orange), the zero-order source with partitioning (red) and the (non-partitioning) zero-order source with saturation threshold (green). The left side plots illustrate the concentration difference (in %) between the zero-order source with partitioning (red) and the (non-partitioning) zero-order source with saturation threshold (green). The solid horizontal blue line indicates the location of the water table.

For both tracers, the one-dimensional profiles show that the standard zero-order source implementation, which doesn't scale the production of tracers in the unsaturated zone, leads to a systematic overestimation of the tracer concentrations in groundwater (Figure 4). This large overestimation (up to one order of magnitude for 37Ar) illustrates the limited capacity of the standard zero-order source implementation to model gas-tracer production in variably saturated conditions, which must usually be considered when studying river-groundwater environments (Brunner et al., 2017).

Alternatively, when the tracer production is completely neglected in the unsaturated zone using the zero-order source with saturation threshold implementation, the concentrations in the groundwater are likely to be underestimated for both 222Rn and 37Ar (Figure 4). These are consistent simulation results since the infiltrated water that reached the water table is less concentrated when the unsaturated zone production of these tracers is ignored, potentially reducing the groundwater concentrations close to the water table. This illustrates the major interest of the proposed zero-order source with partitioning for tracers with high solubility but also for tracers that are heterogeneously emanated in the subsurface and where tracer production in the unsaturated zone cannot be ignored.

Despite a lower solubility of 37Ar compared to 222Rn, production of 37Ar is exponentially decreasing with depth and a high proportion of the concentration (50% at the observation point location) is neglected when the tracer is not simulated in the unsaturated zone (Figure 4). This is less significant for 222Rn (difference at about 15% at the observation point location) because production is considered homogeneous and so the impact of not including the larger production in the unsaturated zone is minimized. However, recent studies have shown that 222Rn emanations rates can be highest within the first few meters below the surface (Mullinger et al., 2009; Peel et al., 2022). High 222Rn production rates in the shallow part of the aquifer might have significant effects on measured 222Rn concentrations in groundwater and the proposed zero-order source with partitioning implementation would significantly improve the simulation of 222Rn in such systems. For both tracers, neglecting the production in the unsaturated zone can significantly influence the groundwater concentrations up to 10 m below the ground surface for our river-groundwater alluvial model. Yet more pronounced differences are expected for thick unsaturated zones and highly dynamic transient conditions.


Limitation of the partitioning zero-order source production

In the preceding sections, it has been shown that a more consistent simulation of environmental gas tracer concentrations can be achieved with single-phase ISSHMs if the production in the unsaturated zone is considered with the proposed instant-equilibration mass-balance formulation (i.e., the zero-order source with partitioning). However, it should be noted that the verification section has underlined some limitations under strongly transient conditions since the gas phase is not explicitly simulated. Thus, while the proposed methodology represents a clear and computationally highly efficient improvement of the current capacity of single-phase ISSHMs to explicitly simulate environmental gas tracers, consistent simulation of environmental gas tracer concentrations in the unsaturated zone, should be done with a multi-phase flow and transport model.


A zero-order source implementation that scales the production of environmental gas tracers according to the water saturation and gas/water partitioning has been proposed. The new implementation, labeled as a zero-order source with partitioning, is based on an instant-equilibration mass-balance and allows realistic simulations of environmental gas tracers in the unsaturated zone with single-phase, variably saturated hydrologic models. The new feature, which has been implemented in HydroGeoSphere, has been verified against a multi-phase fluid flow and transport model (HYDRUS-1D). As illustrated in the verification and illustration examples, the application of the proposed partitioning of tracer production rates can significantly improve the consistency of simulations of environmental gas tracer concentrations under variably saturated conditions. The method is particularly well-suited to enhance the capacity of integrated surface and subsurface hydrologic models (ISSHMs) to explicitly simulate environmental gas tracers in steady-state and dynamic conditions, although some limitations has been shown for extremely dynamic contexts. In addition to the presentation of the new zero-order source with partitioning, we demonstrate for the first time the explicit and simultaneous simulation of 222Rn alongside 37Ar under variably-saturated conditions with an ISSHM. These tracers have so far not been jointly simulated.

Data availability statement

Data used for this study are supplied as SI Dataset P2 available for download from HydroShare (

Author contributions

HD: writing—original draft preparation, investigation, software, conceptualization, methodology. MP: writing—review and editing, conceptualization, and methodology. SM: writing—review and editing, and conceptualization. OS: writing—review and editing. RP: writing—review and editing. PB: writing—review and editing, conceptualization, supervision, and funding acquisition. All authors contributed to the article and approved the submitted version.


HD and MP acknowledge the funding provided by the Swiss National Science Foundation (SNSF; grant number 200021_179017).


The authors acknowledge the Aquanty Inc. scientific team for their support and assistance during the implementation of the zero-order source with partitioning feature into the numerical model code HydroGeoSphere.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher's note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at:


Aquanty Inc. (2022). HydroGeoSphere Theory Manual. Waterloo, ON. p. 101.

Bourke, S. A., Cook, P. G., Shanafield, M., Dogramaci, S., and Clark, J. F. (2014). Characterisation of hyporheic exchange in a losing stream using radon-222. J. Hydrol. 519, 94–105. doi: 10.1016/j.jhydrol.2014.06.057

CrossRef Full Text | Google Scholar

Brunner, P., Cook, P. G., and Simmons, C. T. (2009). Hydrogeologic controls on disconnection between surface water and groundwater. Water Resour. Res. 45:W01422. doi: 10.1029/2008WR006953

CrossRef Full Text | Google Scholar

Brunner, P., and Simmons, C. T. (2012). HydroGeoSphere: a fully integrated, physically based hydrological model. Groundwater 50, 170–176. doi: 10.1111/j.1745-6584.2011.00882.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Brunner, P., Therrien, R., Renard, P., Simmons, C. T., and Franssen, H.-J. H. (2017). Advances in understanding river-groundwater interactions. Rev. Geophys. 55, 818–854. doi: 10.1002/2017RG000556

CrossRef Full Text | Google Scholar

Cartwright, I., and Gilfedder, B. (2015). Mapping and quantifying groundwater inflows to Deep Creek (Maribyrnong catchment, SE Australia) using 222Rn, implications for protecting groundwater-dependant ecosystems. Appl. Geochem. 52, 118–129. doi: 10.1016/j.apgeochem.2014.11.020

CrossRef Full Text | Google Scholar

Cecil, L. D., and Green, J. R. (2000). “Radon-222,” in Environmental Tracers in Subsurface Hydrology, eds P. G. Cook, and A. L. Herczeg (Boston, MA: Springer). doi: 10.1007/978-1-4615-4557-6_6

CrossRef Full Text | Google Scholar

Cook, P. G., and Herczeg, A. L. (2012). Environmental Tracers in Subsurface Hydrology. New York, NY: Springer Science and Business Media. doi: 10.1007/978-1-4615-4557-6

CrossRef Full Text | Google Scholar

Cook, P. G., Wood, C., White, T., Simmons, C. T., Fass, T., and Brunner, P. (2008). Groundwater inflow to a shallow, poorly-mixed wetland estimated from a mass balance of radon. J. Hydrol. 354, 213–226. doi: 10.1016/j.jhydrol.2008.03.016

CrossRef Full Text | Google Scholar

Dann, R., Close, M., Flintoft, M., Hector, R., Barlow, H., Thomas, S., et al. (2009). Characterization and estimation of hydraulic properties in an alluvial gravel vadose zone. Vadose Zone J. 8, 651–663. doi: 10.2136/vzj2008.0174

CrossRef Full Text | Google Scholar

Dwivedi, R., Eastoe, C., Knowles, J. F., McIntosh, J., Meixner, T., Minor, R., et al. (2022). Tandem use of multiple tracers and metrics to identify dynamic and slow hydrological flowpaths. Front. Water. 4:841144. doi: 10.3389/frwa.2022.841144

CrossRef Full Text | Google Scholar

Fabryka-Martin, J. T. (1988). Production of Radionuclides in the Earth and Their Hydrogeologic Significance, with Emphasis on Chlorine-36 and Iodine-129. The University of Arizona.

Google Scholar

Gilfedder, B., Cartwright, I., Hofmann, H., and Frei, S. (2019). Explicit modeling of Radon-222 in HydroGeoSphere during steady state and dynamic transient storage. Groundwater 57, 36–47. doi: 10.1111/gwat.12847

PubMed Abstract | CrossRef Full Text | Google Scholar

Goode, D. J. (1996). Direct simulation of groundwater age. Water Resour. Res. 32, 289–296. doi: 10.1029/95WR03401

CrossRef Full Text | Google Scholar

Gosse J. C. Phillips and, F. M. (2001). Terrestrial in situ cosmogenic nuclides: theory and application. Quatern. Sci. Rev. 20, 1475–1560. doi: 10.1016/S0277-3791(00)00171-2

CrossRef Full Text | Google Scholar

Guillon, S., Sun, Y., Purtschert, R., Raghoo, L., Pili, E., and Carrigan, C. R. (2016). Alteration of natural 37Ar activity concentration in the subsurface by gas transport and water infiltration. J. Environ. Radioact. 155, 89–96. doi: 10.1016/j.jenvrad.2016.02.021

PubMed Abstract | CrossRef Full Text | Google Scholar

Heisinger, B., Lal, D., Jull, A. J. T., Kubik, P., Ivy-Ochs, S., Knie, K., et al. (2002). Production of selected cosmogenic radionuclides by muons: 2. Capture of negative muons. Earth Planet. Sci. Lett. 200, 357–369. doi: 10.1016/S0012-821X(02)00640-4

CrossRef Full Text | Google Scholar

Hoehn, E., and Von Gunten, H. R. (1989). Radon in groundwater: a tool to assess infiltration from surface waters to aquifers. Water Resour. Res. 25, 1795–1803. doi: 10.1029/WR025i008p01795

CrossRef Full Text | Google Scholar

HydroAlgorithmics Pty Ltd (2016). AlgoMesh User Guide. Melbourne, VIC.

Google Scholar

Johnson, C., Lowrey, J. D., Alexander, T., Mace, E., and Prinke, A. (2021). Measurements of the emanation of 37Ar and 39Ar from irradiated rocks and powders. J. Radioanal. Nucl. Chem. 329, 969–974. doi: 10.1007/s10967-021-07827-4

CrossRef Full Text | Google Scholar

Kazemi, G. A., Lehr, J. H., and Perrochet, P. (2006). Groundwater Age. Hoboken, NJ: John Wiley and Sons. doi: 10.1002/0471929514

CrossRef Full Text | Google Scholar

Kolditz, O., Bauer, S., Bilke, L., Böttcher, N., Delfs, J. O., Fischer, T., et al. (2012). OpenGeoSys: an open-source initiative for numerical simulation of thermo-hydro-mechanical/chemical (THM/C) processes in porous media. Environ. Earth Sci. 67, 589–599. doi: 10.1007/s12665-012-1546-x

CrossRef Full Text | Google Scholar

Kollet, S. J., and Maxwell, R. M. (2008). Capturing the influence of groundwater dynamics on land surface processes using an integrated, distributed watershed model. Water Resour. Res. 44. doi: 10.1029/2007WR006004

CrossRef Full Text | Google Scholar

Liggett, J. E., Werner, A. D., and Simmons, C. T. (2012). Influence of the first-order exchange coefficient on simulation of coupled surface–subsurface flow. J. Hydrol. 414–415, 503–515. doi: 10.1016/j.jhydrol.2011.11.028

CrossRef Full Text | Google Scholar

Loosli, H. H., and Purtschert, R. (2005). “Rare gases”, in Isotopes in the Water Cycle: Past, Present and Future of a Developing Science, eds P. K. Aggarwal, J. R. Gat, and K. F. Froehlich (Vienna: IAEA), p. 91–95. doi: 10.1007/1-4020-3023-1_7

CrossRef Full Text | Google Scholar

McCallum, J. L., Cook, P. G., Dogramaci, S., Purtschert, R., Simmons, C. T., and Burk, L. (2017). Identifying modern and historic recharge events from tracer-derived groundwater age distributions. Water Resour. Res. 53, 1039–1056. doi: 10.1002/2016WR019839

CrossRef Full Text | Google Scholar

McCallum, J. L., Cook, P. G., and Simmons, C. T. (2015). Limitations of the use of environmental tracers to infer groundwater age. Groundwater 53, 56–70. doi: 10.1111/gwat.12237

PubMed Abstract | CrossRef Full Text | Google Scholar

Mullinger, N. J., Pates, J. M., Binley, A. M., and Crook, N. P. (2009). Controls on the spatial and temporal variability of 222Rn in riparian groundwater in a lowland chalk catchment. J. Hydrol. 376, 58–69. doi: 10.1016/j.jhydrol.2009.07.015

CrossRef Full Text | Google Scholar

Musy, S., Casolaro, P., Dellepiane, G., Berger, A., Braccini, S., and Purtschert, R. (2022). Quantification of 37Ar emanation fractions from irradiated natural rock samples and field applications. J. Environ. Radioact. 251–252:106966. doi: 10.1016/j.jenvrad.2022.106966

PubMed Abstract | CrossRef Full Text | Google Scholar

Nitao, J. J. (1998). Reference Manual for the NUFT Flow and Transport Code, Version 2.0. Lawrence Livermore National Laboratory.

Google Scholar

Paniconi, C., and Putti, M. (2015). Physically based modeling in catchment hydrology at 50: survey and outlook. Water Resour. Res. 51, 7090–7129. doi: 10.1002/2015WR017780

CrossRef Full Text | Google Scholar

Partington, D., Knowling, M. J., Simmons, C. T., Cook, P. G., Xie, Y., Iwanaga, T., et al. (2020). Worth of hydraulic and water chemistry observation data in terms of the reliability of surface water-groundwater exchange flux predictions under varied flow conditions. J. Hydrol. 590:125441. doi: 10.1016/j.jhydrol.2020.125441

CrossRef Full Text | Google Scholar

Peel, M., Kipfer, R., Hunkeler, D., and Brunner, P. (2022). Variable 222Rn emanation rates in an alluvial aquifer: limits on using 222Rn as a tracer of surface water—Groundwater interactions. Chem. Geol. 599:120829. doi: 10.1016/j.chemgeo.2022.120829

CrossRef Full Text | Google Scholar

Phillips, F., Stone, W. D., and Fabryka-Martin, J. T. (2001). An improved approach to calculating low-energy cosmic-ray neutron fluxes near the land/atmosphere interface. Chem. Geol. 175, 689–701. doi: 10.1016/S0009-2541(00)00329-6

CrossRef Full Text | Google Scholar

Popp, A. L., Pardo-Alvarez, A., Schilling, O. S., Scheidegger, A., Musy, S., Peel, M., et al. (2021). A framework for untangling transient groundwater mixing and travel times. Water Resour. Res. 57:2020WR028362. doi: 10.1029/2020WR028362

CrossRef Full Text | Google Scholar

Renne, P. R., and Norman, E. C. (2001). Determination of the half-life of 37Ar by mass spectrometry. Phys. Rev. C. 63:047302. doi: 10.1103/PhysRevC.63.047302

CrossRef Full Text | Google Scholar

Riedmann, R. A., and Purtschert, R. (2011). Natural 37Ar concentrations in soil air: implications for monitoring underground nuclear explosions. Environ. Sci. Technol. 45, 8656–8664. doi: 10.1021/es201192u

PubMed Abstract | CrossRef Full Text | Google Scholar

Schilling, O. S., Cook, P. G., and Brunner, P. (2019). Beyond classical observations in hydrogeology: the advantages of including exchange flux, temperature, tracer concentration, residence time and soil moisture observations in groundwater model calibration. Rev. Geophys. 57, 146–192. doi: 10.1029/2018RG000619

CrossRef Full Text | Google Scholar

Schilling, O. S., Cook, P. G., Grierson, P. F., Dogramaci, S., and Simmons, C. T. (2020). Controls on interactions between surface water, groundwater and riverine vegetation along intermittent rivers and ephemeral streams in arid regions. Water Resour. Res. 57:e2020WR028429. doi: 10.1029/2020WR028429

CrossRef Full Text | Google Scholar

Schilling, O. S., Gerber, C., Purtschert, R., Kipfer, R., Hunkeler, D., and Brunner, P. (2017a). Advancing physically-based flow simulations of alluvial systems through atmospheric noble gases and the novel 37Ar tracer method. Water Resour. Res. 53, 10465–10490. doi: 10.1002/2017WR020754

CrossRef Full Text | Google Scholar

Schilling, O. S., Irvine, D. J., Hendricks Franssen, H. J., and Brunner, P. (2017b). Estimating the spatial extent of unsaturated zones in heterogeneous river-aquifer systems. Water Resour. Res. 53, 10583-10602. doi: 10.1002/2017WR020409

CrossRef Full Text | Google Scholar

Schilling, O. S., Partington, D. J., Doherty, J., Kipfer, R., Hunkeler, D., and Brunner, P. (2022). Buried paleo-channel detection with a groundwater model, tracer-based observations and spatially varying, preferred anisotropy pilot point calibration. Geophys. Res. Lett. 49, e2022GL098944. doi: 10.1029/2022GL098944

CrossRef Full Text | Google Scholar

Schubert, M., Paschke, A., Lieberman, E., and Burnett, W. C. (2012). Air-water partitioning of 222Rn and its dependence on water temperature and salinity. Environ. Sci. Technol. 46. 3905–3911. doi: 10.1021/es204680n

PubMed Abstract | CrossRef Full Text | Google Scholar

Schulze-Makuch (2005). Longitudinal dispersivity data and implications for scaling behavior. Groundwater. 43, 443–456. doi: 10.1111/j.1745-6584.2005.0051.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Sebben, M. L., Werner, A. D., Liggett, J. E., Partington, D. J., and Simmons, C. T. (2013). On the testing of fully integrated surface-subsurface hydrological models. Hydrol. Process. 27, 1276–1285. doi: 10.1002/hyp.9630

CrossRef Full Text | Google Scholar

Simunek, J., Van Genuchten, M.Th., and Sejna, M. (2016). Recent developments and applications of the HYDRUS computer software packages. Vadose Zone J. 15, 1–25. doi: 10.2136/vzj2016.04.0033

CrossRef Full Text | Google Scholar

Spannagel, G., and Fireman, E. L. (1972). Stopping rate of negative cosmic-ray muons near sea level. J. Geophys. Res. 77, 5351–5359. doi: 10.1029/JA077i028p05351

CrossRef Full Text | Google Scholar

Sramek, O., Stevens, L., McDonough, W. F., Mukhopadhyay, S., and Peterson, R. J. (2017). Subterranean production of neutrons, 39Ar and 21Ne: Rates and uncertainties. Geochim. Cosmochim. Acta 196, 370–387. doi: 10.1016/j.gca.2016.09.040

CrossRef Full Text | Google Scholar

Sun, H., and Furbish, D. J. (1995), Moisture content effect on radon emanation in porous media. J. Contam. Hydrol. 18, 239–255. doi: 10.1016/0169-7722(95)00002-D

CrossRef Full Text | Google Scholar

Tosaka, H., Itoh, K., and Furuno, T. (2000). Fully coupled formulation of surface flow with 2-phase subsurface flow for hydrological simulation. Hydrol. Process. 14, 449–464. doi: 10.1002/(SICI)1099-1085(20000228)14:3<449::AID-HYP948>3.0.CO;2-9

CrossRef Full Text | Google Scholar

Turnadge, C., and Smerdon, B. D. (2014). A review of methods for modelling environmental tracers in groundwater: advantages of tracer concentration simulation. J. Hydrol. 519, 3674–3689. doi: 10.1016/j.jhydrol.2014.10.056

CrossRef Full Text | Google Scholar

Van Huijgevoort, M. H. J., Tetzlaff, D., Sutanudjaja, E. H., and Soulsby, C. (2016). Using high resolution tracer data to constrain water storage, flux and age estimates in a spatially distributed rainfall-runoff model. Hydrol. Process. 30, 4761–4778. doi: 10.1002/hyp.10902

CrossRef Full Text | Google Scholar

Zhuo, W., Lida, T., and Furukawa, M. (2006). Modeling radon flux density from the Earth's Surface. J. Nucl. Sci. Technol. 43, 479–482. doi: 10.1080/18811248.2006.9711127

CrossRef Full Text | Google Scholar

Keywords: environmental gas tracers, integrated surface and subsurface hydrological models, mass-transport, subsurface hydrology, unsaturated zone processes

Citation: Delottier H, Peel M, Musy S, Schilling OS, Purtschert R and Brunner P (2022) Explicit simulation of environmental gas tracers with integrated surface and subsurface hydrological models. Front. Water 4:980030. doi: 10.3389/frwa.2022.980030

Received: 28 June 2022; Accepted: 23 August 2022;
Published: 28 September 2022.

Edited by:

Giacomo Bertoldi, Eurac Research, Italy

Reviewed by:

Reed Maxwell, Princeton University, United States
Majdi Mansour, The Lyell Centre, United Kingdom

Copyright © 2022 Delottier, Peel, Musy, Schilling, Purtschert and Brunner. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

*Correspondence: Hugo Delottier,

Disclaimer: All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article or claim that may be made by its manufacturer is not guaranteed or endorsed by the publisher.