3.1 Young water fraction estimation from seasonal cycles of stable water isotopes in precipitation and stream water: the “direct” input

Kirchner (2016a) designed the young water fraction as the proportion of the transit time distribution younger than a threshold age (τ_yw). By assuming that the transit time distribution mathematical form is the regularized lower incomplete gamma function for all the study catchments, the theoretical young water fraction ($F_{\mathrm{yw}}^{\mathrm{T}}$) can be expressed as

where α and β are the shape and scale factor, respectively.

By using thought experiments, Kirchner (2016a) has demonstrated that for a given shape factor α (ranging from 0.2 to 2) and across a wide range of scale factors β, the theoretical young water fraction can be accurately predicted by the amplitude ratios of seasonal sine curves fitted to stream water and precipitation isotope values by considering a τ_yw of 2–3 months. Operatively, we model seasonal isotope (e.g., δ¹⁸O) cycles in stream water and precipitation as reported in Eqs. (3a) and (3b):

where δ¹⁸O (‰) is the isotopic composition of water sampled at the time t (expressed in decimal years), A (‰) is the amplitude of the seasonal isotope cycle, φ (in radians, with 2π rad = 1 year) is the phase, f (yr⁻¹) is the frequency and k (‰) is the vertical offset of the seasonal isotope cycle. The subscript S refers to stream water, while the subscript P refers to precipitation. The sine wave is fitted to the isotopes measured in precipitation weighted according to the volume of precipitation, reducing the influence of low-precipitation periods and accounting for temporally aggregated rainfall samples (von Freyberg et al., 2018); the sine fit of stream water isotope measurements can be discharge-weighted, using the discharge measured at the moment of sampling as weights, or not (von Freyberg et al., 2018). The sine curves of Eqs. (3a) and (3b) are fitted on the isotope measurements using the iteratively re-weighted least squares (IRLS) method (for reducing the influence of outliers), which leads to estimates of A, φ and k parameters. A function for performing a sine fit using IRLS, based on the IRLS function made available by Kirchner and Knapp (2020), is available in the Supplement.

Accordingly, depending on the unweighted or the flow-weighted fit, an unweighted amplitude (A_S) or a flow-weighted amplitude ($A_{\mathrm{S}}^{*}$) can be obtained, respectively. Such amplitudes can be used to calculate the time-weighted (Eq. 4a) or the flow-weighted (Eq. 4b) young water fractions (F_yw or $F_{\mathrm{yw}}^{*}$, respectively) via the “amplitude ratio approach”:

Gallart et al. (2020a) highlighted the advantages of the flow-weighted analysis (generally yielding $A_{\mathrm{S}}^{*}$ greater than A_S) to compensate for subsampled high-flow periods, which would otherwise lead to a young water fraction underestimate. Accordingly, in this work, we calculate the flow-weighted young water fractions for all the study catchments by applying Eq. (4b). We obtain the standard errors (${\mathrm{SE}}_{F_{\mathrm{yw}}^{*}}$) of the estimated $F_{\mathrm{yw}}^{*}$ using Eq. (5):

where ${\mathrm{SE}}_{A_{\mathrm{S}}^{*}}$ and ${\mathrm{SE}}_{A_{\mathrm{P}}}$ are the standard errors of the regression coefficients $A_{\mathrm{S}}^{*}$ and A_P. The analytical choice of using the amplitude (A_P) fitted to precipitation isotopes, instead of the amplitude (A_Peq) fitted to the equivalent precipitation (i.e., rain plus snowmelt) isotopes, for estimating $F_{\mathrm{yw}}^{*}$ implies that the snowpack (and/or the glacier) is considered as part of the catchment storage. Thus, the damped seasonal cycle observed in the stream is given by the mixing of precipitation with snowpack (and/or the glacier) and subsurface storage (the last two considered as a single entity), as illustrated in Fig. 4.

Figure 4Schematic representation of the “direct input” approach for estimating $F_{\mathrm{yw}}^{*}$. Light blue arrows indicate that meltwater coming from the snowpack preferentially infiltrates. The term C refers to the isotopic composition. The subscript P refers to precipitation, S refers to stream, SP refers to snowpack, Sub refers to subsurface storage and St refers to catchment storage.

Download

Since we have assumed that the transit time distribution belongs to the family of gamma distributions, we can determine the parameter α of such a distribution by solving the following implicit expression for α (Eq. 6):

