Towards a conceptualization of the hydrological processes behind changes of young water fraction with elevation: a focus on mountainous alpine catchments

The young water fraction ($F_{\mathrm{yw}}^{*}$), defined as the fraction of catchment outflow with transit times of less than 2–3 months, is increasingly used in hydrological studies that exploit the potential of isotope tracers. The use of this new metric in catchment intercomparison studies is helpful to understand and conceptualize the relevant processes controlling catchment functioning. Previous studies have shown surprising evidence that mountainous catchments worldwide yield low $F_{\mathrm{yw}}^{*}$. These low values have been partially explained by isolated hydrological processes, including deep vertical infiltration and long groundwater flow paths. However, a thorough framework illustrating the relevant mechanisms leading to a low $F_{\mathrm{yw}}^{*}$ in mountainous catchments is missing.

The main aim of this paper is to give an overview of what drives $F_{\mathrm{yw}}^{*}$ variations according to elevation, thus clarifying why it generally decreases at high elevation. For this purpose, we assembled a data set of 27 study catchments, located in both Switzerland and Italy, for which we calculate$F_{\mathrm{yw}}^{*}$. We assume that this decrease can be explained by the groundwater storage potential, quantified by the areal extent of Quaternary deposits over a catchment (F_qd), and the low-flow duration (LFD) throughout the period of isotope sampling (PoS). In snow-dominated systems, LFD is strictly related to the snowpack persistence, quantified through the mean fractional snow cover area (F_SCA). The drivers are related to the catchment storage contribution to the stream that we quantify by applying a cutting-edge baseflow separation method to the discharge time series of the study sites and by estimating the mean baseflow fraction (F_bf) over the PoS.

Our results suggest that Quaternary deposits could play a role in modulating $F_{\mathrm{yw}}^{*}$ elevation gradients via their capacity to store groundwater, but subsequent confirmation with further, more detailed geological information is necessary. LFD indicates the proportion of PoS in which the stream is sustained and dominated by stored water coming from the catchment storage. Accordingly, our results reveal that the increase of LFD at high elevations, to a large extent driven by the persistence of winter snowpacks and the simultaneous lack of a liquid water input to the catchments, results in lower $F_{\mathrm{yw}}^{*}$. In our data set, F_bf reveals a strong complementarity with $F_{\mathrm{yw}}^{*}$, suggesting that the latter could be estimated as $F_{\mathrm{yw}}^{*}\simeq \mathrm{1}-F_{\mathrm{bf}}$ for catchments without stable water isotope measurements.

As a conclusion, we develop a perceptual model that integrates all the results of our analysis into a framework for how hydrological processes control $F_{\mathrm{yw}}^{*}$ according to elevation. This lays the foundations for an improvement of the theory-driven models.

1 Introduction

Mountainous alpine catchments are often assumed to generate high shares of rapid surface or subsurface runoff due to the presence of exposed bedrock and steep landscapes. Consequently, the role of groundwater storage in high-elevation catchments has been often neglected (Hayashi, 2020). On the contrary, multiple worldwide studies quantified a considerable groundwater input to streamflow in high mountain catchments using tracer or water balance methods (Somers and McKenzie, 2020). Several studies from the Rocky Mountains and Andes show that, on average, about 47 % of groundwater annually sustains the streamflow (Saberi et al., 2019; Somers et al., 2019; Carroll et al., 2018; Harrington et al., 2018; Cowie et al., 2017; Baraer et al., 2009, 2015; Gordon et al., 2015; Frisbee et al., 2011; Liu et al., 2004; Clow et al., 2003). Similar percentages, 49 % and 48 %, are also found in the Himalayas and the Alps, respectively (Chen et al., 2018; Engel et al., 2016; Käser and Hunkeler, 2016; Williams et al., 2016; Wilson et al., 2016; Andermann et al., 2012). It is well known that the water is stored longer than a year or a few years and that stored water plays a key role in streamflow generation processes (McDonnell, 2017; Jasechko, 2019). The study of water age has implications for predicting the timing of nutrient cycles and pollutant transport, since water age and solute dynamics are closely coupled (Li et al., 2021). Nevertheless, water age quantification is not straightforward.

