Journal topic
Hydrol. Earth Syst. Sci., 24, 2895–2920, 2020
https://doi.org/10.5194/hess-24-2895-2020
Hydrol. Earth Syst. Sci., 24, 2895–2920, 2020
https://doi.org/10.5194/hess-24-2895-2020

Research article 03 Jun 2020

Research article | 03 Jun 2020

# On the shape of forward transit time distributions in low-order catchments

On the shape of forward transit time distributions in low-order catchments
Ingo Heidbüchel1, Jie Yang1, Andreas Musolff1, Peter Troch2, Ty Ferré2, and Jan H. Fleckenstein1 Ingo Heidbüchel et al.
• 1Department of Hydrogeology, Helmholtz Centre for Environmental Research – UFZ, 04318 Leipzig, Germany
• 2Department of Hydrology and Atmospheric Sciences, University of Arizona, Tucson, 85721, USA

Correspondence: Ingo Heidbüchel (ingo.heidbuechel@ufz.de)

Abstract

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 low-order catchment and investigate whether it changes systematically with certain catchment and climate properties. To this end, we used a physically based, spatially explicit 3-D model, injected tracer with a precipitation event and recorded the resulting forward TTDs at the outlet of a small (∼6000 m2) 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 best-fit theoretical probability functions could describe the entire TTD shape exactly. Still, we found that generally gamma and log-normal distributions work better for scenarios of low and high soil hydraulic conductivity, respectively.

1 Introduction

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 (pre-event) and young (event) water to streamflow (Heidbüchel et al., 2012). TTDs are time and space-variant 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 first-order 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 piston-flow 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 heavy-tailed 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 well-mixed system is best represented by an exponential distribution, partial mixing can be approximated with gamma distributions and no mixing with the piston-flow model) (van der Velde et al., 2015).

## 1.3 Controls on shape variations

A number of studies reported on the best-fit 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 L-shaped 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 (Peralta-Tapia 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 power-law form and fractal behavior adding macrodispersion and systematic heterogeneity to the domain in the form of depth-decreasing 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; Roa-Garcí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 best-fit 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 flow-weighted 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), flow-weighted 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 transit-time-based catchment transport models (where they compare several time-invariant to time-variable methods), conclude that applying a flow-weighted time approach can indeed yield adequate results for predicting catchment-scale 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; Danesh-Yazdi et al., 2018; Yang et al., 2018) and master transit time distributions (MTTDs) (representing the flow-weighted 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, 3-D 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.

2 Methods

We used HydroGeoSphere (HGS), a 3-D numerical model describing fully coupled surface–subsurface, variably saturated flow and advective–dispersive solute transport (Therrien et al., 2010). Groundwater flow in the 3-D subsurface is simulated with Richards' equation and Darcy's law, and surface runoff is simulated in the 2-D surface domain with Manning's equation and the diffusive-wave approximation of the Saint-Venant 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 real-world 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 zero-order catchment was set up, 100 m long, 75 m wide (∼6000 m2), 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 m3 m−3, van Genuchten parameters alpha αvG=0.5 m−1, beta βvG=1.6, saturated water content θs= 0.39 m3 m−3, residual water content θr= 0.05 m3 m−3, pore-connectivity parameter lp= 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 Schulze-Makuch (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.

Figure 13-D model domain and shape of the virtual catchment from the top (a), front (b) and side (c). The blue square indicates the outflow boundary with constant head condition. The red layer represents the soil, which has a much higher hydraulic conductivity than the underlying bedrock (grey). The orange lines indicate the zone of convergence (but no explicit channel). The two additional catchment shapes (top-heavy and bottom-heavy) we tested in Sect. 2.2.1 are shown in the black box (d).

### 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 follow-up 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.

Figure 2(a) 1-year time series of subsequent precipitation (looped 32 more times for the entire modeling period and rescaled for smaller or larger amounts of subsequent precipitation Psub). Tracer application took place during the first hour of the model runs. (b) Time series of subsequent precipitation for a high-frequency scenario (humid) and a low-frequency scenario (arid). The total precipitation amount is the same for both scenarios.

### 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 Seff≈50 %), an intermediate one (θant=28.8 %; Seff≈70 %) and a wet one (θant=35.6 %; Seff≈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 Seff 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 Seff values (≈30 % in the dry scenarios). That means that although an Seff 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 KS and almost 1500 years for the scenarios with low KS. 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 (Dsoil), saturated soil hydraulic conductivity (KS), antecedent soil moisture content (θant) and subsequent precipitation amount (Psub, essentially a measure of the amount of precipitation that falls after the delivery of the traced event) (Fig. 3).

Figure 3The four properties that were varied to explore their influence on the shape and scale of TTDs: soil depth Dsoil, saturated soil hydraulic conductivity KS, antecedent soil moisture θant and subsequent precipitation amount Psub. The bedrock hydraulic conductivity KBr was kept constant for all of these base-case scenarios.

We tested two soil depths Dsoil, namely depths of 0.5 and 1.0 m, evenly distributed across the entire catchment. Similarly, we chose two saturated soil hydraulic conductivities KS, 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 Psub 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 Dsoil (T = thick; F = flat), the second one KS (H = high; L = low), the third one θant (W = wet; I = intermediate; D = dry) and the fourth one Psub (S = small; M = medium; B = big). For example the name FHIB would indicate a run with a flat (“F”) (shallow) soil, a high (“H”) KS, 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.

Table 1Metrics of the TTDs derived from the modeling of 36 scenarios with different combinations of catchment and climate properties. All times are given in days.

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 m3 m−3) and three had smaller soil porosity (0.24 m3 m−3) than the base-case scenarios (0.39 m3 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 KBr of 10−5 m d−1). In the first run KBr was decreased to 10−7 m d−1, in the following runs it was successively increased to 10−3, 10−2, 10−1, 100 and 2×100 m d−1, matching the KS 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 KS with a shape parameter f= 0.29 m and saturated hydraulic conductivity at the surface KS0= 7 m d−1 (for the high conductivity scenarios) or KS0=0.07 m d−1 (for the low conductivity scenarios) (Fig. 4a):

