Articles | Volume 27, issue 15
Research article
14 Aug 2023
Research article |  | 14 Aug 2023

Uncertainty in water transit time estimation with StorAge Selection functions and tracer data interpolation

Arianna Borriero, Rohini Kumar, Tam V. Nguyen, Jan H. Fleckenstein, and Stefanie R. Lutz

Transit time distributions (TTDs) of streamflow are useful descriptors for understanding flow and solute transport in catchments. Catchment-scale TTDs can be modeled using tracer data (e.g. oxygen isotopes, such as δ18O) in inflow and outflows by employing StorAge Selection (SAS) functions. However, tracer data are often sparse in space and time, so they need to be interpolated to increase their spatiotemporal resolution. Moreover, SAS functions can be parameterized with different forms, but there is no general agreement on which one should be used. Both of these aspects induce uncertainty in the simulated TTDs, and the individual uncertainty sources as well as their combined effect have not been fully investigated. This study provides a comprehensive analysis of the TTD uncertainty resulting from 12 model setups obtained by combining different interpolation schemes for δ18O in precipitation and distinct SAS functions. For each model setup, we found behavioral solutions with satisfactory model performance for in-stream δ18O (KGE > 0.55, where KGE refers to the Kling–Gupta efficiency). Differences in KGE values were statistically significant, thereby showing the relevance of the chosen setup for simulating TTDs. We found a large uncertainty in the simulated TTDs, represented by a large range of variability in the 95 % confidence interval of the median transit time, varying at the most by between 259 and 1009 d across all tested setups. Uncertainty in TTDs was mainly associated with the temporal interpolation of δ18O in precipitation, the choice between time-variant and time-invariant SAS functions, flow conditions, and the use of nonspatially interpolated δ18O in precipitation. We discuss the implications of these results for the SAS framework, uncertainty characterization in TTD-based models, and the influence of the uncertainty for water quality and quantity studies.

1 Introduction

Understanding how catchments store and release water of different ages has significant implications for flow and solute transport, as water ages encapsulate information about flow paths' characteristics (McGuire and McDonnel2006; Botter et al.2011), the contact time of solutes with the soil matrix (Benettin et al.2015a; Hrachowitz et al.2016), and vulnerability assessment (Kumar et al.2020). This plays an important role in water resources protection and management, and it requires a tool that can effectively describe catchment-scale transport processes (Rinaldo and Marani1987). The age of water in outflows is commonly referred to as transit time (TT), i.e. the time that elapses between the entry of a water parcel into the catchment via precipitation and its exit via streamflow or evapotranspiration. Accordingly, the transit time distribution (TTD) describes the whole spectrum of transit times in outflows (Botter et al.2005; Van der Velde et al.2010). Early studies have often assumed simplified steady-state transport models, resulting in time-invariant TTDs (Niemi1977; Rinaldo et al.2006). However, experimental simulations have shown that TTDs are time-variant due to the variability in meteorological forcing (Botter et al.2010; Hrachowitz et al.2010; Heidbüchel et al.2020) and the activation/deactivation of flow paths in response to varying hydrologic conditions (Ambroise2004; Heidbüchel et al.2013). Recent research has introduced new models for representing time-variant TTDs, for example, allowing for the estimation of TTDs without making prior assumptions about their shape (Kirchner2019; Kim and Troch2020) or with the parameterization of the StorAge Selection (SAS) functions (Rinaldo et al.2015; Harman2019). SAS functions describe how catchments selectively remove water of different ages from storage for outflows, and they have led to a new framework of nonstationary transport models based on water age, which have been successfully applied in various studies (Benettin et al.2015b; Queloz et al.2015; Kim et al.2016; Lutz et al.2017; Wilusz et al.2017; Nguyen et al.2021).

Model-based TTDs are subjected to uncertainty, which limits their ability with respect to decision support. In general, model prediction uncertainty stems from model inputs, structure, and parameters (Beven and Freer2001). As TTDs are not directly observable, conservative environmental tracers (e.g. oxygen isotopes, such as δ18O) in inflow and outflows are commonly used to infer water ages (Hrachowitz et al.2013; Birkel and Soulsby2015; Stockinger et al.2015). Long-term, high-frequency tracer data with an appropriate spatial distribution are generally recommended for a sufficient understanding of TTD dynamics across a wide range of fast and heterogeneous hydrological behaviors (Kirchner et al.2004; Danesh-Yazdi et al.2016; von Freyberg et al.2017). Therefore, a lack of appropriate tracer data coverage can hamper our understanding of TTD dynamics at the desired resolution (McGuire and McDonnel2006). Additionally, uncertainty in the driving hydroclimatic fluxes, such as precipitation, discharge, and evapotranspiration, could propagate into the uncertainty in the modeling results. Further uncertainty emerges from the model structure due to the difficulty in representing physical processes because of our incomplete knowledge of complex reality (Ajami et al.2007). Finally, specification of model parameters is also an important source of uncertainty (Beven2006; Kirchner2006), as the best-fit parameters may suffer from equifinality (Schoups et al.2008).

A few studies have investigated the uncertainty in the estimated TTDs with SAS models. Danesh-Yazdi et al. (2018) and Jing et al. (2019) analyzed the effect of interactions between distinct flow domains, external forcing, and recharge rate on the resulting TTDs. Several works (Benettin et al.2017; Wilusz et al.2017; Rodriguez et al.2018, 2021) have explored model parameter uncertainty and suggested that additional types of tracers, data on physical characteristics of the catchment, and parsimonious parameterization may help to further reduce parametric uncertainty in the SAS models. More recently, Buzacott et al. (2020) investigated how gap-filling of the δ18O record in precipitation propagated uncertainty into the simulated mean water transit time (MTT), i.e. the average time it takes for water to leave the catchment (McDonnel et al.2010).

Despite the studies cited above, there are other aspects that are particularly significant for SAS modeling and cause uncertainty in the simulated TTDs that have not yet been thoroughly investigated. First, isotope data are generally sparse globally in space and time (von Freyberg et al.2022), due to laborious and costly sampling campaigns limited to well-equipped areas (Tetzlaff et al.2018). As SAS models require continuous time series of input tracer data, different methods for temporal interpolation could be used to reconstruct isotope values in precipitation; consequently, the interpolated input data are subject to uncertainty. Furthermore, the input data of SAS models are influenced by whether the tracer data in precipitation are collected at a single location within the catchment or at multiple locations. In the latter scenario, there is a need to account for the spatial variability in the tracer composition in precipitation, which is commonly done via spatial interpolation. Choosing data from one approach (i.e. tracer data from a single location) over the other (i.e. tracer data spatially interpolated based on multiple locations, including stations outside the catchment boundaries) can potentially result in different resulting TTDs. Finally, SAS functions, employed to model TTDs, must be parameterized, and their functional forms need to be specified a priori. Commonly used forms are the power law (Benettin et al.2017; Asadollahi et al.2020), beta (van der Velde et al.2012; Drever and Hrachowitz2017), and gamma (Harman2015; Wilusz et al.2017) distributions. However, there is no general agreement on which SAS function should be used, as the hydrological processes that control the patterns and dynamics of the subsurface vary across catchments. Therefore, the most convenient approach is to simply rely on a specific parameterization over another and estimate its parameters (Harman2015). All of these aspects, related to model input, structure, and parameters, induce uncertainty in the simulated TTDs. To date, the role of these individual uncertainty sources and their combined effect on the modeled TTDs have not been adequately discussed.

This study bridges the aforementioned gaps by specifically exploring the combined effect of tracer data interpolation and model parameterizations on the simulated TTDs. We investigated TTD uncertainty using an SAS-based catchment-scale transport model applied to the upper Selke catchment, Germany. We evaluated TTDs resulting from 12 model setups obtained by combining distinct interpolation techniques of δ18O in precipitation and parameterizations of SAS functions. For each model setup, we searched for behavioral parameter sets (i.e. those providing acceptable predictions) based on model performance for in-stream δ18O, and we evaluated the sources of uncertainty and their combined effects in the modeled TTDs. Overall, our results provide new insights into the uncertainty characterization of TTDs, particularly in the absence of high-frequency tracer data, and the use of SAS functions as well as the implications of TTDs' uncertainty on water quantity and quality studies.

2 Study area and data

The upper Selke catchment is located in the Harz Mountains in Saxony-Anhalt, central Germany (Fig. 1). The study site is part of the Bode region, an intensively monitored area within the TERENO (TERrestrial ENvironmental Observatories; Wollschläger et al.2017) network. The catchment has a drainage area of 184 km2, the altitude ranges between 184 and 594 m above mean sea level, and the mean slope is 7.65 %. Land use is dominated by forest (broadleaf, coniferous, and mixed forest) and agricultural land (winter cereals, rapeseed, and maize), representing 72 % and 21 % of the catchment, respectively. The soil is largely composed of Cambisols and the underlying geology consists of schist and claystone, resulting in a predominance of relatively shallow flow paths (Dupas et al.2017; J. Yang et al.2018).

Figure 1The upper Selke catchment, showing precipitation sampling points (purple dots), the river network (blue lines), and the elevation (in meters above sea level) as a colored map. The inset presents the location of the upper Selke catchment in Germany.