Kirchner (2016a, b) proposed a new metric to quantify the share of catchment outflow with transit times lower than roughly 0.2 years or 2–3 months: the young water fraction. This metric can be conveniently inferred from the dampening effect that a catchment has on the seasonal cycle of stable water isotopes in precipitation, i.e., by estimating the ratio of the amplitudes of the seasonal cycles of stable water isotopes in streamflow and in precipitation (Kirchner, 2016a). In this method, the seasonal cycle of stream water isotope measurements is modeled using a sine wave that can be flow weighted, using the discharge measured at the moment of sampling as a weight, or not (von Freyberg et al. 2018). Isotopes measured in precipitation can be modeled with a sine function weighted according to the volume of precipitation to reduce the influence of low-precipitation periods and to account for temporally aggregated rainfall samples (von Freyberg et al., 2018). Flow-weighted fits to the seasonal tracer cycles predict the flow-weighted average young water fraction ($F_{\mathrm{yw}}^{*})$ in streamflow, while unweighted fits to the seasonal tracer cycles predict the unweighted one (F_yw) (Kirchner, 2016b). Gallart et al. (2020a) recently highlighted the advantages of the flow-weighted analysis to compensate for subsampled high-flow periods, thus reducing the underestimation of the young water fraction. Hereafter, we will use the symbol “^*” for referring to a flow-weighted variable, in order to be consistent with previous studies (von Freyberg et al., 2018; Gallart et al., 2020a).

$F_{\mathrm{yw}}^{*}$ is increasingly used in hydrological studies because it has the advantage of being free from the aggregation errors inherent to mean transit time (MTT) estimates obtained through the classical convolution approach (Kirchner, 2016a). Even more so, $F_{\mathrm{yw}}^{*}$ is an informative descriptor of catchment hydrological functions, of nutrients cycles and of pollutant transport (Stockinger et al., 2019; Benettin et al., 2017; Jasechko et al., 2016; Xia et al., 2023). For these reasons, this new metric is useful for catchment intercomparison studies to find what are the main hydroclimatic and landscape characteristics that drive the transit times of water lower than a threshold age, which varies from about 2 to 3 months. Indeed, previous work has tried to study the relationship between $F_{\mathrm{yw}}^{*}$ and catchment characteristics. von Freyberg et al. (2018) found that young water fractions of 22 Swiss catchments are significantly positively correlated with selected hydroclimatic indices and with the fraction of saturated area, suggesting that $F_{\mathrm{yw}}^{*}$ depends on catchment wetness, which promotes rapid flow paths. Interestingly, von Freyberg et al. (2018) found a statistically significant positive correlation with elevation after removing the five snow-dominated catchments, which expressed the smallest $F_{\mathrm{yw}}^{*}$ (von Freyberg et al., 2018). Likewise, Lutz et al. (2018) estimated $F_{\mathrm{yw}}^{*}$ for 24 catchments in Germany and found the smallest values for higher-elevation sites. These results are partially consistent with those of Jasechko et al. (2016), who based on the analysis of 254 watersheds worldwide discovered a reduction of $F_{\mathrm{yw}}^{*}$ in mountainous, steeper terrains. This could be related to deep vertical infiltration caused by fractures generated by high rock stress in complex terrain morphologies or by freely draining soils (i.e., cambisols and luvisols), both associated with high-elevation environments (Lutz et al., 2018; Jasechko et al., 2016; Gleeson et al., 2014). In addition, topographic roughness increases flow path and, correspondingly, transit time (Gleeson and Manning, 2008; Frisbee et al., 2011; Jasechko et al., 2016). Despite these studies, there is still a lack of a unified framework of how the variation among mountainous catchments results in less young water at high elevation.

An early example from the Swiss Alps showed that high celerity originates from massive meltwater infiltration that pushes out groundwater reserves: streamflow following snowmelt is older than meltwater infiltrated in the current year (Martinec, 1975). The resulting effect on water partitioning between the surface and the subsurface should be analyzed considering the temporal concentration of water input during the snowmelt period, but this remains largely unexplored (Rey et al., 2021). Despite this lack of studies on water partitioning during snowmelt, several studies have demonstrated the pivotal role of snowmelt in recharging groundwater during summer in high-elevation environments (Hayashi, 2020; Cochand et al., 2019; Du et al., 2019; Flerchinger et al., 1992).

