The MGC-scale daily precipitation (P) and streamflow (Q) data were calculated between WY 2008 and 2017 (Figures 2A,B; Dwivedi et al., 2019b, 2021). Precipitation was observed at 15-min intervals at eight measurement sites equipped with tipping bucket precipitation gages at seven locations and a heated precipitation gage at the remaining site (Figure 1). From the precipitation time series, Thiessen polygon-derived weights were used to estimate daily catchment-scale mean precipitation (Dwivedi et al., 2019b). Streamflow was measured at 30-min intervals at the MG-Weir site (Figure 1) using a pressure transducer (HOBO by Onset U20-001-01; Onset Computer Corporation) with maximum error of 0.62 kPa and accuracy of 0.02 kPa, and a previously derived stage-discharge relationship (Heidbüchel et al., 2012).

We calculated the amount-weighted time series of Tritium (³H) in Tucson precipitation since 1992 using the following: (i) Eastoe et al. (2004), (ii) unpublished data of the Environmental Isotope Laboratory, The University of Arizona (access date: December 13, 2021); and (iii) C. Eastoe, unpublished data. At low levels of ³H concentration in precipitation, the data had a 1 standard deviation (σ) precision of ±0.5 TU (tritium units) or less. Annual concentration cycles in the ³H data were best captured by the post-2001 period in which semi-annual aggregates represent all precipitation events, and demonstrated a prominent 2–4 year cycle between 2001 and 2009 (Supplementary Figure 1). Between 1992 and 2001, the data represented only large precipitation events. Prior to 1992, mean annual ³H concentrations in Tucson precipitation were obtained using the Doney et al. (1992) model (gray points in Supplementary Figure 2A) as well as log-linearly interpolated values (orange points in Supplementary Figure 2A) between observed or modeled values to obtain amount-weighted precipitation ³H concentrations at the half-yearly time scale (Supplementary Figure 2A). All data were corrected for the elevation difference between the Tucson (747 m asl) and MGC (mean elevation 2,428 m asl) by using the ³H concentrations in three simultaneous precipitation samples (Supplementary Figure 2B) from Tucson and the Palisades Ranger Station (2,422 m asl), Santa Catalina Mountains. Amount-weighted precipitation ³H concentration time series data are shown in Figure 3 and Supplementary Figure 2C.

The ³H time series data for deep groundwater, i.e., the groundwater residing in the fractured bedrock aquifer at MGC, was assembled from observations in streamflow (baseflow conditions; n = 9) and groundwater from Pigeon Spring (n = 5; Supplementary Table 1 and Figure 3 inset). The five discharge measurements from Pigeon Spring (Figure 1) were considered representative of deep groundwater (Dwivedi et al., 2019a). Data are grouped into half-yearly brackets using the following criteria: (i) sampling months 6–10 of year n, and (ii) sampling months 11 of year n-1 to 5 of year n. For groups with three or more measurements, the data are expressed as a mean ± 1σ (Figure 3 inset).

In this study, stable water isotope-based TTD and mTT estimates are based on our previous work (Dwivedi et al., 2021) that estimated TTD and mTT from long-term stable water isotope data in precipitation and streamflow. Using wavelet analysis of tracer flux (i.e., a product of tracer concentration and hydrologic flux for precipitation and stream water), TTD and mTT values were obtained from objectively evaluating the following TTD types: Exponential (Exp), Gamma (Gam), Fixed-path one-dimensional advection-dispersion model [ADE-1x; Maloszewski and Zuber (1982)], Multiple-paths advection-dispersion model [ADE-nx; Kirchner et al. (2001)], and Piston flow (PF). More details are provided in Dwivedi et al. (2021).

The ³H time series data were used in conjunction with the Stewart et al. (2017) method, i.e., Equation (1) of their work, to estimate the mTT and best fitting TTD for deep groundwater. Previous work at MGC determined that deep groundwater was recharged at a time scale of 0.5 years or less (Ajami et al., 2011; Dwivedi et al., 2019a). Because the half-life of ³H is 12.32 years (Lucas and Unterweger, 2000), decay of ³H during recharge was insignificant relative to analytical precision (± 0.5 TU). Therefore, the ³H time series in precipitation was used as the ³H time series in recharge to deep groundwater (Figure 3).