Daily hydroclimatic and monthly tracer data in the upper Selke catchment were available for the period between February 2013 and May 2015. Precipitation (P) was taken from the German Weather Service, whereas discharge (Q) and evapotranspiration (ET) were simulated data obtained from the mesoscale Hydrological Model (mHM; Samaniego et al.2010; Kumar et al.2013), as continuous measurements were not available for the given outlet and period. A thorough evaluation of mHM performance for past measurements has been conducted in previous studies (Zink et al.2017; X. Yang et al.2018; Nguyen et al.2021). The average annual P, Q, and ET are 703, 108, and 596 mm, respectively. The area is characterized by high flow during November–May (average Q=0.88 m 3 s−1) and low flow during June–October (average Q=0.42 m 3 s−1). Evapotranspiration is higher in June (109 mm per month) and lower in December (10 mm per month). The average monthly temperature ranges from −0.7C in January to 17 C in July. The δ18O values in precipitation (δ18OP) and in streamflow (δ18OQ) at a monthly resolution were taken from Lutz et al. (2018) and are displayed in Fig. 2. Values of δ18OP were used in the form of “raw” (i.e. values collected at the catchment outlet) and “processed” (i.e. values collected at multiple locations and spatially interpolated using kriging) data (see Sect. 3.2 for more details). The variability in δ18OP was larger than that in δ18OQ (Fig. 2) because of the damping of the precipitation signal due to mixing and dispersion within the catchment. Temperature dependence caused more depleted (i.e. more negative) δ18OP in winter than in summer (Fig. 2).

Figure 2Data of δ18O in precipitation (kriged values as pink dots and raw values as yellow dots) and streamflow (blue dots).


3 Methods

3.1 Catchment-scale transport model

In this study, we used the tran-SAS model (Benettin and Bertuzzo2018a) for describing the catchment-scale water mixing and solute transport based on SAS functions. The catchment was conceptualized as a single storage S(t) (mm), whose water-age balance can be expressed as follows (Benettin and Bertuzzo2018a):

(1)S(t)=S0+V(t)(2)ST(T,t)t+ST(T,t)T=P(t)-Q(t)ΩQ(ST,t)-ET(t)ΩET(ST,t),(3)with an initial condition of ST(T,t=0)=ST0(T)(4)and a boundary condition of ST(0,t)=0.

Here, S0 (mm) is the initial storage; V(t) (mm) represents the storage variations; P(t) (mm d−1), Q(t) (mm d−1), and ET(t) (mm d−1) are precipitation, discharge, and evapotranspiration, respectively; ST(T,t) (mm) is the age-ranked storage; ST0(T) (mm) is the initial age-ranked storage; and ΩQ(ST,t) (–) and ΩET(ST,t) (–) are the cumulative SAS functions for Q and ET, respectively.

By definition, the TTD of streamflow pQ(T,t) (d−1) is calculated as follows (Benettin and Bertuzzo2018a):

(5) p Q ( T , t ) = Ω Q ( S T , t ) S T S T T .

The isotopic signature in streamflow CQ(t) (‰) can be obtained as follows (Benettin and Bertuzzo2018a):

(6) C Q ( t ) = 0 + C S ( T , t ) p Q ( T , t ) d T ,

where CS(T,t) (‰) is the isotopic signature of a water parcel in storage. Equations (5) and (6) also apply for ET.

In this study, we tested three SAS parameterizations: the power law time-invariant (PLTI; Eq. 7, Queloz et al.2015), power law time-variant (PLTV; Eq. 8, Benettin et al.2017), and beta time-invariant (BETATI; Eq. 9, Drever and Hrachowitz2017) distributions. Here, they are expressed as probability density functions in terms of the normalized age-ranked storage PS(T,t) (–), also known as fractional SAS functions (fSAS):


The parameters k, α, and β determine the catchment's water-age preference for outflows, while B(α,β) is the two-parameter beta function. If k<1, or if α<1 and β>1, the system tends to discharge young water. If k>1, or if α>1 and β<1, the catchment preferably releases old water. The case of k=1 or α=β=1 describes no selection preference (i.e. complete water mixing). PLTV is characterized by k(t) varying linearly over time between two extremes, k1 and k2, as a function of the catchment wetness wi (–), i.e. wi(t)=(S(t)-Smin)/(Smax-Smin), where Smin and Smax are the respective minimum and maximum storage values over the entire period.

3.2 Interpolation techniques for δ18O in precipitation

We tested the model with two spatial representation and two temporal interpolation methods for δ18OP to explore the impact of input tracer data on model performance, results, and uncertainty. To evaluate the effect of spatial representation, we firstly used single-point δ18OP measurements, which we refer to in the following as raw δ18OP. These measurements, obtained from Lutz et al. (2018), were taken at the catchment outlet. The selection of δ18OP at the outlet assumes a precipitation collector close to the stream gauge at the outlet, which is a common occurrence in many catchments for logistical reasons. Indeed, the outlet, where in-stream δ18O is sampled, serves as the location where all precipitation inputs across the catchment are integrated. For convenience, precipitation monitoring is also often conducted at or near the gauging station at the outlet. Secondly, we used spatially interpolated δ18OP with kriging based on multiple locations. The spatial interpolation was conducted in Lutz et al. (2018) using raw δ18OP from 24 precipitation collectors spread over the larger Bode region and using altitude as external drift. In a further step, the kriged δ18OP data were weighted with spatially distributed monthly precipitation to obtain representative estimates for the study catchment. In our study, the kriged δ18OP resulted in slightly more negative values than the raw δ18OP from the catchment outlet (Fig. 2) because of the inclusion of more depleted δ18OP values from locations with higher altitudes during the kriging process. By considering these two options for the spatial representation of δ18OP, we intend to assess the influence of spatial variability and uncertainty in the simulated outputs between two opposing cases i.e. raw isotopes representing the simplest approach and kriged isotopes derived from a more sophisticated method. While there are other possibilities for the spatial representation of δ18OP, our choice allows us to effectively address our research question regarding the effects on SAS models of tracer data in precipitation collected at a single location within the catchment or spatially interpolated from multiple sites.

SAS model results are sensitive to the choice of the temporal resolution of input tracer data, and a finer resolution is generally recommended to achieve a satisfactory level of detail (Benettin and Bertuzzo2018a). Additionally, a forward Euler scheme was employed to solve Eq. (2), whose precision increases with high-frequency time steps. For these reasons, we reconstructed daily δ18OP estimates from monthly values with two interpolation schemes. First, we used a step function in which the values between two consecutive samples assumed the value of the last sample. Second, we used a sine interpolation due to the fact that δ18OP samples typically exhibit pronounced seasonal variations with more depleted values in winter than in summer (Fig. 2). The sine-wave function has been used in several studies to describe temporal variation in δ18OP (McGuire and McDonnel2006; Feng et al.2009; Allen et al.2019). The seasonal pattern of δ18OP values over a period of 1 year can be described as follows (Kirchner2016):

(10) δ 18 O P ( t ) = a P cos ( 2 π f t ) + b P sin ( 2 π f t ) + k P ,

where a and b are regression coefficients (–), t is the time (decimal years), f is the frequency (yr−1), and k is the vertical offset of the isotope signal (‰). The coefficients a and b were estimated by fitting Eq. (10) to monthly δ18OP values using the iteratively re-weighted least-squares (IRLS) estimation (von Freyberg et al.2018). In our study, we reproduced the daily frequency isotopic data through the estimated regression coefficients of Eq. (10). Figure 3 displays the daily kriged and raw δ18OP values simulated via step function and sine interpolation; by employing step function and sine interpolation as techniques to reconstruct tracer data in precipitation, we aim to analyze the effects on SAS-based results from two relatively simple, rather opposing approaches: one focusing on individual measurements and the other on seasonality.

Figure 3Predicted δ18O in precipitation (kriged values as pink lines and raw values as yellow lines) via step function and sine interpolation.


3.3 Experimental design

In this study, different scenarios were used to quantify uncertainty in the modeled results. We tested 12 setups composed of three SAS functions (PLTI, PLTV, and BETATI), two temporal interpolations (step and sine function), and two spatial representations (raw and kriged values) of δ18OP (Table 1). For each setup, we performed a Monte Carlo experiment by running the model with 10 000 parameter sets generated by Latin hypercube sampling (LHS; McKay et al.1979). Model parameters and their search ranges are shown in Table 2. A 5-year warm-up period (i.e. repetition of the input data) from February 2008 to January 2013 was performed to reduce the impact of model initialization. The period from February 2013 to May 2015 was used to infer behavioral parameters (i.e. parameter sets giving acceptable predictions) and, subsequently, to interpret model results. The initial concentration of δ18O in storage was set to 9.2 ‰, coinciding with the mean δ18OQ over the study period.

Table 1List of model setups.

Download Print Version | Download XLSX

Table 2Model parameters and search ranges.

Download Print Version | Download XLSX