From a water modeling perspective and thus from a water age perspective, snowpack storage and groundwater storage can be considered a single entity: they both constitute catchment storage. Therefore, the analytical estimation of $F_{\mathrm{yw}}^{*}$ must reflect this “conceptual” decision of whether to consider the snowpack storage to be part of the catchment storage. This point has been previously addressed by von Freyberg et al. (2018). If total precipitation is considered as catchment input (direct input case), the snowpack is implicitly considered to be part of the catchment storage, and $F_{\mathrm{yw}}^{*}$ results from the combination of snowpack and subsurface storage. In this direct input case, $F_{\mathrm{yw}}^{*}$ is computed from the amplitudes of the seasonal cycles of stable isotopes of water in precipitation (A_P) and streamflow ($A_{\mathrm{s}}^{*}$). If total liquid water input (composed of rainfall and snowmelt before sometimes called equivalent precipitation) is considered to be catchment input, $F_{\mathrm{yw}}^{*}$ is computed based on the amplitudes of the cycles in equivalent precipitation (A_Peq) and in streamflow ($A_{\mathrm{s}}^{*}$). This $F_{\mathrm{yw}}^{*}$ value then results from subsurface storage alone, since snowpack storage is excluded from the catchment storage (von Freyberg et al., 2018). If $F_{\mathrm{yw}}^{*}$ is estimated using a direct input setting (i.e., total precipitation directly), $F_{\mathrm{yw}}^{*}$ is expected to be smaller, since the catchment storage is larger (von Freyberg et al., 2018). Also, Ceperley et al. (2020) investigated the role of water input from snow in $F_{\mathrm{yw}}^{*}$ estimation, concluding that the low values in high alpine snow-dominated catchments result from a combination of snow cover effects and the storage in the subsurface. In the present work, the main aim is not to address how the snowpack affects $F_{\mathrm{yw}}^{*}$ estimation in a single catchment, as this was treated previously (von Freyberg et al., 2018; Ceperley et al., 2020), but to investigate the hydrological processes (also related to the snowpack storage) that lead to variations in $F_{\mathrm{yw}}^{*}$ between catchments located at different elevations with a focus on high-elevation alpine catchments.

Some authors have revealed the possibility of Quaternary deposits (e.g., talus, moraine and alluvium) to store groundwater in high-elevation alpine catchments (Arnoux et al., 2021; Hayashi, 2020; Christensen et al., 2020). The stored water in these deposits can in fact sustain streamflow during the low-flow period (Hayashi, 2020; Arnoux et al., 2021), as supported by the strong positive correlation found by Arnoux et al. (2021) between the fraction of Quaternary deposits and the winter flow index (a low-flow indicator reflecting the groundwater contribution to the stream) for 13 alpine catchments. During winter, the period without liquid water input can last 6 months or more in high-elevation catchments. The occurrence of such long periods of low flows hints towards important amounts of stored water (or old water) that are well connected to the stream network, and it thereby remains accessible throughout the low-flow period (Somers et al., 2019).

To further discuss the role of low flow in $F_{\mathrm{yw}}^{*}$ estimation, let us first consider that $F_{\mathrm{yw}}^{*}$ can be theoretically estimated based on the flow-weighted average of young water fractions (Kirchner, 2016b):

$$\begin{array}{}\text{(1)}& F_{\mathrm{yw}}^{*}={\displaystyle \frac{A_{\mathrm{S}}^{*}}{A_{\mathrm{P}}}}\simeq {\displaystyle \frac{\sum _{i=\mathrm{1}}^{n}Q\left(t_i\right)F_{\mathrm{yw}}\left(t_i\right)}{\sum _{i=\mathrm{1}}^{n}Q\left(t_i\right)}},\end{array}$$

