‘Teﬂon Basin’ or Not? A High-Elevation Catchment Transit Time Modeling Approach

: We determined the streamﬂow transit time and the subsurface water storage volume in the glacierized high-elevation catchment of the Rofenache (Oetztal Alps, Austria) with the lumped parameter transit time model TRANSEP. Therefore we enhanced the surface energy-balance model ESCIMO to simulate the ice melt, snowmelt and rain input to the catchment and associated δ 18 O values for 100 m elevation bands. We then optimized TRANSEP with streamﬂow volume and δ 18 O for a four-year period with input data from the modiﬁed version of ESCIMO at a daily resolution. The median of the 100 best TRANSEP runs revealed a catchment mean transit time of 9.5 years and a mobile storage of 13,846 mm. The interquartile ranges of the best 100 runs were large for both, the mean transit time (8.2–10.5 years) and the mobile storage (11,975–15,382 mm). The young water fraction estimated with the sinusoidal amplitude ratio of input and output δ 18 O values and delayed input of snow and ice melt was 47%. Our results indicate that streamﬂow is dominated by the release of water younger than 56 days. However, tracers also revealed a large water volume in the subsurface with a long transit time resulting to a strongly delayed exchange with streamﬂow and hence also to a certain portion of relatively old water: The median of the best 100 TRANSEP runs for streamﬂow fraction older than ﬁve years is 28%.