$\begin{array}{}\text{(1)}& {K}_{\mathrm{S}}\left(z\right)={K}_{\mathrm{S}\mathrm{0}}{e}^{-\frac{z}{f}}.\end{array}$

This preserved the mean KS values of $\mathrm{2}×{\mathrm{10}}^{-\mathrm{0}}$ (high) and $\mathrm{2}×{\mathrm{10}}^{-\mathrm{2}}$ m d−1 (low) (from the base-case scenarios).

Figure 4(a) Exponential decay in saturated soil hydraulic conductivity with depth for the high (blue) and the low (red) KS scenario. The x axis in the inset has a log scale. The spatial mean KS is indicated by the vertical black lines. (b) Water retention curves (solid) and relative hydraulic conductivities (dotted) for sandy and silty soils. The permanent wilting point (PWP) and the field capacity (FC) are marked as references (dashed).

### 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 high-frequency scenarios (comparable to a humid precipitation distribution and intensity pattern with lower intensities and more frequent events) and low-frequency 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 Psub scenarios.

### 2.2.5 Water retention curve

All the base-case model scenarios were conducted with water retention curves (WRCs) resembling silty soils:

$\begin{array}{}\text{(2)}& \mathit{\theta }={\mathit{\theta }}_{\mathrm{r}}+\frac{{\mathit{\theta }}_{\mathrm{s}}-{\mathit{\theta }}_{\mathrm{r}}}{{\left[\mathrm{1}+{\left({\mathit{\alpha }}_{\mathrm{vG}}\left|\mathit{\psi }\right|\right)}^{{\mathit{\beta }}_{\mathrm{vG}}}\right]}^{\mathit{\nu }}},\end{array}$