where n is the number of time steps (e.g., days) in the period of isotope sampling, PoS, Q(t_i) is the discharge at the time t_i (e.g., daily discharge) and F_yw(t_i) is the young water fraction at the time t_i (e.g., daily young water fraction). As is clear from Eq. (1), $F_{\mathrm{yw}}^{*}$ becomes low if either F_yw(t_i) is low for high flows or if F_yw(t_i) is very low for many time steps or both. The low-flow periods correspond to the recession periods in which there is no new rainfall or meltwater input in the catchments. Thus, during these periods, the catchment storage releases stored water (or old water) to the stream sustaining the streamflow (Hayashi, 2020). Thus, we can anticipate that low Q(t_i) values imply low F_yw(t_i) values. Accordingly, the proportion of the low-flow period during a specified time window should reduce the amount of young water reaching the stream during that time window. Nevertheless, the F_yw(t_i) is higher during high-flow (wet) periods (von Freyberg et al., 2018; Wilusz et al., 2017; Gallart et al., 2020b). Thus, the overall effect of the proportion of low-flow and high-flow periods upon $F_{\mathrm{yw}}^{*}$ remains a priori unclear. It is however tempting to think that the duration of low-flow period or the share of baseflow could explain $F_{\mathrm{yw}}^{*}$ variations at different elevations (since both low-flow duration and the share of baseflow change with elevation). In addition, in high-elevation, snow-dominated catchments, the persistence of the snowpack is the main driver of the long low-flow duration, since the low-flow period at high elevation corresponds to the presence of the seasonal snowpack (corresponding to an absence of liquid water input), while the high-flow period is generally snowmelt driven. Such snowmelt generally occurs in late spring or summer, and it is likely to be older than 2–3 months (because the peak snow fall occurred 3 months earlier). As a result, summer discharge mainly consists of old water, either of current snowmelt that reaches the stream via faster surface or subsurface flow paths, or old snowmelt (main component of groundwater storage) pushed out in the stream by infiltrated rainwater or meltwater. Part of the snowpack can release young water, but this is a minor component in catchments with a seasonal snowpack. In contrast, in catchments with an ephemeral snowpack, it is common to observe intermittent winter snowmelt that is likely younger than 2–3 months: snowmelt is temporally close to snowfall. In this case, streamflow receives relatively more young water from short-lived snowpacks. However, it is still unclear if seasonal or ephemeral snow cover dynamics can affect the $F_{\mathrm{yw}}^{*}$ (Ceperley et al., 2020).

An innovative focus of our work is on variables that were not previously considered for explaining elevation gradients of young water fractions. We specifically exclude catchment size, annual precipitation, bedrock porosity, pasture cover and open water cover that have been discussed and shown to have little correlation in the work of Jasechko et al. (2016).

A special case in terms of explanatory variables is mean annual precipitation: Jasechko et al. (2016) did not observe any significant correlation between the $F_{\mathrm{yw}}^{*}$ and annual precipitation in their worldwide study. Lutz et al. (2018) found that $F_{\mathrm{yw}}^{*}$ decreases with increasing mean annual precipitation, based on 24 catchments in Germany. In contrast, in the relatively wet rainfall-dominated and hybrid catchments studied by von Freyberg et al. (2018), $F_{\mathrm{yw}}^{*}$ was shown to increase with mean monthly precipitation and correspondingly also with elevation. In their study, discharge (unsurprisingly correlated with precipitation) was considered as a proxy of catchment wetness, which favors rapid flow paths and thereby increases $F_{\mathrm{yw}}^{*}$ (von Freyberg et al., 2018). In snow-dominated systems, the use of mean annual precipitation as a proxy for catchment wetness could be misleading because the seasonal snowpack leads to a very dry period of the year despite the high solid water input. In other words, the temporal concentration of the liquid water input is the relevant variable. Indeed, the saturation of the system (i.e., high wetness conditions) can be observed also when the annual precipitation is low if a large volume of water (stored in the snowpack) is released in a relatively concentrated time interval. Indeed, despite the fact that precipitation and, correspondingly, discharge are higher in snow-dominated than in rainfall-dominated catchments, $F_{\mathrm{yw}}^{*}$ is generally lower in snow-dominated systems that are potentially wetter than rainfall-dominated ones. This suggests that the precipitation can only partially explain the variations of $F_{\mathrm{yw}}^{*}$ and that other variables should be put under observation.

Accordingly, we omit here total annual precipitation as an explanatory variable of low $F_{\mathrm{yw}}^{*}$ in snow-dominated catchments (but we consider precipitation for rainfall-dominated and hybrid catchments) and study a new set of hydrological variables to gain new insights into $F_{\mathrm{yw}}^{*}$ along elevation gradients: the fraction of Quaternary deposits (F_qd), the mean fraction of baseflow (F_bf), the low-flow duration (LFD) and the mean fractional snow cover area (F_SCA), defined in detail in Sect. 3.2 and 3.3. We first describe the data set (Sect. 2). Then, we present the $F_{\mathrm{yw}}^{*}$ estimation method (Sect. 3.1) followed by the correlation analysis of the selected hydrological variables with the estimated $F_{\mathrm{yw}}^{*}$, and we bring these results back into the ongoing scientific discussion of $F_{\mathrm{yw}}^{*}$ (Sect. 4.2–4.6).

2 Study sites

We analyze 27 study catchments located both in Switzerland and Italy integrating observations from multiple published data sets (25 catchments) with new additional observations (2 catchments) (Fig. 1). Geomorphological and hydroclimatic characteristics of the study sites are reported in Table 1.