The TTD types evaluated using ³H data had the following fitting parameters: Exp- mTT; Gam- shape parameter (α) and mTT; ADE-1x- Péclet number (Pe) and mTT; ADE-nx- Péclet number (Pe) and mTT; and PF- mTT. The following ranges of TTD model parameters were considered: mTT: 1–50 years when using tritium, which serves as a groundwater age tracer at a time scale of 1–50 years (Aggarwal, 2013; Suckow, 2014; Gleeson et al., 2015); shape parameter (α) for the Gamma TTD: 0.1–15 (Stewart et al., 2017); average catchment-scale Péclet number (Pe): 0.1–100 (Kirchner et al., 2001; Kirchner and Neal, 2013). When evaluating TTD types, the Downhill Simplex method (Nelder and Mead, 1965; Gupta, 2016) was used to search for optimal model parameters and the modified Kling Gupta efficiency or KGE' (Gupta et al., 2009; Kling et al., 2012) was used as the model performance criterion. A perfectly fitting model will have a KGE' value of zero and the worst fitting model will have a KGE' value of ∞. Following Godsey et al. (2010), the KGE' criterion was estimated in a log-transformed space. Both KGE' and the characteristics of the criterion response surface were utilized to search for the optimum model parameters (Dwivedi et al., 2021). Furthermore, additional tests were also conducted to constrain how uncertainty in amount-weighted precipitation and stream water ³H concentrations affected both TTD types and its parameters. In these tests, various combinations of uncertainty in amount-weighted precipitation and stream water ³H concentrations were considered with respect to their means (μ) ± one standard deviation (σ) (see Figure 3), i.e., μ-σ, μ, and μ + σ were specifically related to three possible cases of deep groundwater ³H concentration uncertainty (Supplementary Table 2). Subsequently, both TTD type and TTD parameters were objectively estimated using the previously mentioned TTD types. The same allowable parameter space for all model parameters was used to evaluate the most suitable TTD type and parameters.

F_yw can be estimated from the amplitude ratio of tracer concentrations in outflow and inflow for any tracer (Kirchner, 2016a; von Freyberg et al., 2018). Thus, if the amplitudes of the tracer concentrations in outflow and inflow for period λ are A_Q(λ) and A_P(λ), respectively, then:

This method is preferred when sufficient long-term tracer data are available and can be applied without a priori knowledge of the TTD type (Kirchner, 2016a). However, if a TTD h(τ) is known, F_yw can also be estimated using Equation (2), in which T_yw is the threshold age for young water, defined as the upper limit in Equation (2) for which both Equations (1) and (2) provide equivalent values of F_yw(λ):

For a Gamma type TTD, Equation (2) can be simplified to Equation (3) where β is defined as the ratio of mTT (y) to the shape parameter (α, unitless), κ is the decay constant for a tracer [log_e(2)/half-life; 0 for stable water isotope and 0.056 yr⁻¹ for ³H], and ω is defined as 2 π/λ:

In this study, δ¹⁸O-based F_yw calculations were based on tracer fluxes in precipitation and streamflow, whereas the ³H-based F_yw estimates were based on tracer concentrations in recharging and discharging deep groundwater. For consistency with the literature, we express tracer flux-based fraction of young water values as Fyw* and tracer concentration-based fraction of young water values as F_yw.

The F_yw or Fyw* are typically estimated by fitting sinusoidal curves to long-term tracer data in precipitation and streamflow (e.g., Jasechko, 2016; von Freyberg et al., 2018). In the current study, stable water isotope-based Fyw* was specifically estimated using the iteratively re-weighted least squares (IRLS) method that determines tracer cycle amplitude by fitting sinusoidal functions to tracer flux data in precipitation and streamflow (Equation 1; Kirchner, 2016a; von Freyberg et al., 2018), and the TTD method (Equation 3). The tracer flux is the product of tracer concentration and hydrologic flux, and daily precipitation (aggregated into time intervals corresponding to the availability of isotope data in precipitation) and streamflow data are used to calculate the fraction of precipitation contributing to streamflow. The IRLS method was specifically applied to estimate Fyw* for periods ranging from 2 days to 5 years.

The ³H tracer data in precipitation and deep groundwater at MGC were sparse relative to δ¹⁸O, especially for deep groundwater (Figure 3). As a result, application of sinusoidal curve fitting methods such as IRLS to deep groundwater ³H data was unsatisfactory i.e., fitting a sinusoidal function to the data resulted in significant amplitude uncertainty at a period of 1 year. Therefore, the ³H F_yw estimates were based on the TTD method (Equations 2 and 3 above).