We optimized to find the best solution of Eq. (6) using the best-fitting parameters (${\mathit{\phi }}_{\mathrm{S}}^{*}$, φ_P, $A_{\mathrm{S}}^{*}$, A_P) with α spanning a wide interval between 0.01 and 20. For the relevant math, the reader is referred to Kirchner (2016a). We estimate the uncertainty of α assuming all the fitted parameters having a Gaussian distribution with a standard deviation equal to the regression error (Eqs. 3a and 3b). We generate 1000 random samples of such parameters, and we estimate the optimal α for each parameter set and then compute their standard deviation.

As in past studies (Gallart et al., 2020a; Lutz et al., 2018), we use the parameter α to estimate τ_yw with the following second-order polynomial fit (Kirchner, 2016a):

where T is the period of the tracer cycle (for a seasonal cycle, T=1 year). We estimate the uncertainty of τ_yw by calculating the standard deviation of the threshold ages computed using the 1000 optimal α values previously obtained. A code for estimating α and τ_yw with their uncertainties is made available in the Supplement.

The comparison of $F_{\mathrm{yw}}^{*}$ among different catchments is potentially subject to a bias given by possibly different threshold ages ranging between 2 and 3 months. Accordingly, we couple each $F_{\mathrm{yw}}^{*}$ with the corresponding τ_yw to illustrate what the term “young” means for each study catchment.

3.2 Snow cover persistence quantified through the mean fractional snow cover area (F_SCA)

In this paper, we quantify the snowpack persistence by calculating the mean fractional snow cover area (F_SCA). It is calculated for each catchment over the period 1 October 2017–30 September 2021 (hereafter defined as PoC, i.e., period of calculation) by using the collection of Sentinel-2 L2A satellite images available in Google Earth Engine (Gorelick et al., 2017). Temporally, this relatively recent satellite has increased the visitation frequency to a subweekly temporal resolution and increased the spatial resolution to 20 m for snow cover (Gascoin et al., 2019). High temporal resolution makes Sentinel-2 images preferable to Landsat images, which are available only once every 16 d and whose total number is often further reduced because of cloudiness (Hofmeister et al., 2022). The PoC generally differs from the PoS for the 27 study catchments. This is because Sentinel-2 L2A satellite images are not available before March 2017. For each image available in the PoC, we calculate the normalized difference snow index (NDSI) as suggested in the work of Dozier (1989):

where r_green is the reflectance in the green band (Sentinel-2 band 3) and r_SWIR is the shortwave infrared reflectance band (Sentinel-2 band 11). We classify as snowy pixels those with an NDSI value > 0.4 (Dozier, 1989). Based on the pixel-by-pixel snow classification, we compute the snapshot fractional snow cover area (f_SCA) according to Di Marco et al. (2020) and Hofmeister et al. (2022):

where N_snow is the number of snow cover pixels according to the applied NDSI threshold method, N_tot is the total number of pixels within the catchment area and N_cloudsis the number of pixels classified as clouds and water bodies (Hofmeister et al., 2022). We identify the cloudy pixels directly using the Sentinel-2 band “Scene Classification Map”. We operatively calculate N_snow, N_tot and N_clouds using a Google Earth Engine code.

By using this procedure for calculating f_SCA, we sometimes obtain f_SCA>1. The NDSI threshold method is generally able to distinguish between snow and no-snow pixels (Aalstad et al., 2020). Accordingly, clouds and snow have similar reflectance in the green band, but clouds highly reflect in the shortwave infrared band, while snow reflectance is low in this band. Thus, the N_snow estimation is generally accurate. On the other hand, it is necessary to disregard the false positive pixels deriving from cloud detection (i.e., snow classified as clouds). If f_SCA>1, we calculate f_SCA as $N_{\mathrm{snow}}/N_{\mathrm{tot}}$, since this is the only heuristic solution that guarantees no overestimation. Moreover, by looking at sample Sentinel-2 images during the summer periods for all the catchments, we impose f_SCA=0 during July and August, since when f_SCA≠0, this usually results from clouds falsely identified as snow: imposing f_SCA=0 clearly leads to fewer errors (only missing occasional summer snowfall events of very shallow depth) than falsely accounting for (far more) frequent clouds. The NBPV catchment is an exception: we do not impose f_SCA=0 during July and August, since it generally has snow over the glacier also during summer. Finally, we compute the mean fractional snow cover area (F_SCA) for each catchment by averaging all f_SCA values available for all snow images in the PoC, without interpolation between the time steps.

3.3 Fraction of Quaternary deposits, low-flow duration and the groundwater contribution to the stream

Similarly to Arnoux et al. (2021), we calculate the portion of the catchment area occupied by Quaternary deposits (A_qd) (available from government geological data sets) with respect to the total catchment area (A). Thus, we calculate the fraction of Quaternary deposits (F_qd) as reported by Eq. (10):

