Research article 03 Jun 2020
Research article  03 Jun 2020
On the shape of forward transit time distributions in loworder catchments
 ^{1}Department of Hydrogeology, Helmholtz Centre for Environmental Research – UFZ, 04318 Leipzig, Germany
 ^{2}Department of Hydrology and Atmospheric Sciences, University of Arizona, Tucson, 85721, USA
 ^{1}Department of Hydrogeology, Helmholtz Centre for Environmental Research – UFZ, 04318 Leipzig, Germany
 ^{2}Department of Hydrology and Atmospheric Sciences, University of Arizona, Tucson, 85721, USA
Correspondence: Ingo Heidbüchel (ingo.heidbuechel@ufz.de)
Hide author detailsCorrespondence: Ingo Heidbüchel (ingo.heidbuechel@ufz.de)
Transit time distributions (TTDs) integrate information on timing, amount, storage, mixing and flow paths of water and thus characterize hydrologic and hydrochemical catchment response unlike any other descriptor. Here, we simulate the shape of TTDs in an idealized loworder catchment and investigate whether it changes systematically with certain catchment and climate properties. To this end, we used a physically based, spatially explicit 3D model, injected tracer with a precipitation event and recorded the resulting forward TTDs at the outlet of a small (∼6000 m^{2}) catchment for different scenarios. We found that the TTDs can be subdivided into four parts: (1) early part – controlled by soil hydraulic conductivity and antecedent soil moisture content, (2) middle part – a transition zone with no clear pattern or control, (3) later part – influenced by soil hydraulic conductivity and subsequent precipitation amount, and (4) very late tail of the breakthrough curve – governed by bedrock hydraulic conductivity. The modeled TTD shapes can be predicted using a dimensionless number: higher initial peaks are observed if the inflow of water to a catchment is not equal to its capacity to discharge water via subsurface flow paths, and lower initial peaks are connected to increasing available storage. In most cases the modeled TTDs were humped with nonzero initial values and varying weights of the tails. Therefore, none of the bestfit theoretical probability functions could describe the entire TTD shape exactly. Still, we found that generally gamma and lognormal distributions work better for scenarios of low and high soil hydraulic conductivity, respectively.
 Article
(9222 KB) 
Supplement
(3150 KB)  BibTeX
 EndNote
