Catchment response times – understanding runoff dynamics from catchment distances and celerities

ABSTRACT In this study we explore how varying the river network (RN) density affects the distribution of hillslope to RN distances, the subsurface water celerities and hence the response times. Eleven Norwegian catchments (with areas of 0.007 to 500 km2) were used for the analysis, and the Distance Distribution Dynamics (DDD) hydrological model was calibrated for each catchment and RN. Equally good Kling-Gupta efficiency scores suggest a degree of equifinality in that many constellations of RNs and subsurface celerities have equally good model performance. All catchments display a linear relationship between the calibrated mean subsurface water celerity and mean hillslope to RN distance, consistent with a constant mean response time (MRT). The MRTs vary from 1 to 49 days for the different catchments and agree with MRT estimated from recession analysis and regionalized through regression and catchment characteristics. The latter aids in the estimation of model parameters for ungauged basins.


Introduction
The combined temporal delays of flow in hillslopes and in the river network (RN) constitute the hydrological response of a catchment to a rainfall/snowmelt event. Being able to describe and forecast these temporal delays, or response times (RTs, i.e. the time spent from when rainfall/meltwater enters the catchment until the response can be measured at the catchment's outlet), is really at the heart of the study and modelling of runoff dynamics. When discussing time scales in catchment hydrology, it is important to distinguish between transit time, involving flow velocities, and time of hydrological response, which is the propagation of hydraulic potential and involves flow celerities. These two time scales can differ by orders of magnitude (McGuire and McDonnell 2006, Botter et al. 2010, McDonnell and Beven 2014. In this paper, we discuss RT as a function of distances and flow celerities.
If distances and celerities were readily measurable features, the determination of RT would be simple; flow waves propagate a certain distance, d [m], with a certain celerity, v [m/s], and RT is thus: However, it turns out that characterizing the distances and celerities in play in the catchment response to rainfall/snowmelt events is not at all straightforward. The distribution of distances between points in the hillslope and the river network, hereafter called hillslope to RN distances, is clearly dependent on the configuration of the RN within the catchment. A denser RN will necessarily give shorter hillslope to RN distances. Digitized RNs can be derived from the blue lines of maps or from field observations, or generated automatically from regular grid digital elevation models (DEMs) (Montgomery and Dietrich 1988, 1992, Gandolfi and Bischetti 1997, Tarboton 1997. These methods have issues concerning objectivity, practicality, and realism, respectively, but the latter method is commonly employed for obvious practical reasons and takes advantage of the increasingly available high-resolution, highquality DEMs (Bhowmik et al. 2015, Reddy et al. 2018. One of the problems with automatic methods for deriving an RN from a DEM is the specification of a critical support area (CA), which is the area collecting sufficient water to initiate and maintain a channel (Tarboton et al. 1991, Gandolfi et al. 1999, Lin et al. 2005, Colombo et al. 2007, Reddy et al. 2018. Since the CA clearly expresses hydrological and geomorphological conditions, there have been many studies devoted to developing methods for objectively determining CAs both from hydrological considerations (Lin et al. 2005) and from terrain features, typically the local slope (Montgomery and Dietrich 1988, 1992, Tarboton et al. 1991, Tarboton 1997, Colombo et al. 2007). However, various investigations have found that CA varies in space (Montgomery and Dietrich 1988, Lin et al. 2005, Colombo et al. 2007, Reddy et al. 2018 and in time (Biswal and Marani 2010, Godsey and Kirchner 2014, Papageorgaki and Nalbantis 2017, van Meerveld et al. 2019, and no preferred method for the automatic extraction of CA has, as of yet, been produced. A commonly applied but very coarse simplification is to specify CA as an arbitrary, constant value for a catchment (Lin et al. 2005). Some authors have also investigated the distributions of hillslope distances and of distances within the RN. Tucker et al. (2001) and Skaugen and Onof (2014) describe the hillslope distribution of distances as exponential, which, for constant flow celerity and when ignoring the influence of transport in the RN, explains the exponential recession found for linear reservoirs Mengistu 2016, Skaugen andLawrence 2017). Skaugen and Onof (2014) found that the normal distribution approximated the distribution of distances found within the RN very well.
To summarize the above review, there is no objective method for the determination of the correct, "true," RN. The "true" RN may not even exist since the distributions of both hillslope-RN distances and distances within the RN change for different levels of saturation. Several authors report how the saturation of a catchment influences the density of a catchment's RN, giving a denser RN with increasing saturation (Godsey and Kirchner 2014, van Meerveld et al. 2019, Botter and Durighetto 2020. It is difficult to determine the distribution of catchmentscale celerities through experiments. The hydrological response of a catchment is shaped by the propagation of hydraulic potential in the subsurface, on the surface when there is overland flow, and in the RN. The runoff dynamics observed at the outlet is hence due to a mix of responses which differs from catchment to catchment and as a function of saturation (Robinson et al. 1995, McDonnell andBeven 2014).
Since neither hillslope to RN distances, and distances within the RN, nor catchment-scale flow celerities can be easily observed or determined for a catchment, we choose to study catchment RTs with the use of a rainfall-runoff model. According to the review of catchment transit time modelling by McGuire and McDonnell (2006), this is a common strategy for inferring transit time distributions. We use the Distance Distributions Dynamics (DDD) model (Skaugen and Onof 2014), which is suited for such a task since it parameterizes the distance distributions of both the hillslope to RN distances and distances within the RN and calibrates the subsurface and streamflow celerities. The model has, in total, eight model parameters requiring calibration, of which five are associated with the meteorological input and snowmelt. Only three parameters are hence linked to runoff dynamics and can be well defined through calibration.
In this paper we investigate and characterize the mean response time (MRT) for a selection of Norwegian catchments. Using Geographic Information System (GIS) we derive different RNs for 11 catchments of varying catchment size distributed across southern Norway. A range of CAs are used so that the resulting RNs have varying densities. For each RN and catchment, the DDD model is calibrated, and we investigate how the different RNs affect the distance distributions, the calibrated subsurface and stream celerities, and the RTs. Durighetto et al. (2020) and van Meerveld et al. (2019) performed investigations of a similar nature, in which they related the distance statistics of river lengths and hillslope to RN distances to observed extensions and contractions of RNs for small catchments. The extensions and contractions were caused by varying catchment wetness. In our investigation, we have fixed RNs, and we simulate runoff dynamics by modelling subsurface flow celerity as a function of catchment wetness. A very pertinent and interesting question (which will not be addressed here, but in a future paper) is what role the extensions and contractions of RNs play in runoff dynamics.
We further use recession analysis as an alternative method for quantifying the MRT in our analysis. Previously, when implementing the DDD model, recession characteristics together with distance distributions were used to estimate subsurface celerities under the assumption that subsurface hillslope processes dominate the hydrological response measured at the outlet (Robinson et al. 1995, Botter and Rinaldo 2003, Skaugen and Onof 2014. Since runoff in rain-free periods reflects the drainage of stored water, recession characteristics are directly related to the RTs of a catchment (Stoelzle et al. 2013). RTs derived from recession behaviour are independent of RNs and subsurface celerities and represent, hence, an alternative assessment of MRT.
Finally, we investigate the dependencies between RTs, climate characteristics and catchment characteristics and provide regionalized estimates of MRT. In Norway, the west coast is typically wet with a mean annual precipitation of > 3000 mm and with very shallow glacial and fluvial deposits available for subsurface storage. The inland in Southern Norway, however, can be rich in deposits, but has a relatively dry climate with a mean annual precipitation of <1000 mm. We expect the variation in climatic forcing, soil properties and catchment characteristics to reflect and discriminate the catchment MRT. With regionalized estimates of MRT, model parameters controlling runoff dynamics can be estimated for ungauged catchments.

Study area, DEMs and meteorological and hydrological data
Eleven catchments located in south Norway are used in this study. Five of the catchments are inland with a relatively dry continental climate and with thick glaciofluvial deposits, and six catchments are located on the south and west coast in a wet and maritime climate with poor glaciofluvial deposits (see Fig. 1). See the Norwegian Geological Survey (https:// geo.ngu.no/kart/losmasse_mobil/) for information on deposits and Met.no (http://www.senorge.no/index.html?p=klima) for climatic information. Ten of the catchments vary in size between 1.9 and 432 km 2 , whereas one catchment, 8.88 Muren, is very small, only 0.0075 km 2 . The names and ID numbers of the catchment are given in Table 2. Hereafter, the catchments are only referred to by their ID number from the HYDRAII database of the Norwegian water Resources and Energy Directorate (NVE). The purpose of including such a variety of catchment areas is to investigate the effect of spatial and temporal scales of RN on RTs.
The DDD model calibrated for the nine larger catchments is forced with 24-hourly gridded temperature and precipitation data provided by the Norwegian Meteorological Institute (MET) for the period of 1 September 1999-31 December 2018. The dataset is based on observations from the climate database of MET and has been gridded to 1 × 1 km 2 using a modified optimal interpolation scheme (Lussana et al. 2018). The model is calibrated against 24hourly discharge data provided by the HYDRAII database of the NVE. The two smaller catchments, 8.88 and 28.11, have only short time series of meteorological and hydrological data (April-October 2020 and May-October 2017), with a temporal resolution of 15 and 10 minutes, respectively.
For the larger catchments, we use a DEM with a 10 m resolution, and for the small catchments (28.11 and 8.88), we have used DEMs with a resolution of 1 m and 0.25 m, respectively. The DEMs are downloaded from the Norwegian DEM data repository (https://www.kartverket.no/data/ hoydedata. no/Laserinnsyn). Furthermore, we use the national land resource database (AR50) landcover dataset from the Norwegian Institute of Bioeconomy Research (NIBIO) (https://nibio.no/tema/jord/arealressurser/ar50).
Using a toolkit for drainage networks available in GRASS GIS (Jasiewicz and Metz 2011), we generate GIS layers containing catchment boundaries, flow directions and flow accumulations. RNs are extracted for specified upslope CAs. For seven of the catchments, we extracted a minimum of six RNs using CAs of 0.1, 0.25, 0.5, 1, 1.5 and 2 km 2 . For the small catchments RNs were extracted using CAs of 10, 100, 500, 1000, 2500 and 5000 m 2 for 8.88 and 500, 1000, 2000, 3500, 5000 and 7500 m 2 for 28.11. For catchment 2.279, five RNs are extracted for CAs of 0.25, 0.5, 1, 1.5 and 2 km 2 , and for 12.178 eight RNs are extracted using CAs of 0.005, 0.01, 0.025, 0.05, 0.1, 0.25, 0.5 and 1.0 km 2 . For each GIS-derived RN, we estimate the statistics of hillslope to RN distances and of the distances within the RN. We find that the mean of the hillslope to RN distances, d HS , is linearly related to CA 0.5 with an explained variance ranging from 96% to 100% for the catchments (see Fig. 2). Consequently, the different RNs are, hereafter, quantitively represented by d HS .  Table 2.

The DDD model
The DDD model is a conceptual, semi-distributed, catchmentbased rainfall-runoff model, scripted in the programming language Julia (https://julialang.org/). DDD simulates hydrology at temporal scales of 1 minute to 24 hours and at spatial scales from 50 m 2 to 12 760 km 2 (Viker-Walsøe and Valle 2020, Nazeer et al. 2021). The model was developed at NVE Onof 2014, Skaugen andMengistu 2016) and is in operational use at the Norwegian flood forecasting service. The main objective for developing the model was to reduce the number of model parameters requiring calibration, while maintaining the reliability and required detail of the simulations. Most model parameters are hence derived from digitized maps, but eight parameters are calibrated against runoff (see Table 1). The model input is temperature and precipitation values for each of the 10 elevation zones. The model simulates runoff, the spatial distribution of snow water equivalent (SWE), snow melt, glacial melt, potential and actual evapotranspiration, and saturated and unsaturated soil water. Snowmelt, glacial melt, and evapotranspiration are calculated using an energy-balance (EB) approach (Skaugen and Saloranta 2015). The EB approach applies proxy models to calculate the EB elements (short-and longwave radiation, turbulent fluxes, etc.), but the input data is still only precipitation and temperature. Table 1 shows the parameters calibrated using an automatic calibration routine embedded in a Julia package (https://github.com/robertfeldt/BlackBoxOptim.jl), and Fig. 3 shows the model structure.  In this study, correction factors for liquid and solid precipitation (Pkorr and Skorr in Table 1) are calibrated in addition to those in boldface in Fig. 3. The saturated zone of the subsurface, Soil SZ, consists of five linear reservoirs arranged in parallel. The linear reservoirs are filled from below (which is the slowest) and contribute to runoff according to their saturation state. The linear reservoirs -or, rather, the instantaneous unit hydrographs (IUHs) associated with their recession limb (Pedersen et al. 1980) -are parameterized from the hillslope distance distribution and subsurface flow celerity. The shape of the distribution of hillslope to RN distances determines the shape of the IUH, and the subsurface celerity (one for each IUH) determines the scale (see Onof 2014, Skaugen andMengistu 2016). In the original DDD model, the subsurface celerities were estimated from recession analysis using a recession characteristic Λ sampled as the log of the empirical recession rate (Skaugen and Onof 2014): Table 2. Estimated MRT HS and range of mean distances and mean calibrated celerities using different RNs derived for CAs between 0.1 and 2 km 2 . Catchments 8.88 and 28.11 have separate sets of CAs, 10-5000 m 2 and 500-7500 m 2 , respectively. Catchment 2.279 uses CAs of 0.25-2 km 2 and catchment 12.178 uses CAs of 0.005-1 km 2 . where Q t ð Þ is runoff at time t and Δt is the temporal resolution of the runoff observations. The shape and scale parameter of the gamma distributed Λ are the parameters G Sh and G Sc in Table 1. For estimating the individual celerities of the linear reservoirs from Λ, see Skaugen and Onof (2014). The mean subsurface celerity, v, is estimated as where the mean value of Λ is estimated as Λ ¼ G Sh � G Sc . Values of G Sh and G Sc (and hence ΛÞ can either be estimated from observed recession events Onof 2014, Skaugen andMengistu 2016) or calibrated against runoff. If Λ is estimated from recession, the impact on catchment response from flow in the RN is assumed to be negligible. The distributions of hillslope to RN distances, and of distances within the RN, are parameterized by the mean of the exponential distance distribution, and the mean and standard deviation for the normal distance distribution, respectively (see Skaugen and Onof 2014). Model parameters (Table 1) for each catchment and RN have been calibrated against observed runoff.

Estimating MRT from multiple RNs
When calibrating the DDD model for the different catchments and RNs, the subsurface flow celerity of the hillslopes and Rv (see Table 1) are calibrated. The mean of hillslope to RN distances, d HS , and the mean of the distances within RN, d RN , together with the respective calibrated celerities give us estimates of MRT for the hillslopes (HS) and the RN through: and respectively.

Estimating MRT from recession analysis of observed runoff
By Equations (3) and (4) we can express MRT through the mean of the recession characteristic, Λ: Λ is estimated from the sample of positive Λ s obtained from a runoff time series using Equation (2). Note again that this analysis is independent of assumptions concerning RNs and their densities. Only periods without precipitation are considered, and we filter the runoff series to remove components of the recession that are not dominated by subsurface flows, e.g. short-duration peaks indicative of overland flow early in the recession. We impose a moving average (MA) filter on the runoff time series when estimating Λ. The length of this (central) MA filter is set equal to the decorrelation time, DCT. This gives an autocorrelation of 1=e and is the mean time lag assuming an exponentially decaying autocorrelation function. The same procedure was used by Zawadski (1973) to describe characteristic precipitation patterns. In this context, we interpret the DCT as the time span over which the signal of subsurface drainage processes is captured. Note that MRT Λ includes the RTs of hillslopes and in the RN since Λ is determined from runoff measured at the outlet of the catchment.

Regionalizing MRT
The DDD model has been calibrated for 75 catchments, distributed across Norway, for flood forecasting purposes. The RNs used for the parameterization of these models are those digitized from the blue lines in 1:50 000 maps (N50 Kartdata -Geonorge Register, https://register.geonorge.no/register/versj oner/produktspesifikasjoner/kartverket/n50-kartdata). These RNs are hence derived not from a DEM, but from aerial photography. Through the calibration process, the parameters G Sh and G Sc of the distribution of the recession characteristic Λ are calibrated, which provides an estimate of Λ, and MRT can hence be estimated using (Equation 6). We want to investigate the dependencies of the estimated MRT on catchment characteristics (CCs) to obtain a better understanding of the physical underpinnings of the temporal delays of hydrological response in Norwegian catchments. In addition, such dependencies can be used to estimate MRT for ungauged basins. We investigated correlations between CCs and calibrated MRT to look for descriptor candidates to be used in a multiple regression equation for estimating MRT. Through a stepwise procedure the descriptor candidates were accepted if their inclusion in the equation was significant (p values < .001) and discarded otherwise.

MRT estimated by multiple RNs
A very clear linear relationship between the calibrated mean subsurface celerity and d HS for each RN is apparent for each catchment (Fig. 4). The coefficient of determination for the linear regressions is > 0.99 for all catchments. The mean subsurface celerity is hence proportional to d HS (i.e. the applied RN), and the slope of the linear relationship is, according to (Equation 4), the inverse of the MRT.
It is clear from Table 2 that calibrating the DDD model using different RNs gives very little variation in the resulting skill scores. The KGE values for the different RNs are almost identical and no spatial extent or density of the RN is hence favoured. In addition, no calibration parameters, other than G Sh and G Sc ; are sensitive to changing the RNs (not shown). This pattern was similar for all study catchments. Changing the RNs means changing the distances the hydraulic potential propagates, and this further implies a change in the celerities (see Equation 4) since the RT is unaltered. The equifinality with respect to distances and celerities demonstrates quite clearly how this model exercise cannot provide us with information about the "true" RN or "true" subsurface celerities, and that we need independent information on one of the variables in order to gain further insights into the underlying hydrological dynamics.

MRT estimated from the mean recession characteristic, Λ
Estimating MTT HS using (Equation 4) from multiple RNs and estimating MTT Λ from the recession characteristic Λ (Equation 6) give similar results (Table 3 and Fig. 5). In Table 3 we have also included the DCT and the catchment area. As for the MRT, the DCT seems to be more closely associated with hillslope characteristics than with catchment size (see the discussion section for more comments on the temporal characteristics associated with hillslopes and RNs). One can assume that the MRT estimated from multiple RNs is more robust since it is derived from more datapoints; on the other hand, the MRT estimated from mean recession lacks the uncertainty associated with the model calibration and simulation. What is shown here is that information concerning the temporal characteristics of catchment response can be quantified from recession analysis, independently of RN and  subsurface celerity assumptions. Parameters responsible for the simulation of catchment response in rainfall-runoff models can hence be validated, improving calibration parameter parsimony in such models (see also Kirchner (2009) on runoff response parameterized from recession).

Regionalizing MRT
A multiple regression equation (MRE) for estimating MRT from the CC of 75 calibrated Norwegian catchments has also been established using a stepwise procedure. The significant CCs include lake percentage (Lake%), areal fraction of glaciofluvial deposits classified to have good infiltration capacity (Inf ), catchment steepness (Grad), specific discharge (SpQ), catchment elevation and longitude of the catchment centroid. Lake percentage, specific discharge and elevation were the CCs with the lowest pvalues. Only the CCs with p values equal to or less than 0.001 were retained in order to ensure a robust estimate of MRT. A substantial part of the variability of MRT can be explained by selected CCs with an adjusted coefficient of determination (R 2 ) equal to 0.74. MRT MRE can be estimated by the following equation: We must note that for the 75 catchments used to establish the MRE, the MRT was estimated from a single set of calibrated values of G Sh and G Sc (see section 2.1). The individual estimates of MRT for the 75 catchments can hence be expected to be more uncertain than the MRT HS s in Table 2, which are estimated from several RNs and calibrated celerities. In Fig. 6 we have plotted MRT HS estimated by the RNs, MRT Λ estimated from recession, and MRT MRE estimated using (Equation 7) with 70% of the CCs and resampled (bootstrapped) 50 times. For the MRT MRE based on the bootstrapped samples of CCs, R 2 ranged from 0.67 to 0.84. Figure 6 shows that the regionalized MRT MRE values are, in general, close to MRT HS and MRT Λ .

Discussion
In this study we have focussed on methods for estimating the MRT for catchments in which hillslope subsurface processes dominate the runoff response. Estimating the MRT of catchments is important since it is a unique catchment characteristic and provides information on how the catchment transforms moisture input such as precipitation and snowmelt into runoff. MRT is, moreover, a catchment-scale metric that incorporates the combined effects of different landscape features in the catchment, such as lakes, wetlands, fluvial and glacial deposits, topography, etc., on runoff dynamics. The MRTs found for the 11 catchments in this study appear to be consistently estimated across all RN densities and methods of estimation for individual catchments, and for the different catchment sizes and climatic settings. The values of the estimated MRT can be interpreted as the mean time it takes for the hydraulic potential created by precipitation or snowmelt to propagate through the catchment under mean saturation conditions. The differences in estimated MRT among the catchments are very distinct; there appear to be regional patterns which are linked to regional differences in catchment characteristics, and these are discussed below. The study of runoff dynamics including RTs and recession is, for our set of catchments, basically a study of hillslope   Robinson et al. (1995) for RTs). The MRTs in Table 4 are the mean value of the MRTs estimated by (Equation 4) for hillslopes and (Equation 5) for RNs for the different configurations of RNs. Across the different climates and catchment sizes, the MRT of the hillslopes clearly dominates the total RT. The effect on RT due to transport in RNs is on the order of 1-2% of the total RT. The implication is that the MRT of hillslopes varies by an order of magnitude for the different catchments, but not as a function of catchment size. The consistency of the MRT estimates with varying RNs is, at first glance, somewhat puzzling. The linearity is striking, and we need to investigate whether this is a model artifact or an actual physical feature. Skaugen and Lawrence (2017) showed that reservoir dynamics for flow from hillslope to RN was linear if the hillslopes were subjected to uniform precipitation, constant celerity of flow towards the RN and, most importantly, an exponential distribution function of the hillslope to RN distances (the RT due to the RN itself was ignored). If the exponentially distributed distance d [m], with mean d, is replaced by d=v, where v is the the celerity [m/s] and a constant, then the distance distribution transforms into a RT distribution which is also exponential and with mean d=v [s] (Skaugen and Onof 2014). Figure 2 shows that for average conditions, we have the same RT distribution for all the RNs. This is so if the distance to RN distribution is exponential and for a mean subsurface celerity. Since the hillslope distance to RN distribution is assumed to be exponential -a one-parameter distribution -all distance distributions are proportional to each other, and the mean celerity can be seen as a proportionality constant that gives identical MRTs. To summarize, the observed linearity is a model artefact since the model assumes exponential hillslope to RN distance distributions. One the other hand, such an assumption has been confirmed by observed hillslope to RN distance distributions from blue line maps derived from aerial photography (Tucker et al. 2001, Skaugen andOnof 2014). Skaugen and Onof (2014) suggest that the recession characteristic Λ, estimated from discharge measured at the outlet, gives information primarily related to hillslope subsurface flow dynamics and RTs. They further describe how recession analysis, a catchment-scale analysis, can be used to estimate the mean recession characteristic Λ, which (through Equation 6) gives MRT Λ . It is interesting to note that the relationship of Equation (6) is independent of drainage network since Λ is a constant for all tested RNs. This is, at first glance, in contrast to the findings of Biswal and Marani (2010), who claim that there is a link between RN morphology, i.e. the number of channels, and recession. In Biswal and Marani's (2010) study, however, the extent of the active channel network varies in time and is perceived to contract during recession events. Another point here is that Biswal and Marani (2010) postulated a linear relationship between channel length and runoff. A linear relationship is not supported by observed data, and other authors suggest a power law relationship between river length and runoff (Prancevic and Kirchner 2019, Lapides et al. 2020, Durighetto andBotter 2022). The implication of contracting RNs on runoff dynamics, thinking in terms of distance distributions, is that the distances over which the hydraulic potential is propagated are longer (see also van Meerveld et al. 2019. In our case we apply a fixed RN and let the recession be determined by a decrease in subsurface flow celerity. The effect on runoff dynamics of these two approaches appears to be equivalent. One might consider combining saturation-dependent subsurface celerities and saturation-dependent RNs in a model, and this will be a topic for a future paper. An important question, which is rather difficult to answer, is which process dominates runoff dynamics: dynamic subsurface celerities or dynamic RNs. The literature gives evidence for the presence of both processes, and observations of RN dynamics are increasing, so that models which include both dynamics in RNs and subsurface celerities can now be validated (Botter et al. 2021, Lapides et al. 2020).
An important finding of this study is that optimal, or "true," RNs and subsurface celerity distributions at the catchment scale cannot be determined through this model exercise. Infinitely many combinations of mean hillslope to RN distance and mean subsurface celerity give the same MRT. This equifinality (see Beven 2001, p. 244) is consistent with the dynamic nature of RNs recently reported by many authors (Godsey and Kirchner 2014, Prancevic and Kirchner 2019, Tsegaw et al. 2019, Botter and Durighetto 2020, Durighetto and Botter 2022. Botter et al. (2021) have shown that RNs are systematically and universally characterized by significant expansion/retraction dynamics caused by hydrological and climatic variability. The equifinality with respect to RNs and mean celerities demonstrated in this study supports the notions that RNs cannot be considered static features in catchments and that the distance distributions for both hillslopes and RNs change as a function of catchment wetness. It is, nevertheless, surprising that when considering the skill scores obtained by the DDD model no particular RN density gave better scores than others. This is surprising because, clearly, one RN density accompanied by a characteristic subsurface celerity distribution must be more likely than any other. We would expect that a choice of RN density close to an unknown true mean of this distribution will result in higher KGE values. The failure of the model performance to discriminate between RNs can be caused by a coarse and imperfect model structure, a too-coarse temporal resolution or an insensitive skill score metric. Recently, researchers have initiated programmes to measure the fluctuations of RN extent (see e.g. Godsey and Kirchner 2014, Lapides et al. 2020, Botter et al. 2021, Senatore et al. 2021 and references therein). The measurements are difficult and labour intensive to carry out manually, but technological solutions for the automatic detection of the presence of water in temporary streams are continually improving (Zanetti et al. 2022). Including the RN dynamics in rainfallrunoff models may extend and improve the application of hydrological models in other disciplines such as biology, ecology, water quality studies and carbon cycling (Botter and Durighetto 2020).
The catchment characteristics explaining MRT coincide with our notion of which catchment characteristics are important in controlling hydrological response. The presence of lakes and glaciofluvial deposits for storage and infiltration will dampen the response and increase the MRT, whereas steepness and high annual precipitation values (expressed through the mean specific discharge, SpQ) decrease the MRT. Regionalizing MRT is important for predictions in ungauged basins (PUB) using the DDD model. Since MRT is independent of the RN, we can choose the density of a GIS-derived RN and estimate the mean subsurface celerity using (Equation 4). Since neither RN nor the mean subsurface celerity can be interpreted literally, it is crucial that MRT can be regionalized. This represents an advancement compared to the procedure applied in a previous PUB study using the DDD model , where the parameters G Sh and G Sc of the (gamma) distributed Λ were regionalized. The coefficient of determination for G Sh and G Sc was only 0.45 and 0.35, respectively, so regionalizing only MRT, which is inversely proportional to the product of G Sh and G Sc , is clearly a better option.

Conclusions
In this study we have shown that the mean response time (MRT) of hydrological response is a unique catchment feature. The MRT appears to be more closely associated with hillslope subsurface runoff dynamics than with catchment size. For the 11 catchments investigated, the MRT varied from 1 to 49 days.
The constituents of the MRT -the mean distance from hillslope to RN and mean subsurface celerity -display equifinality properties, and the modelling exercise was not able to suggest RNs or mean subsurface celerities that are more likely than others.
MRT was also estimated from recession analysis, with results that were very close to those obtained using the DDD model for multiple RNs. The method is independent of RNs and subsurface celerities and suggests a potential for investigating the distribution of RTs, not only its mean.
MRT was regionalized with good and robust results. We find that the presence of lakes and thick glaciofluvial deposits increase MRT, whereas high precipitation and steepness decrease MRT. For the use of the DDD model in ungauged basins, a successful regionalization of MRT represents an advancement. Since our modelling has demonstrated the equifinality between distances and celerities, we only really need an estimate of the time scale of the hydrological response, the MRT. We can then use an imprecise RN, whether from blue lines or derived from GIS, and nevertheless estimate a hypothetical mean subsurface celerity using the MRT. If the MRT is estimated with a reasonable precision, the equifinality relationship between distances and celerities ensures a reasonable simulation result.

Disclosure statement
No potential conflict of interest was reported by the authors.