The informal likelihood of the Sequential Uncertainty Fitting (SUFI-2; Abbaspour et al.2004) procedure was applied to account for uncertainty in the parameter sets and resulting modeled estimates. In SUFI-2, the uncertainty is represented by a uniform distribution, which is gradually reduced until a specific criterion is reached. In our study, we calibrated the values of model parameters until the predicted output matched the measured δ18OQ to a satisfactory level, defined by an objective function. As the objective function, we employed the Kling–Gupta efficiency (KGE; Gupta et al.2009), and once the criterion of KGE  0.5 was satisfied, we defined a set of behavioral solutions for each model setup. However, as the aim of this study is to investigate the impact of various sources of uncertainty on simulated outputs, rather than to determine the best model setup, we decided to set a fixed sample size and narrow down those solutions generated by SUFI-2 in the previous step. Setting a fixed sample size ensures the comparability of results across the setups, as different sample sizes could influence the uncertainty analysis (i.e. the greater the number of behavioral solutions, the wider the uncertainty band). By fixing the sample size, we can still meet the requirement of a minimum acceptable KGE value (i.e. KGE  0.5). In this study, we determined the final behavioral solutions by using a fixed sample size that corresponds to the best 5 % parameter sets and modeled results in terms of the KGE.

To assess the range of possible behavioral solutions and understand the level of uncertainty associated with it, we calculated the 95 % confidence interval (CI) derived by computing the 2.5 % and 97.5 % percentile values of the cumulative distribution in the parameters and time series of output variables (Abbaspour et al.2004). These percentile values represent the lower and upper bounds of the CI, respectively. In our experimental setup, the main output variables were the in-stream δ18O signature and backward median transit time (TT50, in days, i.e. the time it takes for half of the water particles to leave the system as streamflow at the catchment outlet). Time series of TT50 were extracted directly from daily TTDs (Eq. 5) and used as a metric for the streamflow age. This was done because TTDs are typically skewed with long tails (Kirchner et al.2001); hence, the median is often a more suitable metric than, for example, the MTT, as it is less impacted by the poor identifiability of the older water components (Benettin et al.2017).

4 Results

4.1 Simulated δ18O in streamflow and model performance

Modeled δ18O values in streamflow (δ18OQ) represented by the 95 % confidence interval (CI) in the ensemble solution are displayed in Fig. 4. The results reveal how the predicted δ18OQ values enveloped the measured isotopic signature by reproducing its seasonal fluctuations, with depleted (i.e. more negative) values in winter and enriched (i.e. less negative) values in summer. Although the behavioral parameter sets were able to capture the seasonal isotopic trend, they poorly reproduced the exact values; therefore, the ensemble simulations are characterized by a non-negligible uncertainty.

Figure 4 shows the distinct effects of the interpolated input tracer data and model parameterization on the simulated δ18OQ values. The step function interpolation generated an erratic isotopic signature in streamflow with flashy fluctuations (Fig. 4a–f). On the other hand, sine interpolation of δ18OP values yielded a smooth response in the simulated δ18OQ values (Fig. 4g–l). No significant visual distinction was found between kriged (Fig. 4a–c) and raw (Fig. 4d–f) δ18OP samples when the step function interpolation was used, except for a slightly larger uncertainty observed with raw δ18OP samples. Furthermore, when employing the sine interpolation and raw δ18OP values (Fig. 4j–l), the simulations overestimated the in-stream measurements in comparison with kriged values (Fig. 4g–i). Finally, distinct SAS parameterizations did not produce remarkable differences in the simulated δ18OQ values, although PLTV generally yielded simulations that better enveloped the measured isotopic signature (Fig. 4b, e, h, k).

Figure 4Predicted δ18O values in streamflow. The dark blue filled circles represent the observed data, and the dashed light blue lines and shaded areas represent the ensemble mean of all possible solutions and their range according to the 95 % CI, respectively.


Despite the differences in the predicted δ18OQ values, all simulations can be considered satisfactory given the KGE values ranging between 0.55 and 0.72 across all tested setups (Fig. 5). This aforementioned model performance can be classified as intermediate (Thiemig et al.2013) to good (Andersson et al.2017; Sutanudjaja et al.2018). When considering the best fit, the combination of step function interpolation and raw δ18OP values performed best. Additionally, PLTV generally yielded slightly better KGE values than PLTI or BETATI when grouping the setups with the same spatiotemporal interpolation of δ18OP. Differences in the mean KGEs were statistically insignificant (t test with p values > 0.05) only between setups a–g, c–i–k, and j–l (Table 1), as the mean KGE values were nearly identical; this largely agrees with the visual analysis (Fig. 5). Contrarily, the differences in the mean KGE values of the remaining setups were statistically significant (p values < 0.05), indicating that the a priori methodological choices (i.e. interpolation techniques of δ18OP values and/or SAS parameterization) strongly impact the overall results. Nonetheless, this does not mean that we can clearly identify the most suitable setup, but there is a need to carefully analyze the multiple potential choices with respect to SAS parameterization and tracer data interpolations as well as to evaluate the uncertainty range in modeled predictions.

Figure 5Box plots of model performance ranges in behavioral solutions. The letters on the x axis refer to the model setup type according to Table 1. Box plots filled with the same colors represent model setups characterized by the same interpolation scheme in space and time. On each box, the central red line indicates the median, and the bottom and top edges of the box indicate the 25th and 75th percentiles, respectively, namely the interquartile range (IQR). The whiskers extend to the most extreme data points not considered outliers, which are the 25th percentile  1.5 × IQR and the 75th percentile + 1.5 × IQR, respectively. The outliers are plotted individually using the red “+” markers.


Ranges of the behavioral SAS parameters for the tested setups are summarized in Table S1 in the Supplement. Parameters for the SAS functions of Q (i.e. kQ, kQ1, kQ2, α, and β) were different across the setups, although they were generally relatively narrow and well identified. However, the behavioral parameters were better constrained when using the step function interpolation, as their 95 % CI was, on average, 34 % narrower than that provided by the sine interpolation, across all the SAS parameterizations. The parameters kQ1 and α were also better identified than kQ2 and β, as their 95 % CI was, on average, 56 % narrower, across all tested setups. Conversely, there was no clear difference in the parameter ranges when using kriged or raw δ18OP values. The evapotranspiration parameter (i.e. kET) was poorly identified in all setups, as any value in the search range provided equally good results. The initial storage (i.e. S0) was only partially constrained, as any value between 335 and 2895 mm was considered acceptable.

4.2 Simulated transit times and model uncertainty

Figure 6 illustrates the 95 % CI of the behavioral solutions for the predicted median transit time (TT50). The results show that the model simulated largely different ranges of TT50 based on the tested setups. When using PLTI and BETATI (Fig. 6a, c, d, f, g, i, j, l), the 95 % CI was relatively stable with smaller fluctuations throughout the simulation period compared with PLTV (Fig. 6b, e, h, k). However, minor differences emerged across the simulated TT50 as a result of the distinct interpolation techniques used for δ18OP. The 95 % CI of TT50 was, on average, 37 % larger, across all tested setups, when using raw δ18OP (Fig. 6d–f and j–l) rather than kriged δ18OP (Fig. 6a–c and g–i). This was especially visible when the step function was used (Fig. 6a–f). Moreover, the sine interpolation generated a 95 % CI of TT50 that was, on average, 62 % narrower across all tested setups (Fig. 6g–l) with respect to the step function (Fig. 6a–f). These differences were more evident for high-flow conditions where the 95 % CI of TT50 showed a significant reduction. In addition, the behavioral solutions obtained with the sine interpolation (Fig. 6g–l) were more skewed towards shorter mean TT50 values, across all tested setups, than those of the step function (Fig. 6a–f).

Figure 6Predicted TT50 of streamflow. The dashed light blue lines and shaded areas represent the ensemble mean of all possible solutions and their range according to the 95 % CI, respectively.


Behavioral solutions obtained with PLTV revealed a similar pattern regardless of the interpolation employed (Fig. 6b, e, h, k). Nonetheless, there was a noticeable difference in the 95 % CI of TT50 under distinct flow regimes. During low flows and dry periods (i.e. summer and autumn), the time series of predicted TT50 showed large uncertainties ranging at most between 259 and 1009 d across the different setups (Fig. 6e). Conversely, during high flows (i.e. winter and spring), the 95 % CI was much narrower and varied at least between 126 and 154 d (Fig. 6h). The large 95 % CI and the notable differences across the tested setups highlight the sensitivity and, in turn, the uncertainty in the predicted TT50 to the model parameterization, temporal interpolation of input data, hydrologic conditions, and nonspatially interpolated δ18OP.

In general, the variability in the predicted TT50 was controlled by the hydrological state of the system (Fig. 6). High-discharge events reduced the TT50 values, whereas low-flow periods were associated with a longer estimated TT50. This is expected, as streamflow during high (low) flows is dominated by near-surface runoff (groundwater) with shallow (deep) flow paths leading to a shorter (longer) TT50. Such differences were particularly visible with PLTV (Fig. 6b, e, h, k), as the exponent kQ(t) shifts the water selection preference over time as a function of the wet/dry conditions. This resulted in the variability in TT50 being more pronounced than that of PLTI and BETATI, whose SAS parameters for Q are constant over time.

4.3 Catchment-scale water release