with van Genuchten parameters αvG (m−1) and βvG (dimensionless), saturated water content θs, residual water content θr (both m3 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 sand-type WRC since it can, in some aspects and to a certain extent, also indicate how a system with the threshold-like initiation of rapid preferential flow behaves. The sand-type WRC causes an increase in hydraulic conductivity at relatively lower soil water contents compared to the silt-type 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 kr was derived with Eq. (3):

$\begin{array}{}\text{(3)}& {k}_{r}={S}_{\mathrm{eff}}^{{l}_{\mathrm{p}}}{\left[\mathrm{1}-{\left(\mathrm{1}-{S}_{\mathrm{eff}}^{{\mathit{\nu }}^{-\mathrm{1}}}\right)}^{\mathit{\nu }}\right]}^{\mathrm{2}},\end{array}$

with effective saturation Seff and pore-connectivity parameter lp (both dimensionless). Other aspects of preferential flow – like bypass flow through macropores in deeper soil layers – are, however, not captured by sand-type WRCs. The van Genuchten parameters for the sand-type 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., Seff 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 Psub 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}×{\mathrm{10}}^{-\mathrm{2}}±\mathrm{7.2}×{\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 Jout (kg d−1) for each time step by the total inflow mass Min (kg):

$\begin{array}{}\text{(4)}& \mathrm{TTD}\left(t\right)={J}_{\mathrm{out}}^{\mathrm{norm}}\left(t\right)=\frac{{J}_{\mathrm{out}}\left(t\right)}{{M}_{\mathrm{in}}}.\end{array}$

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:

$\begin{array}{}\text{(5)}& {\mathrm{TTD}}_{\mathrm{cml}}\left(t\right)={\sum }_{t=\mathrm{0}}^{t}\left({J}_{\mathrm{out}}^{\mathrm{norm}}\left(t\right)\cdot \mathrm{\Delta }t\right).\end{array}$

### 2.4.2 Calculation of TTD metrics

For each TTD we calculated seven metrics to characterize its shape: the first quartile (Q1), the median (Q2), the mean (mTT), the third quartile (Q3), 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 Fyw 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 Fyw 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 log-normal distribution which also has just two adjustable parameters:

1. 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(\frac{\mathrm{4}\mathit{\pi }Dt}{\mathrm{mTT}}\right)}^{-\mathrm{0.5}}\frac{\mathrm{1}}{t}\mathrm{exp}\mathit{\left\{}-\left[{\left(\mathrm{1}-\frac{t}{\mathrm{mTT}}\right)}^{\mathrm{2}}\cdot \frac{\mathrm{mTT}}{\mathrm{4}Dt}\right]\mathit{\right\}}.\end{array}$
2. 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}}\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., sub-exponential) 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.