Transit time distributions (TTDs) characterize hydrologic catchment behavior unlike any other function or descriptor. They integrate information on timing, amount, storage, mixing and flow paths of water and can be modified to predict reactive solute transport (van der Velde et al., 2010; Harman et al., 2011; Musolff et al., 2017; Lutz et al., 2017). If observed in a time series, TTDs bridge the gap between hydrologic response (celerity) and hydrologic transport (velocity) in catchments by linking them via the change in water storage and the varying contributions of old (preevent) and young (event) water to streamflow (Heidbüchel et al., 2012). TTDs are time and spacevariant and hence no TTD of any individual precipitation event completely resembles another one. Therefore, in order to effectively utilize TTDs for the prediction of, for example, the effects of pollution events or water availability, it is necessary to find ways to understand and systematically describe the shape and scale of TTDs so that they are applicable in different locations and at different times. In this paper we look for firstorder principles that describe how the shape and scale of TTDs change, both spatially and temporally. This way we hope to improve our understanding of the dominant factors affecting hydrologic transport and response behavior at the catchment scale.
1.1 Initial use of theoretical probability distributions
Since the concept of TTDs was introduced, many studies have reported on their potential shapes and sought ways to describe them with different mathematical models like, for example, the pistonflow and exponential models (Begemann and Libby, 1957; Eriksson, 1958; Nauman, 1969), the advection–dispersion model (Nir, 1964; Małoszewski and Zuber, 1982), and the two parallel linear reservoirs model (Małoszewski et al., 1983; Stockinger et al., 2014). Dinçer et al. (1970) were the first to combine TTDs for individual precipitation events via the now commonly used convolution integral. Amin and Campana (1996) introduced the gamma distribution to transit time modeling.
Early studies reported that the outflow from entire catchments is characterized best with the exponential model (Rodhe et al., 1996; McGuire et al., 2005). However, neither the advection–dispersion nor the exponential model is able to capture the observed heavier tails of solute signals in streamflow (Kirchner et al., 2000). Instead, the more heavytailed TTDs created by advection and dispersion of spatially distributed rainfall inputs traveling toward the stream can be modeled with TTDs resembling gamma distributions (Kirchner et al., 2001). Likewise, tracer time series from many catchments exhibit fractal 1∕f scaling, which is consistent with gamma TTDs with shape parameter α≈0.5 (Kirchner, 2016).
1.2 General observations on the shape of TTDs
From the application of conceptual and physically based models we know that individual TTDs are highly irregular and that they can rapidly change in time for successive precipitation events (van der Velde et al., 2010; Rinaldo et al., 2011; Heidbüchel et al., 2012; Harman and Kim, 2014). If the early part of TTDs (mainly controlled by unsaturated transport in the soil layer) resembles a power law while the subsoil is responsible for the exponential tailing, the combination of those two parts can result in TTD shapes that are similar to gamma distributions (Fiori et al., 2009). In the field of groundwater hydrology there have been intense discussions on the tailing of breakthrough curves (e.g., on the issue of whether they follow a power law or not) (Haggerty et al., 2000; Becker and Shapiro, 2003; Zhang et al., 2007; Pedretti et al., 2013; Fiori and Becker, 2015; Pedretti and Bianchi, 2018). If disregarded, heavy tails can constitute a significant problem when using TTDs to predict solute transport because the legacy of contamination can be underestimated (not so much from a total mass balance perspective but when providing risk assessments for highly toxic pollutants reaching further into the future). Furthermore, a truncation of heavy TTD tails should be avoided, especially when computing mean transit times (mTTs) since they are highly sensitive to the shape of the chosen transfer function (Seeger and Weiler, 2014). Other complicating matters are special cases of bimodal TTDs that can be caused by varying contributions from fast and slow storages (McMillan et al., 2012) or from urban and rural areas (Soulsby et al., 2015). Apart from individual catchment and event properties, mixing assumptions also affect TTD modeling since certain TTD shapes are inherently linked to specific mixing assumptions (e.g., a wellmixed system is best represented by an exponential distribution, partial mixing can be approximated with gamma distributions and no mixing with the pistonflow model) (van der Velde et al., 2015).
1.3 Controls on shape variations
A number of studies reported on the bestfit shape of gamma distributions generally ranging from α 0.01 to 0.90 (Hrachowitz et al., 2009; Godsey et al., 2010; Berghuijs and Kirchner, 2017; Birkel et al., 2016), which indicates Lshaped distributions with high initial values and heavier tails. Several studies found that α values decrease with increasing wetness conditions (e.g., Birkel et al., 2012; Tetzlaff et al., 2014), causing higher initial values and heavier tails. However, the opposite was observed in a boreal headwater catchment (PeraltaTapia et al., 2016) where α ranged between 0.43 and 0.76 for all years except the wettest year (α=0.98). In the Scottish highlands α showed little temporal variability (and therefore no link to precipitation intensity) but was closely related to catchment landscape organization – especially soil parameters and drainage density – where a high percentage of responsive soils and a high drainage density resulted in small values of α (Hrachowitz et al., 2010).
Conceptual and physically based models have also been used to investigate the (temporally variable) shapes of TTDs. Haitjema (1995) found that the TTD of groundwater can resemble an exponential distribution while Kollet and Maxwell (2008) and Cardenas and Jiang (2010) derived a powerlaw form and fractal behavior adding macrodispersion and systematic heterogeneity to the domain in the form of depthdecreasing poromechanical properties. Increasing the vertical gradient of conductivity decay in the soil decreased the shape parameter α (from 0.95 for homogeneous conditions down to a value of 0.5 for extreme gradients) in a study by Ameli et al. (2016). Somewhat surprisingly, the level of “unstructured” heterogeneity within the soil and the bedrock was found to only have a weak influence on the shape of TTDs (Fiori and Russo, 2008) since the dispersion is predominantly ruled by the distribution of flow path lengths within a catchment. Antecedent moisture conditions and event characteristics influence catchment TTDs at short timescales while land use affects both short and long timescales (Weiler et al., 2003; RoaGarcía and Weiler, 2010); generally TTD shapes appear highly sensitive to catchment wetness history and available storage, mixing mechanisms and flow path connectivity (Hrachowitz et al., 2013). Kim et al. (2016) recorded actual TTDs in a sloping lysimeter and reported that their shapes varied both with storage state and the history of inflows and outflows. They argued that “the observed time variability […] can be decomposed into two parts: [1] “internal” […] – associated with changes in the arrangement of, and partitioning between, flow pathways; and [2] “external” […] – driven by fluctuations in the flow rate along all flow pathways”. From these partly contradictory findings, it is clear that relating bestfit values for the shape parameter α of the gamma distribution to catchment, climate or precipitation event properties does not yield a consistent picture yet. Moreover, the shape of TTDs is also dependent on the resolution of time series data (sampling frequency). While α can decrease with longer sampling intervals – since the nonlinearity of the flow system is overestimated when sampling becomes more infrequent (Hrachowitz et al., 2011) – higher α values can also result from lowering the sampling frequency in both input (precipitation) and output (streamflow) (Timbe et al., 2015).
Replacing transit time with flowweighted time or cumulative outflow (Niemi, 1977; Nyström, 1985) erases a substantial amount of the TTD shape variation associated with the external variability. However, since a change in the inflow often causes both fluctuations along and also a rearrangement between the flow pathways (i.e., internal variability), flowweighted time approaches are not able to completely remove the influence of changes in the inflow rate. Still, Ali et al. (2014), providing a comprehensive assessment of different transittimebased catchment transport models (where they compare several timeinvariant to timevariable methods), conclude that applying a flowweighted time approach can indeed yield adequate results for predicting catchmentscale transport.
1.4 TTD theory
To summarize, soil hydraulic conductivity, antecedent moisture conditions (storage state), soil thickness, and precipitation amount and intensity are amongst the most frequently cited factors that influence the shape of TTDs. Obviously, there is not one single property that controls the TTD shape. Instead, the interplay of several catchment, climate and event characteristics results in the unique shape of every single TTD. One approach to deal with this problem of multicausality is the use of dimensionless numbers. Heidbüchel et al. (2013) introduced the flow path number F which combines several catchment, climate and event properties into one index relating flows in and out of the catchment to the available subsurface storage. It was originally designed to monitor the exceedance of certain storage thresholds for the activation of different dominant flow paths (groundwater flow, interflow, overland flow) at the catchment scale but can also help to categorize and predict TTD shapes. Moreover, from continuous time series of TTDs one can mathematically derive residence time distributions (describing the age distribution of water stored in the catchment), storage selection functions (describing the selection preference of the catchment discharge for younger or older stored water) (Botter et al., 2010, 2011; van der Velde et al., 2012; Benettin et al., 2015; Harman, 2015; Pangle et al., 2017; DaneshYazdi et al., 2018; Yang et al., 2018) and master transit time distributions (MTTDs) (representing the flowweighted average of all TTDs of a catchment) (Heidbüchel et al., 2012; Sprenger et al., 2016; Benettin et al., 2017), which can all take on different shapes depending on climate and catchment properties, just like the individual TTDs. Hence the results presented in this paper can also provide insights into the use of these descriptors of catchment hydrologic processes.
Since McGuire and McDonnell (2006) stated a lack of theoretical work on the actual shapes of TTDs, quite a diverse range of research has been conducted to approach this problem from different angles and has yielded fragments of important knowledge. However, what is still missing is a coherent framework that enables us to structure our understanding of the nature of TTDs so that it eventually becomes applicable to real world hydrologic problems. Already in 2010, McDonnell et al. (2010) asked how the shape of TTDs could be generalized and how it would vary with ambient conditions, from time to time and from place to place. This study sets out to provide such a coherent framework which – although not exhaustive (or entirely correct for that matter) – will provide us with testable hypotheses on how the shape and scale of TTDs change spatially and temporally. As Hrachowitz et al. (2016) put it, “an explicit formulation of transport processes, based on the concept of transit times has the potential to improve the understanding of the integrated system dynamics […] and to provide a stronger link between […] hydrological and water quality models”.
1.5 Our approach
In this study we will make use of a physically based, spatially explicit, 3D model to systematically simulate how different catchment properties and climate characteristics and also their interplay control the shape of forward TTDs. We test which TTD shapes are most appropriate for capturing hydrologic and hydrochemical catchment response at different locations and for specific points in time. Furthermore we will try to interpret the results in the most general way possible, so that the theory can be extended to other potential controls of the TTD shape in the future. Our modeling does not explicitly include preferential flow within the soil and bedrock (like, for example, macropores or fractures), and therefore our TTDs mostly represent systems where water is transported via subsurface matrix flow coupled with overland flow. Still, the exclusion of these components can be considered legitimate and the results meaningful because of the important role that macrodispersion plays in shaping TTDs (Fiori et al., 2009). Hence, we consider our results the basis for further investigations approaching ever more realistic representations of the many hydrological processes taking place at the catchment scale.
We used HydroGeoSphere (HGS), a 3D numerical model describing fully coupled surface–subsurface, variably saturated flow and advective–dispersive solute transport (Therrien et al., 2010). Groundwater flow in the 3D subsurface is simulated with Richards' equation and Darcy's law, and surface runoff is simulated in the 2D surface domain with Manning's equation and the diffusivewave approximation of the SaintVenant equations. The classical advection–dispersion equation for solute transport is solved in all domains. The surface and subsurface domains are numerically coupled using a dual node approach, allowing for the interaction of water and solutes between the surface and subsurface. The general functionality of HGS and its adequacy for solving analytical benchmark tests has been proven in several model intercomparison studies (Maxwell et al., 2014; Kollet et al., 2017), and its solute transport routines have been verified against laboratory (Chapman et al., 2012) and field measurements (Sudicky et al., 2010; Liggett et al., 2015; Gilfedder et al., 2019). Since our modeling approach entails subsurface flow only in porous media (no explicit fractures or macropores are included), the resulting TTDs have to be considered a special subset of distributions lacking some of the dynamics we can expect in realworld catchments while still providing a sound basis for further investigations (like, for example, adding more complex interaction dynamics along the flow pathways).
2.1 Model setup
A small zeroorder catchment was set up, 100 m long, 75 m wide (∼6000 m^{2}), with an average slope of 20 % towards the outlet and elliptical in shape (Fig. 1). The catchment converges slightly towards the center, creating a gradient that concentrates flow. The bedrock is 10 m thick and has a saturated hydraulic conductivity of ${K}_{\mathrm{Br},x}={K}_{\mathrm{Br},y}={\mathrm{10}}^{\mathrm{5}}$ m d^{−1} (horizontal) and ${K}_{\mathrm{Br},z}={\mathrm{10}}^{\mathrm{6}}$ m d^{−1} (vertical). The soil layer is isotropic, of uniform thickness and has a higher hydraulic conductivity. All other parameters are uniform across the entire model domain (based on values typically found in many catchments in central Europe): porosity n=0.39 m^{3} m^{−3}, van Genuchten parameters alpha α_{vG}=0.5 m^{−1}, beta β_{vG}=1.6, saturated water content θ_{s}= 0.39 m^{3} m^{−3}, residual water content θ_{r}= 0.05 m^{3} m^{−3}, poreconnectivity parameter l_{p}= 0.5, and longitudinal and transverse dispersivity α_{L}=5 m and α_{T}=0.5 m, respectively. The magnitude for α_{L} was estimated with regard to the length of the model catchment (100 m) using the relationship described in Gelhar et al. (1992) and SchulzeMakuch (2005). Newer research by Zech et al. (2015) shows that the longitudinal dispersivity is probably up to an order of magnitude smaller than that reported by Gelhar et al. (1992). Still, we do not suspect the value of the local dispersivity to have a large impact on the TTDs (see also Fiori and Russo, 2008) since the local dispersivity is usually negligible compared to the dispersion caused by the spatial distribution of rainfall (the “source zone dispersion” in Fiori et al., 2009). Nevertheless, we tested whether changing α_{L} from 5 to 0.5 m would change our results significantly (see Sect. S1, Fig. S1 and Table S1 in the Supplement). Both bedrock and soil are exclusively porous media without any potential preferential flow paths like macropores or rock fractures.
2.1.1 Boundary conditions
Both the bottom and the sides of the domain were impermeable boundaries. A constant head boundary condition (equal to the surface elevation) was assigned to the lower front edge of the subsurface domain (nodes in the blue square in Fig. 1), allowing outflow from both the bedrock and the soil. A critical depth boundary was assigned to the lower edge of the surface domain (on top of the constant head boundary) to allow for overland flow out of the catchment. The surface of the catchment received spatially uniform precipitation. We used a recorded time series of precipitation from the northeast of Germany (maritime temperate climate: Cfb in the Köppen climate classification) amounting to 690 mm a^{−1} (Fig. 2a). The time series was 1 year long and repeated 32 more times to cover the entire modeling period, which lasted a total of ∼33 years (12 000 d). We made sure that the looping of the precipitation time series would not cause any unwanted artifacts in the resulting TTDs (see Sect. S2 and Fig. S2). Neither evaporation nor transpiration was considered during the simulations. This means that all precipitation we applied was effective precipitation that would eventually discharge at the catchment outlet. The addition of the process of evapotranspiration is planned in a followup modeling study to investigate what influence it exerts on catchment TTDs. The tracer was applied uniformly over the entire catchment during a precipitation event that lasted 1 h and had an intensity of 0.1 mm h^{−1} and a tracer concentration of 1 kg m^{−3}. This resulted in a total applied tracer mass of 0.589 kg.
2.1.2 Initial conditions
The model runs were initialized with three different antecedent soil moisture conditions θ_{ant}: a dry one (θ_{ant}=22.0 %, corresponding to an average effective saturation of the soil layer S_{eff}≈50 %), an intermediate one (θ_{ant}=28.8 %; S_{eff}≈70 %) and a wet one (θ_{ant}=35.6 %; S_{eff}≈90 %). To obtain realistic distributions of soil moisture, we first ran the model starting with full saturation and without any precipitation input and let the soils drain until the average effective saturation reached the states for our initial conditions. We recorded these conditions and used them as initial conditions of the virtual experiment runs. In general, the soil remained wetter close to the outlet in the lower part and became drier in the upper part of the catchment. Note that the process of evapotranspiration was excluded from the modeling so that the lowest achievable saturation was essentially defined by the field capacity. An average effective saturation S_{eff} of approximately 50 % was the lowest that could be achieved by draining the soil layer since the lower part of the catchment stayed highly saturated due to the constant head boundary condition being equal to the surface elevation at the outlet. The upper parts of the catchment, however, were initiated with much lower S_{eff} values (≈30 % in the dry scenarios). That means that although an S_{eff} value of 50 % seems to be quite high, it actually represents an overall dry state of the catchment soil. Throughout the modeling runs the dry initial condition did not occur again as that would have taken 13 years of drainage without any precipitation for the scenarios with high soil hydraulic conductivity K_{S} and almost 1500 years for the scenarios with low K_{S}. The inclusion of evapotranspiration would, however, speed up the drying process of the soil and hence make these initial conditions more realistic.
2.2 Model scenarios
To investigate how different catchment and climate properties influence the shape of forward TTDs we systematically varied four characteristic properties from high to low values and looked at the resulting TTD shapes of all the possible combinations (for a total number of 36 scenarios). The properties we focused on were soil depth (D_{soil}), saturated soil hydraulic conductivity (K_{S}), antecedent soil moisture content (θ_{ant}) and subsequent precipitation amount (P_{sub}, essentially a measure of the amount of precipitation that falls after the delivery of the traced event) (Fig. 3).
We tested two soil depths D_{soil}, namely depths of 0.5 and 1.0 m, evenly distributed across the entire catchment. Similarly, we chose two saturated soil hydraulic conductivities K_{S}, a high one with 2.0 m d^{−1} (similar to fine sand) and a low one with 0.02 m d^{−1} (similar to silt). Three states of antecedent moisture content θ_{ant} were selected to represent initial conditions – 50, 70 and 90 % of effective saturation. Finally the subsequent precipitation amount P_{sub} was varied in three steps from 345 to 690 and up to 1380 mm a^{−1}. The original precipitation time series (690 mm a^{−1}, Fig. 2a) was rescaled to obtain time series with smaller and larger amounts. With two soil depths, two soil hydraulic conductivities, three antecedent moisture conditions and three subsequent precipitation amounts this resulted in 36 model scenarios. Based on these 36 runs we evaluated the differences in the shape of the TTDs. The abbreviated names of the 36 model runs consist of four letters, each representing one of the properties that we varied: the first one describes D_{soil} (T = thick; F = flat), the second one K_{S} (H = high; L = low), the third one θ_{ant} (W = wet; I = intermediate; D = dry) and the fourth one P_{sub} (S = small; M = medium; B = big). For example the name FHIB would indicate a run with a flat (“F”) (shallow) soil, a high (“H”) K_{S}, an intermediate (“I”) antecedent moisture content and a big (“B”) subsequent precipitation amount (see Table 1 for an overview of the names of all 36 scenarios). We are well aware that “thick” and “flat” are technically incorrect descriptions of soil depth. However, in order to have unique identifiers (i.e., individual letters) for all 10 property states we decided to use T and F for describing deep and shallow soils, respectively.
To complement the results obtained from the systematic variation of catchment and climate characteristics we tested the influence of seven other factors: (1) soil porosity, (2) bedrock hydraulic conductivity, (3) exponential decay in hydraulic conductivity with depth in the soil, (4) frequency of precipitation events, (5) soil water retention curve, (6) catchment shape and (7) effect of extreme precipitation after full saturation – conditions during which direct surface runoff may occur. These additional runs with altered soil properties and boundary and initial conditions were performed on the basis of some of the 36 initial runs (in the following sections we always indicate which runs form the basis of the specific scenarios; see also Table S2).
Notable catchment properties we did not test include topography, size, slope and curvature. Apart from investigating the effect of an exponential decay in soil hydraulic conductivity with depth we did not add heterogeneity to the subsurface hydraulic properties. Therefore we cannot make statements about how multiple soil layers or different spatial patterns of hydraulic conductivity would influence TTDs.
2.2.1 Soil porosity
The influence of larger and smaller soil porosity was investigated with six additional runs based on the three scenarios THDM, THIM and THWM. Three of the additional runs had larger (0.54 m^{3} m^{−3}) and three had smaller soil porosity (0.24 m^{3} m^{−3}) than the basecase scenarios (0.39 m^{3} m^{−3}).
2.2.2 Bedrock hydraulic conductivity
Six runs were performed on the basis of the THDB scenario (which had a bedrock hydraulic conductivity K_{Br} of 10^{−5} m d^{−1}). In the first run K_{Br} was decreased to 10^{−7} m d^{−1}, in the following runs it was successively increased to 10^{−3}, 10^{−2}, 10^{−1}, 10^{0} and 2×10^{0} m d^{−1}, matching the K_{S} of the soil layer in the final run.
2.2.3 Decay in saturated hydraulic conductivity with depth
Because all other model scenarios had a constant hydraulic conductivity throughout the soil layer, we wanted to test whether the introduction of an exponential decay in hydraulic conductivity with depth (from high conductivity at the surface to low conductivity at the soil–bedrock interface; see Bishop et al., 2004; Jiang et al., 2009) would have a large influence on the TTD shapes. We based the conductivity decay test on four scenarios (THDB, THWB, TLDB and TLWB), adding relationships of soil depth z and saturated hydraulic conductivity K_{S} with a shape parameter f= 0.29 m and saturated hydraulic conductivity at the surface K_{S0}= 7 m d^{−1} (for the high conductivity scenarios) or K_{S0}=0.07 m d^{−1} (for the low conductivity scenarios) (Fig. 4a):
This preserved the mean K_{S} values of $\mathrm{2}\times {\mathrm{10}}^{\mathrm{0}}$ (high) and $\mathrm{2}\times {\mathrm{10}}^{\mathrm{2}}$ m d^{−1} (low) (from the basecase scenarios).
2.2.4 Precipitation frequency
Five time series with high precipitation event frequency and five time series with low precipitation event frequency were created by means of the rainfall generator used by Musolff et al. (2017) (Fig. 2b). It generates Poisson effective rainfall (Cox and Isham, 1988) which is characterized by exponentially distributed rainfall event amounts and interarrival times. The mean interarrival time was set to 3 d and 15 d for the highfrequency scenarios (comparable to a humid precipitation distribution and intensity pattern with lower intensities and more frequent events) and lowfrequency scenarios (comparable to an arid precipitation distribution and intensity pattern with higher intensities and less frequent events), respectively. The total precipitation for all scenarios (both humid and arid type) was 690 mm so that it matched our medium P_{sub} scenarios.
2.2.5 Water retention curve
All the basecase model scenarios were conducted with water retention curves (WRCs) resembling silty soils:
with van Genuchten parameters α_{vG} (m^{−1}) and β_{vG} (dimensionless), saturated water content θ_{s}, residual water content θ_{r} (both m^{3} m^{−3}), pressure head ψ (m) and $\mathit{\nu}=\mathrm{1}\mathrm{1}/{\mathit{\beta}}_{\mathrm{vG}}$ (see Sect. 2.1 for van Genuchten parameter values). However, we also wanted to investigate how a different WRC in the soil layer (see Fig. 4b) would influence the shape of TTDs. We chose to test a sandtype WRC since it can, in some aspects and to a certain extent, also indicate how a system with the thresholdlike initiation of rapid preferential flow behaves. The sandtype WRC causes an increase in hydraulic conductivity at relatively lower soil water contents compared to the silttype WRC. Hence, for the same precipitation event lateral flow is initiated faster (at lower saturations) in sandy soils since water reaches the soil–bedrock interface more quickly, where it is diverted from vertical to lateral flow. The relative hydraulic conductivity k_{r} was derived with Eq. (3):
with effective saturation S_{eff} and poreconnectivity parameter l_{p} (both dimensionless). Other aspects of preferential flow – like bypass flow through macropores in deeper soil layers – are, however, not captured by sandtype WRCs. The van Genuchten parameters for the sandtype WRC were defined as follows: α_{vG}=14.5 m^{−1} and β_{vG}=2.68. We based the additional eight runs on the scenarios THDB, THWB, THDS, THWS, TLDB, TLWB, TLDS and TLWS.
2.2.6 Catchment shape
In addition to the oval catchment we designed two more shapes to get an idea whether it would have a significant impact on the resulting TTDs (Fig. 1d). One of the catchments had the center of gravity located farther away from the outlet (top; 60 m) the other catchment had the center of gravity located closer to the outlet (bottom; 40 m). This increased the average flow path length from 61 to 70 m for top and decreased it to 55 m for bottom – while catchment length, area and slope stayed the same for all cases. The four additional runs we conducted were based on the scenarios THWM and THDM.
2.2.7 Full saturation and extreme precipitation intensity
We tested these effects for two scenarios (THWB and TLWB) since both of these scenarios were already close to creating overland flow. Full saturation in this case means that the initial condition for these model runs consisted of a fully saturated domain (both in the bedrock and in the soil); i.e., S_{eff} was 100 % (θ_{ant}=39 %). Additionally, we increased the intensity of the input precipitation event (delivering the tracer) from 0.1 mm h^{−1} (normal) to 10 mm h^{−1} (very large, +) and up to 100 mm h^{−1} (extreme, $+++$), in an attempt to create infiltration excess overland flow and record its influence on the shape of TTDs.
2.3 Influence of the sequence of precipitation events
We also tested to what extent the sequences of subsequent precipitation events with different magnitude, intensity and interarrival time influence TTD shapes. This was necessary to assure that our resulting TTD shapes were not primarily a product of the point in time – within the sequence of precipitation events – at which the tracer was applied to the catchment. To this end 15 precipitation event time series were created by means of the rainfall generator used by Musolff et al. (2017). The mean interarrival time was set to 3 d (comparable to a precipitation distribution and intensity pattern found in humid environments with low intensities and more frequent events) and the total precipitation amount for all scenarios was 690 mm, matching our medium P_{sub} scenarios (Fig. S3). The generated precipitation time series resembled our original time series of precipitation, which also had an interarrival time close to 3 d. All other parameters and properties of the 15 model runs were based on the THDM scenario.
2.4 Processing of the output data
The output data from HGS were mainly processed with Microsoft Excel. We summed surface and subsurface flows, computed total tracer outflow from the catchment, created the probability density and cumulative probability density distribution for tracer outflow, calculated the metrics of the forward TTDs, fitted theoretical distributions to our data and smoothed the original TTDs for better visual comparability of the shapes. HGS keeps track of the mass balance of inflow, outflow and storage and calculates the discrepancy (mass balance error) between the three terms (Fig. S4). The mean absolute mass balance error for the 36 runs was negligible ($\mathrm{6.8}\times {\mathrm{10}}^{\mathrm{2}}\pm \mathrm{7.2}\times {\mathrm{10}}^{\mathrm{2}}$ %).
2.4.1 Creation of TTDs
The probability density distributions of transit time (the forward TTDs, d^{−1}) were created by normalizing the mass outflux J_{out} (kg d^{−1}) for each time step by the total inflow mass M_{in} (kg):
The cumulative TTDs (dimensionless) were created by multiplying the normalized mass outflux (d^{−1}) of each time step by the associated time step length Δt (d) before cumulating it:
2.4.2 Calculation of TTD metrics
For each TTD we calculated seven metrics to characterize its shape: the first quartile (Q_{1}), the median (Q_{2}), the mean (mTT), the third quartile (Q_{3}), the standard deviation (σ), the skewness (ν) and the excess kurtosis (γ) (see Sect. S3 and Fig. S5 for details on the calculation and for visual comparison of the metrics). Furthermore we determined the young water fraction F_{yw} as the fraction of water leaving the catchment after 2.3 months (Jasechko et al., 2016; Kirchner, 2016; Wilusz et al., 2017). For more details on how F_{yw} changes with catchment and climate properties, see Sect. S4, Fig. S6 and Table S3.
2.4.3 Fitting
We fitted predefined mathematical probability density functions to the modeled data since condensing the main characteristics of an observed probability distribution into just one to three parameters of a mathematical function is appealing and eases the potential of transferability of the findings. Massoudieh et al. (2014) explored the use of freeform histograms as groundwater age distributions and concluded that mathematical distributions performed better in terms of their ability to capture the observed tracer data relative to their complexity. In order to determine which theoretical probability density function best captures the shape of our modeled TTDs, we chose two probability density functions that are commonly used to describe the transit of water through catchments (inverse Gaussian and gamma), as well as the less common lognormal distribution which also has just two adjustable parameters:

Inverse Gaussian distribution (a particular solution of the advection–dispersion equation) with dispersion parameter D (dimensionless) and mean mTT (d):
$$\begin{array}{}\text{(6)}& \mathrm{InvGau}\left(t\right)={\left({\displaystyle \frac{\mathrm{4}\mathit{\pi}Dt}{\mathrm{mTT}}}\right)}^{\mathrm{0.5}}{\displaystyle \frac{\mathrm{1}}{t}}\mathrm{exp}\mathit{\{}[{\left(\mathrm{1}{\displaystyle \frac{t}{\mathrm{mTT}}}\right)}^{\mathrm{2}}\cdot {\displaystyle \frac{\mathrm{mTT}}{\mathrm{4}Dt}}\left]\mathit{\right\}}.\end{array}$$ 
Gamma distribution with shape parameter α (dimensionless) and scale parameter β (d) (with mean mTT =αβ):
$$\begin{array}{}\text{(7)}& \mathrm{Gamma}\left(t\right)={t}^{\mathit{\alpha}\mathrm{1}}{\displaystyle \frac{{e}^{t/\mathit{\beta}}}{{\mathit{\beta}}^{\mathit{\alpha}}\mathrm{\Gamma}\left(\mathit{\alpha}\right)}}.\end{array}$$Gamma distributions can take on very different shapes when α is changed: α < 1, highly skewed distributions with an initial maximum and heavier (i.e., subexponential) tails; α=1, exponential distribution; α > 1, less skewed, “humped” distributions with an initial value of 0, a mode and lighter tails. They can be stretched or compressed with scale parameter β. Thus when using gamma distributions for the determination of mTTs, it is necessary to choose the correct shape parameter α to avoid problems of equifinality. The same holds true for all multiple parameter distributions.

Lognormal distribution with standard deviation σ and mean μ (both dimensionless) of the natural logarithm of the variable (with mean mTT = exp($\mathit{\mu}+{\mathit{\sigma}}^{\mathrm{2}}/\mathrm{2}$)):
$$\begin{array}{}\text{(8)}& \mathrm{LogN}\left(t\right)={\displaystyle \frac{\mathrm{1}}{t\mathit{\sigma}\sqrt{\mathrm{2}\mathit{\pi}}}}\mathrm{exp}\left[{\displaystyle \frac{{\left(\mathrm{ln}t\mathit{\mu}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma}}^{\mathrm{2}}}}\right].\end{array}$$
We tested two more probability density functions both having three (instead of just two) adjustable parameters:

Threeparameter beta distribution with shape parameters α and β (dimensionless) and upper limit c (d) (with mean mTT $=\mathit{\alpha}c/(\mathit{\alpha}+\mathit{\beta}))$:
$$\begin{array}{}\text{(9)}& \mathrm{beta}\left(t\right)={\displaystyle \frac{{t}^{\mathit{\alpha}\mathrm{1}}(ct{)}^{\mathit{\beta}\mathrm{1}}}{{c}^{\mathit{\alpha}+\mathit{\beta}\mathrm{1}}B(\mathit{\alpha},\mathit{\beta})}}.\end{array}$$The fourth parameter of the beta distribution could be the lower limit a. It is not included in the above definition since in our case it is zero.