Additionally, we use the same winter flow index (WFI) as Arnoux et al. (2021), as indicated by Eq. (11):

where Q_NM7 is the minimum discharge over 7 consecutive days during the winter period (from November to June) and Q_mean is the mean annual discharge. We calculate it for the 27 study catchments during the PoS. To relate WFI to low flow, we apply the recent baseflow separation technique described by Duncan (2019) to the discharge time series of the 27 study catchments (within the PoS indicated in Table 1). In short, this method comprises a single backward pass through the data to fit an exponential master baseflow recession curve (Eq. 12a), followed by the single forward pass (Eqs. 12b and 12c) of the Lyne and Hollick (1979) algorithm. This allows the smoothing of the connection between segments of the master recession by simulating a gradual groundwater recharge during the runoff event (Duncan, 2019):

where M(t_i), Q_q(t_i) and Q_bf(t_i) are the master recession value, the quick recession flow and the baseflow at time t_i, respectively. In this study, we consider daily time steps (i.e., $t_i-t_{i-\mathrm{1}}=\mathrm{1}$ d). This method has two parameters: k is the recession constant, c is a constant flow added to the exponential decay component. We set the recession constant k=0.925 (Nathan and McMahon, 1990): we add no constant flow to the exponential decay (i.e., in terms of the method by Duncan (2019), c=0). A code with the implementation of the Duncan (2019) baseflow filter has been made available in the Supplement.

To express the catchment storage contribution to streamflow in a form that is directly comparable to the $F_{\mathrm{yw}}^{*}$, we define the baseflow fraction (F_bf) as reported in Eq. (13):

where Q_bf(t_i) is the baseflow (mm d⁻¹) at the time t_i (obtained as indicated by Eq. 12c) and Q(t_i) is the discharge (mm d⁻¹) at the time t_i. We tested the uncertainty of k by drawing random samples (10 000) from a normal distribution spanning Nathan and McMahon's (1990) recommended range for k from 0.90 to 0.95, with a mean of 0.925 and a standard deviation equal to 25 % of the range. Thereby, we obtain 10 000 values of F_bf for each catchment of which we compute the standard deviation.

As introduced in Sect. 1, $F_{\mathrm{yw}}^{*}$ can be low if the snapshot young water fraction F_yw(t_i) is very low for many time steps. If we consider the discharge (Q) as a proxy for the catchment wetness, we can reliably assert that F_yw(t_i) is low for low Q(t_i). Thus, another important variable is the duration of the low-flow period. In this study, we define a low-flow period (T_Low) as follows:

Thus, a low-flow period is defined here as a period when 85 % of the total flow is composed of baseflow (i.e., baseflow dominated). Accordingly, we define the low-flow duration (LFD) as the proportion of the time steps (e.g., days) in the PoS that can be considered as a low-flow period according to Eq. (14).

Figure 5Sinusoidal cycles of both precipitation and streamflow fitted to the δ¹⁸O data (using the IRLS method) for six representative study catchments. Amplitudes (‰) and phases (year) are indicated in the figure. Please note that both δ¹⁸O data and sinusoidal cycles of precipitation and stream water are vertically shifted of k_P and k_S, respectively.

Download

4.1 Young water fractions ($F_{\mathrm{yw}}^{*}$) and corresponding threshold ages

Assembling $F_{\mathrm{yw}}^{*}$ values determined by different authors who very likely used different source codes could possibly result in a bias. Indeed, differences in $F_{\mathrm{yw}}^{*}$ among catchments could be driven by the different methods rather than the physical factors. Therefore, the same approach must be applied to all the study catchments to remove the bias introduced by the estimation method of $F_{\mathrm{yw}}^{*}$. For all the study catchments, sinusoidal cycles were fitted to both precipitation and stream water δ¹⁸O data by using the IRLS regression (results for six representative study catchments in Fig. 5; complete results in Fig. S2). We estimate $F_{\mathrm{yw}}^{*}$ via Eq. (4b) by using the best-fitting amplitudes of seasonal cycles. The best-fit amplitudes (A_P, $A_{\mathrm{S}}^{*}$), phases (φ_P, ${\mathit{\phi }}_{\mathrm{S}}^{*}$) and corresponding standard errors are reported in Table 2.

Table 2Summary table with all the relevant quantities estimated for the 27 study catchments.

Download Print Version | Download XLSX