SAS functions provided valuable insights into the catchment-scale water release dynamics. Figure 7 presents the behavioral solutions releasing water of different ages and also shows that the catchment generally experienced a stronger affinity for releasing young water (i.e. kQ<1, or α<1 and β>1), rather than old water (i.e. kQ>1, or α>1 and β<1). These findings are in agreement with other studies in the upper Selke catchment (Winter et al.2020; Nguyen et al.2021). Nonetheless, there were differences in the water release scheme when comparing various combinations of SAS functions and spatiotemporal interpolation techniques of isotopes. The use of PLTV resulted in a substantial number of solutions, approximately 50 % of all behavioral solutions, suggesting a preference for both young and old water. On the other hand, only a few solutions showed affinity for old-water release, and this was more prominent when using the sine interpolation, raw δ18OP values, and PLTI across all tested setups.

Figure 7Percentage of behavioral solutions releasing water of different ages.


5 Discussion

5.1 Uncertainty in TTD modeling

In this study, we characterized the TTD uncertainty arising from some significant and critical aspects for the SAS modeling. These aspects are also the most directly linked to the data interpolation and SAS parameterization that we explored in this work. The uncertainty analysis was carried out across the 12 tested setups corresponding to different combinations of spatiotemporal data interpolation techniques and SAS parameterizations. Our results show that the uncertainty (i.e. 95 % CI) of the simulated TT50 (Fig. 6) was firmly dependent on the choice of model setup, as the 95 % CI was primarily sensitive to the type of SAS function, temporal interpolation, and spatial representation of δ18OP.

Uncertainty in the simulated TT50 differed considerably between time-invariant (i.e. PLTI and BETATI; Fig. 6a, c, d, f, g, i, j, l) and time-variant (i.e. PLTV; Fig. 6b, e, h, k) SAS functions; thus, a large sensitivity is associated with the choice of the SAS parameterization. For example, PLTI and BETATI explicitly assume constant water selection preference over time, as these functions do not consider the temporal variability in the catchment wetness. As a consequence, the resulting TT50 had a moderately stable 95 % CI with smaller fluctuations compared with those of PLTV. On the other hand, including an explicit time dependence in the SAS function strongly affected the 95 % CI of TT50. Notably, PLTV produced a wider 95 % CI during low-flow conditions, which can hinder the TTDs ability to provide robust insights into flow and solute transport behaviors in the study area during low-flow conditions. This highlights the need to further constrain PLTV with additional data, which could involve obtaining tracer data at a finer resolution or additional information on the evapotranspiration and initial storage. In addition, the exceptionally old flow components associated with a very large 95 % CI of TT50 might be a distortion of the actual TT50 values, which can usually be more reliably estimated using radioactive tracers than with stable isotopes (Visser et al.2019). Hence, PLTV-based TT50 greater than the observed period (828 d) should be interpreted carefully. It is important to note that we discussed the fractional SAS (fSAS) functions in this study, but another form of the SAS functions, such as the rank SAS (rSAS) functions, may have different uncertainty. This is mainly due to the difference in how the storage is considered: fSAS functions are expressed as function of the normalized age-ranked storage, which is equal to the cumulative residence time, whereas rSAS functions depend on the age-ranked storage, which is the volume of water in storage ranked from youngest to oldest (Harman2015).

Likewise, the high-frequency reconstruction of δ18OP from monthly samples via interpolation created further uncertainty. The sine interpolation effectively captured the dominant features of the observed δ18OP, such as seasonality. Consequently, sine interpolation successfully reproduced the seasonal trend in in-stream δ18O, although simulations overestimated the measurements (Fig. 4g–l). On the other hand, sine interpolation poorly reproduced rainfall isotopes during short-term flashy events and missed detailed characteristics of the tracer dataset by smoothing the variability in the observed δ18OP (Fig. 3). As a result, high values of δ18OP are underestimated, whereas low values are overestimated. It is critical to recognize these limitations when interpreting modeling results, as uncertainty in the simulated δ18OP may conceal a more pronounced hydrological response of the system (Dunn et al.2008; Birkel et al.2010; Hrachowitz et al.2011). Contrarily, step function interpolation preserved the maxima in the monthly observed δ18OP values by capturing their variation correctly (Fig. 3). Simulations showed a better fit with measured in-stream δ18O (Fig. 4a–f) and higher model performance (Fig. 5). However, combining step function with raw δ18OP resulted in larger uncertainty in the simulated TT50 (Fig. 6d–f). This reflects the need for a comprehensive exploration of the uncertainty range, rather than relying solely on the goodness of fit. Overall, the choice between step function and sine interpolation largely affected the reconstructed input data (Fig. 3), leading to significant differences in the simulated TT50 and associated uncertainty. It is important to note that alternative methods, such as generalized additive models (GAMs; Buzacott et al.2020), offer other options for interpolating tracer data. We conducted further tests with the SAS model using a GAM to reconstruct both kriged and raw δ18OP using smoothing functions; this provides a more sophisticated approach than the intuitive methods used in this study. However, the results, available in the Supplement, show that while a GAM provided more detailed reconstructed input tracer data (Fig. S1), it did not significantly alter the SAS-based results (Figs. S2, S3) or yield any new insights or conclusions about uncertainty with respect to using step function and sine interpolation. Therefore, we conclude that, while highly resolved input data may seem appealing, they do not necessarily lead to substantial benefits for the SAS-based output, supposedly due to the conceptual simplifications in the SAS model.

The spatial representation of δ18OP values had limited impact on the overall pattern of simulated TT50, as the time series were comparable with both kriged (Fig. 6a–c and g–i) and raw (Fig. 6d–f and j–l) δ18OP. Nonetheless, the spatial interpolation of δ18OP from different locations resulted in a reduction in the uncertainty in the TT50, which was particularly evident with the step function (Fig. 6a–f). This difference may be attributed to the fact that the upper Selke is a large (mesoscale) catchment with a substantial gradient in elevation, and, as a consequence, measurement for δ18OP from only one location may be generally overly simplistic. This finding highlights the importance of considering not only the model performance (Fig. 5; raw values with a step function interpolation produced higher KGE values) but also the uncertainty range in predicted TT50.

Finally, we found that the uncertainty was larger under dry conditions when lower flow and longer TT50 were observed. This was especially visible when using the time-variant SAS function (Fig. 6b, e, h, k). It might be due to the fact that, under wet conditions, there is a high level of hydrologic connectivity within the catchment (Ambroise2004; Blume and van Meerveld2015; Hrachowitz et al.2016), which results in nearly all flow paths being active and contributing to the streamflow. This, ultimately, may make TT50 values easier to constrain. Conversely, under dry conditions, when usually only longer flow paths carrying older water are active (Soulsby and Tetzlaff2008; Jasechko et al.2017), water partially flows through a drier soil zone where flow is more erratic (i.e. flow directions and patterns can vary widely) as the conductivity is controlled by soil moisture. As a result, wet areas can be patchy and water flows preferentially at certain locations only, as opposed to spatially uniform flow through the soil matrix; this might make it more challenging to constrain older water ages. Similarly, Benettin et al. (2017) found higher uncertainty in the simulated SAS-based median water ages during drier periods, potentially due to higher uncertainty in the total storage. Moreover, non-SAS function studies have observed major uncertainties and deviations from observations in lumped modeled results during low-flow conditions (Kumar et al.2010). This was primarily due to the lack of spatial variability in the catchment characteristics in lumped models, which is a critical factor controlling low-flow regimes in rivers.

The dissimilarities in the simulated TT50 across the tested setups underline the importance of accounting for uncertainty in model-based TTDs. The uncertainty analysis with SUFI-2 performed in this study was essential to best describe the parameter identifiability and bounds of the behavioral solutions of each output variable. Furthermore, our results highlight the importance of gaining tracer datasets of good quality (i.e. tracer data with a finer temporal resolution), considering the spatial variability in the isotopic composition in precipitation, and, possibly, employing a model parameterization that best describes the catchment-specific storage and release dynamics. The last point can be defined according to a precise conceptual knowledge of the catchment’s functioning and information from previous studies in similar catchments.

5.2 TTD modeling: advantages and limitations

Our results provide visually plausible seasonal fluctuations in the predicted δ18OQ samples (Fig. 4) and satisfactory KGE values (Fig. 5), despite the uncertainty arising from model inputs, structure, and parameters. The good match with observations provides confidence in the simulated TT50 for the upper Selke catchment. The magnitude of the uncertainty resulting from different setups cannot be generalized, but the overall approach for uncertainty assessment presented here could be extended to other areas and TTD studies. However, we recognize some limitations and indicate below possible reasons and, in turn, improvements that future work could achieve.

First, the limited length of the δ18O time series might not describe the system accurately; hence, implementing longer time series could improve the parameter identifiability and provide a more accurate estimation of the TTDs. Second, this study relied on stable water isotopes, which might underestimate the tails of the TTDs (Stewart et al.2010; Seeger and Weiler2014), although recent works have contested this (Wang et al.2022). Possible advancements could be reached by using decaying tracers varying over a longer timescale than stable water isotopes (e.g. tritium; Stewart et al.2012; Morgenstern et al.2015) and imparting more information on old water. Next, future work should retrieve more information on the evapotranspiration ET and initial storage S0, whose parameters were poorly identified. However, this issue is common in transport studies that rely on measurements of in-stream stable water isotopes (Benettin et al.2017; Buzacott et al.2020). As a way forward, information on the ET isotopic compositions might help better constrain ET parameters and assess their affinity for young/old water. Regarding constraining the range of S0, further information can be gained from geophysical surveys in the study areas or groundwater modeling as well as by using decaying isotopes (Visser et al.2019).