Truncated lognormal distribution with the time of truncation λ (d) as the third parameter:
$$\begin{array}{}\text{(10)}& \begin{array}{rl}& \mathrm{Trunc}\left(t\right)=\left\{{\displaystyle \frac{\mathrm{1}}{\left(t+\mathit{\lambda}\right)\mathit{\sigma}\sqrt{\mathrm{2}\mathit{\pi}}}}\mathrm{exp}\left[{\displaystyle \frac{{\left(\mathrm{ln}(t+\mathit{\lambda})\mathit{\mu}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma}}^{\mathrm{2}}}}\right]\right\}/\\ & \left\{\mathrm{1}{\int}_{t=\mathrm{0}}^{\mathit{\lambda}}{\displaystyle \frac{\mathrm{1}}{t\mathit{\sigma}\sqrt{\mathrm{2}\mathit{\pi}}}}\mathrm{exp}\left[{\displaystyle \frac{{\left(\mathrm{ln}t\mathit{\mu}\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma}}^{\mathrm{2}}}}\right]\mathrm{d}t\right\}.\end{array}\end{array}$$
For visual examples of all five types of distributions please refer to Fig. S7.
The method of least squares was used to find the best fit between the modeled TTDs and the theoretical distribution functions (i.e., minimizing the sum of the squared residuals with the Solver function in Excel using one value for each of the 12 000 d of the modeled TTDs).
The fitting was performed on the cumulative probability distributions since their shape is not subject to the more extreme internal variability that probability distributions can experience.
2.4.4 Smoothing
Smoothing was only applied to enhance the visual comparability of the TTDs. All calculations were performed on the unsmoothed TTDs. For details on the smoothing method see Sect. S5 and Fig. S8.
2.5 Flow path number
The flow path number F is a dimensionless number proposed by Heidbüchel et al. (2013) that relates catchment inflow to outflow (in the numerator) while simultaneously assessing available storage space (in the denominator) for each point in time and at the catchment scale. It was introduced to define thresholds for the activation and deactivation of different flow paths that transport water more slowly (e.g., groundwater flow), faster (interflow) or very fast (macropore flow, overland flow). For this paper we modified F slightly so that both numerator and denominator have the dimensions (m^{3}):
where soil depth D_{soil} (m), catchment surface area A_{in} (m^{2}), porosity n (m^{3} m^{−3}) and antecedent moisture content θ_{ant} (m^{3} m^{−3}) are paired with the driving precipitation amount P_{dr} (m^{3}), which is calculated as the average subsequent precipitation amount P_{sub} (m a^{−1}) over the average event duration t_{Ev} (d):
The subsequent precipitation amount P_{sub} (m a^{−1}) is calculated for every time step as the amount of precipitation falling within the year that follows this time step using a moving window. Note that differing from Heidbüchel et al. (2013) we used the event duration t_{Ev} instead of the interevent duration t_{Ie} to compute P_{dr} since it better represents the amount of precipitation falling during an average event filling up the available storage. Furthermore, the subsurface discharge capacity of the soil K_{rem} (m^{3}) consists of the effective saturated soil hydraulic conductivity K_{S} (m d^{−1}), the sum of the average interevent and event duration t_{Ie}+t_{Ev} (d), the porosity n (m^{3} m^{−3}), and the crosssectional area of the soil layer at the outlet of the catchment A_{out} (m^{2}):
The crosssectional area of the soil layer at the outlet of the catchment A_{out} can be considered to represent the connection of the catchment to either a river channel, riparian zone or the alluvial valley fill where medium to rapid subsurface outflow from the catchment can occur. Note that differing from Heidbüchel et al. (2013) we used the sum of the interevent and event duration t_{Ie}+t_{Ev} instead of just the event duration t_{Ev} to compute K_{rem} since it better represents the amount of water that can be removed from the catchment during an average precipitation cycle.
The flow path number F varies in time mainly due to the changes in antecedent moisture content θ_{ant}, since variations in the amount of driving precipitation P_{dr} are damped due to the moving window approach that is used to compute it. That means F can vary quite rapidly (towards either more positive or negative values) during the wetup phase of a catchment and change more slowly (towards 0) during the drydown phase. A positive flow path number F indicates that there is a surplus of water entering the catchment that cannot be removed by subsurface transport at the same rate. Hence, the storage fills up. Conversely, a negative F indicates that the drainage capacity of the catchment exceeds the water inputs and the amount of stored water decreases. Furthermore, values between 0 and 1 signal that the available soil storage space is able to accommodate the net inflow of water, while values larger than 1 mean that the catchment receives more water than it can discharge or store in the subsurface. In turn, the larger the storage capacity in the subsoil, the more F converges towards 0. There is only one notable important exception to this last rule: in highly conductive soils the increase in discharge capacity (caused by an increase in soil depth and the consequential increase in the crosssectional area of the soil layer at the outlet A_{out}) can be larger than the increase in storage capacity itself – leading to F becoming even more negative with increasing storage capacity.
Output from the model runs comprised subsurface discharge, overland discharge and tracer mass outflux in the discharge from which we derived TTDs (for an example see Fig. S9). Additionally, the model provided spatially and temporally resolved tracer concentrations throughout the entire domain. The differences emerging between the individual TTDs can be tracked by looking at the spatiotemporal evolution of the applied tracer impulse throughout the entire catchment. For a detailed example please refer to Sect. S6 and Fig. S10.
3.1 Influence of the sequence of precipitation events
Changing the sequence of precipitation events affects the shape of TTDs to a certain degree. In particular the timing and magnitude of the first precipitation event determines how strong the early response turns out. This can be observed in Fig. 5, where the different TTDs split up into different branches according to the arrival and magnitude of the first event after tracer application. However, following this initial split – with more and more precipitation events taking place – all TTDs tend to converge towards a single line. Examining the cumulative TTDs in Fig. 5 it is obvious that the variability in the TTD shape introduced by different precipitation event sequences is much smaller than the variability introduced by the other catchment and climate properties. While the range of Q_{1} observed for the 15 scenarios with different event sequences is still 14 % of the total range observed for the 36 basecase scenarios, this percentage decreases down to 2 % for Q_{3}. The other distribution metrics describing the shape of the TTDs also vary a lot less between the scenarios with different event sequences compared to the scenarios with different catchment and climate properties (the range of all event sequences is only 1.1 % of the range of all basecase scenarios for the standard deviation, 1.6 % for the skewness and 1.0 % for the excess kurtosis). A table with the distribution metrics for all 15 scenarios can be found in the Supplement (Table S4). Therefore we can assume that the shape of TTDs is not significantly influenced by the precipitation event sequence – at least in environments with a naturally short interarrival time resembling humid climate conditions and an event amount distribution that is exponential.
3.2 Effects on TTD metrics
We found that θ_{ant} affects the young parts of TTDs (the first 10 d) a lot more than the older parts (its influence is hardly discernible after approximately 100 d; see Fig. 6a). By contrast, K_{S} affects the older parts more than the young parts. This difference is due to the fact that θ_{ant} constitutes one of the initial conditions that also directly influences the current soil hydraulic conductivity while the influence of different K_{S} values gains more importance later when the soil moisture conditions become more similar. D_{soil} and P_{sub} influence all parts of the TTDs equally strongly and hence have the smallest influence on the actual shape of the distributions. As can be observed in (b), the influence of K_{S} is a lot stronger in scenarios with wet θ_{ant} while the influence of P_{sub} decreases with increasing θ_{ant}. Panel (c) shows that both θ_{ant} and P_{sub} have a larger influence when K_{S} is high, but for P_{sub} this increased influence is only seen for the longer transit times. The influence of the initial condition θ_{ant} is larger when K_{S} is high because the relative differences in flow through a dry soil and a wet soil are larger for soils with high K_{S} compared to soils with low K_{S}. Panel (d) confirms the impression that D_{soil} only has a minor influence on the shape of TTDs – all parts of the TTDs are equally affected and it does not make a significant difference for the influence of the other factors whether the soils are deeper or shallower. Finally in the lower right panel it is demonstrated that P_{sub} has opposite effects on the influence of θ_{ant} and K_{S}: Larger P_{sub} causes the influence of K_{S} to increase for the longer transit times while the influence of θ_{ant} decreases when P_{sub} becomes larger. The fact that different catchment and climate properties have varying degrees of control on transit times depending on current conditions and the interplay of dominant hydrologic processes has already been observed in the field (Heidbüchel et al., 2013). Table 1 lists all metrics of the 36 TTDs resulting from the basecase scenarios.
3.2.1 Antecedent moisture content
Dry θ_{ant} results in a lower probability for shorter transit times while wet θ_{ant} triggers faster responses and higher initial peaks for TTDs (Fig. 7). When increasing θ_{ant} by 14 % (from S_{eff} 50 % to 90 %), on average Q_{1} decreases by 44 %, Q_{2} by 27 %, the mTT by 19 % and Q_{3} by 15 % (Fig. 6 center, Table 1). The median F_{yw} increases by 16 %. Neither the standard deviation (and hence the width) nor the skewness nor the kurtosis values of the TTDs are affected much by θ_{ant} though. Higher θ_{ant} initially promotes faster lateral transport (both on the surface and in the subsurface) while impeding percolation of tracer towards the bedrock, and therefore more tracer is transported fast towards the outlet and less tracer is entering the deeper soil layers and the bedrock. Longterm trends or interannual shifts in P_{sub} can cause temporal changes in TTDs but substantial shortterm variations are derived mainly from differences in θ_{ant}. Therefore variations in TTD shape and scale can be high even in relatively small catchments. Generally, the influence of θ_{ant} is stronger for catchments with higher K_{S} and for climates with smaller P_{sub} (Fig. 6).
3.2.2 Saturated hydraulic conductivity
High K_{S} values are associated with TTDs that have higher initial values and lighter tails (Fig. 7). Also, a decrease in K_{S} causes more pronounced ups and downs in the TTD, with the effect of individual rainfall events being better discernible even in the later parts of the TTD (Fig. 8b). Increasing K_{S} by 2 orders of magnitude on average shortens Q_{1} by 44 %, Q_{2} by 58 %, the mTT by 59 % and Q_{3} by 62 % (Fig. 6 center, Table 1). The median F_{yw} increases by 13 %. The standard deviation increases with decreasing K_{S}, while the skewness and kurtosis both decrease significantly – TTDs become less skewed and more platykurtic (flatter). The interplay between K_{S} and θ_{ant} is obvious in that the influence of θ_{ant} decreases over time while the influence of K_{S} increases. Initially θ_{ant} controls the soil hydraulic conductivity, the partitioning of the tracer into surface and subsurface flow, and also the spreading within the soil. Later on, as moisture conditions become more similar for scenarios with identical P_{sub} and D_{soil}, K_{S} gains in importance while θ_{ant} becomes less relevant. The influence of K_{S} increases for wet θ_{ant} (especially for short transit times) and for big P_{sub} (especially for long transit times) since both maximize the differences in hydraulic conductivity between catchments – the drier the conditions, the more similar the unsaturated hydraulic conductivities in general (Fig. 6).
3.2.3 Subsequent precipitation amount
Big P_{sub} compresses the TTDs (Fig. 7). Doubling P_{sub} decreases Q_{1} by 63 % on average, Q_{2} by 61 %, the mTT by 57 % and Q_{3} by 58 % (Fig. 6 center, Table 1). The median F_{yw} increases by 22 %. The standard deviation (and hence the width) decreases by 42 %, while the skewness of the TTDs more than doubles. Bigger P_{sub} causes more leptokurtic (peaked) TTDs. Big amounts of P_{sub} increase the total flow through the catchment (both in the soil and bedrock) and hence control how effectively tracer is flushed out of the system. TTDs will have lighter tails and shorter mTTs mainly due to the fact that a bigger P_{sub} flushes the soils faster and only allows a smaller fraction of the precipitation events to infiltrate into the bedrock. The fraction of water entering the bedrock depends strongly on the contact time of that water with the soil–bedrock interface. That means that in regions with small P_{sub} a larger fraction of precipitation has the chance to infiltrate into the bedrock before it is flushed out of the soil layer by subsequent precipitation. Therefore the tails of TTDs in more arid regions tend to be heavier than the TTD tails in humid regions. The influence of P_{sub} is larger for dry θ_{ant} and high K_{S} (especially for the longer transit times) (Fig. 6).
3.2.4 Soil depth
Decreasing D_{soil} causes a larger fraction of tracer to arrive at the outlet faster (Fig. 8a). Halving D_{soil} shortens all the quartiles and the mTT of the TTDs on average by approximately 40 % (Fig. 6a, Table 1), while the median F_{yw} increases by 10 %. The standard deviation (the width of the TTD) is decreased by 19 % and the skewness is increased by about 56 %. Shallower soils cause more leptokurtic (peaked) TTDs, almost doubling the excess kurtosis. Shallower soils saturate faster than deeper soils; they also redirect tracer more quickly from vertical to lateral flow, and therefore the early response in shallower soils is slightly stronger. According to our findings, D_{soil} has only a small amount of influence on TTD shape. In catchments with deeper soils we should, however, expect longer transport times.
3.3 General observations on the shape of TTDs
The simulation results suggest that the TTDs can be visually divided into four distinct parts (Fig. 7), where the shape of three parts is clearly controlled by the catchment and climate properties and the fourth is a transition zone. The shape of the initial part of the TTD (up to ∼10 d) depends strongly on θ_{ant} and K_{S} (in accordance with Fiori et al., 2009) and less strongly on D_{soil}. TTDs in soils with wet θ_{ant} or high K_{S} exhibit higher initial peaks with a larger probability for short transit times. Starting after approximately 10 d a transition period follows where no individual parameter dominates. During this period precipitation drives the emptying of the uppermost soil layers with the presence of faster and/or larger flows (in catchments with higher K_{S}/bigger P_{sub}) being gradually compensated by higher remaining concentrations of tracer (in catchments with lower K_{S}/smaller P_{sub}), so that the tracer mass outflux at the catchment outlet converges towards a very similar value at around 120 days before diverging again. After the transition period, the shape of the TTDs is governed by P_{sub} (i.e., essentially the climate) and K_{S}, with larger P_{sub} and higher K_{S} causing a more rapid decline of outflow and hence a compression of the TTDs. Finally, the shape of the tails of the TTDs is controlled by the hydraulic conductivity of the bedrock K_{Br} (not the soil K_{S}) (see also Fiori et al., 2009). In many cases these tails constitute straight lines in the log–log plots (which is necessary but insufficient for identifying powerlaw functions). Furthermore, all modeled TTDs share one common feature – for every subsequent precipitation event there is a more or less discernible spike. Generally, larger subsequent events cause higher spikes (i.e., a higher proportion of outflow during those events) while the size of the spikes decreases at later times. And although this multitude of local maxima in the probability density curve does invoke a sense of irregularity, the general pattern of shapes of the TTDs is not influenced by the individual subsequent events (Fig. 5 and Table S4), which is why we decided to smooth the TTDs for visual comparison so that the underlying systematic changes in shapes are more clearly visible (see Fig. S8).
Practical implications can be drawn from our results concerning, for example, pollution events. Some catchments are more vulnerable to pollution in the sense that they tend to store pollutants for a longer period of time and hence exhibit long legacy effects. In particular catchments with TTDs with heavy tails belong in that category (i.e., catchments with deeper soils and a moderate hydraulic conductivity difference between soil and bedrock). Also, certain moments in time are worse for pollution events to happen – a spill occurring during dry conditions will stay in the catchment longer than a spill during wet conditions because it is more likely to reach the bedrock and stay in contact with it before it is flushed out of the soils. Accordingly, locations and situations that lead to a longer storage of decaying pollutants will eventually result in the release of fewer solutes downstream.
We also plotted the probability density replacing the actual transit time with the cumulative outflow to check whether this would eradicate the differences between the different distributions (see Fig. S11). We made two interesting observations: (1) For the scenarios with high K_{S}, the differences between the distributions were reduced considerably. For the cumulative probability distributions in particular there were hardly any discernible differences left. The largest discrepancies could still be found in the early part of the distributions where the distributions with high θ_{ant} continued to have larger outflow probabilities. (2) For the scenarios with low K_{S}, the individual distributions did not collapse into a single cumulative probability distribution. They rather split up into three distributions according to their P_{sub} values. That means that for the scenarios with larger P_{sub} a larger amount of cumulative outflow was necessary to flush out the same amount of tracer compared to the scenarios with smaller P_{sub}.
3.4 Distribution fitting
Shape parameters of the bestfit inverse Gaussian (D), gamma (α) and lognormal (σ) distributions as well as flow path numbers (F) for the 36 different scenarios are listed in Table 2. The parameters D, α and σ range from 0.15 to 0.98, from 0.78 to 3.66 and from 0.51 to 1.15, respectively. F ranges from −0.22 to 0.63. First we compared the performances of only these three probability distributions with two parameters. Out of the 36 model scenarios, the inverse Gaussian yielded the best fit 5 times, the gamma 13 times and the lognormal 18 times. In general, the lognormal works a little better for high K_{S}, dry θ_{ant} and small P_{sub} and the gamma for low K_{S}, wet θ_{ant} and big P_{sub}, while the inverse Gaussian is less ideal for capturing the shape of the modeled TTDs (Tables 3 and S5). Contrary to that, the inverse Gaussian represents the mean transit time (mTT) better than the other two distributions. On average, the mTT of the fitted gamma deviates from the observed mean by 24 % (88 d) with a maximum deviation of 423 d for one scenario, underpredicting for dry and overpredicting for wet θ_{ant}, while the inverse Gaussian performs much better in this regard with an average deviation from the mTT of only 5 % (17 d) and a maximum deviation of 102 d. The gamma especially underpredicts the mean when P_{sub} is small. The correct identification of the median transit time works much better for the gamma – here the average deviation of the fitted median from the observed median is only 4 % (12 d) with a maximum deviation of 59 d. The inverse Gaussian and lognormal yield average deviations from the median transit time of 6 % and 5 % (15 and 13 d) with maximum deviations of 50 and 43 d, respectively.
Then we included the two probability distributions with three parameters (beta, truncated lognormal) into the analysis and investigated how they compared to the twoparameter distributions. The performance of the beta was quite similar to the one of the gamma in terms of representing TTD shapes and the median transit times. However, it was able to capture the mTTs a lot better than the gamma, even surpassing the performance of the inverse Gaussian on average (average deviation 4 %, 13 d, maximum deviation 38 d), especially in environments with low K_{S} values. Finally, the truncated lognormal distribution performed best in every regard, capturing TTD shapes, mTTs and median transit times better than all other distributions (mTT average deviation 3 %, 10 d, maximum deviation 91 d; median transit time average deviation 4 %, 11 d, maximum deviation 36 d) (Table 3).
Figure 9 gives an overview of the shape and scale of our modeled TTDs (using the bestfit gamma distribution parameters).
3.5 Predicting the shape of TTDs
Figure 10 shows how the shape and scale of TTDs change with the individual catchment and climate properties. For increasing θ_{ant}, TTDs converge towards Lshaped distributions with shorter mTTs (in highly conductive soils the shape is more affected than the scale, in soils with low K_{S} the scale is more affected than the shape). When K_{S} is increasing, mTT is decreasing (in the case that P_{sub} is big, the shapes of the TTDs also change towards having lighter tails). Quite similar patterns can be observed for increasing D_{soil} and decreasing P_{sub} – with mTTs becoming longer and TTD shapes increasing in tail weight when K_{S} is high and becoming more humped when K_{S} is low.
Nonlinear regression analysis relating the shape and scale parameters of the fitted lognormal and gamma distributions to any single soil, precipitation or storage property (D_{soil}, K_{S}, θ_{ant}, P_{sub}) did not yield satisfying relations that could be used to predict TTD shapes. Here, we would like to present the significant nonlinear relationships we found between the shape parameters of the fitted TTDs and the flow path number F (R^{2}=0.90), mainly because we can draw much more general conclusions on TTD shapes using a dimensionless number (Fig. 11):
Generally, for similar catchments with low K_{S}, gamma distributions are more likely to fit the TTDs. The relatively higher proportion of surface flow within and surface outflow from these catchments seems to favor flow and transport dynamics that are best represented by the shapes of gamma distributions because they are able to capture both rapid response (high initial values) as well as the relatively slow outflow from the soils and the bedrock (long tails). In contrast, similar catchments with high K_{S} and only small proportions of surface flow are more likely to behave according to lognormal distributions with less rapid response from surface flow (low initial values) and faster outflow from the more conductive soils (higher and narrower modes at intermediate transit times). A notable exception is scenarios where catchments with highly conductive soils still experience larger proportions of surface outflow (> 25 %; F > 0.05) due to large amounts of P_{sub} – these dynamics cannot be predicted by the same relationship since they produce distributions with larger contributions of advective transport and lighter tails and hence smaller values of σ (indicated by the black circle in Fig. 11).
3.6 Effects of other factors on the shape of TTDs
3.6.1 Porosity
The influence that soil porosity exerts on the shape of TTDs is quite similar to the influence of D_{soil}. Larger soil porosity causes a dampening of the initial response and increasing transit times in all parts of the TTD (just like deeper soils; see Fig. 12a and Table S6). Increasing porosity also causes larger standard deviations, smaller skewness and smaller kurtosis (i.e., less peaked TTDs).
3.6.2 Hydraulic conductivity of the bedrock
Variations in the saturated hydraulic conductivity of the bedrock K_{Br} affect the shape of TTDs both in the initial part of the distributions but even more so in the tail (Fig. 12b and Table S7). If K_{Br} is increased so that it equals the K_{S} of the soil layer, we basically create one large continuum of homogeneous bedrock (or soil). Hence, the resulting TTD does not contain any abrupt breaks in slope and basically resembles outflow from a larger homogeneous reservoir. For lower K_{Br} breaks in the slope of the TTD tails start to appear indicating that the soil layers have already been emptied while the bedrock still contains water from the traced input precipitation event. For scenarios where K_{Br} is at least 3 orders of magnitude smaller than the soil K_{S}, the tails initially resemble powerlaw distributions with constants (a) around 0.2 and exponents (k) around 1.6 for longer periods of time:
An exponent k smaller than 2 indicates that a mean value of the powerlaw distribution cannot be defined since it is basically infinite; however, in our simulation results, the powerlaw tails eventually break down when the bedrock domain is almost empty. Somewhat counterintuitively, the scenario with the lowest K_{Br} (“very low”) exhibits the shortest quartile and mean transit times. This is clearly an effect of a smaller fraction of water infiltrating into the bedrock and more water being transported laterally in the relatively conductive soil layer. We observe the longest quartile transit times in the scenario where K_{Br} is 1 order of magnitude lower than K_{S} (“high”) and the longest mean transit time when it is 2 orders of magnitude lower (“med. high”). This is due to the fact that for these cases the higher K_{Br} causing faster transport within the bedrock is counterbalanced by the larger fraction of event water that enters into the bedrock, where it is transported more slowly than in the soil. Therefore what seems paradoxical in the first place – longer mTTs when K_{Br} is higher – can be explained by differences in the runoff partitioning between soil and bedrock. This also explains the observation that the standard deviation of the TTDs initially increases with increasing K_{Br} while both skewness and excess kurtosis decrease.
3.6.3 Decay in saturated hydraulic conductivity with depth
For catchments that already have highly conductive soils, adding a decay in K_{S} (with higher K_{S} close to the surface and lower K_{S} close to the soil–bedrock interface) does not change the shape of TTDs to a great extent – all shape metrics remain rather similar and transit times across the entire TTD are moderately shortened (Fig. 12c and Table S8). We observe a larger impact if soil K_{S} is low. In these cases adding a decay reduces the standard deviation and increases the skewness and the kurtosis of the resulting TTDs (i.e., they become narrower, more skewed and more peaked). Additionally, the difference in transit times increases towards the late part of the TTD with mTT and Q_{3} being considerably shorter when there is a decay in K_{S}. This difference between the smaller effects of a K_{S} decay in an already highly conductive soil compared to the larger effects for a low conductivity soil can be explained by the fact that the additional soil zones of higher conductivity are more effectively used for scenarios of generally low conductivity – in soils that are already quite conductive, a larger fraction of the incoming event water will still infiltrate to deeper soil layers before moving laterally, whereas in low conductivity soils the faster lateral transport possible due to the K_{S} decay will be triggered much sooner and for a larger fraction of the incoming event water.
3.6.4 Precipitation frequency
The shape of TTDs is not influenced significantly by precipitation frequency since the mean values of all distribution metrics for the lowfrequency (arid type) and the highfrequency (humid type) scenarios are quite similar to each other (Fig. 12d and Table S9). However, transit times in the highfrequency (humid) environment are shorter (${Q}_{\mathrm{1}}=\mathrm{17}$ %, ${Q}_{\mathrm{2}}=\mathrm{11}$ %, mTT $=\mathrm{9}$ %, ${Q}_{\mathrm{3}}=\mathrm{3}$ %). Additionally, the higher the precipitation frequency, the smaller the variation between individual TTDs. This is mainly due to two facts: when the precipitation frequency is high (1) the interarrival times are shorter, which will more often mobilize event water and avoid longer periods of relative inactivity when the water “just sits” in the soil, and (2) the amounts of precipitation events are on average smaller so that there is a smaller chance of a very big event “flushing” the entire system, creating very short transit times for a preceding event followed by a long period of no or only small precipitation events. These transit time dynamics with regard to different patterns of precipitation have already been observed in the field (Heidbüchel et al., 2013).
3.6.5 Water retention curve
The TTDs from the scenarios with sandtype WRCs have higher initial peaks and lighter tails compared to the ones with silttype WRCs (Fig. 13a and b). Their transit times are consistently shorter over all distributions, and the influence of other parameters (like K_{S} and θ_{ant}) on their shape is reduced. Sandtype TTDs are more skewed and more peaked than silttype TTDs (Table S10). Therefore they more closely resemble TTDs that we would expect in environments where preferential flow is present. Generally, the differences in TTDs between the different WRCs are more pronounced in the scenarios with low K_{S} because the wetting of the upper soil layers and hence the increase in the hydraulic conductivity takes relatively more time such that the differences between the two WRC scenarios are amplified. In the scenarios with silttype WRCs the saturation process causes a slower increase in hydraulic conductivity since soil water potential decreases more gently with increasing soil water content.
3.6.6 Catchment shape
We observe unexpectedly little variation between the TTDs of the differently shaped catchments (Fig. 13c). While Q_{1}, Q_{2} and the mTT are all more or less similar, Q_{3} increases slightly for catchments with a lower center of gravity and on average shorter flow paths (Table S11). The influence of the catchment shape is fractionally larger for dry θ_{ant}. Still, apparently the differences in catchment shape need to be a lot more pronounced than those we explored in order to significantly affect the TTD shape.
3.6.7 Full saturation and extreme precipitation
Starting runs with fully saturated soils increased the fractions of overland flow for both the high and the low K_{S} scenario (THSB and TLSB). For THSB the fraction of outflow during the first 10 d that was overland outflow (SOF_{10}) increased from 1 % to 9 %. For TLSB the increase was even higher, from 76 % to 91 %. The increase had clear effects on the resulting transit times. In particular the very short transit times increased in importance while the longer transit times were less affected. That means the changes we observed in the shape of the TTDs followed the pattern of increasing θ_{ant} (i.e., a higher percentage of increase in the young fraction of the TTD, smaller impact at later times, and shape metrics). Increasing the precipitation amount and intensity of the input event by a factor of 100 (+; from 0.1 to 10 mm h^{−1}) affected only the lowK_{S} scenario (TLSB+), further increasing the fraction of short transit times while the highK_{S} scenario was unaffected (THSB+). We had to increase the precipitation intensity of the input event by a factor of 1000 (to 100 mm h^{−1}) to eventually create substantial amounts of initial overland flow for both scenarios. Once this was triggered, the shape of the TTDs changed considerably. For these scenarios (THSB$+++$ and TLSB$+++$), all quartiles of the TTDs shortened to less than one day and the whole distribution became extremely leptokurtic (Fig. 13d and Table S12).
4.1 Use of theoretical distributions
The fact that TTDs under dry θ_{ant} are better represented by the (humped) lognormal distributions can be explained by the circumstance that the (rather empty) catchment storage has to be filled at least a little bit before faster flow paths are activated and substantial flow out of the system can occur. This means that the early response is much better captured by a distribution that starts with an initial value of close to 0. Furthermore, lognormal distributions also work better in highly conductive soils that produce TTD modes that are higher and narrower than the ones of gamma distributions. Contrary to that, low K_{S} values and wet θ_{ant} favor gamma distributions because initial outflow values are generally higher when the soil is closer to saturation while the TTD modes are lower and wider in soils that are less conductive (Fig. 14).
None of the theoretical distributions we tested captures the shape of all of the observed TTDs adequately over the entire age range. On the one hand, this is due to the misfit after the quite sudden break in slope at the tail end of the distributions; on the other hand – and this is more relevant from a mass balance perspective – it results from a misrepresentation of the initial response. Looking at Fig. 7, 8, 12 and 13, it becomes clear that all TTDs are humped distributions, with none of them exhibiting an initial maximum (with a monotonically decreasing limb afterwards) and none of them possessing a value of 0 after 1.5 min (the first time step reported). Since all inverse Gaussian and lognormal distributions start with a value of 0 and all gamma and beta distributions are either monotonically decreasing or start with a value of 0 they cannot be perfect representations of the modeled TTDs for porous media. Instead, a set of probability distributions – with initial values larger than 0, a rising limb to a maximum probability density and a falling limb with lighter or heavier tails – would theoretically be the best option to represent variable TTDs. We can confirm this expectation since the truncated lognormal distributions we tested do indeed capture the modeled TTD shapes best in most of our scenarios. Still they too are not able to reproduce the break in the TTD tails we observed in the model output after which the tails initially seem to follow a power law. This, however, does not constitute a substantial problem with regard to the correct mass balance since these heavier tails only comprise a very small fraction of the mass that was added to the system as a tracer. Still, if the tailing of the TTDs is relevant to a problem (e.g., when dealing with legacy contamination) one can add the observed breaks in the tails to the distributions (for a description see Sect. S7 and Fig. S7). As for the application of threeparameter distributions, although the beta model performed better than the twoparameter models overall (by a slim margin) we do not recommend using it due to its additional fitting parameter (the upper limit c) which increases equifinality problems (that we set out to eliminate). The same logic applies to the truncated lognormal distribution. It performs best in almost all regards (see Table 3) but is more difficult to parameterize (e.g., we found no good relationships between the parameters σ, λ and F), and no straightforward mathematical expressions exist that define its moments. Therefore we recommend utilizing the twoparameter lognormal distribution for high K_{S} and the gamma distribution for low K_{S} scenarios. When doing that, we have to be careful though and consider the distribution median as a more reliable transit time estimate than the mean (see Table 3).
Further theoretical developments should include the use of TTDs for nonconservative solute transport. This could be achieved by considering the TTD a basic function to which different reaction terms can be added (like “cutting the tail” of solutes that decay after a certain time in the catchment or shifting, damping and extending the TTD for solutes that experience retardation). An example is provided for an exponential decay reaction in Sect. S8 and Fig. S12.
4.2 Connection between the shape of TTDs and the flow path number F
We can pretty accurately predict the general shape of a TTD within the parameter range of our model scenarios using F alone (Fig. 11). Instead of using TTDs with constant shapes for determining variable transit times with transfer functionconvolution models, one can use these relationships to predefine the TTD shapes – reducing the problem of equifinality that stems from the simultaneous determination of shape and scale parameters (Fig. 15). Linked to that, some interesting conclusions can be drawn from the identified relationships between F and the shape parameters α and σ:

A flow path number between −1 and +1 characterizes catchments where the available storage is currently larger than the change in storage caused by the incoming and outgoing flows – over the characteristic timescale of the combined average interevent and event duration t_{Ie}+t_{Ev} (∼5 d).

If the system receives more water than it can remove during t_{Ie}+t_{Ev}, it is inflowdominated, F is positive and the shape of TTDs is generally better represented by gamma distributions.

With increasing F, α decreases to values below 1. This decrease in the shape parameter α is mainly caused by the initial peaks of the TTDs becoming higher. Our simulation results suggest that the tails of the TTDs become lighter with increasing positive F values. Therefore α should increase with increasing positive F values. The circumstance that we find a better relationship between increasing positive F and decreasing α values is due to the fact that – with regard to mass balance – the change in the initial response (higher initial values and peaks) outweighs the change in the tails (that are becoming lighter). Therefore we can conclude that the early response dominates TTD shapes (at least from a mass balance perspective).

If the system has the capacity to remove more water in the subsurface than it receives during t_{Ie}+t_{Ev}, it is outflowdominated, F becomes negative and the shape of the TTDs is generally better represented by lognormal distributions.

When F becomes more negative, σ increases from values around 0.5 to values above 1.0 (although the tails of the modeled TTDs become lighter), indicating higher peaks.