The $F_{\mathrm{yw}}^{*}$ values achieved in this study are consistent with published $F_{\mathrm{yw}}^{*}$ values of Ceperley et al. (2020) and with direct-input $F_{\mathrm{yw}}^{*}$ of von Freyberg et al. (2018). Accordingly, a two-sample Kolmogorov–Smirnov test accepts the null hypothesis that the new and past $F_{\mathrm{yw}}^{*}$ estimates are from the same continuous distribution at the 5 % significance level. The $F_{\mathrm{yw}}^{*}$ estimates are reported in Fig. 6a against the mean catchment elevation and are also listed in Table 2. $F_{\mathrm{yw}}^{*}$ increases with mean catchment elevation until 1500 m a.s.l., which corresponds to the elevation above which all catchments are snow dominated (with NBPV detected as an outlier as will be discussed in Sect. 4.5). This pattern is also reflected by the median $F_{\mathrm{yw}}^{*}$ within each flow regime: the median $F_{\mathrm{yw}}^{*}$ is 0.13 for rainfall-dominated catchments, 0.29 for hybrid catchments and 0.12 for snow-dominated catchments. Such results are consistent with previous studies that have shown the tendency toward low $F_{\mathrm{yw}}^{*}$ in mountainous catchments (Ceperley et al., 2020; von Freyberg et al., 2018; Lutz et al., 2018; Jasechko et al., 2016).

Figure 6(a) $F_{\mathrm{yw}}^{*}$ as function of mean catchment elevation. Points dimension is proportional to τ_yw (b), obtained with Eq. (14) of Kirchner (2016a), as function of the shape factor α.

Download

Even though we remove the bias introduced by the estimation method of $F_{\mathrm{yw}}^{*}$, the application of Eq. (4b) implicitly introduces another bias if we want to use $F_{\mathrm{yw}}^{*}$ for intercomparison purposes (which is the goal of this work). By computing the amplitude ratio, we estimate $F_{\mathrm{yw}}^{\mathrm{T}}$ (Eq. 2), without defining a corresponding τ_yw (which can be estimated via Eq. 7). Therefore, part of the scatter of $F_{\mathrm{yw}}^{*}$ between catchments might be because of different τ_yw rather than physical factors, also if the τ_yw is expected to vary modestly between 2 and 3 months (Kirchner, 2016a). Accordingly, past studies estimated $F_{\mathrm{yw}}^{*}$ using the amplitude ratio approach without information about the corresponding τ_yw (Stockinger et al., 2019; von Freyberg et al., 2018; Jasechko et al., 2016). Nevertheless, we estimate τ_yw and α for each study catchment: the resulting estimates are reported in Fig. 6b and are listed in Table 2. Our estimates of the α parameter ($\mathrm{0.19}\le \mathit{\alpha }\le \mathrm{2.1}$) are consistent with the shape factor range ($\mathrm{0.2}\le \mathit{\alpha }\le \mathrm{2}$) investigated by Kirchner (2016a). Consequently, the τ_yw obtained by applying Eq. (7) falls between 1.38 and 3.16 months. As expected, τ_ywvaries in a narrow range consistent with the explanation provided by Kirchner (2016a). However, in order to have a fully coherent metric for all the catchments, the optimal procedure would be to set τ_yw and to calculate the young water fraction corresponding to this τ_yw. Nevertheless, establishing a constant τ_yw for all the catchments could be a tricky choice. Indeed, by setting a τ_yw higher than that obtained via Eq. (7), we are improperly using the TTD to estimate the young water fraction. From this reasoning, the only solution is to set τ_yw equal to the overall minimum τ_yw. This choice ensures that the TTD is used properly to estimate the young water fraction in all the sites. What we could expect is that all the young water fractions would be lower with respect to those obtained with the amplitude ratio approach. However, changes in young water fraction depending on changes in τ_yw are unintuitive, since they vary with the TTD shape. Accordingly, a constant reduction of τ_yw would change the area under the transit time pdf differently based on the α value. Thus, the overall effect of the reduced τ_yw upon the young water fraction remains a priori unknown.

Since in this study we are using the amplitude ratio approach, our $F_{\mathrm{yw}}^{*}$ estimates refer to the proportion of runoff younger than an inconsistent threshold age. This variation by catchment (albeit limited) is the main limitation of this work.

4.2 The role of Quaternary deposits

In line with the results of Arnoux et al. (2021), we find a negative statistically significant correlation between $F_{\mathrm{yw}}^{*}$ and WFI (${\mathit{\rho }}_{\mathrm{Spearman}}=-\mathrm{0.5}$, p value < 0.01; see Fig. S6), suggesting (unsurprisingly) that more groundwater contribution to streamflow increases the water age. WFI and F_qd values for all the study catchments are reported in Table 2. To analyze the relationship of $F_{\mathrm{yw}}^{*}$ with Quaternary deposits, we exclude the SOU and BCC catchments, since they show F_qd=0 (see Table 2). The inclusion of catchments with F_qd=0 would bias the analysis, since an absent parameter cannot modulate the share of groundwater and thus the young water fraction in the stream.