Table 1 Catchment geomorphological and hydroclimatic characteristics. The catchment area and mean slope are directly calculated in Google Earth Engine. For the slope calculation, we use the Shuttle Radar Topography Mission (SRTM) DEM (Farr et al., 2007). We obtained mean elevation and precipitation information of the existing data set directly from published papers (von Freyberg et al., 2018; Ceperley et al., 2020). Discharge (Q), precipitation (P) and isotopic composition (δ¹⁸O) data are all referred to the period of sampling (PoS) indicated in this table. The letter in brackets in the first column indicates the hydroclimatic regime: (R) is rainfall dominated, (H) is hybrid and (S) is snow dominated.

Figure 1 Location of the 27 study catchments with indication of the hydroclimatic regime.

Specifically, we assembled the 22 Swiss catchments studied by von Freyberg et al. (2018) with the three alpine catchments investigated by Ceperley et al. (2020) (Vallon de Nant, Noce Bianco at Pian Venezia and Bridge Creek catchment) into a single data set. Hereafter, we refer to these catchments with the ID reported in the above-mentioned published papers (Table 1). We also consider two additional high-elevation catchments located near the Nivolet Pass (Valsavaranche, Aosta Valley, Italy) (Gisolo et al., 2022). In this alpine environment, we monitor the mainstream, called “Dora del Nivolet”, and a secondary river called “Source”. Hereafter we refer to these catchments with the IDs DOR and SOU, respectively. A detailed description of the DOR and SOU catchments is reported in the Supplement.

The von Freyberg et al. (2018) data set includes catchments with areas between 0.7 and 351 km² and mean elevations between 472 and 2369 m a.s.l. With the five catchments added here, the complete data set includes catchment areas between 0.14 and 359 km² and spans mean elevation between 472 and 3049 m a.s.l. The mean monthly precipitation ranges between 61.3 and 168.7 mm per month, while mean discharge ranges between 28.6 and 138.9 mm per month. The mean slope ranges from 4 to 34^∘, and our study sites reveal an increase of steepness with elevation (Fig. 2a). Precipitation increases with elevation until 1500 m a.s.l., above which it decreases (Fig. 2c), highlighting a change of precipitation regime as described by previous studies (Santos et al., 2018). The five catchments added to the initial data set of von Freyberg et al. (2018) allow the analysis to explore the high-elevation regions (mean elevation > 1500 m a.s.l.) that were previously poorly represented. Most of the study catchments reveal a sedimentary bedrock but dolomitic and metamorphic bedrocks, characteristic of high-elevation sites, are also included in our data set. Moreover, the presence of unconsolidated Quaternary deposits is widespread among our study catchments: only two catchments (BCC and SOU) do not reveal this type of geology. The complete data set now explores case studies from the Swiss plateau and pre-alpine area; from the Jura; and from five different alpine regions, including the northern part of the Swiss Alps, the southern Swiss Alps (Alpi Ticinesi), the western Italian Alps (Alpi Graie), the Rätische Alps and the Dolomites. Overall, this represents a good range of geologies as well as of climatic conditions.

Figure 2 (a) Mean slope against mean elevation. Vertical bars represent the mean slope standard deviation, horizontal bars represent min–max elevation range, (b) boxplot of the mean slope values, (c) mean precipitation and discharge against elevation, and (d) boxplots of the mean precipitation and discharge values. Here and later: the boxplots show the median and the interquartile range (IQR), the whiskers are defined as the IQR multiplied by 1.5 and outliers are plotted with red “+” markers.

Figure 3 Boxplots of mean monthly flow for all the study catchments grouped according to their flow regime (rainfall dominated, hybrid, snow dominated). Colored areas represent the monthly flow of each study catchment belonging to the relative regime.