5.3 Implications of TTD uncertainties

This study characterized the uncertainty in TTDs, which summarize the catchment's hydrologic transport behavior and, therefore, comprise decisive information for water managers. The value of TT50 has relevant implications for both water quantity and quality, as does its uncertainty. The larger the 95 % CI in the simulated TT50, the greater the difference in the TT50 values, which, ultimately, implies distinct hydrological processes, water availability, groundwater recharge, and solute export dynamics (McDonnel et al.2010).

For example, knowing the TTD and its uncertainty may be crucial for characterizing the catchment's response to climatic change (Wilusz et al.2017). Considering the increasing severity of droughts in the past decades (Dai2013), a catchment with a shorter TT50 and a dominant release of young water might be more affected by droughts than a catchment with a longer TT50, which means that its stream is fed by relatively old water sources. Therefore, a short TT50 reveals a low drought resilience of the catchment and limited water availability, which could limit streamflow generation processes and change the in-stream water quality status during drought periods (Winter et al.2023). Likewise, TTD uncertainty may affect the understanding of the water infiltration rate, hydrological processes, and aquifer recharge, as a shorter TT50 suggests that water is quickly routed to the catchment outlet rather than infiltrating deeply into the groundwater. Finally, TTD uncertainty can have an impact on the quantification of the modern groundwater age, i.e. groundwater younger than 50 years (Bethke and Johnson2008). According to Jasechko (2019), the correct identification of modern groundwater abundance and distribution can help determine its renewal (Le Gal La Salle et al.2001; Huang et al.2017), groundwater wells and depths most likely to contain contaminants (Visser et al.2013; Opazo et al.2016), and the part of the aquifer flushed more rapidly.

Uncertainty in TTDs also impacts on assessing the fate of dissolved solutes, such as nitrates (X. Yang et al.2018; Nguyen et al.2021, 2022; Lutz et al.2022), pesticides (Holvoet et al.2007; Lutz et al.2017), and chlorides (Kirchner et al.2000; Benettin et al.2013). These solutes constitute a crucial source of diffuse water pollution in agricultural areas (Jiang et al.2014; Kumar et al.2020), as they are spread on the soil in large quantities, especially during the growing season. The exposure time of solutes with the soil matrix has strong consequences for biogeochemical reactions, such as denitrification in the case of nitrates (Kolbe et al.2019; Kumar et al.2020). A short TT50 suggests that water can be rapidly conveyed to the stream network (Kirchner et al.2001), with limited time for denitrification. This explains the elevated in-stream concentration and short-term impact of nitrate export compared with that of a longer TT50, which is typically associated with old-water release and a low nitrate concentration (Nguyen et al.2021). Similarly, pesticide transport is highly affected by the TTD uncertainty, as a long TT50 suggests little pesticide degradation due to decreased microbial activity along deeper flow paths (Rodríguez-Cruz et al.2006). In other cases, a shorter TT50 may limit the time for degradation, causing a peak in the in-stream concentration (Leu et al.2004). Overall, a longer TT50 can delay or buffer the catchment's reactive solute response at the outlet (Dupas et al.2016; Van Meter et al.2017). This creates a long-term effect of hydrological legacies and a continuous problem with diffuse pollution of nitrates (Ehrhardt et al.2019; Winter et al.2020) and pesticides (Lutz et al.2013), which can persist in the catchment for several years. Finally, TTD uncertainties also play an important role in chloride transport, although chlorides are commonly known to be conservative (Svensson et al.2012). A short TT50 may indicate rapid chloride mobilization, whereas a long TT50 implies chloride persistence in groundwater; therefore, chloride accumulates and is released at lower rates, with impacts on the ecosystem functions, vegetation uptake, and metabolism (Xu et al.1999).

Understanding the uncertainty in TTDs is crucial for the aforementioned implications. While previous studies have used only a specific SAS function and/or specific tracer data interpolation technique in time and space, here we show that there could be a wide range of different results in terms of water ages, model performance, and parameter uncertainty. This is due to the specific choice regarding SAS parameterization and tracer data interpolation. With this, we want to convey that uncertainty is omnipresent in TTD-based models, and we need to recognize it, especially when dealing with sparse tracer data and multiple choices for model parameterization. Therefore, we want to encourage future studies to explore these uncertainties in other catchments and different geophysical settings, with the final aim to investigate whether these uncertainties may affect the conclusions of water quantity and quality studies for management purposes.

6 Conclusions

This study explored the uncertainty in TTDs of streamflow, resulting from 12 model setups obtained from different SAS parameterizations (i.e. PLTI, PLTV, and BETATI), and reconstruction of the precipitation isotopic signature in time and space via interpolation (step function vs. sine fit and raw vs. kriged values).

We found satisfactory KGE values, whose differences across the tested setups were statistically significant, meaning that the choice of the setup matters. As a consequence, distinct setups led to considerably different simulated TT50 values. The choice between using time-variant or time-invariant SAS functions was crucial, as the time-invariant functions generated moderate fluctuations in the 95 % CI of the estimated TT50 because of the constant water selection preference over time. On the other hand, the time-variant SAS function captured the dynamics of the catchment wetness, resulting in more pronounced fluctuations in TT50. However, the time-variant SAS function also produced a larger 95 % CI in TT50, notably during drier periods, which might indicate the need to constrain the function with additional data (e.g. finer tracer data resolution and/or information on evapotranspiration and storage). Significant differences in TT50 were observed depending on the employed temporal interpolations. Results from the sine interpolation produced a smaller uncertainty in TT50, with the time series skewed towards smaller values. However, such results must be interpreted carefully, as the sine interpolation poorly reproduced flashy events in precipitation, thus indicating that some more dynamic transport processes were not fully considered. Conversely, the step function interpolation resulted in a larger uncertainty in the TT50, but it was able to better reproduce the measured δ18OP data by capturing the peak values, as opposed to the sine interpolation. Dry conditions were another reason for uncertainty, as indicated by the high variance in the simulated TT50 values, which is potentially attributed to the water preferentially moving at certain locations, making wet areas patchy, so it may be more challenging to accurately constrain older water ages. Finally, there was comparable pattern in the modeled results when using kriged vs. raw isotopes, but the kriged values yielded an uncertainty reduction in TT50. This highlights the potential advantage of spatially interpolated values over a single measurement representative of the entire area, particularly in a mesoscale catchment that varies with respect to elevation.

These findings provide new insights into TTD uncertainty when high-frequency tracer data are missing and the SAS framework is used. Regardless of the degree of efficiency or uncertainty, the decision on which setup is more plausible depends on the best conceptual knowledge of the catchment functioning. We consider the presented approach to be potentially applicable to other studies to enable a better characterization of TTD uncertainty, improve TTD simulations and, ultimately, inform water management. These aspects are particularly crucial in view of increasingly extreme climatic conditions and worsening water pollution under global change.

Code and data availability

The version of the model used in this study (v1.0) and its documentation are available at (last access: August 2020) and (Benettin and Bertuzzo2018b). The iteratively re-weighted least-squares (IRLS) method used to obtain the modeled daily kriged and raw isotope (δ18O) in precipitation information with the sine interpolation is presented at (von Freyberg et al.2018). The hydroclimatic time series, δ18O data, and interpolated δ18O time series can be accessed at (Borriero2022).


The supplement related to this article is available online at:

Author contributions

AB conducted the model simulations, carried out the analysis, interpreted the results, and wrote the original draft of the paper. SRL and RK designed and conceptualized the study and provided data for model simulations. TVN provided technical support for modeling and helped organize the structure and content of the paper. AB, SRL, RK, and TVN conceived the methodology and experimental design. All co-authors helped AB interpret the results. All authors contributed to the review, final writing, and finalization of this work.

Competing interests

At least one of the (co-)authors is a member of the editorial board of Hydrology and Earth System Sciences. The peer-review process was guided by an independent editor, and the authors also have no other competing interests to declare.


Publisher’s note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


The authors thank the Helmholtz Centre for Environmental Research – UFZ of the Helmholtz Association for funding this study; the German Weather Service and the mHM model team for providing the necessary input data; and the developers of the tran-SAS model and IRLS code for making them publicly available.

Financial support

The article processing charges for this open-access publication were covered by the Helmholtz Centre for Environmental Research – UFZ.

Review statement

This paper was edited by Laurent Pfister and reviewed by Ciaran Harman and one anonymous referee.


Abbaspour, K. C., Johnson, C. A., and van Genuchten, M. T.: Estimating uncertain flow and transport parameters using a sequential uncertainty fitting procedure, Vadose Zone J., 3, 1340–1352,, 2004. a, b

Ajami, N. K., Duan, Q., and Sorooshian, S.: An integrated hydrologic Bayesian multimodel combination framework: Confronting input, parameter, and model structural uncertainty in hydrologic prediction, Water Resour. Res., 43, W01403,, 2007. a