By focusing first only on the snow-dominated catchments, a linear fit on the data returns a negative slope of −0.36 (R²=0.52), indicating a reduction of $F_{\mathrm{yw}}^{*}$ with increasing F_qd (Fig. 7a). Moreover, we find a Spearman rank correlation coefficient of −0.6 with a p value of 0.13, meaning a negative but not statistically significant correlation between $F_{\mathrm{yw}}^{*}$ and F_qd. This result can be explained by considering several factors. First, water storage in Quaternary deposits is not the only groundwater storage contribution to the stream in such environments: additional storage is provided by the bedrock fractures (Gleeson et al., 2014; Jasechko et al., 2016; Martin et al., 2021), possibly caused by rock stress and high erosion rates and by the bedrock geology, which has influence on groundwater retention capacity (Hayashi, 2020). Second, the area covered by Quaternary deposits could be an insufficient proxy of the groundwater storage potential: the knowledge of the thickness of these deposits (i.e., their volume) and the bedrock topography are crucial factors for controlling groundwater storage (Arnoux et al., 2021; Hayashi, 2020), but corresponding data are not available to date.

Figure 7(a) Young water fraction against fraction of Quaternary deposits. Points dimension is proportional to τ_yw. (b) Fraction of Quaternary deposits against mean catchments elevation.

Download

$F_{\mathrm{yw}}^{*}$ values of the hybrid catchments reveal a weak positive correlation with Quaternary deposits (ρ_Spearman=0.13, p value = 0.74), while for rainfall-dominated catchments they show a negative correlation (${\mathit{\rho }}_{\mathrm{Spearman}}=-\mathrm{0.52}$, p value = 0.2); however, both correlations are not statistically significant. These weak correlations suggest that F_qd represents only a limited part of the catchment geology responsible for groundwater flow and that it can only be considered as a first-order measure of geological groundwater storage.

We furthermore observe that F_qd decreases with mean catchment elevation in our data set (Fig. 7b), revealing a negative statistically significant correlation (${\mathit{\rho }}_{\mathrm{Spearman}}=-\mathrm{0.5}$, p value < 0.01). This negative correlation reflects the fact that F_qd decreases when the mean slope increases (Arnoux et al., 2021).

To conclude, we stress that more catchments and more geological information would be required to statistically validate these observations about the role of the groundwater storage potential for explaining young water fraction variations.

4.3 The role of groundwater flow (baseflow) in $F_{\mathrm{yw}}^{*}$

4.3.1 Baseflow under different hydroclimatic regimes

The baseflow time series resulting from the baseflow separation of Duncan (2019) for six representative study catchments (two of each regime) are reported in Fig. 8 (complete results in Fig. S3). This figure shows the effect of groundwater recharge from rain and snowmelt through the “smoothed” baseflow proposed by Duncan (2019). This “smoothing” simulates a delayed storage contribution to the stream following the recharge phase during an input event. This recharge phase promotes the system wetness, thus favoring increasing quick flow. The increasing quick flow during events also leads to an increase of F_yw(t_i), as found previously (von Freyberg et al., 2018). However, the relative amount of baseflow remains high during events: the mean baseflow fraction during the high-flow period is 0.49 and 0.52 for hybrid and rainfall-dominated catchments, while it is 0.63 for snow-dominated catchments. In agreement with worldwide stable-isotope-based hydrograph separation results (Jasechko, 2019), this outcome underlines the mobilization of stored water (i.e., old water) during rainfall and snowmelt events, and this process seems to be particularly relevant in high-elevation catchments.

Figure 8Baseflow separation for six representative study catchments using the Duncan (2019) filter. The black area represents the daily discharge, while the colored area represents the estimated daily baseflow. The darker color represents a time step in which at least 85 % of the daily discharge is composed by baseflow.

Download

Accordingly, in snow-dominated systems, the snowmelt largely transits through the groundwater store (Hayashi, 2020; Cochand et al., 2019; Du et al., 2019; Flerchinger et al., 1992; Martinec, 1975), as schematized in Fig. 4, and the very high baseflow in high mountain catchments during summer is a direct sign of meltwater infiltration and percolation to groundwater that pushes old snowmelt (the main groundwater storage component) out to the stream network, as also found by Martinec (1975). This is also supported by the fact that groundwater, in such catchments, often has the isotopic signature of snowmelt (Michelon et al., 2023; Pavlovskii et al., 2018).