Introduction
High-elevation catchments were often seen as 'Teflon'-like basins [1], meaning that water input to the catchment immediately runs off as overland flow or fast subsurface flow with very limited subsurface interaction [2] and minimal infiltration, due to the permeability of bedrock typically assumed to be negligible [3]. The concept of how catchments retain water is becoming almost equally recognized as how catchments release water [4][5][6][7]. However, very recent studies mentioned the formerly unrecognized subsurface storage potential of mountain catchments [4,8,9], i.e., water is not only temporarily stored as snow and ice, but also in the soil, fractured bedrock, moraines, talus, alluvium, alluvial fans, permafrost, rock glaciers and rock slides [4,[10][11][12][13]. This water contributes to streamflow by shallow to deep flow paths [3,14], both in periods with rain, snowmelt and ice melt input and in periods without water input to the catchment (i.e., when the catchment is in a frozen state [8]), and hence is an important contributor to streamflow. This subsurface water may play an important role in the alpine hydrologic cycle as shown by numerous studies in glacier-free and glacierized catchments [12,[15][16][17][18][19][20]. Snowmelt, ice melt and rain recharge the subsurface storage, with snowmelt being the most important contributor in many mountain catchments worldwide [21,22]. The water volume and flow (the quotient is the transit time) thereby affect both, downstream water supply and quality [4,23].
Streamflow age is a crucial descriptor of catchment functioning, affecting the control of runoff generation, biogeochemical cycling and contaminant transport. It is defined as the elapsed time for water transmitted through the catchment from input to the stream at the outlet [23,24]. The transit time distribution (TTD) describes the age composition of a streamflow water sample, and can be inferred from seasonal tracer cycles with the convolution integral method by relating past tracer input (e.g., precipitation) to tracer output (e.g., streamflow) [23]. The mean transit time (MTT) thereby describes the average travel time for a water parcel from entering to leaving a catchment [23]. Dynamic storage controls streamflow dynamics (i.e., the runoff response), whereas mobile storage, also connected to streamflow, shows large differences in water age and controls transport in a catchment (i.e., the tracer response) [4]. An overview of perceptual storage terms in catchment hydrology and respective estimation methods can be found in Staudinger et al. [4]. The young water fraction (F yw ) is the proportion of streamflow below a certain threshold age (typically less than 3 months) and can be estimated by comparing fitted seasonal sine wave tracer cycles of catchment input and output [25,26].
Streamflow age (including F yw ) and subsurface storage studies were already conducted in catchments in a variety of environmental settings: in low mountain catchments, e.g., [27,28], in northern high-latitude catchments, e.g., [29,30], in high mountain catchments, whereby only a few were covered by a small percentage of glacier ice (<5%), e.g., [4,31,32], or in permafrost regions, e.g., [33]. Seeger and Weiler [32] combined a surface energy-balance model with a lumped parameter transit time model and reported MTTs of 0.6 to 10.5 years for five snow-dominated catchments in Switzerland. Staudinger et al. [4] revealed dynamic and mobile storages of approximately 100 to 500 mm and 1000 to 13,000 mm for four snow-dominated catchments in Switzerland (these were included in the analyses of [32]), respectively. Freyberg et al. [31] estimated a flow-weighted F yw of 0.1 to 0.22 for the same Swiss catchments analyzed in [32]. A global study of 254 river systems revealed a flow-weighted mean F yw of 0.34 with a threshold age of 2.3 ± 0.8 months, whereas almost all of the mountain rivers had a lower F yw than the mean [11]. Freyberg et al. [31] tested the delayed effect of snowmelt during their application of the F yw approach and found a significant difference in the amplitude of fitted seasonal tracer cycle of input water for only one out of 21 catchments compared to the method with a direct input of precipitation. They concluded that in their study catchments, F yw values were not sensitive to the use of delayed input. Tetzlaff et al. [29] extended the traditional transit time estimation approach with a snow model that provides snow melt volume and isotope ratio and applied it in an Arctic catchment. They concluded that considering all available inflows (i.e., snow melt, soil ice thaw and rainfall) led to the best estimated transit time. Hence we hypothesize that the spatio-temporal variability of the water input in glacierized high-elevation catchments (i.e., especially the snow and ice melt volume and their isotope ratios) is important, which challenges the traditional transit time and F yw approach (where rain is considered to contribute dominantly).
It can be concluded that studies addressing the F yw , catchment TTD or storage estimates in catchments with >5% glacierized area are missing and estimates of the above-mentioned metrics are unknown, but high-elevation catchment storage (not snow and ice) may become even more important in terms of a changing climate and shrinking glaciers. The question arises how the occurrence of ice melt as a contributor and significant part of the water cycle affects the estimation of F yw , catchment TTD and storage in glacierized high-elevation catchments. We make use of a model coupling approach in which we combine a surface energy-balance model to simulate the snow and ice melt volume and isotope ratio with a lumped parameter transit time model to estimate the catchment TTD and storage. Overall, we aim to identify the role of subsurface water (and storage) in a glacierized high-elevation catchment with the use of δ 18 O. Specifically, the objectives are: 1.
The estimation of catchment transit time and mobile storage by coupling of a surface energy-balance snow and ice melt model with a lumped parameter transit time model, 2.
The estimation of the F yw with delayed input of snow and ice melt using the sine wave approach and, 3.
The comparison of the TTD and the F yw .

Study Area
We conducted our study in the 35% glacierized high-elevation catchment Rofental (Figure 1), a LTSER site in the Austrian Alps (https://www.lter-austria.at/rofental/). The catchment (98 km 2 ) ranges from approximately 1890 to 3770 m a.s.l. (mean elevation is 2930 m a.s.l.). The Rofenache is draining the catchment and has a distinct glacial regime with a mean flow of 4.6 m 3 s −1  at the outlet in Vent (1891 m a.s.l., 46.85722 • N, 10.91083 • E). Mean 7-day annual minimum flow, mean flow and mean annual peak flow for the period 2014 to 2017 are 0.43, 3.93 and 24.59 mm d −1 , respectively.
Detailed site descriptions, including various climatic and hydrologic characteristics, can be found in Schmieder et al. [34] or in Strasser et al. [35]. Several cryospheric, hydrologic and geodetic studies have been conducted in the catchment and data sets from more than 150 years ago up to very recently are freely available [35]. Hydrologic tracer studies have been conducted intensively during the International Hydrological Decade (1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974) [36] and since 2014 [37]. Strasser et al. [35] provide a detailed overview about the research history in the Rofental. Figure 1 shows the catchment, as well as the water sampling and 'Austrian Network of Isotopes in Precipitation and Surface Waters' (ANIP) sites used in this publication.

Study Area
We conducted our study in the 35% glacierized high-elevation catchment Rofental (Figure 1), a LTSER site in the Austrian Alps (https://www.lter-austria.at/rofental/). The catchment (98 km 2 ) ranges from approximately 1890 to 3770 m a.s.l. (mean elevation is 2930 m a.s.l.). The Rofenache is draining the catchment and has a distinct glacial regime with a mean flow of 4.6 m 3 s −1  at the outlet in Vent (1891 m a.s.l., 46.85722° N, 10.91083° E). Mean 7-day annual minimum flow, mean flow and mean annual peak flow for the period 2014 to 2017 are 0.43, 3.93 and 24.59 mm d −1 , respectively.
Detailed site descriptions, including various climatic and hydrologic characteristics, can be found in Schmieder et al. [34] or in Strasser et al. [35]. Several cryospheric, hydrologic and geodetic studies have been conducted in the catchment and data sets from more than 150 years ago up to very recently are freely available [35]. Hydrologic tracer studies have been conducted intensively during the International Hydrological Decade (1965)(1966)(1967)(1968)(1969)(1970)(1971)(1972)(1973)(1974) [36] and since 2014 [37]. Strasser et al. [35] provide a detailed overview about the research history in the Rofental. Figure 1 shows the catchment, as well as the water sampling and 'Austrian Network of Isotopes in Precipitation and Surface Waters' (ANIP) sites used in this publication.

Data
We used daily runoff data (streamflow was measured at the gauging station in Vent) within this study. Hourly meteorological variables (air temperature, precipitation amount, relative humidity, wind speed and global radiation) were measured at various stations (see [38,39] for more details) and were interpolated with the hydroclimatological modeling system AMUNDSEN [38,40] to 50 m

Data
We used daily runoff data (streamflow was measured at the gauging station in Vent) within this study. Hourly meteorological variables (air temperature, precipitation amount, relative humidity, wind speed and global radiation) were measured at various stations (see [38,39] for more details) and were interpolated with the hydroclimatological modeling system AMUNDSEN [38,40] to 50 m pixel resolution grids. For further model input in this study, these grids were aggregated to 100 m elevation bands (n = 19) in daily resolution.
Streamflow isotope samples (n = 98) were collected in non-equitemporal intervals for 2014 to 2016 (grab samples) and in daily intervals for 2017 (automatic water sampler) at the gauging station at Vent. The water samples were mostly collected during daytime in the ablation season. For days with multiple samples a flow-weighted mean was calculated. The sampling time is documented in Table S3.
Monthly long-term precipitation isotope data of the surrounding ANIP stations (Figure 1) were used to calculate monthly lapse rates of δ 18 O. These lapse rates were applied to the ANIP station in Obergurgl (1942 m a.s.l., 46.866667 • N, 11.024444 • E, see Figure 1) to regionalize the δ 18 O of precipitation registered there (100 m elevation bands). The Obergurgl δ 18 O values of precipitation were assumed to be transferable to the catchment since regression of precipitation δ 2 H values from Vent and Obergurgl revealed a very high correlation (R 2 = 0.89) for the years 1972 to 1975 [41]. Lapse rate data (Tables S1 and S2) and the regression plot ( Figure S1) can be found in the supplement.
Snowmelt isotopic data was not derived from sampling, but was modeled from precipitation isotopic data (as described in Section 2.3). Supraglacial meltwater (hereafter ice melt) was sampled on Hochjochferner (n = 72) at an elevation of approximately 2700 m a.s.l. at various temporal (including sub-daily and month-to month) and spatial resolutions (spatial extent: 100s of meter). Samples have been collected when most of the glacier was snow-free.
The isotope data (not ANIP data) has been partially described already [10,34,37,42]. Water samples were stored cold and dark before analyses with cavity ring-down spectroscopy (Picarro Inc., Santa Clara, CA, USA). The accuracy, indicated by the standard deviation of at least six repeated measurements, was ≤0.1% . 18 O/ 16 O isotope ratios are reported as δ 18 O (in % ) relative to the Vienna standard mean ocean water (VSMOW). The modeling workflow is described in the following sections and a schematic overview is displayed in Figure S2.

The Surface Energy-Balance Model ESCIMO
We used the modified version [32] of the surface energy-balance model ESCIMO [43] to simulate the rain, snow and ice melt input fluxes and the respective δ 18 O values for 100 m elevation bands. ESCIMO has been successfully applied at low and high mountain sites, e.g., [44,45]. A similar setup like in our study has been already applied successfully [4,31,32], where snow melt δ 18 O was simulated as the weighted average of individual snowfall events without addressing isotopic fractionation explicitly. δ 18 O of ice melt was not considered in those studies since the glacierized area of the test catchments was <5%. ESCIMO inputs are meteorological variables as described in Section 2.2 and δ 18 O in precipitation. We enhanced ESCIMO to simulate ice melt volumes. Therefore ESCIMO is enabled to calculate melt rates when snow water equivalent is zero. Potential ice melt rates are then multiplied with the glacierized fraction of the respective elevation band to calculate semi-distributed ice melt input fluxes. We ascribed the simulated ice melt volume in each elevation band a δ 18 O value to account for spatial variability as done in similar studies, e.g., [46]. This value is the median of the measured ice melt δ 18 O values at the elevation of 2700 m a.s.l. (−15.4% ) and changes according to the mean lapse rate (−0.221% /100 m) of δ 18 O in precipitation between November and May (as described in Section 2.2). The temporal variability of δ 18 O in ice melt was not considered.

The Lumped Parameter Transit Time Model TRANSEP
The lumped parameter transit time model TRANSEP [47] was used to estimate the time-invariant TTD and response time distribution (RTD), as well as the dynamic and mobile storage. TRANSEP has been applied successfully in catchments with various environmental conditions, e.g., [28,48,49], and also in high mountain catchments, e.g., [32]. ESCIMO output data (i.e., catchment water input P E : the sum of rain, snow and ice melt volumes; C E : the amount-weighted δ 18 O value of rain, snow and ice melt) were spatially aggregated to the catchment scale and used as input data for TRANSEP. The model was run at a daily resolution. TRANSEP simulates discharge for each time step t with Equation (1), where g(τ R ) is the transfer function for discharge (Equation (2)) with τ R , ϕ, τ f and τ s being the response time, fraction of fast reservoir and the MTT of the fast and the slow reservoir, respectively. Equation (2) is calibrated with measured streamflow data. (1) δ 18 O values of streamflow are simulated with Equation (3), where h(τ T ) is the transfer function for streamflow δ 18 O (Equation (4)), with τ T being the transit time, and α and β being the shape and scale parameter, respectively. Equation (4) is calibrated with measured streamflow δ 18 O.
After initial testing, we chose the two parallel linear reservoirs (TPLR) for the discharge transfer function (i.e., Equation (2)) and the gamma distribution (GM) for the tracer transfer function (i.e., Equation (4)). These gave the best results with respect to the objective function.
To optimize the simulated discharge, we used a scaled variant [50] of the Kling-Gupta efficiency (KGE) [51]. It emphasizes the flow variability term (α): A two times higher weight compared to the bias (β) and the correlation term (r) was used. Therefore we maximized KGE Q (Equation (5)) for the period 2014 to 2017.
We optimized the simulated δ 18 O in streamflow by minimizing OF C (Equation (6)), using the KGE (with the same weighting as used for optimizing discharge). We split the time series into two periods according to the temporal sampling intervals (as described in Section 2.2): (1) low frequency period (lf) 2014-2016 and (2) high-frequency period (hf) 2017. For a simpler readability we reported the KGE with standard scaling hereafter.
For the optimization procedure we used a Monte Carlo simulation with 10,000 runs and uniformly distributed parameters. The initial parameter ranges for Equations (2) and (4) (5)) and the lowest value of OF C (Equation (6)) were considered for further analyses. The mean response time (MRT) and MTT were calculated with Equations (7) and (8).
Dynamic and mobile storage has been calculated by multiplying mean flow with the MRT and the MTT, respectively. Table S3 includes the observed runoff and δ 18 O value of streamflow, as well as the catchment water input and its respective δ 18 O value as modeled with ESCIMO.

Estimating the Young Water Fraction with the Sine Wave Approach
We calculated the F yw following Kirchner [25] and Stockinger et al. [28]. First, seasonal δ 18 O cycles in input water and streamflow (Equations (9) and (10)) were fitted by multiple linear regression with an iteratively reweighted least squares algorithm (provided by Freyberg et al. [31]). C E and C S are the δ 18 O values of input water as modeled with ESCIMO and measured streamflow at time t, and were weighted with P E and observed runoff, respectively. a E , b E , k E , a S , b S and k S are the cosine and sine coefficient parameter and the vertical shift for input water and streamflow, respectively. f is the frequency (1/365. 25).
The amplitudes (A) of the seasonal δ 18 O cycles of input water and streamflow (Equations (11) and (12)) are calculated as: and the phases φ E and φ S (in rad) are calculated as: We iteratively solved Equation (15) to calculate the shape parameter α of a gamma distribution.
The scale parameter β (Equation (16)) was calculated as: By solving Equation (17) we calculated the threshold age τ yw that defines F yw with water younger than this value. F yw is calculated with Equation (18) using the lower incomplete gamma function Γ(τ yw ,α,β).
We call this F yw estimation approach 'delayed daily F yw ' since it accounts for the delayed contribution of snow and ice melt as simulated with ESCIMO. We also calculated the non-delayed monthly F yw and the delayed monthly F yw for comparison reasons. For the non-delayed monthly F yw we used the input precipitation data as used for the ESCIMO simulation (see in Section 2.2) and aggregated it to the monthly resolution (because the ANIP data used for the isotope interpolation is in monthly resolution). This approach does neither account for any temporal storage as snow and ice, nor the delayed release from snow and ice. Since resolution may play a role in estimating the F yw [28], we also estimated the delayed monthly F yw to compare the results related to the effect of the delayed and the non-delayed input (both then at the same temporal resolution). For the delayed monthly F yw we simply aggregated P E (sum) and C E (volume-weighted average) and performed the analyses in monthly resolution.
We calculated the uncertainty in all three F yw approaches by fitting each sine wave 10,000 times. For this, we used randomly sampled sine and cosine coefficient parameter (a and b in Equations (9) and (10)) out of a normal distribution with the mean being the respective parameter as estimated with Equations (9)- (18), and a standard deviation equal to the standard error of the multiple linear regression analyses. Fit statistics, as well as regression parameter and their standard errors can be found in Table S4.  Figure 2. The median (−15.4% ) was slightly higher and the interquartile range is damped compared to snowmelt. Streamflow represents the measured data (temporal variability) as used for the optimization of TRANSEP, had a median of −14.9% and was damped compared to the rain and snow melt input.

Flux and δ 18 O of Water Input into the Rofenache Catchment
Monthly isotopic lapse rates were estimated using surrounding long-term data from ANIP stations. R 2 for monthly δ 18 O-elevation relationships ranged from 0.46 to 0.94 (Table S2). The lapse rates showed seasonality with steeper gradients in winter (March: −0.265‰/100 m) and flatter gradients in summer (July: −0.112‰/100 m). Further information is provided in the supplement (Tables S1 and S2).
The interpolated precipitation input (non-delayed) into the catchment (black line in Figure 3a), i.e., the sum of snowfall and rain, reveals a slight seasonal cycle with values of about 100 mm month −1 . The delayed liquid water input into the catchment (light blue line in Figure 3a), i.e., the sum of rain, snow and ice melt modeled with ESCIMO, had a strong seasonal characteristic with close to or zero input during the winter months (December to March) and very high inputs >350 mm month −1 during the peak ablation period (June to August).

Flux and δ 18 O of Water Input into the Rofenache Catchment
Monthly isotopic lapse rates were estimated using surrounding long-term data from ANIP stations. R 2 for monthly δ 18 O-elevation relationships ranged from 0.46 to 0.94 (Table S2). The lapse rates showed seasonality with steeper gradients in winter (March: −0.265% /100 m) and flatter gradients in summer (July: −0.112% /100 m). Further information is provided in the supplement (Tables S1 and S2).
The interpolated precipitation input (non-delayed) into the catchment (black line in Figure 3a), i.e., the sum of snowfall and rain, reveals a slight seasonal cycle with values of about 100 mm month −1 . The delayed liquid water input into the catchment (light blue line in Figure 3a), i.e., the sum of rain, snow and ice melt modeled with ESCIMO, had a strong seasonal characteristic with close to or zero input during the winter months (December to March) and very high inputs >350 mm month −1 during the peak ablation period (June to August).   The non-delayed input represents precipitation, i.e., the sum of snow and rainfall, and the delayed input is the sum of rain, snow and ice melt as modeled by ESCIMO. The delayed input is aggregated to monthly resolution to make both approaches comparable (description in Section 2.5), i.e., the monthly sum P E (a) and the volume-weighted monthly mean C E (b).

Runoff and δ 18 O in Streamflow of the Rofenache at the Gauging Station in Vent
The δ 18 O value of catchment input shows a strong seasonality for both, the delayed ESCIMO catchment input (P E ) and the interpolated non-delayed precipitation catchment input (light blue and black lines in Figure 3b). The delayed monthly δ 18 O values (C E ) ranged from −18.7% to −10.2% and exhibited more damped amplitude than the non-delayed δ 18 O values (−22.0% to −5.5% ).

Young Water Fraction Estimated with the Sine Wave Approach
We fitted sine waves to the input and output data with an iteratively reweighted least squares algorithm for calculating Fyw using three different input characterizations, i.e., the delayed daily, the delayed monthly, and the non-delayed monthly input (cf. Section 2.5). The p values were below 0.1 for all regression coefficients of Equations 9 and 10, except for aE when fitting the delayed monthly input. R 2 was 0.38, 0.58, 0.82 and 0.71 for the fits of the delayed daily, delayed monthly and non-delayed monthly input and streamflow, respectively (Table S4). The delayed daily Fyw was 0.47 with an interquartile range of 0.45 to 0.50 and the delayed monthly Fyw and the non-delayed monthly Fyw were 0.54 (0.47-0.60) and 0.28 (0.26-0.30), respectively ( Table 1). The young water threshold ages (τyw) were 56, 44 and 67 days for the delayed daily, delayed monthly and non-delayed monthly Fyw, respectively. The phase shift (ϕS-ϕE) was largest for the non-delayed monthly Fyw (1.24 rad), followed by the delayed daily (0.76 rad) and delayed monthly Fyw (0.37 rad). The iterative solution of the algorithm (Equation 15) did not work for each of the 10,000 runs when applied for the delayed

Young Water Fraction Estimated with the Sine Wave Approach
We fitted sine waves to the input and output data with an iteratively reweighted least squares algorithm for calculating F yw using three different input characterizations, i.e., the delayed daily, the delayed monthly, and the non-delayed monthly input (cf. Section 2.5). The p values were below 0.1 for all regression coefficients of Equations 9 and 10, except for a E when fitting the delayed monthly input. R 2 was 0.38, 0.58, 0.82 and 0.71 for the fits of the delayed daily, delayed monthly and non-delayed monthly input and streamflow, respectively (Table S4). The delayed daily F yw was 0.47 with an interquartile range of 0.45 to 0.50 and the delayed monthly F yw and the non-delayed monthly F yw were 0.54 (0.47-0.60) and 0.28 (0.26-0.30), respectively ( Table 1). The young water threshold ages (τ yw ) were 56, 44 and 67 days for the delayed daily, delayed monthly and non-delayed monthly F yw , respectively. The phase shift (φ S -φ E ) was largest for the non-delayed monthly F yw (1.24 rad), followed by the delayed daily (0.76 rad) and delayed monthly F yw (0.37 rad). The iterative solution of the algorithm (Equation (15)) did not work for each of the 10,000 runs when applied for the delayed monthly F yw , likely due to negative values of φ S -φ E (i.e., streamflow runs off before water input), which may have resulted from the high p value and large standard error of a E (Table S4). Hence we calculated the interquartile range of the delayed monthly F yw by using the ratio A S /A E and were not able to calculate the interquartile range for α and τ yw (see Table 1). The shape parameter α for the three approaches ranges from 0.25 to 0.93. Table 1 lists the described metrics and their interquartile ranges. Table 1. Amplitude, threshold age, shape parameter and phase shift for the three applied F yw approaches (with delayed daily, delayed monthly and non-delayed monthly input). The interquartile range (in brackets) is estimated from 10,000 calculations.  Figure 5 shows the cumulated RTD and TTD using the convolution integral method (TRANSEP), as well as the gamma distribution derived during the delayed daily F yw approach (sine wave approach). The delayed daily F yw had a narrow band (interquartile range of 10,000 runs), as well as the RTD and TTD estimated with TRANSEP (interquartile range of best 100 runs). The cumulated TRANSEP TTD at the threshold age of 56 days ranged between 0.43 and 0.45 (25th−75th percentile) with a median of 0.44 for the best 100 runs ( Table 2). Table 2 reveals the transfer function parameters for the TPLR and the GM. Median (interquartile range) MRT and MTT were 0.13 (0.08 to 0.19) and 9.5 (8.2-10.5) years, respectively. The median of the cumulated TTD at the threshold age of 5 years (F ow ) was 0.28 (i.e., 28% of water was older than 5 years). The median (interquartile range) of the dynamic and mobile storages estimated with TRANSEP was 195 (110-284) mm and 13,846 (11975-15382) mm, respectively.

Age Distribution and Subsurface Storage Potential
monthly Fyw, likely due to negative values of ϕS-ϕE (i.e., streamflow runs off before water input), which may have resulted from the high p value and large standard error of aE (Table S4). Hence we calculated the interquartile range of the delayed monthly Fyw by using the ratio AS/AE and were not able to calculate the interquartile range for α and τyw (see Table 1). The shape parameter α for the three approaches ranges from 0.25 to 0.93. Table 1 lists the described metrics and their interquartile ranges. Table 1. Amplitude, threshold age, shape parameter and phase shift for the three applied Fyw approaches (with delayed daily, delayed monthly and non-delayed monthly input). The interquartile range (in brackets) is estimated from 10,000 calculations.

Method AE (‰) As (‰) τyw (d) α (-) ϕS-ϕE (rad) Fyw (-)
Delayed daily input 3.  Figure 5 shows the cumulated RTD and TTD using the convolution integral method (TRANSEP), as well as the gamma distribution derived during the delayed daily Fyw approach (sine wave approach). The delayed daily Fyw had a narrow band (interquartile range of 10,000 runs), as well as the RTD and TTD estimated with TRANSEP (interquartile range of best 100 runs). The cumulated TRANSEP TTD at the threshold age of 56 days ranged between 0.43 and 0.45 (25th−75th percentile) with a median of 0.44 for the best 100 runs ( Table 2). Table 2 reveals the transfer function parameters for the TPLR and the GM. Median (interquartile range) MRT and MTT were 0.13 (0.08 to 0.19) and 9.5 (8.2-10.5) years, respectively. The median of the cumulated TTD at the threshold age of 5 years (Fow) was 0.28 (i.e., 28% of water was older than 5 years). The median (interquartile range) of the dynamic and mobile storages estimated with TRANSEP was 195 (110-284) mm and 13,846 (11975-15382) mm, respectively.

How Large is the Subsurface Storage Potential?
Using TRANSEP, we found a relatively large dynamic (195 mm) and a very large mobile storage (13,846 mm), which were, though, in the range of values reported for other high-mountain catchments (up to 13,000 mm; see Section 1; e.g., [4,32]). Overburden material (such as glacial deposits, talus or alluvium) and soil/bedrock may serve as a physical storage for this water. Since subsurface information from wells or boreholes were not available we could not yet evaluate this hypothesis, but overburden material covers a high percentage of the catchment (42%, [34]) and is assumed to be many meters deep in some parts of the catchment [52]. Soil is usually <1 m thin, providing a relatively limited storage potential, hence bedrock is attributed to store and transmit larger portions of water, most likely in fractures. The bedrock storage might be obvious due to the large elevation difference of almost 2000 m from the highest peak to the catchment outlet, indicating a large volume of material, which provides potential storage [53]. The total catchment water storage (not snow and ice) hence provides a huge tank for tracer mixing. The large range in dynamic and mobile storage in the TRANSEP results (Table 2) indicates a high uncertainty, as known also for other transit time models (e.g., [23]). This uncertainty is triggered by the high variability of the parameter τ s and β (Table 2), which were not well-constrained during the calibration procedure (see Section 4.3.3).

How Old is Streamflow?
Contrary to the large storage volumes estimated with TRANSEP, we also found a high F yw of 0.47 as estimated with the sine wave approach. Our estimate was markedly higher than the flow-weighted mean of 254 catchments (0.34) in [11], or estimates from other high-elevation catchments (e.g., maximum flow-weighted F yw was 0.22 [31]). This can be attributed to mainly two effects. Firstly, we used the delayed input of snow and ice melt simulated with ESCIMO. This flattened out the seasonal cycle of δ 18 O (see Figure 3 and parameter A E in Table 1). The non-delayed F yw of 0.28 was much smaller and closer to the other values reported in the literature. Secondly, and related to the first effect, snow and ice melt had a disproportional high impact on runoff and streamflow δ 18 O because snow and ice melt input volumes were relatively large and relatively fast transmitted (especially large portions of ice melt) to the stream network during a limited period of the hydrological year (i.e., the ablation period). Up to now, there are no studies available that address this topic, but we thought this was intuitively plausible and should be investigated in future studies with more detail.
Interestingly, we also found a marked F ow : 28% of streamflow water was considered to be older than five years. TRANSEP F yw was 0.44, so the resulting concomitant occurrence of high fractions of younger and older water in streamflow might be a combination of two facts: (1) the high F yw triggered by the fast transmittal of ice and snow melt as outlined above and (2) the high F ow provided by the large storage described in Section 4.1. Longer transit times (>5 years) are seen as critical when applying the convolution integral method with stable isotopes [54][55][56], but method restrictions were discussed in Section 4.3.3.
The estimated MTT of 9.5 years was potentially biased, as recently postulated by Kirchner [25,26]. Although it is known that "stable isotopes are effectively blind to the long tails of transit time distributions" [25], we still used δ 18 O (as also suggested by Seeger and Weiler [32]), since no other tracers (e.g., 3 H) were available. Our estimate was at the top end compared to MTTs (0.6-10.5 years) of other snow-melt dominated high-elevation catchments in Switzerland [32]. Behrens et al. [57] estimated the winter baseflow MTT (4 years) for the period 1972 to 1975 at approximately four years for the same catchment, but they used a different tracer ( 3 H), a different transfer function (exponential model), non-delayed input and they focused on winter baseflow transit time (low flow regime). Unfortunately, any other direct comparison to other glacierized catchments was still lacking.

Modeling Catchment Water Input and its δ 18 O Value with ESCIMO
Using monthly δ 18 O values in precipitation to simulate daily δ 18 O water input values is a restriction and hinders the determination of highly variable δ 18 O values of rain and snowmelt, as frequently observed in empirical studies (e.g., [34,58]). Additionally, the use of the δ 18 O lapse rate in precipitation is a source of uncertainty, since local effects (e.g., windward vs. leeward) may not be reproduced by the regional pattern of the lapse rate. Still, our approach represents a stable approach in that it uses long-term data. Our lapse rates compared well to the lapse rate in Austria (−0.19% /100 m) [59] and Switzerland (−0.27% /100 m) [60]. Due to the absence of high-frequency and high-density input δ 18 O data in our study area, we were bound to the design used in this study, but future work should include intense sampling and/or disaggregation methods.
Although we used δ 18 O ice melt data measured in the field to derive δ 18 O ice melt input data for ESCIMO, the spatio-temporal variability, as observed in empirical studies [10,61,62], provides an alternative source of uncertainty. An isotopic lapse rate was not directly estimated from the ice melt samples, rather it was estimated from the mean winter precipitation lapse rate. We used temporally invariant δ 18 O values of ice melt due to missing temporally distributed data throughout the four-year period, but temporal variability has also not yet been described intensively [10,61,62]. Our estimate (0.221% /100 m) compared well to values presented in other studies; e.g., He et al. [46] directly estimated an isotopic glacier melt lapse rate from δ 18 O measurements at −0.226% per 100 m in a 17% glacierized high-elevation catchment in the Tien Shan and used it for their modeling study.
The resolution of ESCIMO (100 m elevation bands and one day) seems a reliable balance in terms of data availability and sufficient process representation. Recent studies used fully distributed approaches to model the δ 18 O snowmelt input into the catchment [29,63,64], but the advantage of this more data-intensive approach over the semi-distributed one needs to be evaluated. It remains a critical topic how much detail is needed when considering that the distributed data is aggregated to the catchment scale for further input in simpler models, such as lumped parameter transit time models or for the calculation of F yw with the sine wave approach.
A shortcoming of ESCIMO is the missing implementation of isotopic fractionation during the transition from snowfall to snowmelt. We used the volume-weighted mean of snowfall as snowmelt δ 18 O similar to Seeger and Weiler [32]. According to the hypothesis of Beria et al. [60], the median of individual snowfall events is close to the median of outflowing snowmelt water over one season, which makes this approach plausible. On the other side, the short-term temporal variability (e.g., depleted signal during first melt impulse) is not directly incorporated, and the data-intensive accounting for isotopic fractionation due to top layer sublimation and percolation through the snowpack should be accounted for in detailed studies, as initially presented by [65]. Validation data, especially with appropriate temporal and spatial resolution (δ 18 O and areal snowmelt) is often absent. In our study, liquid water input by ice and snow melt into the catchment seemed plausible and was comparable to the results of other studies in the catchment [38,39].
Although the approach was comprised of some restrictions, mostly due to data availability, we still thought that considering snow and ice melt-at least at the semi-distributed scale-was mandatory when characterizing catchment water inputs for the further estimation of the F yw or as input in transit time models.

The Young Water Fraction of a Glacierized High-Elevation Catchment
We could not support the hypothesis of Jasechko et al. [11] that mountain catchments have less young streamflow, since we found a higher F yw compared to most other mountain catchments (see Section 4.2). We attributed this result to the use of the delayed input, as outlined in Section 4.2, and the occurrence of high snow and ice melt contributions to streamflow with typical residence times of hours (glacier) to weeks (snowpack) [66,67], which might not have been included in the dataset of [11]. The delayed monthly F yw (0.54) was twice the non-delayed monthly F yw (0.28). In a pre-study, a mean F yw of 0.16 was estimated using non-delayed monthly precipitation δ 18 O data from Vent for the periods 1972-1977 and 2015-2016 [41]. Contrary to the marked differences we observed for our dataset when computing F yw with non-delayed input data (direct precipitation input) or delayed input data (unretained precipitation + snow/ice melt from ESCIMO), Freyberg et al. [31] found no significant difference between the results of the two ways to compute F yw for 22 Swiss catchments, likely due to the absence of marked ice melt inputs from glaciers, as they did not include glacierized catchments in their study.
The temporal resolution played a role when applying the F yw approach, since we estimated a seven percent higher F yw when estimating the delayed monthly F yw compared to the delayed daily F yw , but the difference was not as large as the one presented by Stockinger et al. [28], who estimated F yw with daily data as almost twice that of F yw when using weekly data.
Although all regression coefficient parameters of the delayed daily and the non-delayed monthly sine wave fits had p values below 0.1, the delayed daily fit had a markedly lower variance explained (R 2 = 0.38) compared to the non-delayed fit (R 2 = 0.82). This may be due to the missing input volume during the winter months and the extreme high input volumes during the ablation period. Such variable input is a serious challenge for the volume-weighted sine wave fit. Unfortunately, F yw estimates from other highly glacierized catchments were absent for comparison and hence we call for future studies.

TRANSEP Applied in Glacierized High-Elevation Catchment
We estimated the catchment TTD (not baseflow TTD) with δ 18 O in a highly glacierized catchment for the period 2014 to 2017. Therefore we used streamflow δ 18 O samples collected during various flow regime periods, from winter baseflow to annual peak flow. This time-scale should provide both, short-term (daily) to mid-term (several years) transit time information [24]. Due to the limited number of samples, all 98 samples were used for the optimization procedure, similar to [68]. The TTD may be biased in that it does not capture very short (sub-daily) and long (decadal) transit time information due the sampling interval and due to the tracer suitability issue, respectively [24,[54][55][56].
We also used a time-invariant TTD modeling approach due to data limitations, since time-variant TTD modeling approaches typically require longer tracer records and/or higher sampling frequencies [68][69][70][71][72][73]. Although headwater catchments-often characterized by high responsiveness and high seasonal variability in the climatic conditions-are known to be non-steady systems [74,75], this approach is still feasible [32,68], but time-variant approaches are becoming popular, e.g., [63]. The MTT estimate (8.2-10.5 years) must be carefully interpreted since δ 18 O is increasingly recognized to be blind to the tail of the TTD and transit times longer than five years [25,26,32,[54][55][56]68].
The cumulated TRANSEP TTD at the threshold age of 56 days compared well to the F yw estimated with the sine wave approach, which was already reported elsewhere (e.g., [28]). This may result from the facts that the GM underlied both approaches and that the same input was used, even though the shape parameter α differs (compare α values in Tables 1 and 2), resulting in different TTDs. This led to the conclusion, that steady-state transit time models (like TRANSEP) likely provide unbiased mid-term transit time information (up to several months). Two main problems arose during the modeling procedure: i) runoff in 2014 has a markedly lower goodness of fit; ii) during the streamflow and tracer recession period in late 2015 the δ 18 O of streamflow was not well emulated by TRANSEP. Up to now we are not able to explain these features.
Four degrees of uncertainty are worth being highlighted within the scope of our study: • Structural uncertainty: We used a top-down modeling approach including a pre-defined model structure and tested it in our complex study catchment. The transfer functions were chosen by prior modeling experiments and the TPLR for runoff and the flexible GM for streamflow δ 18 O were used. The latter allows for both, fast tracer throughputs and relatively long transit times [76], and did not assume a well-mixed reservoir, which was suitable for application in our catchment. Both transfer functions were proved best regarding the objective functions.

•
Parameter uncertainty: Three out of five TRANSEP parameters were well-constrained during the Monte Carlo simulation (ϕ, τ f and α). τ s and β were less well-constrained, leading to the observed uncertainty of the TTD, RTD and storage estimates (see Figure 5 and Table 2). The physical interpretation of the low α suggests a highly non-linear streamflow tracer response ( Table 2, [77]).

•
Input uncertainty: We used a single ESCIMO run as input data for TRANSEP and focused on the applicability of and the uncertainty within TRANSEP, but input uncertainty should be investigated in future studies as this represents an important source of uncertainty in TTD modeling, e.g., [68]. • Uncertainty due to the optimization procedure: The choice of the objective function was arbitrary, but after initial model experiments it became apparent that the splitting of the streamflow δ 18 O time series, as well as the higher weighted flow variability term (as used for the KGE), were relevant.
Future improvements can include (i) the calibration of lumped parameter transit time models against F yw (e.g., [78]), making the approach not independent, or (ii) using a bottom-up modeling chain, starting with a perceptual model (and not a priori determined TTDs) or using time-variant TTDs.

Conclusions
We estimated the dynamic and mobile storage, the F yw (with the sine wave approach), and the TTD (with the convolution integral method using TRANSEP) in a highly glacierized catchment for the period 2014 to 2017. For input in the latter, we used lumped δ 18 O values of rain, snow and ice melt, as simulated with ESCIMO. Being the first study carried out in such a simulation setup and environment, we found large storages and a relatively high F ow , but contrarily also a high F yw . This led to the conclusion that the basin behaves to some degree like a 'Teflon basin', especially because of the large contribution of fast transmitted ice melt, and to some degree like a huge sponge with a very much delayed release of water, especially due to the large potential subsurface water storage volume. The two behaviors are probably distributed more in space and less in time. We also found a highly non-linear streamflow tracer response. An appropriate input characterization is obligatory for both, the estimation of the F yw and the application of lumped parameter transit time models, i.e., accounting for the delayed contribution of snow and ice melt in glacierized catchments. Regarding mid-term transit time information (up to a few months), F yw estimated with TRANSEP and the sine wave approach were well comparable. Unfortunately, very short (sub-daily) and long (decadal) transit time information is not covered in the modeling approach, but we assumed that this might play an important role in the biogeochemical cycle of the catchment, especially if one considers very fast catchment responses to rain on ice events or the potential of deeper groundwater contribution to winter baseflow. Further studies to enhance the methodology for glacierized catchments may focus on time-variant TTDs and more detailed input characterization (e.g., higher sampling frequency and density, disaggregation/regionalization methods and accounting for isotopic fractionation of rain, snow and ice melt).
Author Contributions: U.S. and M.W. conceptualized the study. U.S., M.W. and S.S. reviewed and edited the manuscript. S.S. and J.S. developed the model code and performed the simulation with the coupled ESCIMO-TRANSEP approach. J.S. conducted the field work, calculated the F yw with the sine wave approach and prepared the original draft.
Funding: This work is associated to the research project HydroGeM 3 , which was funded by the Austrian Academy of Sciences. We also very much appreciate the University of Innsbruck for support in the financial aspect of the publication.