Allen, S. T., Jasechko, S., Berghuijs, W. R., Welker, J. M., Goldsmith, G. R., and Kirchner, J. W.: Global sinusoidal seasonality in precipitation isotopes, Hydrol. Earth Syst. Sci., 23, 3423–3436,, 2019. a

Ambroise, B.: Variable “active” versus “contributing” areas or periods: a necessary distinction, Hydrol. Process., 18, 1149–1155,, 2004. a, b

Andersson, J. C. M., Arheimer, B., Traoré, F., Gustafsson, D., and Ali, A.: Process refinements improve a hydrological model concept applied to the Niger River basin, Hydrol. Process., 31, 4540–4554,, 2017. a

Asadollahi, M., Stumpp, C., Rinaldo, A., and Benettin, P.: Transport and water age dynamics in soils: A comparative study of spatially integrated and spatially explicit models, Water Resour. Res., 56, e2019WR025539,, 2020. a

Benettin, P. and Bertuzzo, E.: tran-SAS v1.0: a numerical model to compute catchment-scale hydrologic transport using StorAge Selection functions, Geosci. Model Dev., 11, 1627–1639,, 2018a. a, b, c, d, e

Benettin, P. and Bertuzzo, E.: tran-SAS v1.0 (1.0), Zenodo [code],, 2018b. a

Benettin, P., van der Velde, Y., van der Zee, S. E. A. T. M., Rinaldo, A., and Botter, G.: Chloride circulation in a lowland catchment and the formulation of transport by travel time distributions, Water Resour. Res., 49, 4619–4632,, 2013. a

Benettin, P., Bailey, S. W., Campbell, J. L., Green, M. B., Rinaldo, A., Likens, G. E., McGuire, J. K., and Botter, G.: Linking water age and solute dynamics in streamflow at the Hubbard Brook Experimental Forest, NH, USA, Water Resour. Res., 51, 9256–9272,, 2015a. a

Benettin, P., Kirchner, J. W., Rinaldo, A., and Botter, G.: Modeling chloride transport using travel time distributions at Plynlimon, Wales, Water Resour. Res., 51, 3259–3276,, 2015b. a

Benettin, P., Soulsby, C., Birkel, C., Tetzlaff, D., Botter, G., and Rinaldo, A.: Using SAS functions and high-resolution isotope data to unravel travel time distributions in headwater catchments, Water Resour. Res., 53, 1864–1878,, 2017. a, b, c, d, e, f

Bethke, C. M. and Johnson, T. M.: Groundwater age and groundwater age dating, Annu. Rev. Earth Planet. Sci., 36, 121–152,, 2008. a

Beven, K.: A manifesto for the equifinality thesis, J. Hydrol., 320, 18–36,, 2006. a

Beven, K. and Freer, J.: Equifinality, data assimilation, and uncertainty estimation in mechanistic modelling of complex environmental systems using the GLUE methodology, J. Hydrol., 249, 11–29,, 2001. a

Birkel, C. and Soulsby, C.: Advancing tracer-aided rainfall-runoff modelling: A review of progress, problems and unrealised potential, Hydrol. Process., 29, 5227–5240,, 2015. a

Birkel, C., Dunn, S. M., Tetzlaff, D., and Soulsby, C.: Assessing the value of high-resolution isotope tracer data in the stepwise development of a lumped conceptual rainfall–runoff model, Hydrol. Process., 24, 2335–2348,, 2010. a

Blume, T. and van Meerveld, H. J.: From hillslope to stream: methods to investigate subsurface connectivity, WIREs Water, 2, 177–198,, 2015. a

Borriero, A.: Hydroclimatic and isotope data – Upper Selke, Zenodo [data set],, 2022. a

Botter, G., Bertuzzo, E., Bellin, A., and Rinaldo, A.: On the Lagrangian formulations of reactive solute transport in the hydrologic response, Water Resour. Res., 41, W04008,, 2005. a

Botter, G., Bertuzzo, E., and Rinaldo, A.: Transport in the hydrologic response: Travel time distributions, soil moisture dynamics, and the old water paradox, Water Resour. Res., 46, W03514,, 2010. a

Botter, G., Bertuzzo, E., and Rinaldo, A.: Catchment residence and travel time distributions: The master equation, Geophys. Res. Lett., 38, L11403,, 2011. a

Buzacott, A. J. V., van der Velde, Y., Keitel, C., and Vervoort, R. W.: Constraining water age dynamics in a south-eastern Australian catchment using an age-ranked storage and stable isotope approach, Hydrol. Process., 34, 4384–4403,, 2020. a, b, c

Dai, A.: Erratum: Increasing drought under global warming in observations and models, Nat. Clim. Change, 3, 171,, 2013. a

Danesh-Yazdi, M., Foufoula-Georgiou, E., Karwan, D. L., and Botter, G.: Inferring changes in water cycle dynamics of intensively managed landscapes via the theory of time-variant travel time distributions, Water Resour. Res., 52, 7593–7614,, 2016. a

Danesh-Yazdi, M., Klaus, J., Condon, L. E., and Maxwell, R. M.: Bridging the gap between numerical solutions of travel time distributions and analytical storage selection functions, Hydrol. Process., 32, 1063–1076,, 2018. a

Drever, M. C. and Hrachowitz, M.: Migration as flow: using hydrological concepts to estimate the residence time of migrating birds from the daily counts, Methods Ecol. Evol., 8, 1146–1157,, 2017. a, b

Dunn, S. M., Bacon, J. R., Soulsby, C., Tetzlaff, D., Stutter, M. I., Waldron, S., and Malcolm, I. A.: Interpretation of homogeneity in δ18O signatures of stream water in a nested sub-catchment system in north-east Scotland, Hydrol. Process., 22, 4767–4782,, 2008. a

Dupas, R., Jomaa, S., Musolff, A., Borchardt, D., and Rode, M.: Disentangling the influence of hydroclimatic patterns and agricultural management on river nitrate dynamics from sub-hourly to decadal time scales, Sci. Total Environ., 571, 791–800,, 2016. a

Dupas, R., Musolff, A., Jawitz, J. W., Rao, P. S. C., Jäger, C. G., Fleckenstein, J. H., Rode, M., and Borchardt, D.: Carbon and nutrient export regimes from headwater catchments to downstream reaches, Biogeosciences, 14, 4391–4407,, 2017. a

Ehrhardt, S., Kumar, R., Fleckenstein, J. H., Attinger, S., and Musolff, A.: Trajectories of nitrate input and output in three nested catchments along a land use gradient, Hydrol. Earth Syst. Sci., 23, 3503–3524,, 2019. a

Feng, X., Faiia, A. M., and Posmentier, E. S.: Seasonality of isotopes in precipitation: A global perspective, J. Geophys. Res, 114, D08116,, 2009. a

Gupta, H. V., Kling, H., Yilmaz, K. K., and Martinez, G. F.: Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377, 80–91,, 2009. a

Harman, C. J.: Time-variable transit time distributions and transport: Theory and application to storage-dependent transport of chloride in a watershed, Water Resour. Res., 51, 1–30,, 2015. a, b, c

Harman, C. J.: Age-Ranked Storage-Discharge Relations: A Unified Description of Spatially Lumped Flow and Water Age in Hydrologic Systems, Water Resour. Res., 55, 7143–7165,, 2019. a

Heidbüchel, I., Troch, P. A., and Lyon, S. W.: Separating physical and meteorological controls of variable transit times in zero-order catchments, Water Resour. Res., 49, 7644–7657,, 2013. a

Heidbüchel, I., Yang, J., Musolff, A., Troch, P., Ferré, T., and Fleckenstein, J. H.: On the shape of forward transit time distributions in low-order catchments, Hydrol. Earth Syst. Sci., 24, 2895–2920,, 2020. a

Holvoet, K. M. A., Seuntjens, P., and Vanrolleghem, P. A.: Monitoring and modeling pesticide fate in surface waters at the catchment scale, Ecol. Model., 209, 53–64,, 2007. a

Hrachowitz, M., Soulsby, C., Tetzlaff, D., Malcolm, I. A., and Schoups, G.: Gamma distribution models for transit time estimation in catchments: Physical interpretation of parameters and implications for time-variant transit time assessment, Water Resour. Res., 46, W10536,, 2010. a

Hrachowitz, M., Soulsby, C., Tetzlaff, D., and Malcolm, I. A.: Sensitivity of mean transit time estimates to model conditioning and data availability, Hydrol. Process., 25, 980–990,, 2011. a

Hrachowitz, M., Savenije, H., Bogaard, T. A., Tetzlaff, D., and Soulsby, C.: What can flux tracking teach us about water age distribution patterns and their temporal dynamics?, Hydrol. Earth Syst. Sci., 17, 533–564,, 2013. a

Hrachowitz, M., Benettin, P., van Breukelen, B. M., Fovet, O., Howden, N. J. K., Ruiz, L van der Velde, Y., and Wade, A. J.: Transit times – the link between hydrology and water quality at the catchment scale, WIREs Water, 3, 629–657,, 2016. a, b

Huang, T., Pang, Z., Li, J., Xiang, Y., and Zhao, Z.: Mapping groundwater renewability using age data in the Baiyang alluvial fan, NW, China, Hydrogeol. J., 25, 743–755,, 2017. a