In order to be consistent with previous studies (von Freyberg et al., 2018; Staudinger et al., 2017), we classify the 23 Swiss catchments according to the hydroclimatic regimes proposed by Staudinger et al. (2017) which group the regimes defined by Weingartner and Aschwanden (1992) in three categories: rainfall dominated (R), hybrid (H) and snow dominated (S). For the four Italian catchments, where the aforementioned classification schemes cannot be rigorously applied, we use that proposed by Stoelzle et al. (2020). This classification scheme is based on mean and maximum catchment elevation, periods of typical low flow, snow onset and beginning of snowmelt and was already used by Stoelzle et al. (2020) to classify catchments outside the Swiss borders (e.g., German catchments). According to this classification scheme, the four Italian catchments (DOR, SOU, BCC and NBPV) are all categorized as snow dominated (S). The classification of BCC is also consistent with the one given in a previous study without considering the application of a formal classification scheme (Penna et al., 2016). Across the three considered streamflow regimes, a shift of the monthly hydrograph peak (computed using discharge data in the PoS) from winter to summer months is observed (Fig. 3): this flow peak shifting is a clear sign of the increasing predominance of snowmelt in the streamflow generation processes. Our data set includes NBPV, whose area is 42 % glacier covered and consequently exhibits a characteristic glacier-dominated streamflow regime with a monthly peak in late summer (Zuecco et al., 2019; Carturan, 2016). NBPV has been classified as snow dominated following the Stoelzle et al. (2020) classification scheme. Nevertheless, its characteristics suggest it may belong to a fourth category of glacier-dominated catchments. Unfortunately, this category has not been considered by the aforementioned classification scheme, and the definition of the classifiers for a new category is outside the scope of this work. In this catchment, the effect of glacier melt on $F_{\mathrm{yw}}^{*}$ cannot be neglected, and this was partially discussed by Ceperley et al. (2020). In our data set, also the Dischmabach (DIS) and the Vallon de Nant (VdN) catchments are 2 % and 3 % glacier covered, but we assume that the effect on $F_{\mathrm{yw}}^{*}$ is negligible when compared with that of the seasonal snowpack.

3 Material and methods

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

$$\begin{array}{}\text{(2)}& F_{\mathrm{yw}}^{\mathrm{T}}=P\left(\mathit{\tau }<{\mathit{\tau }}_{\mathrm{yw}}\right)=\underset{\mathrm{0}}{\overset{{\mathit{\tau }}_{\mathrm{yw}}}{\int }}{\displaystyle \frac{{\mathit{\tau }}^{\mathit{\alpha }-\mathrm{1}}}{{\mathit{\beta }}^{\mathit{\alpha }}\mathrm{\Gamma }\left(\mathit{\alpha }\right)}}{e}^{-\frac{\mathit{\tau }}{\mathit{\beta }}}\mathrm{d}\mathit{\tau },\end{array}$$

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):

$$\begin{array}{}\text{(3a)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{S}}\left(t\right)=A_{\mathrm{S}}\sin \left(\mathrm{2}\mathit{\pi }ft-{\mathit{\phi }}_{\mathrm{S}}\right)+k_{\mathrm{S}},\text{(3b)}& {\displaystyle }& {\displaystyle }{\mathit{\delta }}^{\mathrm{18}}{\mathrm{O}}_{\mathrm{P}}\left(t\right)=A_{\mathrm{P}}\sin \left(\mathrm{2}\mathit{\pi }ft-{\mathit{\phi }}_{\mathrm{P}}\right)+k_{\mathrm{P}},\end{array}$$

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”:

$$\begin{array}{}\text{(4a)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{yw}}={\displaystyle \frac{A_{\mathrm{S}}}{A_{\mathrm{P}}}},\text{(4b)}& {\displaystyle }& {\displaystyle }{F}_{\mathrm{yw}}^{*}={\displaystyle \frac{A_{\mathrm{S}}^{*}}{A_{\mathrm{P}}}}.\end{array}$$

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):

$$\begin{array}{}\text{(5)}& {\mathrm{SE}}_{F_{\mathrm{yw}}^{*}}={\displaystyle \frac{\mathrm{1}}{A_{\mathrm{P}}}}\sqrt{{\mathrm{SE}}_{A_{\mathrm{S}}^{*}}^{\mathrm{2}}+{\left({\displaystyle \frac{A_{\mathrm{S}}^{*}}{A_{\mathrm{P}}}}\right)}^{\mathrm{2}}\cdot {\mathrm{SE}}_{A_{\mathrm{P}}}^{\mathrm{2}}},\end{array}$$

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 4 Schematic 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.

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):

$$\begin{array}{}\text{(6)}& {\mathit{\phi }}_{\mathrm{S}}^{*}-{\mathit{\phi }}_{\mathrm{P}}=\mathit{\alpha }\arctan \left(\sqrt{{\left(A_{\mathrm{S}}^{*}/A_{\mathrm{P}}\right)}^{-\frac{\mathrm{2}}{\mathit{\alpha }}}-\mathrm{1}}\right).\end{array}$$

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):