When examining the overall flow (and not only at the high-flow periods), F_bf is generally lower for hybrid catchments (mean of F_bf=0.67) than for rainfall-dominated (mean of F_bf=0.74) and snow-dominated catchments (mean of F_bf=0.83). The values of F_bf for all the study sites are reported in Table 2. In the BCC catchment, the F_bf (0.87) is consistent with the previous findings of Penna et al. (2016) who found, using stable water isotopes, that on average between 80 % and 98 % of the discharge in BCC is composed of pre-event water (assumed to represent groundwater). On average, the F_bf computed over the entire PoS is higher than that computed during the high-flow periods. This result suggests, unsurprisingly, that the largest percentage of base flow is released during low-flow periods. Accordingly, the variations of F_bf with elevation among different catchments (Fig. 9b) can be explained considering the changes in low-flow duration (LFD) with elevation, as will be discussed in Sect. 4.4.

Figure 9(a) Young water fraction plotted against fraction of baseflow: vertical and horizontal bars represent ± standard deviation. Gray area represents the 95 % prediction bounds of a linear regression of $F_{\mathrm{yw}}^{*}$ on F_bf. Points dimension is proportional to τ_yw. (b) Fraction of baseflow and young water fraction against mean elevation. Bars with black edge indicate F_bf (left axis), while bars with gray edge indicate $F_{\mathrm{yw}}^{*}$(right axis). Vertical bars represent ± standard deviation.

Download

Baseflow filters were already applied in previous studies, and their results were correlated with $F_{\mathrm{yw}}^{*}$. For example, von Freyberg et al. (2018) found a strong positive correlation (ρ_Spearman=0.73, p value < 0.001) between $F_{\mathrm{yw}}^{*}$ and the quick-flow index (QFI), calculated as the mean ratio between (Q−Q_bf) and Q, where Q is the daily discharge and Q_bf is the daily baseflow calculated in their paper with the Lyne and Hollick (1979) baseflow filter. By relating the F_bf to $F_{\mathrm{yw}}^{*}$, we have found a strong negative correlation (${\mathit{\rho }}_{\mathrm{Spearman}}=-\mathrm{0.73}$, p value < 0.001), as shown in Fig. 9a, consistent with the results of von Freyberg et al. (2018).

In snow-free systems, $F_{\mathrm{yw}}^{*}$ is by definition related to F_bf: baseflow is composed of groundwater and groundwater is the dominant source of old water in such systems (in absence of large lakes). In snow-influenced systems, through the “direct input” approach for estimating $F_{\mathrm{yw}}^{*},$ we consider the snowpack (i.e., a temporarily old water storage) as part of the catchment storage. However, the share of snowmelt (with age > 3 months) that flows off quickly as surface or fast subsurface runoff will not show up in F_bf. In other words, F_bf is not able to take into account all the snowmelt but only the part of meltwater that infiltrates and recharges the groundwater storage, which is a large portion of the overall snowmelt.

4.4 Low-flow duration (LFD) and $F_{\mathrm{yw}}^{*}$

The values of LFD for all the study sites are reported in Table 2. Specifically, LFD is lower for hybrid catchments (median of LFD = 0.39), and it is increasingly higher for rainfall (median of LFD = 0.50) and snow-dominated catchments (median of LFD = 0.62). In hybrid catchments, the presence of rain and snowmelt events spanning large parts of the year and the relatively low evapotranspiration (compared to rainfall-dominated catchments) due to reduced temperatures (Goulden et al., 2012) dramatically reduces the duration of low-flow periods, and this is also visible from the recurring discharge peaks (Fig. 8). In low-lying, rain-dominated catchments, evapotranspiration and precipitation are respectively higher and lower than in hybrid catchments, leading to longer low-flow periods (usually during summer and autumn). Under current climate and according to our data set, in snow-dominated catchments, we observe longer winter low-flow periods (streamflow decreasing below 0.5 to 1 mm d⁻¹ for the highest locations; see Fig. S7) on an annual scale than in hybrid catchments. To gain additional insights into the high LFD in snow-dominated catchments and the low LFD in hybrid catchments, it is necessary to further consider the role of snowpack persistence, discussed in the following section. The variations of LFD with elevation are shown in Fig. 11b.