Jasechko, S.: Global isotope hydrogeology – review, Rev. Geophys., 57, 835–965,, 2019. a

Jasechko, S., Wassenaar, L. I., and Mayer, B.: Isotopic evidence for widespread cold-season-biased groundwater recharge and young streamflow across central Canada, Hydrol. Process., 31, 2196–2209,, 2017. a

Jiang, S., Jomaa, S., and Rode, M.: Modelling inorganic nitrogen leaching in nested mesoscale catchments in central Germany, Ecohydrol., 7, 1345–1362,, 2014. a

Jing, M., Heße, F., Kumar, R., Kolditz, O., Kalbacher, T., and Attinger, S.: Influence of input and parameter uncertainty on the prediction of catchment-scale groundwater travel time distributions, Hydrol. Earth Syst. Sci., 23, 171–190,, 2019. a

Kim, M. and Troch, P. A.: Transit time distributions estimation exploiting flow-weighted time: Theory and proof-of-concept, Water Resour. Res., 56, e2020WR027186,, 2020. a

Kim, M., Pangle, L. A., Cardoso, C., Lora, M., Volkmann, T. H. M., Wang, Y., Harman, C. J., and Troch, P. A.: Transit time distributions and StorAge Selection functions in a sloping soil lysimeter with time-varying flow paths: Direct observation of internal and external transport variability, Water Resour. Res., 52, 7105–7129,, 2016. a

Kirchner, J. W.: Getting the right answers for the right reasons: Linking measurements, analyses, and models to advance the science of hydrology, Water Resour. Res., 42, W03S04,, 2006. a

Kirchner, J. W.: Aggregation in environmental systems – Part 1: Seasonal tracer cycles quantify young water fractions, but not mean transit times, in spatially heterogeneous catchments, Hydrol. Earth Syst. Sci., 20, 279–297,, 2016. a

Kirchner, J. W.: Quantifying new water fractions and transit time distributions using ensemble hydrograph separation: theory and benchmark tests, Hydrol. Earth Syst. Sci., 23, 303–349,, 2019. a

Kirchner, J. W., Feng, X., and Neal, C.: Fractal stream chemistry and its implications for contaminant transport in catchments, Nature, 403, 524–527,, 2000. a

Kirchner, J. W., Feng, X., and Neal, C.: Catchment-scale advection and dispersion as a mechanism for fractal scaling in stream tracer concentrations, J. Hydrol., 254, 82–101,, 2001. a, b

Kirchner, J. W., Feng, X., Neal, C., and Robson, A. J.: The fine structure of water-quality dynamics: the (high-frequency) wave of the future, Hydrol. Process., 18, 1353–1359,, 2004. a

Kolbe, T., de Dreuzy, J. R., Abbot, B. W., Aquilina, L., Babey, T., Green, C. T., Fleckenstein, J. H., Labasque, T., Laverman, A. M., Marçais, J., Peiffer, S., Thomas, Z., and Pinay, G.: Stratification of reactivity determines nitrate removal in groundwater, P. Natl. Acad. Sci. USA, 116, 2494–2499,, 2019. a

Kumar, R., Samaniego, L., and Attinger, S.: The effects of spatial discretization and model parameterization on the prediction of extreme runoff characteristics, J. Hydrol., 392, 54–69,, 2010. a

Kumar, R., Samaniego, L., and Attinger, S.: Implications of distributed hydrologic model parameterization on water fluxes at multiple scales and locations, Water Resour. Res., 49, 360–379,, 2013. a

Kumar, R., Heße, F., Rao, P. S. C., Musolff, A., Jawitz, J. W., Sarrazin, F., Samaniego, L., Fleckenstein, J. H., Rakovec, O., Thober, S., and Attinger, S.: Strong hydroclimatic controls on vulnerability to subsurface nitrate contamination across Europe, Nat. Commun., 11, 6302,, 2020. a, b, c

Le Gal La Salle, C., Marlin, C., Leduc, C., Taupin, J. D., Massault, M., and Favreau, G.: Renewal rate estimation of groundwater based on radioactive tracers (3H, 14C) in an unconfined aquifer in a semi-arid area, Iullemeden Basin, Niger, J. Hydrol., 254, 145–156,, 2001. a

Leu, C., Singer, H., Stamm, C., Muller, S. R., and Schwarzenbach, R. P.: Simultaneous Assessment of Sources, Processes, and Factors Influencing Herbicide Losses to Surface Waters in a Small Agricultural Catchment, Environ. Sci. Technol., 38, 3827–3834,, 2004. a

Lutz, S. R., van Meerveld, H. J., Waterloo, M. J., Broers, H. P., and van Breukelen, B. M.: A model-based assessment of the potential use of compound-specific stable isotope analysis in river monitoring of diffuse pesticide pollution, Hydrol. Earth Syst. Sci., 17, 4505–4524,, 2013. a

Lutz, S. R., Velde, Y. V. D., Elsayed, O. F., Imfeld, G., Lefrancq, M., Payraudeau, S., and van Breukelen, B. M.: Pesticide fate on catchment scale: conceptual modelling of stream CSIA data, Hydrol. Earth Syst. Sci., 21, 5243–5261,, 2017. a, b

Lutz, S. R., Krieg, R., Müller, C., Zink, M., Knöller, K., Samaniego, L., and Merz, R.: Spatial patterns of water age: using young water fractions to improve the characterization of transit times in contrasting catchments, Water Resour. Res., 54, 4767–4784,, 2018. a, b, c

Lutz, S. R., Ebeling, P., Musolff, A., Nguyen, T. V., Sarrazin, F. J., Van Meter, K. J., Basu, N. B., Fleckenstein, J. H., Attinger, S., and Kumar, R.: Pulling the rabbit out of the hat: Unravelling hidden nitrogen legacies in catchment-scale water quality models, Hydrol. Processes, 36, e14682,, 2022. a

McDonnell, J. J., McGuire, K., Aggarwal, P., Beven, K. J., Biondi, D., Destouni, G., Dunn, S., James, A., Kirchner, J., Kraft, P., Lyon, S., Maloszewski, P., Newman, B., Pfister, L., Rinaldo, A., Rodhe, A., Sayama, T., Seibert, J., Solomon, K., Soulsby, C., Stewart, M., Tetzlaff, D., Tobin, C., Troch, P., Weiler, M., Western, A., Wörman, A., and Wrede, S.: How old is streamwater? Open questions in catchment transit time conceptualization, modelling and analysis, Hydrol. Process., 24, 1745–1754,, 2010. a, b

McGuire, K. J. and McDonnel, J. J.: A review and evaluation of catchment transit time modeling, J. Hydrol., 330, 543–563,, 2006. a, b, c

McKay, M. D., Beckman, R. J., and Conover, W. J.: A Comparison of Three Methods for Selecting Values of Input Variables in the Analysis of Output from a Computer Code, Technometrics, 21, 239–245,, 1979. a

Morgenstern, U., Daughney, C. J., Leonard, G., Gordon, D., Donath, F. M., and Reeves, R.: Using groundwater age and hydrochemistry to understand sources and dynamics of nutrient contamination through the catchment into Lake Rotorua, New Zealand, Hydrol. Earth Syst. Sci., 19, 803–822,, 2015. a

Nguyen, T. V., Kumar, R., Lutz, S. R., Musolff, A., Yang, J., and Fleckenstein, J. H.: Modeling Nitrate Export From a Mesoscale Catchment Using StorAge Selection Functions, Water Resour. Res., 57, e2020WR028490,, 2021. a, b, c, d, e

Nguyen, T. V., Sarrazin, F. J., Ebeling, P., Musolff, A., Fleckenstein, J. H., and Kumar, R.: Toward Understanding of Long-Term Nitrogen Transport and Retention Dynamics Across German Catchments, Geophys. Res. Lett., 49, e2022GL100278,, 2022. a

Niemi, A. J.: Residence time distributions of variable flow processes, Int. J. Appl. Radiat. Is., 28, 855–860,, 1977. a

Opazo, T., Aravena, R., and Parker, B.: Nitrate distribution and potential attenuation mechanisms of a municipal water supply bedrock aquifer, Appl. Geochem., 73, 157–168,, 2016. a

Queloz, P., Carraro, L., Benettin, P., Botter, G., Rinaldo, A., and Bertuzzo, E.: Transport of fluorobenzoate tracers in a vegetated hydrologic control volume: 2. Theoretical inferences and modeling, Water Resour. Res, 51, 2793–2806,, 2015. a, b

Rinaldo, A. and Marani, M.: Basin Scale Model of Solute Transport, Water Resour. Res., 23, 2107–2118,, 1987. a

Rinaldo, A., Botter, G., Bertuzzo, E., Uccelli, A., Settin, T., and Marani, M.: Transport at basin scales: 1. Theoretical framework, Hydrol. Earth Syst. Sci., 10, 19–29,, 2006. a

Rinaldo, A., Benettin, P., Harman, C. J., Hrachowitz, M., McGuire, K. J., ven der Velde, Y., Bertuzzo, E., and Botter, G.: Storage selection functions: A coherent framework for quantifying how catchments store and release water and solutes, Water Resour. Res., 51, 4840–4847,, 2015. a