$$\begin{array}{}\text{(7)}& {\displaystyle \frac{{\mathit{\tau }}_{\mathrm{yw}}}{T}}\approx \mathrm{0.0949}+\mathrm{0.1065}\mathit{\alpha }-\mathrm{0.0126}{\mathit{\alpha }}^{\mathrm{2}},\end{array}$$

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):

$$\begin{array}{}\text{(8)}& \mathrm{NDSI}={\displaystyle \frac{r_{\mathrm{green}}-r_{\mathrm{SWIR}}}{r_{\mathrm{green}}+r_{\mathrm{SWIR}}}},\end{array}$$

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):

$$\begin{array}{}\text{(9)}& f_{\mathrm{SCA}}={\displaystyle \frac{N_{\mathrm{snow}}}{N_{\mathrm{tot}}-N_{\mathrm{clouds}}}},\end{array}$$

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):

$$\begin{array}{}\text{(10)}& F_{\mathrm{qd}}={\displaystyle \frac{A_{\mathrm{qd}}}{A}}.\end{array}$$

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

$$\begin{array}{}\text{(11)}& \mathrm{WFI}={\displaystyle \frac{Q_{\mathrm{NM}\mathrm{7}}}{Q_{\mathrm{mean}}}},\end{array}$$

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):

$$\begin{array}{}\text{(12a)}& {\displaystyle }& {\displaystyle }M\left(t_{i-\mathrm{1}}\right)={\displaystyle \frac{M\left(t_i\right)-c}{k}}+c\text{(12b)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{q}}\left(t_i\right)=k{Q}_{\mathrm{q}}\left(t_{i-\mathrm{1}}\right)+\left(M\left(t_i\right)-M\left(t_{i-\mathrm{1}}\right)\right){\displaystyle \frac{\mathrm{1}+k}{\mathrm{2}}}\text{(12c)}& {\displaystyle }& {\displaystyle }{Q}_{\mathrm{bf}}\left(t_i\right)=M\left(t_i\right)-Q_{\mathrm{q}}\left(t_i\right),\end{array}$$

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):

$$\begin{array}{}\text{(13)}& F_{\mathrm{bf}}={\displaystyle \frac{\mathrm{1}}{n}}\sum _{i=\mathrm{1}}^n{\displaystyle \frac{Q_{\mathrm{bf}}\left(t_i\right)}{Q\left(t_i\right)}},\end{array}$$

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:

$$\begin{array}{}\text{(14)}& T_{\mathrm{Low}}=\forall t_i:{\displaystyle \frac{Q_{\mathrm{bf}}\left(t_i\right)}{Q\left(t_i\right)}}\ge \mathrm{0.85}.\end{array}$$

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 5 Sinusoidal 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.

4 Results and discussion – towards a harmonious and exhaustive framework of the hydrological processes that drive the young water fraction variations with elevation

We present and discuss hereafter the $F_{\mathrm{yw}}^{*}$ and τ_yw estimates (Sect. 4.1) and the identified relations between $F_{\mathrm{yw}}^{*}$ and the studied explanatory variables (Sect. 4.2–4.5), followed by the perceptual model that describes the main processes driving the $F_{\mathrm{yw}}^{*}$ variations with mean catchment elevation and that harmonizes our results with previous studies (Sect. 4.6).

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 2 Summary table with all the relevant quantities estimated for the 27 study catchments.

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 α.

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.

$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 8 Baseflow 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.

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.

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.3.2 The complementarity between the fraction of baseflow (F_bf) and the young water fraction ($F_{\mathrm{yw}}^{*}$)

A by-product of this work is that the F_bf, estimated with the Duncan (2019) baseflow filter, is roughly the complementary term of $F_{\mathrm{yw}}^{*}$ (Fig. 9a and b), which is an important result for catchments where isotope measurements are missing. In such catchments, $F_{\mathrm{yw}}^{*}$ could potentially be estimated without the application of the amplitude ratio approach as follows:

$$\begin{array}{}\text{(15)}& F_{\mathrm{yw}}^{*}\simeq \mathrm{1}-F_{\mathrm{bf}}.\end{array}$$