Low-flow periods are typically baseflow dominated (or old water dominated). Accordingly, as anticipated in Sect. 4.3, the variation of F_bf between catchments reflects the proportion of the low-flow duration during the PoS. We observe that the higher the LFD is, the higher the F_bf is: in fact, they are strongly positively correlated (ρ_Spearman=0.97, p value < 0.001) as shown in Fig. 10. The negative correlation between LFD and $F_{\mathrm{yw}}^{*}$ is lower (${\mathit{\rho }}_{\mathrm{Spearman}}=-\mathrm{0.74}$, p value < 0.001; Fig. 11a) but nevertheless suggests that LFD is an important predictor for $F_{\mathrm{yw}}^{*}$.

Figure 10F_bf against the low-flow duration (LFD).

Download

Figure 11(a) $F_{\mathrm{yw}}^{*}$ against the low-flow duration, LFD. Boxplots of $F_{\mathrm{yw}}^{*}$ for catchments belonging to the same regime are plotted in correspondence to the median LFD. Points dimension is proportional to τ_yw. (b) LFD against mean elevation.

Download

Figure 12Time series of f_SCA for six representative study catchments (two for each hydroclimatic regime), illustrating the gradual increase of the F_SCA passing from rainfall-dominated to snow-dominated catchments.

Download

4.5 The role of snowpack persistence

We explore next the presence of an ephemeral or seasonal snowpack as a relevant factor for the time concentration of liquid water input and for LFD. We consider the F_SCA, calculated as reported in Sect. 3.2, as a proxy of the snowpack persistence. The F_SCA values for all the study sites are reported in Table 2. The underlying f_SCA time series for six representative study catchments (two for each hydroclimatic regime) are reported in Fig. 12 (for complete results, see Fig. S4). All the catchments characterized by a seasonal snow cover (i.e., snow dominated) reveal a high F_SCA (>0.40, median of F_SCA=0.51). Gradually smaller F_SCA values correspond to increasingly more ephemeral snowpacks with intermittent snowmelt events during the winter season (Petersky and Harpold, 2018), as reflected by the spiky f_SCA time series of hybrid and rainfall-dominated catchments (Fig. 12).

Figure 13(a) Young water fraction against F_SCA. The gray area represents the perceptual bell-shaped behavior of $F_{\mathrm{yw}}^{*}$ with increasing F_SCA. Points dimension is proportional to τ_yw. (b) F_SCA against mean elevation.

Download

Our results exhibit a bell-shaped behavior of $F_{\mathrm{yw}}^{*}$ with varying F_SCA (Fig. 13a). Specifically, we observe a general increase of $F_{\mathrm{yw}}^{*}$ for F_SCA values roughly below 0.3. This result can be explained considering that especially in hybrid catchments (median of F_SCA=0.28), but partially also in rain-dominated catchments (median of F_SCA=0.13), streamflow receives relatively more young water from ephemeral snowpacks. These short-lived snowpacks melt during the winter season resulting in only a short delay between precipitation input and melt (i.e., no water aging in the snowpack), and correspondingly meltwater flows off quickly into the stream (reducing LFD, Fig. 14), e.g., in presence of a frozen surface soil layer. In fact, ephemeral and slightly thick snowpacks do not protect the underlying soil from freezing (Harrison et al., 2021; Rey et al., 2021). Even for low-elevation locations ($<\approx \mathrm{1500}$ m a.s.l.), freezing conditions are regularly observed during winter (Keller et al., 2017).

Figure 14Low-flow duration (LFD) against F_SCA.

Download

For F_SCA values roughly higher than 0.3, we observe a decrease of $F_{\mathrm{yw}}^{*}$ with F_SCA; here all the catchments of our data set are snow dominated. The mechanisms at play here are as follows: (i) in catchments with seasonal snowpacks, streamflow receives snowmelt in spring and summer that is at least partly older than 2–3 months (because part of the snow fell more than 3 months before the melt occurs); and (ii) the building up of a persistent, deep snowpack can promote deep vertical infiltration during the main melt period, either by insulating the soil and thereby preventing/reducing freezing (Harrison et al., 2021; Rey et al., 2021; Jasechko et al., 2016) or by gradual soil thawing during the melt period (Rey et al., 2021; Scherler et al., 2010). The temporal dynamic of snow accumulation and melt supports the pivotal role of snowmelt in recharging groundwater during summer in high-elevation environments (Cochand et al., 2019; Du et al., 2019; Flerchinger et al., 1992). A similar result was also found for dolomitic catchments (such as BCC and OVA) by Lucianetti et al. (2020), who discovered that different proportions of rain and snow contribute to the recharge of springs in the Dolomites, with a gradually higher meltwater contribution in springs with increasing elevation. This role of snowmelt supports our analytical choice of computing $F_{\mathrm{yw}}^{*}$ through the “direct input” approach, thus considering the snowpack as part of the catchment storage. In addition, the potentially large shares of meltwater that recharge groundwater via deep vertical infiltration also result in old water sustaining winter baseflow (Fig. S5): the persistent snowpack and the absence of a liquid water input favor a groundwater storage release that creates a longer winter low-flow period that increases LFD (Fig. 14), thus reducing $F_{\mathrm{yw}}^{*}$, as discussed in Sect. 4.4.