3. Log-normal 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)=\frac{\mathrm{1}}{t\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\mathrm{exp}\left[-\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:

1. Three-parameter beta distribution with shape parameters α and β (dimensionless) and upper limit c (d) (with mean mTT $=\mathit{\alpha }c/\left(\mathit{\alpha }+\mathit{\beta }\right)\right)$:

$\begin{array}{}\text{(9)}& \mathrm{beta}\left(t\right)=\frac{{t}^{\mathit{\alpha }-\mathrm{1}}\left(c-t{\right)}^{\mathit{\beta }-\mathrm{1}}}{{c}^{\mathit{\alpha }+\mathit{\beta }-\mathrm{1}}B\left(\mathit{\alpha },\mathit{\beta }\right)}.\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.

2. Truncated log-normal distribution with the time of truncation λ (d) as the third parameter:

$\begin{array}{}\text{(10)}& \begin{array}{rl}& \mathrm{Trunc}\left(t\right)=\left\{\frac{\mathrm{1}}{\left(t+\mathit{\lambda }\right)\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\mathrm{exp}\left[-\frac{{\left(\mathrm{ln}\left(t+\mathit{\lambda }\right)-\mathit{\mu }\right)}^{\mathrm{2}}}{\mathrm{2}{\mathit{\sigma }}^{\mathrm{2}}}\right]\right\}/\\ & \left\{\mathrm{1}-{\int }_{t=\mathrm{0}}^{\mathit{\lambda }}\frac{\mathrm{1}}{t\mathit{\sigma }\sqrt{\mathrm{2}\mathit{\pi }}}\mathrm{exp}\left[-\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 (m3):

$\begin{array}{}\text{(11)}& F\left(t\right)=\frac{{P}_{\mathrm{dr}}\left(t\right)-{K}_{\mathrm{rem}}}{{D}_{\mathrm{soil}}\left(n-{\mathit{\theta }}_{\mathrm{ant}}\left(t\right)\right){A}_{\mathrm{in}}},\end{array}$

where soil depth Dsoil (m), catchment surface area Ain (m2), porosity n (m3 m−3) and antecedent moisture content θant (m3 m−3) are paired with the driving precipitation amount Pdr (m3), which is calculated as the average subsequent precipitation amount Psub (m a−1) over the average event duration tEv (d):

$\begin{array}{}\text{(12)}& {P}_{\mathrm{dr}}\left(t\right)=\frac{{t}_{\mathrm{Ev}}{P}_{\mathrm{sub}}\left(t\right){A}_{\mathrm{in}}}{\mathrm{365.25}}.\end{array}$

The subsequent precipitation amount Psub (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 tEv instead of the interevent duration tIe to compute Pdr 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 Krem (m3) consists of the effective saturated soil hydraulic conductivity KS (m d−1), the sum of the average interevent and event duration tIe+tEv (d), the porosity n (m3 m−3), and the cross-sectional area of the soil layer at the outlet of the catchment Aout (m2):

$\begin{array}{}\text{(13)}& {K}_{\mathrm{rem}}=\left({t}_{\mathrm{Ie}}+{t}_{\mathrm{Ev}}\right){K}_{\mathrm{S}}n{A}_{\mathrm{out}}.\end{array}$

The cross-sectional area of the soil layer at the outlet of the catchment Aout 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 tIe+tEv instead of just the event duration tEv to compute Krem 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 Pdr 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 wet-up phase of a catchment and change more slowly (towards 0) during the dry-down 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 cross-sectional area of the soil layer at the outlet Aout) can be larger than the increase in storage capacity itself – leading to F becoming even more negative with increasing storage capacity.

3 Results

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 Q1 observed for the 15 scenarios with different event sequences is still 14 % of the total range observed for the 36 base-case scenarios, this percentage decreases down to 2 % for Q3. 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 base-case 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.

Figure 515 TTDs resulting from 15 different precipitation time series with all other catchment and climate properties being equal. The first few events have the largest influence on the TTD shapes, while subsequent events gradually even out the differences. Inset shows cumulative distributions.

## 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, KS 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 KS values gains more importance later when the soil moisture conditions become more similar. Dsoil and Psub 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 KS is a lot stronger in scenarios with wet θant while the influence of Psub decreases with increasing θant. Panel (c) shows that both θant and Psub have a larger influence when KS is high, but for Psub this increased influence is only seen for the longer transit times. The influence of the initial condition θant is larger when KS is high because the relative differences in flow through a dry soil and a wet soil are larger for soils with high KS compared to soils with low KS. Panel (d) confirms the impression that Dsoil 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 Psub has opposite effects on the influence of θant and KS: Larger Psub causes the influence of KS to increase for the longer transit times while the influence of θant decreases when Psub 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 base-case scenarios.

Figure 6Influence of different properties on different parts of the TTD. Shown is the average percentage decrease in transit time for each quartile (Q1, Q2, Q3) and the mean (μ) of the TTD caused by a decrease in Dsoil from 1 to 0.5 m (green), an increase in KS from 0.02 to 2 m d−1 (purple), an increase in θant from 50 % to 90 % effective saturation Seff (red) and an increase in Psub from 0.3 to 1.4 m a−1 (blue). Panel (a) shows the average decrease in transit time for changing each of the four properties, and panels (b–e) show the decrease in transit time conditional on the variation of one of the four properties (θant, KS, Dsoil and Psub). Two examples are illustrated by the black circles: (1) The dashed blue line in (c) shows that the increase in Psub has a larger influence on the third quartile transit time (Q3) – a decrease of ∼75 % instead of just ∼50 % – for a catchment with a high KS compared to a catchment with a low KS. (2) The thick red line in (e) shows that the increase in θant from 50 % to 90 % Seff has a smaller influence on the second quartile transit time (Q2) – a decrease of just ∼15 % instead of ∼35 % – for a catchment with a big Psub compared to a catchment with a small Psub.

### 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 Seff 50 % to 90 %), on average Q1 decreases by 44 %, Q2 by 27 %, the mTT by 19 % and Q3 by 15 % (Fig. 6 center, Table 1). The median Fyw 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. Long-term trends or interannual shifts in Psub can cause temporal changes in TTDs but substantial short-term 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 KS and for climates with smaller Psub (Fig. 6).

Figure 7Results of the 36 model runs. TTDs are grouped by soil depth (a and b= deep (thick); c and d= shallow (flat)) and soil hydraulic conductivity (a and c = high; b and d = low). Yellow colors indicate dry, green intermediate and blue wet antecedent moisture conditions; thick lines indicate large, mid-sized lines medium and thin lines small amounts of subsequent precipitation amounts. Dashed black lines divide TTDs into four parts, each part controlled by different properties. Note the log–log axes. Insets show cumulative TTDs.

### 3.2.2 Saturated hydraulic conductivity

High KS values are associated with TTDs that have higher initial values and lighter tails (Fig. 7). Also, a decrease in KS 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 KS by 2 orders of magnitude on average shortens Q1 by 44 %, Q2 by 58 %, the mTT by 59 % and Q3 by 62 % (Fig. 6 center, Table 1). The median Fyw increases by 13 %. The standard deviation increases with decreasing KS, while the skewness and kurtosis both decrease significantly – TTDs become less skewed and more platykurtic (flatter). The interplay between KS and θant is obvious in that the influence of θant decreases over time while the influence of KS 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 Psub and Dsoil, KS gains in importance while θant becomes less relevant. The influence of KS increases for wet θant (especially for short transit times) and for big Psub (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 Psub compresses the TTDs (Fig. 7). Doubling Psub decreases Q1 by 63 % on average, Q2 by 61 %, the mTT by 57 % and Q3 by 58 % (Fig. 6 center, Table 1). The median Fyw increases by 22 %. The standard deviation (and hence the width) decreases by 42 %, while the skewness of the TTDs more than doubles. Bigger Psub causes more leptokurtic (peaked) TTDs. Big amounts of Psub 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 Psub 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 Psub 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 Psub is larger for dry θant and high KS (especially for the longer transit times) (Fig. 6).

### 3.2.4 Soil depth

Decreasing Dsoil causes a larger fraction of tracer to arrive at the outlet faster (Fig. 8a). Halving Dsoil shortens all the quartiles and the mTT of the TTDs on average by approximately 40 % (Fig. 6a, Table 1), while the median Fyw 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, Dsoil has only a small amount of influence on TTD shape. In catchments with deeper soils we should, however, expect longer transport times.

Figure 8Influence of soil depth (a) and saturated soil hydraulic conductivity (b) on the shape of TTDs. Lighter shades of one color indicate shallower soils; dashed lines indicate higher hydraulic conductivity. Insets show cumulative TTDs.

## 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 KS (in accordance with Fiori et al., 2009) and less strongly on Dsoil. TTDs in soils with wet θant or high KS 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 KS/bigger Psub) being gradually compensated by higher remaining concentrations of tracer (in catchments with lower KS/smaller Psub), 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 Psub (i.e., essentially the climate) and KS, with larger Psub and higher KS 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 KBr (not the soil KS) (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 power-law 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 KS, 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 KS, the individual distributions did not collapse into a single cumulative probability distribution. They rather split up into three distributions according to their Psub values. That means that for the scenarios with larger Psub a larger amount of cumulative outflow was necessary to flush out the same amount of tracer compared to the scenarios with smaller Psub.

## 3.4 Distribution fitting

Shape parameters of the best-fit inverse Gaussian (D), gamma (α) and log-normal (σ) 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 log-normal 18 times. In general, the log-normal works a little better for high KS, dry θant and small Psub and the gamma for low KS, wet θant and big Psub, 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 Psub 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 log-normal 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.

Table 2Shape parameters of the best-fit inverse Gaussian (D), gamma (α) and log-normal (σ) distributions and associated flow path numbers (F) for the 36 different scenarios.

Then we included the two probability distributions with three parameters (beta, truncated log-normal) into the analysis and investigated how they compared to the two-parameter 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 KS values. Finally, the truncated log-normal 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).

Table 3Average and maximum deviations of mean and median transit times between the best-fit theoretical probability distributions and the modeled TTDs (given as the ratio of average deviation of the fitted distributions to the average modeled mean and median transit times as well as the average deviation in days). The sum of the squared residuals indicates the goodness of fit between the shape of theoretical probability distributions and modeled TTDs.

Figure 9Gamma shape parameters (α) and mean transit times (mTTs) for individual scenarios with different combinations of catchment and climate properties. The red boxes contain exemplary gamma distributions with shape and scale corresponding to the red dot location. The dotted black line marks the shape parameter value of 1 that corresponds to an exponential distribution.

Figure 9 gives an overview of the shape and scale of our modeled TTDs (using the best-fit 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 L-shaped distributions with shorter mTTs (in highly conductive soils the shape is more affected than the scale, in soils with low KS the scale is more affected than the shape). When KS is increasing, mTT is decreasing (in the case that Psub is big, the shapes of the TTDs also change towards having lighter tails). Quite similar patterns can be observed for increasing Dsoil and decreasing Psub – with mTTs becoming longer and TTD shapes increasing in tail weight when KS is high and becoming more humped when KS is low.

Figure 10Change of gamma shape parameters (α) and mean transit times (mTTs) for four catchment and climate properties: yellow colors indicate dry, green intermediate and blue wet θant; thick marker lines indicate large, mid-sized lines medium and thin lines small Psub; solid lines indicate low, dashed lines high KS; lighter shades of a color indicate shallow, darker shades deep Dsoil. The dotted black line marks the shape parameter value of 1 that corresponds to an exponential distribution.

Nonlinear regression analysis relating the shape and scale parameters of the fitted log-normal and gamma distributions to any single soil, precipitation or storage property (Dsoil, KS, θant, Psub) 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 (R2=0.90), mainly because we can draw much more general conclusions on TTD shapes using a dimensionless number (Fig. 11):

$\begin{array}{}\text{(14)}& \begin{array}{rl}& \mathrm{shape}\phantom{\rule{0.25em}{0ex}}\mathrm{parameter}\phantom{\rule{0.25em}{0ex}}\mathit{\alpha }\left(F\right)=\mathrm{0.64}{\left|F\right|}^{-\mathrm{0.20}},\\ & \mathrm{if}\phantom{\rule{0.25em}{0ex}}{K}_{\mathrm{S}}<\mathrm{0.2}\phantom{\rule{0.125em}{0ex}}m{d}^{-\mathrm{1}},\end{array}\text{(15)}& \begin{array}{rl}& \mathrm{shape}\phantom{\rule{0.25em}{0ex}}\mathrm{parameter}\phantom{\rule{0.25em}{0ex}}\mathit{\sigma }\left(F\right)=\mathrm{0.12}\mathrm{ln}\left|F\right|+\mathrm{1.19},\\ & \mathrm{if}\phantom{\rule{0.25em}{0ex}}{K}_{\mathrm{S}}\ge \mathrm{0.2}\phantom{\rule{0.125em}{0ex}}m{d}^{-\mathrm{1}}.\end{array}\end{array}$

Generally, for similar catchments with low KS, 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 KS and only small proportions of surface flow are more likely to behave according to log-normal 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 Psub – 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).

Figure 11Relationship between the dimensionless flow path number F and the shape parameters α (upper panel, scenarios with low KS) and σ (lower panel, scenarios with high KS) of the gamma and the log-normal distribution, respectively. The dotted trend lines are the best-fit regressions for the relationship between the flow path number and the shape parameters α (light blue) and σ (orange). The points in the black circle are excluded from the regression analysis since they are associated with scenarios of excessive surface outflow.

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

Figure 12Overview of how certain catchment and climate characteristics influence the shape of TTDs. (a) Porosity – solid lines indicate small, dotted lines large porosity. (b) Hydraulic conductivity of the bedrock – being equal or lower than the KS of the soil layer. (c) Decay in saturated soil hydraulic conductivity with depth – darker shades of one color represent scenarios with decay, lighter shades scenarios without decay. (d) Precipitation frequency – orange TTDs are low-frequency (“arid type”) scenarios, blue TTDs are high-frequency (“humid type”) scenarios. The shaded areas between the lines illustrate the higher shape variability for the low-frequency TTDs. Insets show cumulative TTDs.

### 3.6.2 Hydraulic conductivity of the bedrock

Variations in the saturated hydraulic conductivity of the bedrock KBr 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 KBr is increased so that it equals the KS 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 KBr 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 KBr is at least 3 orders of magnitude smaller than the soil KS, the tails initially resemble power-law distributions with constants (a) around 0.2 and exponents (k) around 1.6 for longer periods of time:

$\begin{array}{}\text{(16)}& \mathrm{TTD}\left(t\right)=a{t}^{-k}.\end{array}$

An exponent k smaller than 2 indicates that a mean value of the power-law distribution cannot be defined since it is basically infinite; however, in our simulation results, the power-law tails eventually break down when the bedrock domain is almost empty. Somewhat counterintuitively, the scenario with the lowest KBr (“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 KBr is 1 order of magnitude lower than KS (“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 KBr 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 KBr 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 KBr 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 KS (with higher KS close to the surface and lower KS 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 KS 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 Q3 being considerably shorter when there is a decay in KS. This difference between the smaller effects of a KS 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 KS 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 low-frequency (arid type) and the high-frequency (humid type) scenarios are quite similar to each other (Fig. 12d and Table S9). However, transit times in the high-frequency (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 sand-type WRCs have higher initial peaks and lighter tails compared to the ones with silt-type WRCs (Fig. 13a and b). Their transit times are consistently shorter over all distributions, and the influence of other parameters (like KS and θant) on their shape is reduced. Sand-type TTDs are more skewed and more peaked than silt-type 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 KS 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 silt-type WRCs the saturation process causes a slower increase in hydraulic conductivity since soil water potential decreases more gently with increasing soil water content.

Figure 13Overview of how certain catchment and climate characteristics influence the shape of TTDs (continued). Panels (a–b): Water retention curves (WRCs) – light blue and yellow lines indicate silt-type soil WRCs, dark blue and orange lines indicate sand-type soil WRCs. (a) Scenarios with high KS, (b) scenarios with low KS. (c) Catchment shape – lighter shades of a color indicate top-heavy, darker shades bottom-heavy catchments. (d) Full saturation and extreme precipitation – black lines indicate fully saturated initial conditions, pink lines fully saturated initial conditions and very large event precipitation (+), red lines fully saturated initial conditions and extreme event precipitation ($+++$). The horizontal lines in the box above the diagram indicate periods where actual overland flow was recorded during the respective runs. The insets show the cumulative TTDs.

### 3.6.6 Catchment shape

We observe unexpectedly little variation between the TTDs of the differently shaped catchments (Fig. 13c). While Q1, Q2 and the mTT are all more or less similar, Q3 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 KS scenario (THSB and TLSB). For THSB the fraction of outflow during the first 10 d that was overland outflow (SOF10) 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 low-KS scenario (TLSB+), further increasing the fraction of short transit times while the high-KS 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 Discussion

## 4.1 Use of theoretical distributions

The fact that TTDs under dry θant are better represented by the (humped) log-normal 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, log-normal 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 KS 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).

Figure 14Modeled TTDs for low KS with high θant (blue) and high KS with low θant (yellow). Best-fit theoretical distributions (dotted lines) for the individual scenarios for the log-normal model (a, c) and the gamma model (b, d). Breaks in the tails of the modeled TTDs are marked by the solid black lines. Small panels show cumulative TTDs.

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 log-normal 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 log-normal 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 three-parameter distributions, although the beta model performed better than the two-parameter 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 log-normal 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 two-parameter log-normal distribution for high KS and the gamma distribution for low KS 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.

Figure 15Predicted TTD shapes based on their relationship to the flow path number F, resulting from different antecedent moisture conditions θant (from blue – wet – on the left to yellow – dry – on the right) and subsequent precipitation amounts Psub. TTDs for low KS are gamma distributions (b); for high KS they are log-normal distributions (c). Individual TTDs start with time shifts so that they do not overlap (individual start times correspond to the Psub markers in a).

## 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 function-convolution models, one can use these relationships to pre-define 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 σ:

1. 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 tIe+tEv (∼5 d).

2. If the system receives more water than it can remove during tIe+tEv, it is inflow-dominated, F is positive and the shape of TTDs is generally better represented by gamma distributions.

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

4. If the system has the capacity to remove more water in the subsurface than it receives during tIe+tEv, it is outflow-dominated, F becomes negative and the shape of the TTDs is generally better represented by log-normal distributions.

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

6. 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 log-normal 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 cross-sectional area of the outflow boundary Aout changed disproportionately. If, for example, the catchment area Ain increased but the cross-sectional area Aout remained the same, then the subsurface outflow capacity Krem would decrease and hence F would change.

This research can also contribute to the field of catchment evolution. One could argue that in low-order catchments positive flow path numbers are not sustainable over longer periods of time because that would mean that the subsurface outflow capacity of the (zero-order) 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 KS or Dsoil). That, in turn, could also mean that L-shaped (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 KS scenarios (where differences are reduced) we only generate external flow variability while for the low KS scenarios (where differences remain) we also cause internal flow variability (Kim et al, 2016). That means that in the high KS scenarios an increase in Psub 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 KS scenarios an increase in Psub 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 Psub 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 Krem) 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, small-scale variations) would also be a worthwhile exercise that is, however, outside of the scope of our study. Therefore, since some of the potential shape-controlling parameters are still excluded from the analysis (like, for example, KBr or the precipitation event amount PEv), 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.

Figure 16Gamma distributions (solid lines) capture the middle part of the modeled TTDs (dashed lines; thickness corresponds to Psub amount) quite well but do not correctly represent the initial parts, breaks in the tails and heavier tails. Inset: gamma distributions (thick and thin black solid lines) combine either high initial values with heavier tails or zero initial values with lighter tails while modeled TTDs often are best described by high initial values and lighter tails (blue dashed line) or low (albeit nonzero) initial values with heavier tails (yellow dashed line).

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.

5 Conclusion

In our simulations for a virtual low-order catchment we observed that the shape of TTDs changes systematically with the four investigated catchment and climate properties (Dsoil, KS, θant and Psub) 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):

1. The shape of TTDs converges towards L-shaped 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.

2. 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 Psub 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 Psub 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 KBr and the soil KS.

3. Gamma functions are able to capture the time variance of TTDs in an appropriate way, especially for low KS and wet θant scenarios, while log-normal distributions work well for high KS and dry θant scenarios.

However, neither gamma nor log-normal 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 better-fitting 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.

Data availability
Data availability.

All data used in this study are presented either in the main paper or in the Supplement.

Supplement
Supplement.

Author contributions
Author contributions.

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.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

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.

Financial support
Financial support.

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 open-access publication were covered by a Research Centre of the Helmholtz Association.

Review statement
Review statement.

This paper was edited by Nunzio Romano and reviewed by three anonymous referees.

References

Ali, M., Fiori, A., and Russo, D.: A comparison of travel-time 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/0022-1694(95)02880-3, 1996.

Becker, M. W. and Shapiro, A. M.: Interpreting tracer breakthrough tailing from different forced-gradient 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 world-wide water circulation patterns from cosmic-ray and bomb tritium, Geochim. Cosmochim. Ac., 12, 277–296, https://doi.org/10.1016/0016-7037(57)90040-6, 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, 3259-3276, https://doi.org/10.1002/2014WR016600, 2015.

Benettin, P., Soulsby, C., Birkel, C., Tetzlaff, D., Botter, G., and Rinaldo, A.: Using SAS functions and high-resolution isotope data to unravel travel time distributions in headwater catchments, Water Resour. Res., 53, 1864–1878, 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 non-stationary 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.: High-frequency storm event isotope sampling reveals time-variant 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 topography-driven 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 spatial-temporal model of rainfall, Philos. T. R. Soc. Lond., 415, 317–328, 1988.

Danesh-Yazdi, M., Klaus, J., Condon, L. E., and Maxwell, R. M.: Bridging the gap between numerical solutions of travel time distributions and analytical storage selection functions, Hydrol. Process., 32, 1063–1076, 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 oxygen-18, 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 field-scale 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 Radon-222 in HydroGeoSphere During Steady State and Dynamic Transient Storage, Groundwater, 57, 36-47, 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 late-time 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/0022-1694(95)02732-5, 1995.

Harman, C. J.: Time-variable transit time distributions and transport: Theory and application to storage-dependent transport of chloride in a watershed, Water Resour. Res., 51, 1–30, https://doi.org/10.1002/2014WR015707, 2015.

Harman, C. J. and Kim, M.: An efficient tracer test for time-variable 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 zero-order 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 long-term 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 time-variant 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/hess-17-533-2013, 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 time-varying 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/hess-20-279-2016, 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.: Catchment-scale advection and dispersion as a mechanism for fractal scaling in stream tracer concentrations, J. Hydrol., 254, 82–101, https://doi.org/10.1016/S0022-1694(01)00487-5, 2001.

Kollet, S. J. and Maxwell, R. M.: Demonstrating fractal scaling of baseflow residence time distributions using a fully-coupled 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, IH-MIP2: 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/hess-21-5243-2017, 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/0022-1694(83)90193-2, 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/0022-1694(82)90147-0, 1982.

Massoudieh, A., Visser, A., Sharifi, S., and Broers, H. P.: A Bayesian modeling approach for estimation of a shape-free 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.: Surface-subsurface 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 catchment-scale 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 time-variable 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/0009-2509(69)85074-8, 1969.

Niemi, A. J.: Residence time distributions of variable flow processes, Int. J. Appl. Radiat. Is., 28, 855–860, https://doi.org/10.1016/0020-708X(77)90026-6, 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 storage-dependent 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àndez-Garcia, D., Bolster, D., and Sanchez-Vila, X.: On the formation of breakthrough curves tailing during convergent flow tracer tests in three-dimensional heterogeneous aquifers, Water Resour. Res., 49, 4157–4173, https://doi.org/10.1002/wrcr.20330, 2013.

Peralta-Tapia, A., Soulsby, C., Tetzlaff, D., Sponseller, R., Bishop, K., and Laudon, H.: Hydroclimatic influences on non-stationary 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 Spatio-temporal 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.

Roa-García, M. C. and Weiler, M.: Integrated response and transit time distributions of watersheds by combining hydrograph separation and long-term transit time modeling, Hydrol. Earth Syst. Sci., 14, 1537–1549, https://doi.org/10.5194/hess-14-1537-2010, 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.

Schulze-Makuch, D.: Longitudinal dispersivity data and implications for scaling behavior, Groundwater, 43, 443–456, https://doi.org/10.1111/j.1745-6584.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/hess-18-4751-2014, 2014.

Soulsby, C., Birkel, C., Geris, J., and Tetzlaff, D.: Spatial aggregation of time-variant 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 three-dimensional numerical model describing fully-integrated 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 trade-offs in the assessment of mean transit times of tropical montane catchment waters under semi-steady-state conditions, Hydrol. Earth Syst. Sci., 19, 1153–1168, https://doi.org/10.5194/hess-19-1153-2015, 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 catchment-scale mixing and its effect on time-varying 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 time-variable 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 regional-scale alluvial systems, Groundwater, 45, 473–484, https://doi.org/10.1111/j.1745-6584.2007.00320.x, 2007.