F converges towards 0 for systems with increasing available storage (because the denominator keeps increasing) or if inflows and outflow capacity are evenly balanced. For these cases both gamma and lognormal distributions become more and more dominated by smaller initial and early values as well as the later arrival of the peak concentration, which is illustrated by α becoming larger and by σ becoming smaller. This should not be interpreted as growing dominance of advective over dispersive transport because the TTD tails still become heavier in these situations.
The theoretical framework around the flow path number F can also be used to assess the impact that other catchment and climate properties have on TTD shapes. For example catchment size would only have an impact on TTD shape if the crosssectional area of the outflow boundary A_{out} changed disproportionately. If, for example, the catchment area A_{in} increased but the crosssectional area A_{out} remained the same, then the subsurface outflow capacity K_{rem} would decrease and hence F would change.
This research can also contribute to the field of catchment evolution. One could argue that in loworder catchments positive flow path numbers are not sustainable over longer periods of time because that would mean that the subsurface outflow capacity of the (zeroorder) catchment is permanently insufficient and the catchment is not capable of efficiently discharging all of the incoming precipitation via the subsurface. Consequently, the catchment storage would be filled up completely and overland flow would be occurring on a regular basis. Since widespread overland flow is rarely observed in most catchments it could be argued that most catchments have already evolved towards negative flow path numbers (e.g., by increasing K_{S} or D_{soil}). That, in turn, could also mean that Lshaped (or initially slightly humped) TTDs with heavier tails and gamma shape parameters α around 0.5 are the natural endpoint of catchment evolution.
4.3 Replacing transit time with cumulative outflow
For certain scenarios we still see differences in the probability distributions if we replace transit time with cumulative outflow (see Fig. S11). This observation can be explained by the fact that for the high K_{S} scenarios (where differences are reduced) we only generate external flow variability while for the low K_{S} scenarios (where differences remain) we also cause internal flow variability (Kim et al, 2016). That means that in the high K_{S} scenarios an increase in P_{sub} increases the flow in all of the available flow paths proportionally (without changing the flow path partitioning or activating previously unused flow paths) while for the low K_{S} scenarios an increase in P_{sub} causes pronounced shifts in the flow path partitioning where the additional amount of precipitation can bypass the subsurface by predominantly utilizing overland flow paths (leading to the observation that a larger amount of P_{sub} is necessary to flush out an equal amount of tracer). This can serve as direct proof that replacing transit time with cumulative outflow does not erase all differences between TTDs; however, it also shows that it may be adequate for many applications where large shifts in flow path partitioning are not expected.
4.4 Limitations and outlook
Our results can be considered valid for systems that do not experience a large fraction of preferential flow in the soil and bedrock since we only model flow taking place in the porous matrix of the subsurface domain. This is the likely reason that we also encounter α values that are larger than 1 – although such high α values were not found in previous studies (Hrachowitz et al., 2009; Godsey et al., 2010; Berghuijs and Kirchner, 2017; Birkel et al., 2016). Therefore, in terms of expanding the modeling effort, it would be very beneficial to include both evapotranspiration and macropore flow into the simulations. An inclusion of these processes will shift the flow path number F towards more negative values. On the one hand, evapotranspiration will provide an additional way to remove water from the subsurface (representing another sink term similar to K_{rem}) and macropore flow will enhance the subsurface outflow capacity of the catchment resulting in a shift towards TTDs with higher initial peaks. On the other hand, evapotranspiration also has the potential to reduce θ_{ant} below moisture levels obtainable with free drainage alone. This more extreme dryness could lead to even more humped TTDs with initial values closer to 0. The inclusion of additional heterogeneity in soil properties (layering, smallscale variations) would also be a worthwhile exercise that is, however, outside of the scope of our study. Therefore, since some of the potential shapecontrolling parameters are still excluded from the analysis (like, for example, K_{Br} or the precipitation event amount P_{Ev}), this study is not meant to represent the full and complete truth about TTD shapes. It is rather an attempt to find some structure in the way TTD shapes change with certain parameters and boundary conditions, an attempt to illuminate essential dynamics and to explore overarching principles in catchment hydrology. Therefore, the next important step is to verify the generality of these model findings and the resulting theory on catchment response with field observations. In particular since under many circumstances, e.g., in areas where soils are characterized by macropores and preferential flow pathways, traditional hydrological modeling (i.e., the applicability of the Richards equation) may not be suitable.
An interesting question that remains is whether backward TTDs can be linked to catchment and climate properties in a similar fashion to the one we used, since backward TTDs are comprised of many individual water inputs that entered the catchment over a very long period of time with potentially greatly varying initial conditions. That leads to the question of whether it is more important to know the conditions at the time of entry to the catchment or the conditions at the time of exit from the catchment (or both) in order to make predictions about TTD shapes and mTTs. Remondi et al. (2018) were among the first to tackle this problem by water flux tracking with a distributed model. They found that mainly soil saturation and groundwater storage affected backward TTDs.
In our simulations for a virtual loworder catchment we observed that the shape of TTDs changes systematically with the four investigated catchment and climate properties (D_{soil}, K_{S}, θ_{ant} and P_{sub}) so that it is possible to predict the change using the dimensionless flow path number F. The results can be summarized in three main conclusions (see also Fig. 11):