F_SCA is strongly correlated with the mean catchment elevation in our data set (ρ_Spearman=0.97, p value < 0.01; Fig. 13b). A posteriori, we could have considered mean elevation instead of F_SCA as a proxy for snowpack persistence. However, a priori, it could be approximative to describe the snow cover persistence only with the increasing elevation: the persistence of snow in a catchment also depends on overall topographic and climatic characteristics, specifically relating to snow and aspect (Painter et al., 2023). In fact, catchments with very different characteristics (e.g., different elevation ranges and different areas) can reveal a similar mean elevation, but the snowpack persistence could considerably change. This is the reason why we focused on F_SCA that integrates these physical factors.

Figure 15Perceptual model of the hydrological processes that drive the young water fraction variations with elevation. This model emerges from our analysis and harmonizes these results with those of previous studies. For snow-dominated and hybrid catchments, we indicate the dominant processes, occurring during summer and during winter, that lead to low and high $F_{\mathrm{yw}}^{*}$, respectively.

Download

The above mechanisms are unable to explain the hydrological function of the glacier-dominated NBPV catchment, which has a very high $F_{\mathrm{yw}}^{*}$ and is an outlier among the snow-dominated catchments (Fig. 13a). The high $F_{\mathrm{yw}}^{*}$ of the high-elevation glacier-covered (42 %) catchment can be explained considering that the glacier melt produces high amounts of streamflow that transit the glacier system very quickly during the summer, given generally fast englacial and subglacial flow paths and the often limited water storage capacity in the glacier forefield (Müller et al., 2022; Saberi et al., 2019; Jansson et al., 2003). Schmieder et al. (2019) also found a high young water fraction in an Austrian glacier-covered (35 %) catchment, leading them to the conclusion that the basin behaves locally like a “Teflon basin” with quickly transmitted ice melt.

4.6 Process interplay along elevation: perceptual model

The identified key drivers of young water fractions for rainfall-dominated, hybrid and snow-dominated catchments can conveniently be summarized into a perceptual model of the involved hydrological processes and their seasonal interplay (Fig. 15).

High-elevation catchments are characterized by long winter low-flow periods, resulting from the build-up of a seasonal snowpack, and are sustained by the emptying of groundwater (or old water) storage. Accordingly, such storage releases stored water, mainly old meltwater, for prolonged periods where the snowpack can last for several months (typically from December to early April) before releasing water during the melting season. Such seasonal snowpack can protect the underlying soils from freezing, thus promoting meltwater infiltration and groundwater recharge. From this viewpoint, snowpack is considered as part of the catchment storage, and there is a thin line between groundwater and meltwater in snow-dominated catchments. Snowmelt or rain events push out old meltwater to the stream during summer, suggested by the relatively high amount of daily baseflow during the melting season. During this period, the high catchment wetness might even lead to saturation and thereby favor fast flow paths of meltwater or rainwater, which in turn can temporarily increase the young water fraction. Despite this increase during high-flow periods, the prevailing winter low-flow periods in such systems lead to a reduction of the average annual young water fraction.

In catchments with an ephemeral snowpack, at lower elevations, snowmelt events occur regularly during winter such that water released from the corresponding short-lived snowpack is likely younger than 3 months. Moreover, ephemeral snowpacks do not protect the underlying soils from freezing, and rapid flow paths can emerge during episodic or long-term soil surface freezing, by increasing the young water fraction. The high $F_{\mathrm{yw}}^{*}$ of such systems is also explained by the simultaneity of snowmelt and rain events during extended parts of the year (leading to large volume of annual precipitation) and the relatively low (compared to rainfall-dominated catchments) evapotranspiration. Both processes increase the catchment's wetness and reduce the low-flow period's length.

Finally, at the lowest elevations, lower amounts of precipitation and higher evapotranspiration favor longer low-flow periods, mainly sustained by old groundwater from alluvial aquifers, which lead to both a $F_{\mathrm{yw}}^{*}$ and a catchment wetness reduction. Further, the relatively flat topographies at the lowest elevations favor slow flow paths increasing the transit times of water.

How well current hydrological models can represent the interplay of these processes along elevation gradients is left for future research, but our perceptual model builds a solid basis for an improvement of theory-driven models (Clark et al., 2016).