Some of our case studies show considerable “residuals” of $\mathrm{1}-(F_{\mathrm{bf}}+F_{\mathrm{yw}}^{*})$ (Fig. 9b). This is partially due to the uncertainty of the parameters used for estimating F_bf. In this regard, Duncan (2019) suggests some calibration guidelines to obtain optimal parameters set for baseflow estimation per catchment. In this work, we did not use the calibration guidelines, but we simply used the recession parameter proposed by Nathan and McMahon (1990) in order to achieve factual and reproducible results. In addition, the estimation of baseflow during an event is generally less rigorous than during the recession phase (Duncan, 2019), affecting the F_bf estimation. Moreover, $F_{\mathrm{yw}}^{*}$ values are influenced by the sampling rate: the higher the frequency of sampling is, the higher the young water fraction is (Gallart et al., 2020a; Stockinger et al., 2016). Thus, the young water fraction calculated with the amplitude ratio approach generally underestimates the “theoretical” young water fraction, and we simply compensate by computing the flow-weighted young water fraction ($F_{\mathrm{yw}}^{*}$). In hybrid and snow-dominated catchments, these “residuals” can also be explained by considering that the F_bf does not include surface runoff or fast lateral subsurface flow of meltwater, likely older than the estimated threshold ages, following a snowmelt event. On the other hand, these residuals might also be related to the non-linear recession behavior of catchments, which was shown by Santos et al. (2018) to dominate Swiss low-elevation (i.e., rain dominated) catchments, when the exponential recession assumption of the baseflow filter necessarily leads to less reliable results (Duncan, 2019).

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 10 F_bf against the low-flow duration (LFD).

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.

Figure 12 Time 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.

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.

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 14 Low-flow duration (LFD) against F_SCA.

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 15 Perceptual 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.

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).

5 Conclusion

This study proposes a conceptualization of the processes behind changes in young water fraction ($F_{\mathrm{yw}}^{*}$) with elevation, defined here following Kirchner (2016a) as stream water that is younger than a threshold age of about 2–3 months. The analysis is focused on amplitude-ratio-based young water fractions for a set of 27 study catchments located in Switzerland and Italy, which span a wide range of geological and hydroclimatic conditions. The young water fraction estimates ($F_{\mathrm{yw}}^{*}$) obtained from the phase and amplitude information of seasonal isotope cycles correspond to different young water threshold ages for the different catchments, which is a limitation of this work. However, the threshold ages vary only modestly from about 1.5 to 3 months.

Our analysis focuses on mountainous catchments to fill the knowledge gap, referring to the surprisingly low young water fractions at high elevations (>1500 m a.s.l.), but we have also considered catchments at lower elevations to obtain a complete picture of the dominant hydrological processes at different elevations.

We have focused on variables and processes that were not previously considered for explaining elevation gradients of young water fraction. We have investigated the role of (i) groundwater storage potential, (ii) catchment storage contribution to the stream, (iii) low-flow duration and (iv) snowpack persistence. Our results suggest that (ii), (iii) and (iv) are interconnected: low-flow periods are generally sustained by old water deriving from the catchment storage, and the length of such periods is driven by the snowpack persistence at high elevations. The proportion of low-flow periods during the period of isotope sampling strongly influences the amount of old water contributing to the stream, thus reducing the estimated $F_{\mathrm{yw}}^{*}$. Consequently, the low-flow duration, which varies with elevation, can be retained as a driver of the $F_{\mathrm{yw}}^{*}$ changes with elevation. Given the importance of low-flow periods, we have also investigated the role of groundwater storage potential, represented here by the portion of catchment area covered by Quaternary deposits. Our results suggest that an exhaustive description of the groundwater storage potential should be completed with more detailed geological information, e.g., the geology and topography of bedrock, the fraction of fractured bedrock and the deposits' thickness, which is challenging to retrieve from a geological data set. We have brought together the results of this analysis in a perceptual model that describes a framework for how hydrological processes control the $F_{\mathrm{yw}}^{*}$ according to elevation, laying the foundations for an improvement of theory-driven models.

The strong complementarity between $F_{\mathrm{yw}}^{*}$ and the mean fraction of baseflow obtained for our data set suggests that $F_{\mathrm{yw}}^{*}$ could be estimated starting from automated baseflow separation techniques for catchments in which stable water isotope measurements are not available. This complementarity should however be validated in future work, by considering, e.g., alternative baseflow separation techniques and different hydroclimatic conditions.

Finally, the conceptualization of the hydrological processes described in this paper do not fit the high young water fraction of the single glacier-dominated catchment of our data set. In conclusion, we encourage future studies to compare and to collect isotopic data from glacier-dominated catchments to better understand the processes in such systems that, under glacier retreat due to climate change, will see a gradual transition to purely snow-dominated systems.