Rodriguez, N. B., McGuire, K. J., and Klaus, J.: Time-Varying Storage–Water Age Relationships in a Catchment With a Mediterranean Climate, Water Resour. Res., 54, 3988–4008,, 2018. a

Rodriguez, N. B., Pfister, L., Zehe, E., and Klaus, J.: A comparison of catchment travel times and storage deduced from deuterium and tritium tracers using StorAge Selection functions, Hydrol. Earth Syst. Sci., 25, 401–428,, 2021. a

Rodríguez-Cruz, M. S., Jones, J. E., and Bending, G. D.: Field-scale study of the variability in pesticide biodegradation with soil depth and its relationship with soil characteristics, Soil Biol. Biochem., 38, 2910–2918,, 2006. a

Samaniego, L., Kumar, R., and Attinger, S.: Multiscale parameter regionalization of a grid-based hydrologic model at the mesoscale, Water Resour. Res., 46, W05523,, 2010. a

Schoups, G., van de Giesen, N. C., and Savenije, H. H. G.: Model complexity control for hydrologic prediction, Water Resour. Res., 44, W00B03,, 2008. a

Seeger, S. and Weiler, M.: Reevaluation of transit time distributions, mean transit times and their relation to catchment topography, Hydrol. Earth Syst. Sci., 18, 4751–4771,, 2014. a

Soulsby, C. and Tetzlaff, D.: Towards simple approaches for mean residence time estimation in ungauged basins using tracers and soil distributions, J. Hydrol., 363, 60–74,, 2008. a

Stewart, M. K., Morgenstern, U., and McDonnell, J. J.: Truncation of stream residence time: How the use of stable isotopes has skewed our concept of streamwater age and origin, Hydrol. Process., 24, 1646–1659,, 2010. a

Stewart, M. K., Morgenstern, U., McDonnell, J. J., and Pfister, L.: The “hidden streamflow” challenge in catchment hydrology: A call to action for stream water transit time analysis, Hydrol. Process., 26, 2061–2066,, 2012. a

Stockinger, M. P., Lücke, A., McDonnell, J. J., Diekkrüger, B., Vereecken, H., and Bogena, H. R.: Interception effects on stable isotope driven streamwater transit time estimates, Geophys. Res. Lett., 42, 5299–5308,, 2015. a

Sutanudjaja, E. H., van Beek, R., Wanders, N., Wada, Y., Bosmans, J. H. C., Drost, N., van der Ent, R. J., de Graaf, I. E. M., Hoch, J. M., de Jong, K., Karssenberg, D., López López, P., Peßenteiner, S., Schmitz, O., Straatsma, M. W., Vannametee, E., Wisser, D., and Bierkens, M. F. P.: PCR-GLOBWB 2: a 5 arcmin global hydrological and water resources model, Geosci. Model Dev., 11, 2429–2453,, 2018. a

Svensson, T., Lovett, G. M., and Likens, G. E.: Is chloride a conservative ion in forest ecosystems?, Biogeochemistry, 107, 125–134,, 2012. a

Tetzlaff, D., Piovano, T., Ala-Aho, P., Smith, A., Carey, S. K., Marsh, P., Wookey, P. A., Street, L. E., and Soulsby, C.: Using stable isotopes to estimate travel times in a data-sparse Arctic catchment: Challenges and possible solutions, Hydrol. Process., 32, 1936–1952,, 2018. a

Thiemig, V., Rojas, R., Zombrano-Bigiarini, M., and De Roo, A.: Hydrological evaluation of satellite-based rainfall estimates over the Volta and Baro-Akobo Basin, J. Hydrol., 499, 324–338,, 2013. a

Van der Velde, Y., De Rooij, G. H., Rozemeijer, J. C., Van Geer, F. C., and Broers, H. P.: Nitrate response of a lowland catchment: On the relation between stream concentration and travel time distribution dynamics, Water Resour. Res., 46, W11534,, 2010. a

van der Velde, Y., Torfs, P. J. J. F., van der Zee, S. E. A. T. M., and Uijlenhoet, R.: Quantifying catchment-scale mixing and its effect on time-varying travel time distributions, Water Resour. Res., 48, W06536,, 2012. a

Van Meter, K. J., Basu, N. B., and Van Cappellen, P.: Two centuries of nitrogen dynamics: legacy sources and sinks in the Mississippi and Susquehanna River Basins, Global Biogeochem. Cy., 31, 2–23,, 2017. a

Visser, A., Broers, H. P., Purtschert, R., Sültenfuß, J., and de Jonge, M.: Groundwater age distributions at a public drinking water supply well field derived from multiple age tracers (85Kr, 3H /3He, and 39Ar), Water Resour. Res., 49, 7778–7796,, 2013. a

Visser, A., Thaw, M., Deinhart, A., Bibby, R., Safeeq, M., Conklin, M., Esser, B., and van der Velde, Y.: Cosmogenic isotopes unravel the hydrochronology and water storage dynamics of the Southern Sierra critical zone, Water Resour. Res., 55, 1429–1450,, 2019. a, b

von Freyberg, J., Studer, B., and Kirchner, J. W.: A lab in the field: high-frequency analysis of water quality and stable isotopes in stream water and precipitation, Hydrol. Earth Syst. Sci., 21, 1721–1739,, 2017. a

von Freyberg, J., Allen, S. T., Seeger, S., Weiler, M., and Kirchner, J. W.: Sensitivity of young water fractions to hydro-climatic forcing and landscape properties across 22 Swiss catchments, Hydrol. Earth Syst. Sci., 22, 3841–3861,, 2018. a, b

von Freyberg, J., Rücker, A., Zappa, M., Schlumpf, A., Studer, B., and Kirchner, J. W.: Four years of daily stable water isotope data in stream water and precipitation from three Swiss catchments, Sci. Data, 9, 46,, 2022. a

Wang, S., Hrachowitz, M., Schoups, G., and Stumpp, C.: Stable water isotopes and tritium tracers tell the same tale: No evidence for underestimation of catchment transit times inferred by stable isotopes in SAS function models, Hydrol. Earth Syst. Sci. Discuss. [preprint],, in review, 2022. a

Wilusz, D. C., Harman, C. J., and Ball, W. P.: Sensitivity of Catchment Transit Times to Rainfall Variability Under Present and Future Climates, Water Resour. Res., 53, 10231–10256,, 2017. a, b, c, d

Winter, C., Lutz, S. R., Musolff, A., Kumar, R. Weber, M., and Fleckenstein, J. H.: Disentangling the Impact of Catchment Heterogeneity on Nitrate Export Dynamics From Event to Long-Term Time Scales, Water Resour. Res., 57, e2020WR027992,, 2020. a, b

Winter, C., Nguyen, T. V., Musolff, A., Lutz, S. R., Rode, M., Kumar, R., and Fleckenstein, J. H.: Droughts can reduce the nitrogen retention capacity of catchments, Hydrol. Earth Syst. Sci., 27, 303–318,, 2023.  a

Wollschläger, U., Attinger, S., Borchardt, D., Brauns, M., Cuntz, M., Dietrich, P., Fleckenstein, J. H., Friese, K., Friesen, J., Harpke, A., Hildebrandt, A., Jäckel, G., Kamjunke, N., Knöller, K., Kögler, S., Kolditz, O., Krieg, R., Kumar, R., Lausch, A., Liess, M., Marx, A., Merz, R., Mueller, C., Musolff, A., Norf, H., Oswald, S. E., Rebmann, C., Reinstorf, F., Rode, M., Rink, K., Rinke, K., Samaniego, L., Vieweg, M., Vogel, H. J., Weitere, M., Werban, U., Zink, M., and Zacharias, S.: The Bode hydrological observatory: a platform for integrated, interdisciplinary hydro-ecological research within the TERENO Harz/Central German Lowland Observatory, Environ. Earth Sci., 76, 29,, 2017. a

Xu, G., Magen, H., Tarchitzky, J., and Kafkafi, U.: Advances in Chloride Nutrition of Plants, Adv. Agron., 68, 97–150,, 1999. a

Yang, J., Heidbüchel, I., Musolff, A., Reinstorf, F., and Fleckenstein, J. H.: Exploring the Dynamics of Transit Times and Subsurface Mixing in a Small Agricultural Catchment, Water Resour. Res., 54, 2317–2335,, 2018. a

Yang, X., Seifeddine, J., Zink, M., Fleckenstein, J. H., Borchardt, D., and Rode, M.: A New Fully Distributed Model of Nitrate Transport and Removal at Catchment Scale, Water Resour. Res., 54, 5856–5877,, 2018. a, b

Zink, M., Kumar, R., Cuntz, M., and Samaniego, L.: A high-resolution dataset of water fluxes and states for Germany accounting for parametric uncertainty, Hydrol. Earth Syst. Sci., 21, 1769–1790,, 2017. a

Short summary
We analyzed the uncertainty of the water transit time distribution (TTD) arising from model input (interpolated tracer data) and structure (StorAge Selection, SAS, functions). We found that uncertainty was mainly associated with temporal interpolation, choice of SAS function, nonspatial interpolation, and low-flow conditions. It is important to characterize the specific uncertainty sources and their combined effects on TTD, as this has relevant implications for both water quantity and quality.