The shape of TTDs converges towards Lshaped distributions with high initial values if a catchment's capacity to store inflow decreases or if the actual inflow to a catchment does not equal its subsurface outflow capacity.

Heavier tails are produced when the system is in a more “relaxed” state, where all potential flow paths (deep and shallow, slower and faster) are equally used for transport. This is generally the case if P_{sub} is relatively small. Lighter tails appear when the system is in a more “stressed” state, where the shallow and faster flow paths are disproportionally used for transport. This can be associated with larger P_{sub} values. In addition, we observe a distinct break in the TTD tails if there is a sufficiently large difference in hydraulic conductivity between the bedrock K_{Br} and the soil K_{S}.

Gamma functions are able to capture the time variance of TTDs in an appropriate way, especially for low K_{S} and wet θ_{ant} scenarios, while lognormal distributions work well for high K_{S} and dry θ_{ant} scenarios.
However, neither gamma nor lognormal distributions are able to correctly represent the early part of the simulated distributions with nonzero initial values combined with a mode shortly after (i.e., the humped form) that we observe in most cases. Moreover, we noticed the general pattern that TTDs with high initial values tend to have lighter tails than TTDs with low initial values. Gamma distributions, unfortunately, exhibit the opposite behavior (with high initial values being associated with heavier tails than low initial values; see Fig. 16). Based on the results from our modeling efforts, we therefore encourage the exploration of betterfitting theoretical distributions. These distributions should be able to (a) represent high initial values paired with lighter tails as well as low initial values paired with heavier tails and (b) take on a “humped” form with nonzero initial values. We found that truncated distributions fulfil these requirements a lot better but have more degrees of freedom and are harder to parameterize.
Ideally, this work will help to generate new or to expand existing hypotheses on hydrologic and hydrochemical catchment response that can be tested in future field experiments.
All data used in this study are presented either in the main paper or in the Supplement.
The supplement related to this article is available online at: https://doi.org/10.5194/hess2428952020supplement.
IH, PT and TF conceptualized the study. Formal analysis was carried out by IH. Funding acquisition was organized by JHF. The investigation was carried out by IH, AM, JY and JHF. JY edited the software. IH wrote the original draft of the paper, and further writing, reviewing and editing was performed by IH, AM, JHF, JY, PT and TF.
The authors declare that they have no conflict of interest.
We would like to thank Carlotta Scudeler for her guidance on hydrologic modeling and her contribution to a previous version of this paper. Thanks also to René Therrien for his help with the HGS modeling and to Ilja van Meerveld and Stefanie Lutz for excellent discussions of the manuscript. Finally, we would like to acknowledge the work of at least six anonymous reviewers that provided necessary criticism and valuable suggestions for improvement.
This research was supported by the Helmholtz Research Programme
“Terrestrial Environment”, topic 3: “Sustainable Water Resources
Management”, with the integrated project: “Water and Matter Flux Dynamics
in Catchments”.
The article processing charges for this openaccess publication were covered by a Research Centre of the Helmholtz Association.
This paper was edited by Nunzio Romano and reviewed by three anonymous referees.
Ali, M., Fiori, A., and Russo, D.: A comparison of traveltime based catchment transport models, with application to numerical experiments, J. Hydrol., 511, 605–618, https://doi.org/10.1016/j.jhydrol.2014.02.010, 2014.
Ameli, A. A., Amvrosiadi, N., Grabs, T., Laudon, H., Creed, I. F., McDonnell, J. J., and Bishop, K.: Hillslope permeability architecture controls on subsurface transit time distribution and flow paths, J. Hydrol., 543, 17–30, https://doi.org/10.1016/j.jhydrol.2016.04.071, 2016.
Amin, I. E. and Campana, M. E.: A general lumped parameter model for the interpretation of tracer data and transit time calculation in hydrologic systems, J. Hydrol., 179, 1–21, https://doi.org/10.1016/00221694(95)028803, 1996.
Becker, M. W. and Shapiro, A. M.: Interpreting tracer breakthrough tailing from different forcedgradient tracer experiment configurations in fractured bedrock, Water Resour. Res., 39, 1024, https://doi.org/10.1029/2001WR001190, 2003.
Begemann, F. and Libby, W. F.: Continental water balance, ground water inventory and storage times, surface ocean mixing rates and worldwide water circulation patterns from cosmicray and bomb tritium, Geochim. Cosmochim. Ac., 12, 277–296, https://doi.org/10.1016/00167037(57)900406, 1957.
Benettin, P., Kirchner, J. W., Rinaldo, A., and Botter, G.: Modeling chloride transport using travel time distributions at Plynlimon, Wales, Water Resour. Res., 51, 32593276, https://doi.org/10.1002/2014WR016600, 2015.
Benettin, P., Soulsby, C., Birkel, C., Tetzlaff, D., Botter, G., and Rinaldo, A.: Using SAS functions and highresolution isotope data to unravel travel time distributions in headwater catchments, Water Resour. Res., 53, 1864–1878, https://doi.org/10.1002/2016WR020117, 2017.
Berghuijs, W. R. and Kirchner, J. W.: The relationship between contrasting ages of groundwater and streamflow, Geophys. Res. Lett., 44, 8925–8935, https://doi.org/10.1002/2017GL074962, 2017.
Birkel, C., Geris, J., Molina, M. J., Mendez, C., Arce, R., Dick, J., Tetzlaff, D., and Soulsby, C.: Hydroclimatic controls on nonstationary stream water ages in humid tropical catchments, J. Hydrol., 542, 231–240, https://doi.org/10.1016/j.jhydrol.2016.09.006, 2016.
Birkel, C., Soulsby, C., Tetzlaff, D., Dunn, S., and Spezia, L.: Highfrequency storm event isotope sampling reveals timevariant transit time distributions and influence of diurnal cycles, Hydrol. Process., 26, 308–316, https://doi.org/10.1002/hyp.8210, 2012.
Bishop, K., Seibert, J., Köhler, S., and Laudon, H.: Resolving the double paradox of rapidly mobilized old water with highly variable responses in runoff chemistry, Hydrol. Process., 18, 185–189, https://doi.org/10.1002/hyp.5209, 2004.
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, https://doi.org/10.1029/2009WR008371, 2010.
Botter, G., Bertuzzo, E., and Rinaldo, A.: Catchment residence and travel time distributions: The master equation, Geophys. Res. Lett., 38, L11403, https://doi.org/10.1029/2011GL047666, 2011.
Cardenas, M. B. and Jiang, X. W.: Groundwater flow, transport, and residence times through topographydriven basins with exponentially decreasing permeability and porosity, Water Resour. Res., 46, W11538, https://doi.org/10.1029/2010WR009370, 2010.
Chapman, S. W., Parker, B. L., Sale, T. C., and Doner, L. A.: Testing high resolution numerical models for analysis of contaminant storage and release from low permeability zones, J. Contam. Hydrol., 136, 106–116, https://doi.org/10.1016/j.jconhyd.2012.04.006, 2012.
Cox, D. R. and Isham, V.: A simple spatialtemporal model of rainfall, Philos. T. R. Soc. Lond., 415, 317–328, 1988.
DaneshYazdi, 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, https://doi.org/10.1002/hyp.11481, 2018.
Dinçer, T., Payne, B. R., Florkowski, T., Martinec, J., and Tongiorgi, E.: Snowmelt runoff from measurements of tritium and oxygen18, Water Resour. Res., 6, 110–124, https://doi.org/10.1029/WR006i001p00110, 1970.
Eriksson, E.: The Possible Use of Tritium for Estimating Groundwater Storage, Tellus, 10, 472–478, https://doi.org/10.3402/tellusa.v10i4.9265, 1958.
Fiori, A. and Becker, M. W.: Power law breakthrough curve tailing in a fracture: The role of advection, J. Hydrol., 525, 706–710, https://doi.org/10.1016/j.jhydrol.2015.04.029, 2015.
Fiori, A. and Russo, D.: Travel time distribution in a hillslope: Insight from numerical simulations, Water Resour. Res., 44, W12426, https://doi.org/10.1029/2008WR007135, 2008.
Fiori, A., Russo, D., and Di Lazzaro, M.: Stochastic analysis of transport in hillslopes: Travel time distribution and source zone dispersion, Water Resour. Res., 45, W08435, https://doi.org/10.1029/2008WR007668, 2009.
Gelhar, L. W., Welty, C., and Rehfeldt, K. R.: A critical review of data on fieldscale dispersion in aquifers, Water Resour. Res., 28, 1955–1974, https://doi.org/10.1029/92WR00607, 1992.
Gilfedder, B. S., Cartwright, I., Hofmann, H., and Frei, S.: Explicit Modeling of Radon222 in HydroGeoSphere During Steady State and Dynamic Transient Storage, Groundwater, 57, 3647, https://doi.org/10.1111/gwat.12847, 2019.
Godsey, S. E., Aas, W., Clair, T. A., de Wit, H. A., Fernandez, I. J., Kahl, J. S., Malcolm, I. A., Neal, C., Neal, M., Nelson, S. J., Norton, S. A., Palucis, M. C., Skjelkvåle, B. L., Soulsby, C., Tetzlaff, D., and Kirchner, J. W.: Generality of fractal 1∕f scaling in catchment tracer time series, and its implications for catchment travel time distributions, Hydrol. Process., 24, 1660–1671, https://doi.org/10.1002/hyp.7677, 2010.
Haggerty, R., McKenna, S. A., and Meigs, L. C.: On the latetime behavior of tracer test breakthrough curves. Water Resour. Res., 36, 3467–3479, https://doi.org/10.1029/2000WR900214, 2000.
Haitjema, H. M.: On the residence time distribution in idealized groundwatersheds, J. Hydrol., 172, 127–146, https://doi.org/10.1016/00221694(95)027325, 1995.
Harman, C. J.: Timevariable transit time distributions and transport: Theory and application to storagedependent transport of chloride in a watershed, Water Resour. Res., 51, 1–30, https://doi.org/10.1002/2014WR015707, 2015.
Harman, C. J. and Kim, M.: An efficient tracer test for timevariable transit time distributions in periodic hydrodynamic systems, Geophys. Res. Lett., 41, 1567–1575, https://doi.org/10.1002/2013GL058980, 2014.
Harman, C. J., Rao, P. S. C., Basu, N. B., McGrath, G. S., Kumar, P., and Sivapalan, M.: Climate, soil, and vegetation controls on the temporal variability of vadose zone transport, Water Resour. Res., 47, W00J13, https://doi.org/10.1029/2010WR010194, 2011.
Heidbüchel, I., Troch, P. A., and Lyon, S. W.: Separating physical and meteorological controls of variable transit times in zeroorder catchments, Water Resour. Res., 49, 7644–7657, https://doi.org/10.1002/2012WR013149, 2013.
Heidbüchel, I., Troch, P. A., Lyon, S. W., and Weiler, M.: The master transit time distribution of variable flow systems, Water Resour. Res., 48, W06520, https://doi.org/10.1029/2011WR011293, 2012.
Hrachowitz, M., Soulsby, C., Tetzlaff, D., Dawson, J. J. C., Dunn, S. M., and Malcolm, I. A.: Using longterm data sets to understand transit times in contrasting headwater catchments, J. Hydrol., 367, 237–248, https://doi.org/10.1016/j.jhydrol.2009.01.001, 2009.
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 timevariant transit time assessment, Water Resour. Res., 46, W10536, https://doi.org/10.1029/2010WR009148, 2010.
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, https://doi.org/10.1002/hyp.7922, 2011.
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, https://doi.org/10.5194/hess175332013, 2013.
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, Wiley Interdisciplinary Reviews: Water, 3, 629–657, https://doi.org/10.1002/wat2.1155, 2016.
Jasechko, S., Kirchner, J. W., Welker, J. M., and McDonnell, J. J.: Substantial proportion of global streamflow less than three months old, Nat. Geosci., 9, 126–129, https://doi.org/10.1038/NGEO2636, 2016.
Jiang, X. W., Wan, L., Wang, X. S., Ge, S., and Liu, J.: Effect of exponential decay in hydraulic conductivity with depth on regional groundwater flow, Geophys. Res. Lett., 36, L24402, https://doi.org/10.1029/2009GL041251, 2009.
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 timevarying flow paths: Direct observation of internal and external transport variability, Water Resour. Res., 52, 7105–7129, https://doi.org/10.1002/2016WR018620, 2016.
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, https://doi.org/10.5194/hess202792016, 2016.
Kirchner, J. W., Feng, X., and Neal, C.: Fractal stream chemistry and its implications for contaminant transport in catchments, Nature, 403, 524–527, https://doi.org/10.1038/35000537, 2000.
Kirchner, J. W., Feng, X., and Neal, C.: Catchmentscale advection and dispersion as a mechanism for fractal scaling in stream tracer concentrations, J. Hydrol., 254, 82–101, https://doi.org/10.1016/S00221694(01)004875, 2001.
Kollet, S. J. and Maxwell, R. M.: Demonstrating fractal scaling of baseflow residence time distributions using a fullycoupled groundwater and land surface model, Geophys. Res. Lett., 35, L07402, https://doi.org/10.1029/2008GL033215, 2008.
Kollet, S., Sulis, M., Maxwell, R. M., Paniconi, C., Putti, M., Bertoldi, G., Coon, E. T., Cordano, E., Endrizzi, S., Kikinzon, E., Mouche, E., Mügler, C., Park, Y.J., Refsgaard, J. C., Stisen, S., and Sudicky, E.: The integrated hydrologic model intercomparison project, IHMIP2: A second set of benchmark results to diagnose integrated hydrology and feedbacks, Water Resour. Res., 53, 867–890, https://doi.org/10.1002/2016WR019191, 2017.
Liggett, J. E., Partington, D., Frei, S., Werner, A. D., Simmons, C. T., and Fleckenstein, J. H.: An exploration of coupled surface–subsurface solute transport in a fully integrated catchment model, J. Hydrol., 529, 969–979, https://doi.org/10.1016/j.jhydrol.2015.09.006, 2015.
Lutz, S. R., Velde, Y. V. D., Elsayed, O. F., Imfeld, G., Lefrancq, M., Payraudeau, S., and Breukelen, B. M. V.: Pesticide fate on catchment scale: conceptual modelling of stream CSIA data, Hydrol. Earth Syst. Sc., 21, 5243–5261, https://doi.org/10.5194/hess2152432017, 2017.
Małoszewski, P., Rauert, W., Stichler, W., and Herrmann, A.: Application of flow models in an alpine catchment area using tritium and deuterium data, J. Hydrol., 66, 319–330, https://doi.org/10.1016/00221694(83)901932, 1983.
Małoszewski, P. and Zuber, A.: Determining the turnover time of groundwater systems with the aid of environmental tracers: 1. Models and their applicability, J. Hydrol., 57, 207–231, https://doi.org/10.1016/00221694(82)901470, 1982.
Massoudieh, A., Visser, A., Sharifi, S., and Broers, H. P.: A Bayesian modeling approach for estimation of a shapefree groundwater age distribution using multiple tracers, Appl. Geochem., 50, 252–264, https://doi.org/10.1016/j.apgeochem.2013.10.004, 2014.
Maxwell, R. M., Putti, M., Meyerhoff, S., Delfs, J.O., Ferguson, I. M., Ivanov, V., Kim, J., Kolditz, O., Kollet, S. J., Kumar, M., Lopez, S., Niu, J., Paniconi, C., Park, Y.J., Phanikumar, M. S., Shen, C., Sudicky, E. A., and Sulis, M.: Surfacesubsurface model intercomparison: A first set of benchmark results to diagnose integrated hydrology and feedbacks, Water Resour. Res., 50, 1531–1549, https://doi.org/10.1002/2013WR013725, 2014.
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, https://doi.org/10.1002/hyp.7796, 2010.
McGuire, K. J., McDonnell, J. J., Weiler, M., Kendall, C., McGlynn, B. L., Welker, J. M., and Seibert, J.: The role of topography on catchmentscale water residence time, Water Resour. Res., 41, https://doi.org/10.1029/2004WR003657, 2005.
McGuire, K. J. and McDonnell, J. J.: A review and evaluation of catchment transit time modeling, J. Hydrol., 330, 543–563, https://doi.org/10.1016/j.jhydrol.2006.04.020, 2006.
McMillan, H., Tetzlaff, D., Clark, M., and Soulsby, C.: Do timevariable tracers aid the evaluation of hydrological model structure? A multimodel approach, Water Resour. Res., 48, W05501, https://doi.org/10.1029/2011WR011688, 2012.
Musolff, A., Fleckenstein, J. H., Rao, P. S. C., and Jawitz, J. W.: Emergent archetype patterns of coupled hydrologic and biogeochemical responses in catchments, Geophys. Res. Lett., 44, 4143–4151, https://doi.org/10.1002/2017GL072630, 2017.
Nauman, E. B.: Residence time distribution theory for unsteady stirred tank reactors, Chem. Eng. Sci., 24, 1461–1470, https://doi.org/10.1016/00092509(69)850748, 1969.
Niemi, A. J.: Residence time distributions of variable flow processes, Int. J. Appl. Radiat. Is., 28, 855–860, https://doi.org/10.1016/0020708X(77)900266, 1977.
Nir, A.: On the interpretation of tritium “age” measurements of groundwater, J. Geophys. Res., 69, 2589–2595, https://doi.org/10.1029/JZ069i012p02589, 1964.
Nyström, U.: Transit time distributions of water in two small forested catchments, Ecol. Bull., 37, 97–100, 1985.
Pangle, L. A., Kim, M., Cardoso, C., Lora, M., Meira Neto, A. A., Volkmann, Wang, Y., Troch, P. A., and Harman, C. J.: The mechanistic basis for storagedependent age distributions of water discharged from an experimental hillslope, Water Resour. Res., 53, 2733–2754, https://doi.org/10.1002/2016WR019901, 2017.
Pedretti, D. and Bianchi, M.: Reproducing tailing in breakthrough curves: Are statistical models equally representative and predictive?, Adv. Water Resour., 113, 236–248, https://doi.org/10.1016/j.advwatres.2018.01.023, 2018.
Pedretti, D., FernàndezGarcia, D., Bolster, D., and SanchezVila, X.: On the formation of breakthrough curves tailing during convergent flow tracer tests in threedimensional heterogeneous aquifers, Water Resour. Res., 49, 4157–4173, https://doi.org/10.1002/wrcr.20330, 2013.
PeraltaTapia, A., Soulsby, C., Tetzlaff, D., Sponseller, R., Bishop, K., and Laudon, H.: Hydroclimatic influences on nonstationary transit time distributions in a boreal headwater catchment, J. Hydrol., 543, 7–16, https://doi.org/10.1016/j.jhydrol.2016.01.079, 2016.
Remondi, F., Kirchner, J. W., Burlando, P., and Fatichi, S.: Water Flux Tracking With a Distributed Hydrological Model to Quantify Controls on the Spatiotemporal Variability of Transit Time Distributions, Water Resour. Res., 54, 3081–3099, https://doi.org/10.1002/2017WR021689, 2018.
Rinaldo, A., Beven, K. J., Bertuzzo, E., Nicotina, L., Davies, J., Fiori, A., Russo, D., and Botter, G.: Catchment travel time distributions and water flow in soils, Water Resour. Res., 47, W07537, https://doi.org/10.1029/2011WR010478, 2011.
RoaGarcía, M. C. and Weiler, M.: Integrated response and transit time distributions of watersheds by combining hydrograph separation and longterm transit time modeling, Hydrol. Earth Syst. Sci., 14, 1537–1549, https://doi.org/10.5194/hess1415372010, 2010.
Rodhe, A., Nyberg, L., and Bishop, K.: Transit times for water in a small till catchment from a step shift in the oxygen 18 content of the water input, Water Resour. Res., 32, 3497–3511, https://doi.org/10.1029/95WR01806, 1996.
SchulzeMakuch, D.: Longitudinal dispersivity data and implications for scaling behavior, Groundwater, 43, 443–456, https://doi.org/10.1111/j.17456584.2005.0051.x, 2005.
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, https://doi.org/10.5194/hess1847512014, 2014.
Soulsby, C., Birkel, C., Geris, J., and Tetzlaff, D.: Spatial aggregation of timevariant stream water ages in urbanizing catchments, Hydrol. Process., 29, 3038–3050, https://doi.org/10.1002/hyp.10500, 2015.
Sprenger, M., Seeger, S., Blume, T., and Weiler, M.: Travel times in the vadose zone: Variability in space and time, Water Resour. Res., 52, 5727–5754, https://doi.org/10.1002/2015WR018077, 2016.
Stockinger, M. P., Bogena, H. R., Lücke, A., Diekkrueger, B., Weiler, M., and Vereecken, H.: Seasonal soil moisture patterns: Controlling transit time distributions in a forested headwater catchment, Water Resour. Res., 50, 5270–5289, https://doi.org/10.1002/2013WR014815, 2014.
Sudicky, E. A., Illman, W. A., Goltz, I. K., Adams, J. J., and McLaren, R. G.: Heterogeneity in hydraulic conductivity and its role on the macroscale transport of a solute plume: From measurements to a practical application of stochastic flow and transport theory, Water Resour. Res., 46, W01508, https://doi.org/10.1029/2008WR007558, 2010.
Tetzlaff, D., Birkel, C., Dick, J., Geris, J., and Soulsby, C.: Storage dynamics in hydropedological units control hillslope connectivity, runoff generation, and the evolution of catchment transit time distributions, Water Resour. Res., 50, 969–985, https://doi.org/10.1002/2013WR014147, 2014.
Therrien, R., McLaren, R. G., Sudicky, E. A., and Panday, S. M.: HydroGeoSphere: a threedimensional numerical model describing fullyintegrated subsurface and surface flow and solute transport, Groundwater Simulations Group, University of Waterloo, Waterloo, ON, 2010.
Timbe, E., Windhorst, D., Celleri, R., Timbe, L., Crespo, P., Frede, H.G., Feyen, J., and Breuer, L.: Sampling frequency tradeoffs in the assessment of mean transit times of tropical montane catchment waters under semisteadystate conditions, Hydrol. Earth Syst. Sci., 19, 1153–1168, https://doi.org/10.5194/hess1911532015, 2015.
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, https://doi.org/10.1029/2010WR009105, 2010.
van der Velde, Y., Torfs, P. J. J. F., van der Zee, S. E. A. T. M., and Uijlenhoet, R.: Quantifying catchmentscale mixing and its effect on timevarying travel time distributions, Water Resour. Res., 48, W06536, https://doi.org/10.1029/2011WR011310, 2012.
van der Velde, Y., Heidbüchel, I., Lyon, S. W., Nyberg, L., Rodhe, A., Bishop, K., and Troch, P. A.: Consequences of mixing assumptions for timevariable travel time distributions, Hydrol. Process., 29, 3460–3474, https://doi.org/10.1002/hyp.10372, 2015.
Weiler, M., McGlynn, B. L., McGuire, K. J., and McDonnell, J. J.: How does rainfall become runoff? A combined tracer and runoff transfer function approach, Water Resour. Res., 39, 1315, https://doi.org/10.1029/2003WR002331, 2003.
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, https://doi.org/10.1002/2017WR020894, 2017.
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, https://doi.org/10.1002/2017WR021896, 2018.
Zech, A., Attinger, S., Cvetkovic, V., Dagan, G., Dietrich, P., Fiori, A., Rubin, J., and Teutsch, G.: Is unique scaling of aquifer macrodispersivity supported by field data?, Water Resour. Res., 51, 7662–7679, https://doi.org/10.1002/2015WR017220, 2015.
Zhang, Y., Benson, D. A., and Baeumer, B.: Predicting the tails of breakthrough curves in regionalscale alluvial systems, Groundwater, 45, 473–484, https://doi.org/10.1111/j.17456584.2007.00320.x, 2007.