When using stable water isotope data, the temporal variability of δ¹⁸O in precipitation was addressed by means of Fyw* uncertainty analyses where the temporal variability of δ¹⁸O was expressed as three statistics: daily mean, mean + 1σ, and mean −1σ calculated for both precipitation (P) and stream water (Q); consideration of all pair combinations resulted in nine scenarios. For each period, the minimum, mean (referred to as the ensemble mean below), maximum, and 1σ of the Fyw* results were computed for all nine scenarios. Finally, the ensemble means for F_yw estimates from ³H data were calculated similarly, using mean, mean + 1σ, and mean −1σ of the Gamma TTD parameters α and mTT.

Piston-flow (PF) and Gamma (Gam) TTD types performed adequately and yielded TTD parameters within permissible parameter spaces (Table 1, Figure 4), whereas the other TTD types (Exp, ADE-1x, ADE-nx) did not. In addition, the response surfaces for certain TTD types were found to be rough (e.g., Figure 4E), and therefore, the model parameters for each TTD type were estimated from three separate model runs (Figure 4 and Supplementary Table 3). The model performance expressed in terms of KGE' was slightly better for the PF relative to the Gam TTD type (PF KGE' was 29% lower). Comparison of the PF and Gam TTDs further suggested “approximate equifinality” (Kirchner, 2016b) in the PF TTD results (Figures 4B–E For example, three separate PF TTD model runs yielded mTTs of 32.5, 29.5, and 35.5 years (Figure 4E; Case 5 in Supplementary Table 3) with very similar KGE' values (~0.4; Case 5 in Supplementary Table 3). In contrast, three separate Gam TTD runs yielded similar mTTs (mTT ~ 26 years; minimum: 25.9 years and maximum: 30.4 years) and α parameters (5.23; minimum: 2.2 and maximum: 14.6). Overall, the PF mTT varied between 4 and 33 years with a coefficient of variation of 0.57 (Table 1). The Gam mTT and α parameters varied between 26 and 30 years (mean = 27 years; coefficient of variation = 0.05) and 2.17 to 14.58 (unitless) (mean = 6.53; coefficient of variation = 0.64), respectively.

Modeled ³H concentrations in deep groundwater were generally within their observed ranges, and remained within 1σ of their simulated means for both PF and Gam TTDs, but modeled ³H concentration variability was lower for the Gam TTD (Figure 5). The error bars in Figure 5 were determined by considering the analytical uncertainty in ³H concentrations and were based on the first set of model runs. Although the ADE-1x TTD-based model produced ³H concentrations that were within their observed ranges, the estimated TTD parameters were sometimes at the edge of the allowable parameter space (section 3H-Based TTD Type and mTT); simulated ³H concentrations from the ADE-nx TTD-based model were far from observed concentrations and indicated that this TTD type is not applicable at MGC.

Considering Fyw* variability due to δ¹⁸O variability in precipitation and stream water, the IRLS method yielded highly variable results, particularly for periods of <1 year (Figure 6A). For a seasonal tracer cycle, i.e., period = 0.5 years, the Fyw* estimate was 52 ± 2.1 % (mean ± one standard deviation). However, for periods close to 0.5 years, Fyw* estimates were 17 ± 0.3 % (for λ = 0.45 years) and ~99% (for λ = 0.56 years), respectively. Similarly, for the annual tracer cycle, the Fyw* estimate was 11 ± 0.7 %, but for periods close to 1 year, Fyw* estimates were 22 ± 1.3% (λ = 0.83 years) and 99 ± 0.6% (λ = 1.25 years). Note that the field site experiences two dominant wet seasons: a winter wet season due to Pacific cold fronts and a summer wet season due to the North American Monsoon (Eastoe and Dettman, 2016) such that both seasonal and annual periods may occur in isotopic time series data for precipitation. Thus, the Fyw* estimates using sinusoidal curve fitting to the tracer cycle data showed significant variability with period, which is not yet recognized in the F_yw literature.

The time series of ³H in groundwater and streamflow (Figure 3, inset) were too sparse and coarse for reliable estimation of A_Q/A_P using the IRLS method (Equation 1; see also Supplementary Section Sources of uncertainty for F_yw when considering annual cycles of ³H in precipitation and deep groundwater). Instead, the TTD method (Equation 3) was used to calculate F_yw with model parameters drawn from Section 3H-based TTD type and mTT for a Gamma TTD. The pattern of F_yw with period was characterized by a gradual increase with period (Figure 6B). For an annual tracer cycle, the ensemble mean-based F_yw was 1.6 ± 2.4 × 10⁻³ % (blue triangle in Figure 6B). For the highest period considered in our analysis, 27 years, F_yw was 6.4 ± 1.1%. Between periods 1 and 27 years, F_yw aand its slope with respect to period increased gradually. For example, dF_yw/dλ changed from 0.06%/year (5 < λ < 10 years), to 0.19%/year (10 < λ < 15 years), to 0.36%/year (15 < λ < 20 years), and finally 0.47%/year (20 < λ < 25 years).

The tritium-based F_yw values reported in the current study are significantly lower than previously published estimates. The ensemble-mean F_yw determined from annual ³H cycles at MGC was 1.6 x 10⁻³ %, or effectively 0% (Figure 6B; Section F_yw based on ³H). For comparison, the lowest tritium-based F_yw value in Stewart et al. (2017) for an annual cycle was ~8% for a system composed of two homogeneous sub-systems, each having an mTT of 25 years with a Gamma TTD shape parameter α = 10. We attribute this difference i.e., 8 vs. ~0%, to either the constant T_yw value that Stewart et al. (2017) used for both α and tracer cycle period, or the use of multiple lumped parameter models by Stewart et al. (2017), as opposed to the F_yw values in the current work that are based on fitting a single TTD to the whole ³H dataset for deep groundwater. Using TTD parameter estimates from rSAS (rank StorAge Selection) functions, Rodriguez et al. (2021) reported a F_yw estimate of 1.8% for a forested headwater catchment that is closer to the near-zero F_yw estimate at MGC. This difference may result from differences in sampling protocols and/or calculation methods; Rodriguez et al. (2021) sampled stream water under varying flow conditions, in contrast to deep groundwater and base flow sampling in the current work.

A near-zero ³H-based F_yw at MGC calls into question the suitability of the ³H-based F_yw approach for deeper groundwater. The deep groundwater residing in fractured bedrock at MGC has a large mTT (~27 years), and the annual tracer cycle will be highly damped as a result (F_yw~ 0; Section F_yw based on ³H). From the standpoint of the IRLS method (Equation 1), useful ³H data in precipitation or deep groundwater should have an amplitude A_P or A_Q greater than the ³H measurement precision (0.5 TU, Section Tritium in Precipitation and Deep Groundwater) for some period between 1 and 27 years. At MGC, this is true for A_P at periods <19 years, but not for A_Q at any period up to 27 years (Supplementary Section Sources of uncertainty for F_yw when considering annual cycles of ³H in precipitation and deep groundwater); consequently, the available data are inadequate for calculating F_yw. This is apparent in the lack of consistent annual periodicity in the Tucson Basin precipitation data (Figure 3), which may not be possible to overcome even with a much larger ³H dataset.

Deep groundwater storage within the fractured bedrock aquifer is important because it recharges the surrounding valley fill aquifer through mountain front and/or mountain block recharge pathways (Wilson and Guan, 2004; Ajami et al., 2011). In the traditional approach [see Rodriguez et al. (2021) and references therein], subsurface storage is a function of the tracer-based mTT multiplied by the long-term mean discharge. Application of this approach to MGC for deep groundwater results in a storage estimate of 6.7 m using a ³H-based mTT of 27 years and the observed long-term streamflow (WY 2008 to WY 2017). In contrast, Gleeson et al. (2015) determined that global present-day groundwater is equivalent to a depth of only 3 m on the land surface. Acknowledging such issues, Kirchner (2016b)—using stable water isotopes and virtual experiments—suggested the use of volume-weighted rather than time-weighted mTTs. However, Peters et al. (2013) and Dwivedi et al. (2021) have shown that volume-weighted mTTs and time-weighted mTTs differ by only a factor of ~2 as opposed to orders of magnitude. Volume-weighted mTTs are also difficult to obtain via ³H data due to the lack of multi-decade observations of both streamflow and tracer concentrations in outflow. Further, the use of mean long-term discharge is problematic as the contribution of deep groundwater to streamflow is likely to be lower than contributions from soil water or other near surface storages. Using end-member mixing analysis, Dwivedi et al. (2019a) reported that deep groundwater contributes 4.5% of the long-term streamflow at MGC. If this fraction is included in the storage calculation, then the storage estimate decreases to a more plausible 0.3 m. Taken together, the results of the current work demonstrate that the tracer-based mTT can be used to estimate storage volumes, provided they are interpreted with appropriate caution.