Articles | Volume 25, issue 9
Research article
07 Sep 2021
Research article |  | 07 Sep 2021

Reduction of vegetation-accessible water storage capacity after deforestation affects catchment travel time distributions and increases young water fractions in a headwater catchment

Markus Hrachowitz, Michael Stockinger, Miriam Coenders-Gerrits, Ruud van der Ent, Heye Bogena, Andreas Lücke, and Christine Stumpp

Deforestation can considerably affect transpiration dynamics and magnitudes at the catchment scale and thereby alter the partitioning between drainage and evaporative water fluxes released from terrestrial hydrological systems. However, it has so far remained problematic to directly link reductions in transpiration to changes in the physical properties of the system and to quantify these changes in system properties at the catchment scale. As a consequence, it is difficult to quantify the effect of deforestation on parameters of catchment-scale hydrological models. This in turn leads to substantial uncertainties in predictions of the hydrological response after deforestation but also to a poor understanding of how deforestation affects principal descriptors of catchment-scale transport, such as travel time distributions and young water fractions. The objectives of this study in the Wüstebach experimental catchment are therefore to provide a mechanistic explanation of why changes in the partitioning of water fluxes can be observed after deforestation and how this further affects the storage and release dynamics of water. More specifically, we test the hypotheses that (1) post-deforestation changes in water storage dynamics and partitioning of water fluxes are largely a direct consequence of a reduction of the catchment-scale effective vegetation-accessible water storage capacity in the unsaturated root zone (SU, max) after deforestation and that (2) the deforestation-induced reduction of SU, max affects the shape of travel time distributions and results in shifts towards higher fractions of young water in the stream. Simultaneously modelling streamflow and stable water isotope dynamics using meaningfully adjusted model parameters both for the pre- and post-deforestation periods, respectively, a hydrological model with an integrated tracer routine based on the concept of storage-age selection functions is used to track fluxes through the system and to estimate the effects of deforestation on catchment travel time distributions and young water fractions Fyw.

It was found that deforestation led to a significant increase in streamflow accompanied by corresponding reductions of evaporative fluxes. This is reflected by an increase in the runoff ratio from CR=0.55 to 0.68 in the post-deforestation period despite similar climatic conditions. This reduction of evaporative fluxes could be linked to a reduction of the catchment-scale water storage volume in the unsaturated soil (SU, max) that is within the reach of active roots and thus accessible for vegetation transpiration from ∼258mm in the pre-deforestation period to ∼101mm in the post-deforestation period. The hydrological model, reflecting the changes in the parameter SU, max, indicated that in the post-deforestation period stream water was characterized by slightly yet statistically not significantly higher mean fractions of young water (Fyw∼0.13) than in the pre-deforestation period (Fyw∼0.12). In spite of these limited effects on the overall Fyw, changes were found for wet periods, during which post-deforestation fractions of young water increased to values Fyw∼0.37 for individual storms. Deforestation also caused a significantly increased sensitivity of young water fractions to discharge under wet conditions from dFyw/dQ=0.25 to 0.36.

Overall, this study provides quantitative evidence that deforestation resulted in changes in vegetation-accessible storage volumes SU, max and that these changes are not only responsible for changes in the partitioning between drainage and evaporation and thus the fundamental hydrological response characteristics of the Wüstebach catchment, but also for changes in catchment-scale tracer circulation dynamics. In particular for wet conditions, deforestation caused higher proportions of younger water to reach the stream, implying faster routing of stable isotopes and plausibly also solutes through the sub-surface.

1 Introduction

Plant transpiration is, globally, the largest continental water flux (Jasechko, 2018). Notwithstanding considerable uncertainties (Coenders-Gerrits, 2014), its magnitude depends on the interplay between canopy water demand and sub-surface water supply (Eagleson, 1982; Milly and Dunne, 1994; Donohue et al., 2007; Yang et al., 2016; Jaramillo et al., 2018; Mianabadi et al., 2019). The latter is regulated by water volumes that are within the reach of roots and can be taken up by plants. Many plant species across humid climate zones develop only rather shallow root systems (Schenk, 2005) that do not directly tap the groundwater (Fan et al., 2017). In regions that are dominated by such shallow-rooting vegetation, the pore volume between field capacity and permanent wilting point that is within the reach of active roots becomes a core property of many terrestrial hydrological systems (Rodriguez-Iturbe et al., 2007). This maximum vegetation-accessible water storage volume in the unsaturated root zone of soils, hereafter referred to as vegetation-accessible water storage capacity SU, max (mm), constitutes a major partitioning point of water fluxes. It regulates the temporally varying ratio between drainage, such as groundwater recharge or shallow lateral flow on the one hand and transpiration fluxes on the other hand (Savenije and Hrachowitz, 2017), which can in turn generate considerable feedback effects on downwind precipitation and drought generation (e.g. Seneviratne et al., 2013; Ellison et al., 2017; Teuling, 2018; Wang-Erlandsson et al., 2018; Wehrli et al., 2019).

Traditionally, SU, max is determined as the product of root depths or root distributions and porewater content between field capacity and permanent wilting point. Although correct in principle, this method has several weaknesses for applications at the catchment scale as much of the required data are typically not available at sufficient levels of detail. While soil maps and the associated soil water retention curves have become globally available at resolutions <1km (Arrouays et al., 2017; Hengl et al., 2017), they are characterized by considerable uncertainties. Similarly, direct and detailed observations of root systems are very scarce. They are, globally, limited to a few thousand individual plants only (e.g. Schenk and Jackson, 2002; Fan et al., 2017), and many of the observations are based on biomass extrapolations after excavating only the first metre of soil or less (Schenk and Jackson, 2003). Consequently, soil and root data largely remain inaccurate snapshots in space. As such, they are likely to be inadequate reflections of the spatial heterogeneity of soils and roots. In addition, these available data are also mostly snapshots in time and therefore disregard the adaptive behaviour of plant communities, whose compositions, and thus characteristics, at ecosystem level continuously evolve over multiple scales in space and time in response to changes in ambient conditions (e.g. Laio et al., 2006; Brunner et al., 2015; Tron et al., 2015).

There is increasing evidence that vegetation does not only actively adapt to its (changing) environment, but that it also does so in a way that allows the most efficient use of available energy and resources (e.g. Guswa, 2008; Schymanski et al., 2008). The vegetation, i.e. a collective of individual different plants within an area of interest that is present at any given moment at any given location, has survived past conditions. This in itself is a manifestation of the successful adaption of individual plants to their environment in the past. They have optimally allocated resources to balance sub- and above-surface growth to simultaneously meet water, nutrient and light requirements. This implies that these plants developed root systems that, amongst other factors, ensure continuous access to sufficient water – but not more – to bridge dry periods. An individual plant that is not adapted to meet its water and nutrient requirements through its root system as well as its light requirements through its foliage system in competition with other plants will disappear and be replaced by a better adapted plant. The root system of vegetation at ecosystem level and the associated vegetation-accessible water storage capacity SU, max are therefore at a dynamic equilibrium with and responding to the ever-changing conditions of its environment. Similarly, any type of direct human interference with vegetation, such as deforestation, has an impact on transpiration water demand, the extent and structure of active root systems and consequently on SU, max (Nijzink et al., 2016a).

For a meaningful quantification of SU, max at larger scales, such as the catchment scale, it is therefore necessary to adopt a Darwinian perspective (Harman and Troch, 2014) and to estimate effective values of SU, max, reflecting the collective and adaptive behaviour of all individual plants within a catchment. Results from many previous studies suggest, broadly speaking, three methods to do so. The first is the use of inverse approaches that treat SU, max as a model calibration parameter (Fenicia et al., 2008; Speich et al., 2018; Bouaziz et al., 2020; Knighton et al., 2020). Alternatively, the second type of method is based on optimality principles that maximize variables such as net primary production or carbon gain (Kleidon, 2004; Guswa, 2008; Hwang et al., 2009; Yang, et al., 2016; Speich et al., 2018, 2020), nitrogen uptake (McMurtrie et al., 2012) or transpiration rates (Collins and Bras, 2007; Sivandran and Bras, 2012). Lastly, SU, max and its evolution over time can be directly estimated through magnitudes of annual water deficits as determined from observed water balance data (Gentine et al., 2012; Donohue et al., 2012; Gao et al., 2014; DeBoer-Euser et al., 2016; van Oorschot et al., 2021).

For transpiration, shallow-rooting plants extract porewater of unsaturated soils that is held against gravity, i.e. between field capacity and permanent wilting point, and within the reach of roots. Significant vertical or lateral drainage only occurs at water contents above field capacity. By extracting soil water below that, transpiration therefore generates a root-zone water storage reservoir between field capacity and permanent wilting point that is characterized by a storage capacity SU, max, i.e. a maximum vegetation-accessible storage volume, and that is at any given moment filled with a specific water volume SU(t), depending on the past sequence of water inflow and release.

Storage reservoirs such as SU, max or others such as groundwater bodies are key for hydrological functioning (Sprenger et al., 2019b) as they provide a buffer against hydrological extremes such as floods and droughts. With larger storage reservoirs, the hydrological memory of a system can increase as more water can be stored and held over longer periods of time (e.g. Hrachowitz et al., 2015; Sprenger et al., 2019b). This also implies that while increased actual volumes of water stored in and thus the degree of filling of storage reservoirs, e.g. SU(t), can reduce water ages (Harman, 2015), increased sizes of storage reservoirs, e.g. SU, max, can increase water ages, both thereby controlling catchment travel time distributions (TTDs; Soulsby et al., 2010). As fundamental descriptors of hydrological functioning, TTDs describe the age structure of water held in and released from catchments (Birkel et al., 2015; Rinaldo et al., 2015), which is critical for regulating solute transport and thus nutrient and contaminant dynamics (Hrachowitz et al., 2016).

However, neither the effects of land cover change (Blöschl et al., 2019) nor the individual roles of different storage compartments in terrestrial hydrological systems are well understood (McDonnell et al., 2010; Penna et al., 2018, 2020). This is mostly a consequence of the lack of suitable observational technology to directly observe their respective volumes at larger scales. It remains therefore also unclear how deforestation affects SU, max (e.g. due to a less developed and complex rooting system for subsequent younger vegetation) and how changes in SU, max may propagate to affect both the partitioning of water fluxes as well as the age structure of water stored in and released from catchments as described by residence and travel time distributions.

For the study site of this paper, the Wüstebach experimental catchment (Germany), a previous study quantified the effects of deforestation on the partitioning of water fluxes (Wiekenkamp et al., 2016). It was found that forest removal significantly reduced evaporative fluxes. This led to more persistent higher soil moisture levels and eventually to increases in streamflow. Similarly, in the same catchment, Wiekenkamp et al. (2020) found evidence for increased post-deforestation occurrence of preferential flows, while Stockinger et al. (2019) reported minor post-deforestation reductions in travel times.

To establish a quantitative mechanistic link between these studies, we here aim to trace back and attribute the above-reported post-deforestation changes in the hydrological response of the Wüstebach to deforestation-induced changes in (sub-surface) system properties. The overall objective of this study is thus to analyse whether changes in these (sub-surface) properties can explain why deforestation affects water flux partitioning and reduces travel times in the Wüstebach in an attempt to improve our quantitative understanding of critical zone processes (Brooks et al., 2015). Specifically, we test the hypotheses that (1) post-deforestation changes in water storage dynamics and partitioning of water fluxes are largely a direct consequence of a reduction of the catchment-scale effective vegetation-accessible water storage capacity in the unsaturated root zone (SU, max) after deforestation and that (2) the deforestation-induced reduction of SU, max affects the shape of travel time distributions and results in shifts towards higher fractions of young water in the stream.

Figure 1Map of the Wüstebach study catchment showing the spatial distribution of soil types. The riparian zone is defined by the parts of the catchment covered by Gleysols, Planosols and Halfbogs. The red line indicates the outline of the deforested part of the catchment, as can also be seen in the aerial images (© Google Earth, Maxar Technologies 2020) from 2013 and 2016.

Figure 2(a) Time series of observed weekly precipitation P; (b) daily cumulative evaporative fluxes for the pre- and post-deforestation periods, where the dark brown line indicates potential evaporation EP and the orange lines and the yellow shaded areas show the actual evaporation EA modelled using the best fit parameter sets and the associated 5th/95th percentiles of all feasible solutions of the pre- and post-deforestation periods, respectively. The dashed red line indicates the modelled EA in the post-deforestation period using the best fit pre-deforestation parameter set; (c) observed (dark blue line) and modelled daily streamflow Q; light blue line indicates the best fit model and the shaded area the 5th/95th percentile of all feasible solutions for the pre- and post-deforestation periods, respectively. The dashed red line indicates the modelled streamflow in the post-deforestation period using the best fit pre-deforestation parameter set; (d) zoom-in to the observed and modelled streamflow for the October 2012–October 2014 period. The grey shaded area indicates the deforestation period.


2 Study site

The experimental Wüstebach headwater catchment (0.39 km2; Fig. 1a) is part of the Lower Rhine/Eifel Observatory of the Terrestrial Environmental Observatories network (TERENO; Bogena et al., 2018) located in the Eifel National Park in Germany (5030'16” N, 0620'00” E). The catchment is characterized by a humid, temperate climate with warm summers, mild winters and a mean annual temperature of around 7 C (Zacharias et al., 2011). Mean annual precipitation is about 1200 mm yr−1 and mean annual runoff about 700 mm yr−1 (Fig. 2). Although most of the precipitation occurs in the winter months, the fraction that falls as snow is typically less than 10 % of the annual precipitation, and snow cover is present for no more than 3–4 weeks per year.

The catchment is drained by a perennial second-order stream and extends from 595 to 630 ma.s.l. The landscape is characterized by the gentle slopes of the surrounding hills and a flatter riparian area close to the stream, covering approximately 10 % of the catchment (Fig. 1a). The underlying bedrock is largely Devonian shales with sandstone inclusions (Richter, 2008) covered by periglacial layers (Borchardt, 2012). While Cambisols dominate the hillslopes, Gleysols and Histosols characterize much of the riparian area (Bogena et al., 2015). The average soil depth in the catchment reaches about 1.6 m with a maximum of 2 m (Graf et al., 2014). In 1946, after the Second World War, the catchment was homogeneously and completely afforested (Fig. 1) with Sitka spruce (Picea sitchensis) and Norway spruce (Picea abies; Etmann, 2009). The maximum observed rooting depth of these spruce trees in the catchment is 50 cm, and no roots were observed below this depth. In the course of the development of the area into a national park, approximately 21 % of the catchment, including the entire riparian zone, was deforested in September 2013 and has been kept largely vegetation free since (Wiekenkamp et al., 2016; Fig. 1).

3 Data

3.1 Hydro-meteorological data

Daily hydro-meteorological data were available for the period 1 October 2009–30 September 2016 (Fig. 2). Precipitation P (mm d−1) and mean daily temperature T (C) were available from the Monschau–Kalterherberg meteorological station operated by the German Weather Service (Deutscher Wetterdienst DWD station 3339), located 9 km north-west of the Wüstebach catchment. The precipitation data were corrected for evaporation and wind drift losses according to Richter (1995) and as described in detail by Graf et al. (2014). Stream discharge Q (mm d−1) at the outlet of the Wüstebach was observed with a V-notch weir for low-flow measurements and a Parshall flume for medium to high flows (Bogena et al., 2015). Daily potential evaporation EP (mm d−1) was estimated using the Penman–Monteith equation. Daily depth-weighted average soil water content for the study period was estimated from a network of soil moisture sensors placed at 5, 20 and 50 cm depths at >100 locations across the study catchment as described by Graf et al. (2014) and Bogena et al. (2015). In addition, throughfall rates PE (mm d−1) were measured at one continuously forested location in the study catchment (Fig. 1) with an array of samplers as described in detail by Stockinger et al. (2015) over irregular intervals over the period 1 October 2012–30 September 2016.

Figure 3(a) Observed volume-weighted monthly δ18O signals in precipitation (grey dots; size of dots indicates the precipitation volume) and streamflow (green dots) as well as the best fit modelled δ18O signal in the stream (green line) and the 5th/95th percentile of all feasible solutions from pre- and post-deforestation calibration (green shaded area); (b) zoom-in of observed and modelled δ18O signals in the stream for the October 2012–October 2014 period.


Figure 4Model structure used in this study. The light blue boxes indicate the hydrologically active individual storage volumes in the hillslope and riparian zones, respectively. The darker blue box SS,p indicates a hydrologically passive, i.e. dSS,p/dt=0, mixing volume. The blue lines indicate liquid water fluxes, and the green lines indicate vapour fluxes. Model parameters are shown in red adjacent to the model component they are associated with. All symbols are defined in Table 1.


Table 1Water balance, state and flux equations used in the hydrological model. Symbols shown in bold are model parameters. Subscripts H and R indicate hillslope and riparian zone, respectively. Model variables: P is total precipitation (mm d−1), PS is solid precipitation (snow) (mm d−1), PM is snowmelt (mm d−1), PR is rain (mm d−1), PE is throughfall (mm d−1), EP is potential evaporation (mm d−1), EI is interception evaporation (mm d−1), RF is preferential recharge (mm d−1), RS is slow recharge (mm d−1), ET is transpiration (mm d−1), QS is flow from the slow-responding reservoir (mm d−1), QR is flow from the fast-responding riparian reservoir (mm d−1), Q is the total flow (mm d−1) and EA is the total actual evaporation (mm d−1). Model parameters: TT is the threshold temperature (C), FM is a melt factor (mmd-1C−1), Imax is the interception capacity (mm), SU, max is the root-zone storage capacity (mm), γ is a shape factor (–), RS,max is the maximum percolation rate (mm d−1), Lp is a transpiration water stress factor (–), fQS is a factor determining the fraction of groundwater flow that is upwelling into the riparian zone (–), kS is the storage coefficient of the slow-responding reservoir (d−1), kR is the storage coefficient for the fast-responding riparian reservoir (d−1) and f is the areal fraction of the riparian zone (–).

Download Print Version | Download XLSX

3.2 Stable isotope data

Regular weekly δ18O data from bulk precipitation samples collected in a cooled wet deposition gauge at the meteorological station Schleiden–Schöneseiffen (Meteomedia station) 3 km north-east of the catchment were available for the period 1 October 2010–24 September 2012. After that, precipitation was sampled at half-daily intervals until 30 September 2016 using an automatic, cooled sampler (Eigenbrodt GmbH, Germany). The half-daily samples were precipitation volume weighed to daily sampling intervals (Stockinger et al., 2016, 2017). Weekly stream water grab samples for stable water isotope analysis were taken at the outlet of the Wüstebach catchment in the 1 October 2010–30 September 2016 period (Fig. 3a; Bogena et al., 2020).

Isotope analysis was carried out using laser-based cavity ring-down spectrometers (L2120-i/L2130-i, Picarro Inc.). Internal standards calibrated against VSMOW, Greenland Ice Sheet Precipitation (GISP) and Standard Light Antarctic Precipitation (SLAP2) were used for calibration and to ensure long-term stability of analyses (Brand et al., 2014). The long-term precision of the analytical system was ≤0.1% for δ18O.

4 Methods

To quantify effects of deforestation on SU, max and, due to the role of SU, max as a mixing volume also in the age structure of water as described by TTDs and the associated young water fractions Fyw, the following stepwise experiment was designed. (1) Quantify changes in the partitioning of annual water fluxes between the pre- and post-deforestation periods based on observed water balance data. (2) Estimate the effect of these changes on the magnitudes of pre- and post-deforestation SU, max, respectively, using the same data. (3) Calibrate a hydrological model to simultaneously reproduce streamflow and stream δ18O dynamics for the pre-deforestation period. (4) Use the calibrated parameter sets to run the model in the post-deforestation period and evaluate the model's post-deforestation performance without further calibration. (5) Re-calibrate the model for the post-deforestation period and evaluate whether changes in calibrated SU, max (and other parameters) are plausible and reflect changes in SU, max directly estimated from water balance data in step (2). Finally, (6) use the calibrated pre- and post-deforestation parameter sets, respectively, to track modelled water fluxes through the system and quantify changes in TTDs and Fyw between the pre- and post-deforestation periods.

4.1 Water balance-based estimation of SU, max

To survive, plants need continuous access to water to satisfy canopy water demand. The root systems of vegetation are therefore adapted to provide access to water volumes that correspond to annual water deficits that result from the combination of (1) the phase lag between and (2) the difference in the respective magnitudes of seasonal precipitation and solar radiation signals (Donohue et al., 2012; Gentine et al., 2012; Gao et al., 2014). On a daily basis, these water deficits SD,j(t) can be estimated as the cumulative sum of daily throughfall PE (mm d−1) minus transpiration ET (mm d−1). The maximum deficit SD,j for a specific year j is then equivalent to the soil water volume that was accessible to and actually accessed by vegetation through its root system for transpiration during the dry season over that period when ET exceeded PE (deBoer-Euser et al., 2016; Nijzink et al., 2016a; van Oorschot et al., 2021):

(1)SD,jt=t0tPEt-ETtdt,if SD,jt0,0,if SD,jt>0,(2)SD,j=maxSD,jt,

where t is the time step (d) and t0 is the last preceding time step for which the storage deficit SD,j(t)=0. As an approximation, Eq. (1) implies that if SD,j(t)=0, the water content in the root-accessible pore space at day t is at field capacity and cannot hold additional water. If water supply then exceeds canopy water demand on that day, i.e. PE(t)ET(t)>0, this water surplus is drained from the root zone, e.g. to recharge groundwater or directly to the stream, and cannot be used for transpiration.

Daily throughfall PE, i.e. precipitation that actually reaches the soil, was estimated on the basis of the water balance of a canopy interception storage (Nijzink et al., 2016a):

(3) d S I t d t = P t - E I t - P E t ,

where EI (mm d−1) is daily interception evaporation and SI (mm) the canopy interception storage. For each time step, EI can then be computed as

(4) E I t = E P t , if  E p t d t < S I t . S I t d t , if  E P t d t S I t .

This then further allows us to estimate PE according to

(5) P E t = 0 , if  S I t < I max , S I t - I max d t , if  S I t I max ,

where Imax (mm) is the canopy interception capacity. In the absence of more detailed information PE was estimated with a range of different interception capacities, i.e. Imax=0, 1, 2, 3, and 4 mm, in a sensitivity analysis approach.

Note that the catchment average PE after deforestation was estimated as the area-weighted mean of PE in the deforested area (21 % of the catchment area) computed with an assumed Imax=0 mm and PE from the remaining area computed based on the above range of Imax between 0 and 4 mm. In a next step, assuming negligible groundwater imports or exports (cf. Bouaziz et al., 2018), data errors and storage changes, long-term mean transpiration ET¯ was estimated according to the water balance:

(6) E T = P E - Q ,

where PE (mm d−1) is the long-term mean throughfall and Q (mm d−1) is the long-term mean observed stream discharge. Daily transpiration ET (mm d−1) for use in Eq. (1) is then estimated by scaling the long-term mean transpiration to the signal of daily potential evaporation to approximate the seasonal fluctuation of energy input (Bouaziz et al., 2020):

(7) E T t = E P t - E I t E T E P - E I .

A range of previous studies provided evidence that mature forests develop root systems that allow access to sufficiently large porewater storage volumes SU, max to bridge droughts with return periods TR∼40 years (Gao et al., 2014; deBoer-Euser et al., 2016; Nijzink et al., 2016a; Wang-Erlandsson et al., 2016). The maximum annual water deficits SD,j (Eq. 2) for all j years in the pre-deforestation study period were therefore used to fit a Gumbel extreme value distribution (Gumbel, 1941). This subsequently allowed the estimation of a water deficit with a 40-year return period, which is for this study defined as vegetation-accessible water storage SU, max so that SU, max=SD,40yr.

Note that due to the limited length of the data series, the SU, max estimates are rather uncertain and need to be understood as merely indicative approximations. This is in particular true for the post-deforestation period, where attempts to explicitly link SU, max to a specific return period are subject to additional uncertainty: as the catchment was not reforested and natural recovery of vegetation is negligible (see aerial images in Fig. 1), it is not implausible to assume that the development of the root system after the disturbance is far from equilibrium and likely to be actively evolving over time. Also note that although ET is, for brevity, referred to as transpiration throughout this paper, it also contains soil evaporation. However, no explicit and quantitative distinction could be made between these two fluxes with the available data. A further critical assumption of the above method required that roots do not tap the groundwater and that water for transpiration is exclusively extracted from the unsaturated soil. In contrast to other landscapes (Fan et al., 2017; Roebroek et al., 2020), it is likely that this assumption largely holds in the Wüstebach as throughout the catchment the groundwater levels, also in the riparian zone, remain largely below a depth of 50 cm during the relatively dry growing season (Bogena et al., 2015), when storage deficits SD typically accumulate (May to October) and no roots have so far been observed for the dominant picea species below that depth in the Wüstebach catchment. This is also broadly consistent with the results of Evaristo and McDonnell (2017), who show rather limited groundwater use by picea species.

4.2 Model architecture

A semi-distributed, process-based catchment model, iteratively customized and tested within the previously developed DYNAMITE modular modelling framework (Hrachowitz et al., 2014; Fovet et al., 2015), was adapted with additional, hydrologically passive storage volumes to allow for simultaneous representation of water fluxes and tracer transport (Hrachowitz et al., 2013) based on the general concept of storage-age selection functions (SAS; Rinaldo et al., 2015). This model type was chosen over simpler, more data-based methods (e.g. McGuire and McDonnell, 2006; Kirchner, 2016) as it did not only allow a simultaneous representation of water and tracer fluxes, but also allowed attribution of an observed pattern to specific process hypotheses and the associated model parameters that represent (sub-surface) system properties, thereby providing potential quantitative mechanistic explanations of why deforestation affects the hydrology in the Wüstebach. As an intermediate model type between purely data-driven (e.g. Kirchner, 2016) and spatially explicit physically based models (e.g. Maxwell et al., 2016), it requires assumptions about underlying processes and effective parameters and does not allow a detailed spatial analysis. Yet this model type provides the possibility of testing these process hypotheses at the scale of the semi-distributed model units, thereby integrating and accounting for the natural heterogeneity of system properties across the model domain (Hrachowitz and Clark, 2017).

4.2.1 Hydrological model

The model domain of the Wüstebach catchment was spatially discretized into two functionally distinct response units, i.e. hillslopes and riparian areas. These are represented in the model as two parallel suites of storage components, linked by a common groundwater body as shown in Fig. 4 (e.g. Euser et al., 2015; Nijzink et al., 2016b). According to elevation data and distribution of soil types (Fig. 1), 90 % of the catchment area was classified as hillslope and the remaining 10 % as riparian area. Below a threshold temperature TT (C), precipitation P (mm d−1) accumulates as snow PS (mm d−1) in SSnow (mm). Above that temperature precipitation is falling as rain PR (mm d−1) and snowmelt PM (mm d−1) is released from SSnow according to a melt factor FM (mmd-1C-1) using a simple degree-day method (e.g. Arsenault et al., 2015; Ala-aho et al., 2017; Gao et al., 2017). The total liquid water input PR+PM (mm d−1) entering the hillslope is routed through the canopy interception storage SI,H (mm). Water that is not evaporated as EI,H (mm d−1) enters the unsaturated root zone SU,H (mm), whose storage capacity is defined by the calibration parameter SU,max,H (mm). Water can be released from SU,H as combined root-zone transpiration and soil evaporation flux ET,H (mm d−1) or eventually recharge the groundwater SS,a (mm) over a fast, preferential recharge pathway as RF,H (mm d−1) and a slower percolation flux RS,H (mm d−1). Similarly, water entering the riparian zone, i.e. PR+PM (mm d−1), is routed through SI,R (mm). Excess water PE,R (mm d−1) that is not evaporated infiltrates into the unsaturated root-zone SU,R (mm), defined by calibration parameter SU,max,R (mm). In addition, a fraction of the upwelling groundwater RS,R (mm d−1) replenishes SU,R and, thus, in addition to precipitation, sustains soil moisture levels in the riparian zone (e.g. Hulsman et al., 2021a), while the remainder QS (mm d−1) drains directly into the stream. While water stored in SU,R is available for transpiration (and soil evaporation) ET,R (mm d−1), water that cannot be held is released as RF,R (mm d−1) to a fast-responding reservoir SF,R (mm), from where it reaches the stream as QR (mm d−1). The relevant model equations can be found in Table 1.

4.2.2 Tracer transport model

The δ18O composition of water fluxes and storages was tracked through the model using the SAS approach (Rinaldo et al., 2015), which allows a catchment-scale description of conservative transport based on time-variant travel time distributions. The method builds on the fact that a water volume S (mm) stored in any storage component can, at any moment t (d), consist of parcels of water of different ages T (d). The composition of ages in the stored volume at t depends on the history of water inflows and outflows. Consequently, it evolves over time as new inputs enter into and outflows are released from the storage component, whereby each inflow I (mm d−1) and outflow volume O (mm d−1) can have a different age composition. A convenient way to implement the SAS approach is the use of age-ranked storage ST(T,t) (mm), which represents “at any time t the cumulative volumes of water in a storage component as ranked by their age T” (Benettin et al., 2017). Similarly, decomposing each inflow and outflow of a storage component into their respective cumulative, age-ranked volumes IT(T,t) and OT(T,t) (mm d−1), respectively, then allows us to update the age-ranked storage ST(T,t) at each time step according to the general water age balance (Botter et al., 2011; van der Velde et al., 2012; Benettin et al., 2015a, 2017; Harman, 2015):


where the term ST/T represents the aging of water in storage. Reflecting the slightly more abstract approach by Rodriguez and Klaus (2019) and similar to previous studies based on the functionally equivalent mixing coefficient approach (e.g. Fenicia et al., 2010; McMillan et al., 2012; Birkel and Soulsby, 2016; Hrachowitz et al., 2015), the water age balance is here individually formulated for each storage reservoir j (e.g. SI,H, SU,H), which each can have varying numbers N and M of inflows I (e.g. PR, PM, RS,H) and outflows O (e.g. PM, RS,H, QS), respectively (see Fig. 4). It is assumed that the entire volume of a precipitation signal P(t) entering the system at t has an age T of zero so that the associated IT,P,j(T,t)=PT(T,t)=P(t) for all T. As all other inflows to any following storage component in the system are outflows of storage components prior in the sequence (see Fig. 4), the corresponding IT,n,j(T,t) entering a storage component are identical to the OT,m,j(T,t) released from the storage component above.

Each age-ranked outflow OT,m,j(T,t) of a specific storage component j depends on the outflow volume Om,j(t) along this outflow pathway and the cumulative age distribution Po,m,j(T,t) of that outflow:

(9) O T , m , j T , t = O m , j t P O , m , j T , t .

The outflow volume Om,j(t) is estimated via the hydrological model (see Sect. 4.2.1; Fig. 4) and thus assumed to be known. In contrast, the cumulative age distribution Po,m,j(T,t) can in general not be directly parameterized, as it depends on the temporally varying age distribution of water in the storage component j represented by ST,j(T,t) and thus on the history of past inflows and outflows (Botter et al., 2011; Harman, 2015). Instead, it is possible to define a SAS function ωo,m,j (or Ωo,m,j in its cumulative form) for each outflow m from each storage component j that describes how outflow is sampled (or selected) from the temporally varying water volumes of different ages present in the age-ranked storage ST,j(T,t) at any time t:

(10) P O , m , j T , t = Ω O , m , j S T , j T , t , t .

From the cumulative age distribution Po,m,j(T,t) the associated probability density function, which represents the outflow age distribution po,m,j(T,t), frequently also referred to as backward travel time distribution of that outflow (TTD; e.g. Benettin et al., 2015a; Wilusz et al., 2017), can be obtained according to

(11) p O , m , j T , t = ϖ O , m , j S T , j T , t , t S T , j T .

Note that conservation of mass requires that any SAS function ωO,m,j integrates to the total storage volume Sj(t) present in j at any time t. To avoid the resulting need for rescaling ωO,m,j at each time step, it is helpful to normalize the age-ranked storage to ST,norm,j(T,t)=ST,j(T,t)/Sj(t) so that it remains bounded to the interval [0,1] and defines a residence time distribution (RTD).

For this study, beta distributions, which are conveniently bound between the limits [0, 1] and defined by two shape parameters α and β, were used as SAS functions ωo,m,j to sample water of different ages for outflows from storage components. The parameters β were fixed at a value of 1 for all SAS functions ωo,m,j used here. However, there is substantial evidence for preferential flow through macropores in the shallow sub-surface (e.g. Weiler and Naef, 2003; Zehe et al., 2006, 2007; Weiler and McDonnell, 2007; Beven, 2010; Beven and Germann, 2013; Klaus et al., 2013; Angermann et al., 2017; Loritz et al., 2017). Such preferential flow can, with increasing wetness, increasingly bypass water volumes stored in small pores with little exchange (Sprenger et al., 2016, 2018, 2019a; Cain et al., 2019; Evaristo et al., 2019; Knighton et al., 2019). This then leads to an increasing preferential release of younger water as the system becomes wetter (Brooks et al., 2010). To mimic this, the shape parameters α of the preferential fluxes RF,H and RF,R released from the two unsaturated root-zone storage components Sj=SU,H and SU,R (Fig. 4) were allowed to vary as a function of the water volumes stored in SU,H and SU,R, respectively (Hrachowitz et al., 2013; van der Velde et al., 2015):

(12) α m , j t = 1 - S j t S U,max , j 1 - α 0 ,

where α0 is a calibration parameter representing a lower bound so that αm,j(t) can vary between α0 and 1. A value of αm,j=1 indicates complete mixing in dry conditions. Any value below that entails incomplete mixing and thus increases the preference towards releasing younger water in wet conditions (Benettin et al., 2017). Although there is evidence for the presence of preferential flow in other components of the system, such as in the groundwater (e.g. Berkowitz and Zehe, 2020), initial model testing suggested that the inclusion of the additional calibration parameters is not warranted by the available data. For simplicity and following the principle of model parsimony, we assumed complete mixing for all other outflows from all other storage components (Fig. 4; cf. Fenicia et al., 2010; Kuppel et al., 2018a; Rodriguez et al., 2018). Parameter α was therefore fixed to a value of 1 for these SAS functions.

The δ18O precipitation input signals are damped to the level of fluctuation observed in the stream by sub-surface storage volumes that remain to some extent hydrologically passive (e.g. Birkel et al., 2011b). While the hydrologically active storage volumes are represented by the individual storage components of the model (Fig. 4; Eqs. 8–14), an additional hydrologically passive storage volume SS,p (mm) was added as a calibration parameter to the active groundwater storage SS,a (Zuber, 1986; Hrachowitz et al., 2015, 2016), so that SS,tot=SS,a+SS,p (Fig. 4). While dSS,p/dt=0, the age-ranked groundwater storage was computed as ST,S,tot and the outflows from the groundwater component consequently sampled from the entire storage volume SS,tot, thereby representing the combined contributions from SS,a and SS,p to the age structure of the outflow QS according to Eq. (39). Note that the effects of the hydrologically passive water volume stored in the unsaturated soil below the wilting point are assumed to be negligible due to the small size of that storage volume and the low diffusive exchange rates with the hydrologically active storage volume in the unsaturated zone.

Each individual volume with a different age in IT,n,j(T,t) and, as a consequence, also in ST,j(T,t) is also characterized by a different tracer concentration CI,n,j(IT,n,j(T,t),t) and CS,j(ST,j(T,t),t), respectively. For a conservative tracer such as δ18O that is not significantly affected by decay, evapoconcentration, retention or any other biogeochemical transformation (e.g. Bertuzzo et al., 2013; Benettin et al., 2015b; Hrachowitz et al., 2015), the concentration CO,m,j(t) in any outflow at any time t can then be obtained from


Due to data availability, age tracking was here limited to 4 years in the pre-deforestation period and 3 years in the post-deforestation period. For ages beyond that it can only be said that water is older than these 4 and 3 years, respectively. The TTDs reported hereafter are thus truncated at these ages. The model generates TTDs for all fluxes and storage components (Fig. 4) for each time step. As a summary metric, we will here use the fraction of young water Fyw as a robust descriptor of the left tail of TTDs. Following the definition of Kirchner (2016), Fyw is here the fraction of water that is younger than 3 months, which can be extracted directly from any TTD generated by the model. Note that we here only analyse water ages in streamflow as these are the only ones that are directly constrained by available data, while for all other model components, such as transpiration ET, such direct data support was not available, and the resulting age estimates may thus be characterized by considerable additional uncertainty.

Table 2Parameter prior distributions and 5th/95th percentiles of the posterior distributions. Note that * parameter f, characterizing the areal proportion of the riparian zone, was fixed according to soil and elevation data and ** the interception capacity Imax was assumed to be identical on the hillslopes and the riparian zone in the pre-deforestation period.

Download Print Version | Download XLSX

Table 3Signatures of flow and δ18O and the associated performance metrics used for model calibration and evaluation. The performance metrics used include the Nash–Sutcliffe efficiency (ENS), the volume error (EV) and the relative error (ER).


Download Print Version | Download XLSX

4.3 Model calibration and post-calibration evaluation

The model was run with a daily time step and has a total of 14 free calibration parameters, which were calibrated for the model to simultaneously reproduce flow and δ18O dynamics in the stream. The uniform prior parameter distributions (Table 2) were sampled using a Monte Carlo approach with 3×106 realizations. To limit equifinality (Beven, 2006) and to ensure robust posterior parameter distributions for a meaningful process representation (e.g. Kuppel et al., 2018b), an extensive multi-objective calibration strategy was applied. Briefly, this was done using a total of 14 performance metrics that describe the model's skill in reproducing different signatures associated with streamflow (EQ) and δ18O dynamics (Eδ18O) as shown in Table 3.

To be accepted as feasible, solutions had to exceed a threshold value of 0.5 for all performance metrics, with the exception of ENS,δ18O, for which a threshold of 0.2 was used. To further constrain the model, we only accepted solutions that could reproduce the dynamics in observed soil moisture as well as the average observed magnitudes of canopy throughfall. To do so we used a simplified limits-of-acceptability approach (e.g. Coxon et al., 2014) with a rectangular step function so that all solutions that fall within the limits of the step function receive a weight of 1 while all others are assigned a weight of 0 and are thus rejected (Bouaziz et al., 2021). More specifically, we rejected solutions whose modelled normalized relative soil moisture fell outside the acceptable limits, here defined as ±0.15 of the observed relative soil moisture, in more than 75 % of the time steps in the calibration periods. Similarly, we rejected solutions for which the modelled mean ratio PE/P in the continuously forested part of the catchment was outside ±0.15 of the observed mean ratio PE/P=0.71. This strategy was chosen instead of directly calibrating the time series of associated model variables SU (Eq. 20) and PE (Eq. 18) to explicitly account for commensurability errors between the point scale and the scale of the model application (Bouaziz et al., 2021). Subsequently, the 14 metrics of the solutions retained as feasible were combined into two equally weighted classes, describing streamflow (Q) and tracer (δ18O) dynamics, respectively. This then allowed us to obtain solutions with balanced overall model performances using the mean Euclidean distance DE (–) from the “perfect” model (i.e. DE=1; Hrachowitz et al., 2014; Hulsman et al., 2020):


where N=12 is the number of different performance metrics describing streamflow and M=2 the number of different performance metrics for δ18O. To construct the posterior parameter distributions and the corresponding model uncertainty intervals, the retained parameter sets where then weighted according to a likelihood measure L=DEp (cf. Freer et al., 1996), where the exponent p was set to a value of 10 to emphasize models with good overall calibration performance.

In a first step, the model was calibrated for the pre-deforestation period 1 October 2009–31 August 2013. Note that due to a lack of regular and weekly δ18O precipitation data before 1 October 2010, the performance metric Eδ18O describing the δ18O dynamics was computed from that date onwards only. The feasible parameter sets were then used to test the model without further calibration in the post-deforestation period. In a second step, the model was re-calibrated for the 1 September 2013–30 September 2016 post-deforestation period and the changes in the resulting model performance and posterior distributions compared to those from the pre-deforestation calibration. The estimation of the effects of deforestation on TTDs is based on model parameter sets obtained from calibration in the pre-deforestation and post-deforestation periods, respectively.

Figure 5(a) Positions of the individual years of the study period in the Budyko framework. The x axis shows the aridity index IA=EP/P; the y axis indicates the evaporative ratio EA/P and the runoff ratio CR=1-EA/P. Pre-deforestation years are shown with blueish shades, post-deforestation years with greenish shades. The bold black lines indicate the energy and water limits, respectively. The dashed grey line is the theoretical–analytical Turc–Mezentsev relationship (Turc, 1954; Mezentsev, 1955). (b) The range of time series of storage deficits as computed according to Eq. (2), using values of Imax from 0 to 4 mm. The maximum annual storage deficits SD,j are indicated by the arrows. The grey shaded area indicates the deforestation period. (c) Estimation of SU, max as the storage deficit associated with a 40-year return period SD,40 yr using the Gumbel extreme value distribution for the pre-deforestation period. The blueish dots indicate the range of maximum annual storage deficits SD,j for each year in the 4-year pre-deforestation period. The dark-grey shaded area indicates the envelope of least-square fits for the individual values of Imax. The light-grey shaded area indicates the envelope of the 5th/95th confidence intervals. The red line shows the plausible range for SU, max.


Figure 6(a) Model performance metrics for all variables and signatures. DE is the Euclidean distance to the perfect model. It combines all other performance metrics (Table 3) into one number (Eq. 42). All performance metrics are formulated in a way that a value of 1 indicates a perfect fit. The boxplots summarize the performances of all parameter sets retained as feasible. The circle symbols indicate the performance of the best-performing model in terms of DE. The dark-red shades indicate pre-deforestation model performance based on calibration in the pre-deforestation period. Orange shades indicate post-deforestation performance using the pre-deforestation parameter sets without further re-calibration. Yellow shades show the post-deforestation performance after model re-calibration in the post-deforestation period. (b, c) show flow duration curves, (d, e) show the peak distributions and (f, g) show the autocorrelation functions for the pre- (red) and post-deforestation periods (orange and yellow), respectively. The black lines indicate the observed values, the dashed lines indicate the best fits and the shaded areas indicate the 5th/95th uncertainty interval of all solutions retained as feasible. The dark-red shades indicate pre-deforestation model results based on calibration in the pre-deforestation period. Orange shades indicate post-deforestation model results using the pre-deforestation parameter sets without further re-calibration. Yellow shades show the post-deforestation model results after model re-calibration in the post-deforestation period.


5 Results

5.1 Observed deforestation effects on the hydrological system

Initial analysis of water balance data suggests that the hydro-meteorological conditions as expressed by the aridity index IA=EP/P do not show significant differences between the pre-deforestation (IA=0.50±0.02) and post-deforestation periods (IA=0.51±0.03), respectively (Fig. 5a). However, and in spite of these comparable climatic conditions, the results show a shift in the partitioning of water fluxes between runoff Q and actual evaporation EA (note that EA=EI+ET). While the fraction of precipitation that was released into the atmosphere as vapour was reduced (EA/P; Fig. 5a), the mean runoff ratio (CR=1-EA/P) increased correspondingly from CR=0.55±0.04 to CR=0.68±0.03 after deforestation of 21 % of the catchment with p=0.049 based on a Wilcoxon rank sum test. In absolute terms this entails that, notwithstanding rather stable mean annual precipitation P=1269±24mm yr−1 and potential evaporation EP=632±9mm yr−1 over the entire study period, the annual actual evaporation EA decreased from 576±11 to 401±6mm yr−1, whereas annual runoff Q increased by ∼25% from 694±47 to 870±63mm yr−1.

In spite of similar climatic conditions, the above is reflected in a significantly higher (p=0.047) mean annual maximum storage deficit in the pre-deforestation period than in the post-deforestation period. In the pre-deforestation period values between 105±23mm for Imax=0 mm and 95±21 mm for Imax=4 mm, respectively, were found (Fig. 5b), whereas in the post-deforestation period the mean storage deficit only reached between 49±10 and 33±7mm for the same values of Imax (Fig. 5b). Note that in both periods, SD,j is relatively insensitive to the magnitude of Imax (cf. Gerrits et al., 2009). From the above maximum annual storage deficits SD,j, the corresponding catchment-scale vegetation-accessible water storage capacity, assuming vegetation adaptation to dry conditions with 40-year return periods (see Sect. 4.1), was estimated at values of SU, max=258±125mm for the pre-deforestation (R2=0.91, p=0.04; Fig. 5c) and SU, max=101±149mm for the post-deforestation period (R2=0.83, p=0.27; not shown).

Figure 7Observed mean PE/P (dashed line), the range around observed mean PE/P defined as acceptable (grey shaded area), the distribution of modelled mean PE/P from all solutions that satisfy the behavioural thresholds for all performance metrics (Table 3) as well as the mean PE/ of the best solution in terms of DE (orange symbol) for (a) the 2009–2013 pre-deforestation period and (b) the 2013–2016 post-deforestation period. Note that only modelled solutions (yellow) that fall into the acceptable observed range (grey shaded) are kept as feasible. The fractions of time steps in the pre-deforestation (c) and post-deforestation (d) periods in which the modelled relative soil moisture SU,rel falls within the pre-defined acceptable range around the observed relative catchment-average soil moisture. The blue symbols indicate the best solution in terms of DE, and the distributions indicate the set of solutions that satisfy the behavioural thresholds for all performance metrics (Table 3). The grey shaded areas indicate the region of acceptable solutions, i.e. solutions that fall at least 75 % of the time steps into the acceptable interval. Note that only modelled solutions (light blue) that fall into the acceptable observed range (grey shaded) are kept as feasible. Pre-deforestation (e) and post-deforestation (f) time series of the acceptable range around the observed normalized, relative soil moisture (light-grey shade) and range of modelled normalized relative soil moisture for all solutions that satisfy all performance metrics (“unconstrained”; light blue) and for the set of feasible solutions that satisfy both PE/P and soil moisture constraints as shown in (a–d) (“constrained”; blue). The dark-blue line indicates the modelled normalized, relative soil moisture of the best solution in terms of DE.


Figure 8Posterior distributions of selected parameters shown as empirical cumulative distribution functions (lines) and the associated relative frequency distributions (bars). Red shades indicate calibration in the pre-deforestation period. Yellow shades indicate post-deforestation calibration. The dots indicate the parameter values associated with the respective best fit models.


5.2 Modelled deforestation effects on the hydrological system

5.2.1 Model calibration for the pre-deforestation period

The model parameter sets retained as feasible after calibration in the 2009–2013 pre-deforestation period reproduce the general features of the hydrograph in that period rather well (Fig. 2c, d), similar to a previous modelling study (Cornelissen et al., 2014). This is true for both the timing and magnitudes of high flows, with an associated Nash–Sutcliffe efficiency ENS,Q=0.83 for the best-performing model in terms of DE (Fig. 6a) but also for low flows (ENS,log(Q)=0.70), with the exception of some overestimation in summer 2011. The modelled runoff ratio comes with CR=0.54 (5th/95th interquantile range (IQR): 0.52–0.58) very close to the observed runoff ratio of CR=0.55 (ER,CR=0.98). In addition, the model could also simultaneously mimic most other observed flow signatures reasonably well (Fig. 6a), in particular the flow duration curve (ENS,FDC=0.79; Fig. 6b), the peak distribution (ENS,PD=0.85; Fig. 6d) and the autocorrelation function (ENS,AC=0.98; Fig. 6f). The limits-of-acceptability constraints for PE/P allowed the identification and removal of a few additional parameter sets (∼5%) that likely overestimate throughfall PE (Fig. 7a). The soil moisture constraint was more effective as it allowed us to reject a considerable additional proportion of solutions (>90%) that did not sufficiently well match the observed soil moisture dynamics according to the pre-defined limits of acceptability (Fig. 7c). With the parameter sets eventually retained as feasible, the modelled temporal dynamics of relative soil moisture broadly reflect the observed ones (Fig. 7e). Similarly, the model captures the substantial attenuation of the precipitation δ18O variability (ER,RD=0.98; Fig. 6a) while at the same time largely preserving the limited but visible low-frequency temporal fluctuations in the stream δ18O composition (Fig. 3a, b). In comparison to the flow performance metrics, the Nash–Sutcliffe efficiency of the δ18O composition for the best model is somewhat lower (ENS,δ18O=0.37; Fig. 6a), which mostly results from the low variability of such a damped signal, where even very small absolute errors (MAE = 0.11 ) and a few scattered outliers can lead to very low Nash–Sutcliffe efficiencies (cf. Hrachowitz et al., 2009).

The posterior distributions (Table 2, Fig. 8) show that most model parameters are reasonably well identified. Individually calibrated for their respective landscape class, i.e. hillslope and riparian zone, SU,max,H=242 mm (5th/95th IQR: 213–311 mm) and SU,max,R=213 mm (186–280 mm) showed similar optimal values and distributions (Fig. 7a, b), reflecting the catchment-wide relatively homogenous forest cover in the pre-deforestation period (Fig. 1). Remarkably, these calibrated values also come close to catchment-scale estimates of SU, max=258±125mm that were directly derived from water balance data without any calibration, as described in Sect. 5.1 (Fig. 5c).

5.2.2 Application of the pre-deforestation model to the post-deforestation period

In a next step, the parameter sets obtained from the above calibration in the pre-deforestation period were used to run the model without further re-calibration in the post-deforestation period. This entails the implicit and clearly wrong assumption that the physical characteristics of the system remained unaffected by deforestation. The consequence of that can be seen in Fig. 2c and d (red line). While the low flows remain well reproduced, the post-deforestation application of the model substantially and systematically underestimates high flows, partly by 50 % or more, such as in November 2013 or August 2014. The inability of the model to reproduce several aspects of post-deforestation high-flow dynamics of the system is also evident in the lower model performance metrics associated with high flows (Fig. 6a). Besides the time series of flow (ENS,Q=0.65), notably the model's skill in capturing the rising limb density (ENS,PD=0.78), the autocorrelation function (ENS,AC=0.58; Fig. 6g) and the runoff ratio (ER,CR=0.81) were negatively affected. In contrast to the pre-deforestation period, the modelled runoff ratio CR=0.55 (0.54–0.58) in the post-deforestation period considerably underestimates the observed CR=0.68±0.03 (Fig. 5a). The problems of describing the high-flow periods are accompanied by the model's reduced ability to describe the post-deforestation δ18O dynamics in stream water (ENS,δ18O=0.11), although the observed general degree of damping of the δ18O signal (ER,RD=0.98) remains well reproduced as shown in Figs. 3 and 6a. While the low ENS,δ18O values are partly an effect of the above-explained low signal-to-noise ratio of such a damped signal and thus of the chosen performance metric, the model also struggles to adequately reproduce the lower-frequency fluctuations, such as between February and July 2014, when the model indicated rather stable δ18O values, while the observed values show a slight yet clear increasing trend over the same period (Fig. 3b). Together with the lower overall model performance metric DE (Fig. 6a), these results illustrate that the pre-deforestation model parameter sets provide an unsuitable characterization of the system characteristics in the post-deforestation period.

5.2.3 Recalibrate the model for the post-deforestation period

To estimate the effect of forest removal on the characteristics of the hydrological system and thus on the model parameters, the model was in a next step recalibrated for the post-deforestation period. This led to a slight improvement of the overall model performance from DE=0.77 to 0.80 (Fig. 6a). Most notably, it can be observed that the recalibrated model can much better reproduce the increased high flows in that period (Fig. 2c, d), as reflected by improvements in the performance metrics associated with high flows (Fig. 6a), but most notably ENS,Q=0.70, ENS,FDC=0.95 (Fig. 6c) or ENS,AC=0.92 (Fig. 6g). Similarly, the limits-of-acceptability constraints ensured a choice of solutions that broadly reflect the observed throughfall ratios PE/P (Fig. 7b) as well as the observed soil moisture dynamics (Fig. 7d, f). In addition, and perhaps most importantly, the runoff ratio also increased and was with a modelled value of CR=0.62 (0.56–0.63) closer to the observed CR=0.68 (ER,CR=0.91). This further implies that, in contrast to the initial model, the recalibrated model also features expected reductions of evaporative fluxes EA by about 10 %, which can be seen in Fig. 2b. Mirroring the improvements in the reproduction of flows, re-calibration also allowed the model to better capture the stream water δ18O dynamics (ENS,δ18O=0.24; MAE = 0.10 ; Fig. 6a). While there is little change in the model's ability to mimic the general level of damping of the δ18O signal (ER,RD=0.99) and its low-frequency fluctuations, the more pronounced, albeit in absolute terms still small, high-frequency fluctuations as short-term responses to individual storms are better described (Fig. 3a, b).

Inspection of the posterior parameter distributions reveals that the catchment-scale SU, max experienced considerable reductions after re-calibration. While in the hillslope parts of the catchment, which were less affected by deforestation (∼10% of the hillslope area; Fig. 1), an average decrease by ∼50mm to SU,max,H=212mm (137–270 mm) can be seen (Fig. 8a), the completely deforested riparian area exhibits an average decrease by ∼100mm to SU,max,R=93 mm (92–190 mm; Fig. 8b). As an indicative value, the area-weighted catchment average SU, max=199 mm of the best-performing parameter set falls into the plausible range of SU, max=101±149mm as described in Sect. 5.1. While there is little evidence for reductions of Imax on the less deforested hillslopes (Fig. 8d), a clear decrease in interception capacities by on average ∼2mm to Imax,R=1.1 mm (0.1–1.3 mm; Fig. 8e) can be observed in the fully deforested riparian zone.

Figure 9Panels in the left-hand-side column show pre-deforestation (a) discharge, and the coloured dots indicate to which period (dry, wet-up, wet, drying) the individual selected time steps belong; (b) the 5th/95th percentiles of the empirical cumulative TTDs for wet (blue) and dry (red) periods, respectively; (c) the ensemble of the individual TTDs at the time steps indicated in (a). Panels in the middle column (d–g) compare the 5th/95th percentiles of empirical cumulative TTDs between pre-deforestation (dark shades) and post-deforestation (light shades) periods for dry, wet-up, wet and drying conditions, respectively. Panels in the right-hand-side column show post-deforestation (h) discharge, and the coloured dots indicate to which period (dry, wet-up, wet, drying) the individual selected time steps belong; (i) the 5th/95th percentiles of the empirical cumulative TTDs for wet (blue) and dry (red) periods, respectively; (j) the ensemble of the individual TTDs at the time steps indicated in (h). All distributions shown are truncated at 3 (post-deforestation) or 4 years (pre-deforestation), which coincides with the tracked periods. For the remaining fractions, i.e. the difference to 1, it can only be said that they are older than 3 or 4 years but nothing more than that. The grey shaded areas indicate regions with ages >3 months, thereby exceeding Fyw.


Figure 10(a, b) Pre- and post-deforestation time series of young water fractions Fyw in discharge. The colour code indicates the transition between dry, wetting-up, wet and drying conditions. The bold black line shows the mean Fyw of the best model fit, and the grey shaded area shows the 5th/95th percentile of Fyw for all feasible model solutions; (c, d) pre- and post-deforestation sensitivity of Fyw to discharge, using the same colour code as above to indicate dry, wetting-up, wet and drying conditions. The arrows in (d) indicate whether there are statistically significant (; p<0.05) changes or not () in the sensitivities between the post-deforestation period and the pre-deforestation period.


Figure 11Individual catchment overall SAS ω functions for individual time steps under different wetness conditions in the (a) pre-deforestation period and (b) post-deforestation period. The insets show the relative water content in SU,rel,mod=SU/SU, max at the individual time steps.


5.3 Deforestation effects on travel time distributions, SAS functions and young water fractions

While the volume-weighted mean δ18O compositions of observed precipitation with 7.9  and stream water with 8.2  are comparable, a substantial difference in their fluctuations, with standard deviations of 3.6  and 0.2 , respectively, is evident (Fig. 3a, b). This difference suggests a remarkably elevated degree of damping rarely found elsewhere (e.g. Speed et al., 2010), indicative of the importance of old water contributions to the stream in the study catchment. No significant difference in damping ratios was observed between the pre- and post-deforestation periods, which further corroborates the prevalence of old water.

Tracking the δ18O signals through the model then allowed us to estimate TTDs. Note that any results reported hereafter are necessarily conditional on the assumptions made in and the uncertainties arising from the modelling process.

In general and consistent with the observed high degree of damping, it was found that pre-deforestation of the system was characterized by rather old water. The range of truncated TTDs of stream water exhibits considerable variability in response to changing wetness conditions, with on average about 27 % of the discharge younger than 3 years (Fig. 9b, c). In spite of the low mean Fyw∼0.12 (Fig. 10a), stream water can contain up to 34 % water younger than 3 months (i.e. Fyw∼0.34) for individual storm events in the wet period while frequently dropping to <1% during elongated summer dry periods (Figs. 8c, 9a), similar to what has been reported elsewhere (e.g. Gallart et al., 2020b). It can also be observed that the age composition of stream water (Fig. 9c) and the associated Fyw (Fig. 10a) do considerably vary throughout wet periods. Dry periods are characterized by considerably less variability and more stable stream water TTDs. This is further corroborated by the significantly higher sensitivity of Fyw to changes in streamflow in wet-up and wet periods (dFyw/dQn0.35 and 0.25, respectively) as compared to dry periods (dFyw/dQn0.05; Fig. 10c). In spite of the low mean Fyw∼0.12 (Fig. 10a), the above also entails that very fast switches towards higher young water fractions can be observed when the system is wetting up after dry periods as well as for storm events throughout the wet season. In general, the above observations are also encapsulated in the catchment-overall storage-age selection functions ω that represent the ratio of stream water TTD over the combined RTD of all model storage elements (Benettin et al., 2015a). While for dry periods undersampling of young water ages with relatively little variability is evident, it can also be seen that in particular during wet-up and wet periods a considerable yet highly variable preference for very young water can be seen (Fig. 11a), similar to what has been reported previously in other environments (e.g. Benettin et al. 2015a; Remondi et al., 2018).

The overall picture did not change in the post-deforestation period. Similarly to the pre-deforestation period, the TTDs can exhibit considerable variability. However, in contrast to the pre-deforestation period, and depending on the wetness conditions, considerable shifts towards younger water can be observed for the TTDs (Fig. 9d–g). There are little discernible changes in Fyw during the dry summer months (Fig. 9d). However, storms in wet-up periods, mostly during autumn, led to considerable increases in the fractions of water younger than 10–20 d (Fig. 9e). During wet periods clear shifts towards younger water can be observed throughout the entire spectrum of tracked ages (Fig. 9f). During the wet period ∼36% of the stream water is on average younger than the tracked 3 years (Fig. 9i). The mean Fyw only slightly increased to 0.13 (Fig. 10b) compared to 0.12 in the pre-deforestation period (Fig. 10a), which corroborates earlier results by Stockinger et al. (2019) that suggested only minor fluctuations in mean Fyw over multiple moving time windows. For individual winter storm events, Fyw slightly increased to up to ∼0.37 (Figs. 9j, 10b) compared to Fyw of up to ∼0.34 in the pre-deforestation period (Figs. 9c, 10a). Besides the generally higher Fyw during wet periods, the Fyw became more sensitive to flow during wet conditions, with dFyw/dQn0.36 (Fig. 10d), similar to what has been previously reported by von Freyberg et al. (2018) and Gallart et al. (2020a). The above-described post-deforestation changes are also manifest in the corresponding storage-age selection function ω (Fig. 11b) for that period. While the degree of undersampling of young water during dry periods significantly decreased, a substantially higher preference for young water during wet-up and wet periods can be observed than during the pre-deforestation period, with a clear overall shift towards younger water for all wetness conditions.

6 Discussion

6.1 Observed deforestation effects on the hydrological system

The observed post-deforestation changes to the hydrological response, in particular the increase in CR from ∼0.55 to ∼0.68, correspond well to the findings of an earlier study in the Wüstebach, based on a shorter study period (2011–2015; Wiekenkamp et al., 2016), which estimated an increase in CR from ∼0.58 to ∼0.66 during that period using eddy-covariance measurements. The overall pattern found here also broadly reflects the effects of land cover/use change in many different environments (Creed et al., 2014; Jaramillo and Destouni, 2014; Renner et al., 2014; van der Velde et al., 2014; Moran-Tejada et al., 2015; Nijzink et al., 2016a; Zhang et al., 2017; Jaramillo et al., 2018). The vast majority of these studies suggest that forest removal leads to an increase in the runoff ratio CR at the cost of reduced evaporation EA, although the magnitudes of these changes do substantially vary between individual catchments and studies, which is consistent with our physical understanding of the importance of forest for transpiration in hydrological systems.

Under the assumption that reduction of EA is largely a direct consequence of forest removal in the Wüstebach, a plausible hypothesis to directly attribute this shift in water partitioning from EA to Q to a physical process can be formulated as follows: the roots of harvested trees stopped extracting water for transpiration from the sub-surface. In addition, the limited turbulent exchange of vapour at depth effectively limits soil evaporation to the first few centimetres of the soil (e.g. Brutsaert, 2014). Thus, the felling of trees led to a situation where under comparable atmospheric water demand EP, water volumes held at depths below that and previously within the reach of active roots became largely unavailable for transpiration and evaporation after deforestation. This implies that the water volumes accessible to satisfy atmospheric water demand, i.e. SU, max and Imax, are drastically reduced. Most notably, the available water balance data suggest that catchment-scale SU, max decreased from pre-deforestation SU, max=258±125mm to post-deforestation SU, max=101±149mm.

Note, however, that in particular the estimates for the post-deforestation period are characterized by considerable uncertainty and therefore need to be understood as merely indicative, as they are inferred from only 3 years of data and a system that is likely to be far from equilibrium, because the deforested part cannot have adapted yet (e.g. Nijzink et al., 2016a; Teuling and Hoek van Dijke, 2020). These considerable uncertainties are also reflected in the surprisingly low post-deforestation SU, max.

Notwithstanding these limitations, the above results illustrate that here the reduction of transpiration due to deforestation is likely a direct consequence of the considerable reduction of SU, max and thus the catchment-scale sub-surface pore volume between field capacity and permanent wilting point that can be actively accessed by vegetation to satisfy the evaporative demand. These post-deforestation decreases in transpiration due to reductions in accessible water volumes SU, max further lead to reduced soil water storage deficits SD,j (Eq. 2) in dry seasons, which is consistent with observed post-deforestation increases in soil moisture (Wiekenkamp et al., 2016).

6.2 Modelled deforestation effects on the hydrological system

The model application provided further evidence for the central role of SU, max as a dominant control on the hydrological response as well as for the direct effects of deforestation on SU, max. The model calibration in the pre-deforestation period resulted in a set of solutions that could simultaneously reproduce multiple signatures, as expressed by 14 individual performance metrics, while also satisfying two additional limits-of-acceptability constraints. Overall this suggests a rather robust representation of the system.

In a next step, the parameter sets obtained from the calibration in the pre-deforestation period were used to run the model without further re-calibration in the post-deforestation period. This entails the implicit and clearly wrong assumption that the physical characteristics of the system remained unaffected by deforestation. As a consequence, that model exhibited a considerably reduced ability to reproduce the hydrological response in the post-deforestation period, in particular high flows as well as the runoff ratio CR. The latter implies that the model also overestimates post-deforestation evaporative fluxes EA. Therefore, it can, without re-calibration, not deal with the observed changes in the partitioning between drainage and evaporative fluxes (Fig. 5a). A likely explanation for the pattern produced by the model is that, in contrast to the real world, no reduction in EA due to the reduced forest cover is achieved because the model still relies on the catchment-scale vegetation-accessible storage volume SU, max that characterizes the extent of the catchment-scale active root system before deforestation. This SU, max falsely provides sufficient water supply to sustain EA at high levels comparatively close to EP throughout the year (see red line in Fig. 2b), although, in the parts of the catchment where trees were removed, water stored at depths below a few centimetres is not available for significant evaporation anymore in reality. Such an overestimation of SU, max implies also that in the model a more pronounced water storage deficit can and does develop throughout dry periods. The model therefore assumes that soils dry out to deeper depths. Consequently, to establish connectivity and to eventually generate flow during and after rainstorms, more water needs to be stored in the model than in the real-world system to overcome this deficit. This water is then in the model held against gravity and thus only available for evaporation but not for drainage, thereby underestimating in particular the magnitude of high flows. Although it is reasonable to assume that groundwater recharge is affected in a similar way, the model can better reproduce low flows. The reason for this is that the draining groundwater body, which sustains summer low flows, is, due to limited recharge during these drier periods, largely disconnected from and thus largely unaffected by sub-surface–vegetation interaction in shallower parts of the sub-surface. In the parts of the catchment where trees were removed, a similar reasoning also holds for the interception capacity Imax and the associated likely overestimation of interception evaporation EI, yet, due to the smaller magnitude of Imax, to a lesser extent than for SU, max.

Re-calibration of the model in the post-deforestation period led to a considerably improved representation of the hydrological response and in particular of the high flows as well as the runoff ratio CR. The latter implies that the modelled partitioning of water fluxes and in particular EA (see orange line in Fig. 2b) is more consistent with the observed post-deforestation reductions in EA. In addition, analysis of the modelled fluxes indicates that a higher proportion of flows, mostly during wet-up periods, is rapidly released from the root zones as fluxes RF,H and RF,R (Fig. 4; Table 1), representing preferential flows. Such a post-deforestation increase in preferential flow occurrence is supported by observations recently reported by Wiekenkamp et al. (2020).

It is of course unsurprising that re-calibration leads to an improved model performance in the post-deforestation period. Without further analysis, such a mere model fitting exercise allows in the presence of model equifinality only little insight into the underlying processes (Beven, 2006; Kirchner, 2006). To gain more confidence that the improvements in the recalibrated model are at least partly due to the right reasons (Kirchner, 2006), the changes in the posterior parameter distributions resulting from the two calibration runs were thus analysed. In the pre-deforestation period, the range of the posterior distributions of SU,max,H and SU,max,R (Fig. 8a, b) as well as the modelled catchment average SU, max=240 mm, estimated as an area-weighted average of SU,max,H and SU,max,R, come close to the catchment-scale estimate of SU, max=258±125mm that was directly derived from water balance data without any calibration (Fig. 5c). The modelled post-deforestation reductions of SU,max,H and SU,max,R are evident in the shifts of their respective posterior distributions (Fig. 8a, b) and the lower catchment average SU, max=199 mm of the best-performing parameter set, falling into the plausible range of SU, max=101±149mm as estimated from water balance data. In addition and quite remarkably, the re-calibrated model is able to broadly represent the differences in forest removal on the hillslopes and in the riparian zone. While in the fully deforested riparian area SU,max,R decreased by ∼100mm (Fig. 8b), SU,max,H on the only partly deforested hillslopes decreased by merely ∼50mm (Fig. 8a). Similarly, there is little evidence for reductions of Imax,H on the less deforested hillslopes (Fig. 8d). However, a clear decrease in interception capacities by on average ∼2mm to Imax,R=1.1 mm (0.1–1.3 mm; Fig. 8e) can be observed in the riparian zone. Comparing to the posterior distributions of other parameters, the results illustrate that the storage parameters SU, max and Imax of the completely deforested riparian zone, and to a lesser extent of the hillslope, were subject to the most pronounced changes. For most other parameters, the pre- and post-deforestation posterior distributions exhibit much less pronounced differences (Fig. 8). Together, these results suggest that deforestation mostly affects SU, max and Imax, while there is less evidence for systematic changes in other parameters. However, it can also be observed that the individual parameter values associated with the best model solutions in the pre- and post-deforestation periods, respectively, do vary to a stronger degree for most parameters. Notwithstanding the distinct overall effects of forest removal on the individual posterior distributions, this clearly highlights the influence of parameter compensation effects and related uncertainties. This is also illustrated by a few parameters, such as RS,max (Fig. 8c, Eq. 22), that remain poorly constrained.

It was hypothesized above that reductions in evaporative fluxes are directly and exclusively linked to reduced water volumes SU, max and Imax, respectively, which are accessible and available for evaporation and transpiration at the catchment scale. In the theoretical ideal case, the representations of the associated storage capacities in the model, i.e. the parameters SU, max and Imax, should thus be the only ones to significantly change after deforestation. However, note that this is unlikely for two reasons. First, while it is plausible to assume that these storage capacities are significantly affected by forest removal, it is not unlikely that other system characteristics and their mutual interactions, so far unknown and not considered, are similarly influenced, potentially causing considerable ontological uncertainty. Second, model parameter interactions that arise as artefacts to compensate overly simplistic process representations and/or data uncertainty are also likely to affect parameters seemingly unrelated to deforestation. Note that in spite of these uncertainties and the associated compensation effects, in particular SU, max remains rather well constrained. However, after preliminary unsuccessful testing, no further attempts were made to re-calibrate only the above discussed four storage parameters, i.e. SU,max,H, SU,max,R, Imax,H and Imax,R, acknowledging the limitations introduced by parameter compensation effects.

Overall these results suggest that the model formulation together with the multi-objective calibration strategy ensured the identification of solutions that provide a robust description of the system and allow a simultaneous representation of flow and isotope dynamics in the stream. There are indications that at least some processes and parameters can be directly linked to real world quantities. In particular, the results provide strong evidence that the parameters SU,max,H and SU,max,R are not merely abstract quantities, but that it is plausible to assume that they, taken together, provide a catchment-scale representation of vegetation-accessible and vegetation-accessed water volumes, which can be estimated based on water balance data without calibration as defined by Eq. (2), thereby providing an alternative to small-scale in situ observations. As such, the parameter SU, max is also a means to directly and independently estimate the catchment-scale effects of deforestation, and plausibly other types of land cover disturbances, on sub-surface system properties, which underlie and control the changes in the post-disturbance partitioning of water fluxes into drainage and evaporative fluxes.

6.3 Deforestation effects on travel time distributions, SAS functions and young water fractions

Tracking water fluxes through the system, it was observed that wet periods are characterized by substantially more variability and more stable stream water TTDs than dry periods. This is largely a consequence of increased bypass flow that has little interaction with resident water as the system gets wetter and that may reach the stream over preferential flow paths and increased contributions from the riparian zone with its shorter flow paths. In other words, in a wet system where little additional water can be stored, the precipitation volumes of individual storm events control the shape of TTDs, resulting in considerable variability (Heidbüchel et al. 2020). In the summer dry season, however, precipitation is to a higher degree buffered in the root zone and used for transpiration (Stockinger et al., 2014). Conversely, streamflow is then mostly sustained by groundwater which is characterized by large volumes of older water. This effectively attenuates fluctuations by the proportionally much lower volumes of younger precipitation water that cannot be stored and is thus quickly released to the stream.

In particular, at the beginning of the wet period, elsewhere also referred to as “autumn flush” (e.g. Dawson et al., 2011), the switches towards younger water at given flow levels occur considerably faster in the post-deforestation period than in the pre-deforestation period. Therefore, where, at the same discharge, previously relatively little young water reached the stream, a much higher fraction of young water can now be observed in the stream. Underlining the role of transpiration (e.g. Douinot et al., 2019; Kuppel et al., 2020), this is a direct effect of the reduced evaporative removal of relatively young near-surface water (Maxwell et al., 2019) in the post-deforestation period, which in turn is intimately linked to the reduced water supply for evaporative fluxes, i.e. smaller storage volumes SU, max and Imax. This modelled relatively young, surface-near water, not taken up by vegetation anymore is thus to a higher degree flushed from the system mostly via preferential flow paths to the stream (i.e. RF,H, RF,R) and thus bypassing older resident water with little exchange, which is consistent with recent observations of more frequent activation of preferential flow paths (Wiekenkamp et al., 2020). Once connectivity and the associated higher degree of bypass flow are established in the wet period, the post-deforestation peak sensitivity of dFyw/dQ to flow increased to ∼0.36, as under these conditions when little additional water can be stored in the shallow sub-surface, Fyw is largely controlled by magnitude of the individual precipitation signals and to a lesser extent by the footprint of the pre-storm history of evaporative fluxes in the shallow sub-surface storage. In contrast, no significant post-deforestation changes could be observed for the sensitivity of Fyw to discharge during dry periods, as during that period, the composition of water ages is controlled by large volumes of old water.

Altogether these results suggest that even in systems dominated by old water, such as the Wüstebach, the removal of forest has the potential to increase the importance of bypass flow through fast flow paths and thus increase the risk of fast, often underestimated propagation of contaminant pulses into groundwater and stream water (e.g. Hartmann et al., 2021).

6.4 Uncertainties, unresolved questions and limitations

As emphasized above, all results are conditional on the assumptions made throughout the modelling process. These assumptions, present in model structure, parameterization and parameters, can lead to uncertainties. However, notwithstanding these potential uncertainties, extensive preliminary model testing together with the use of multiple model calibration and evaluation criteria suggests that there is relatively strong evidence to support the main results in this study: the post-deforestation reduction of evaporative fluxes can, at least partially, be linked to a relatively clear reduction in the catchment-scale storage capacities SU, max and Imax, which in turn triggered a shift towards younger water ages in the stream, particularly during wet-up and wet conditions.

This is further corroborated when comparing the estimates of SU, max to estimates of physically plausible upper limits of SU, max. By definition, SU, max is physically bound by the depth of the groundwater table. Although fluctuating, the groundwater table in the Wüstebach remains at depths below 1 m for much of the year even in the riparian zone (Bogena et al., 2015) and can be expected to be considerably deeper on the hillslopes. Thus assuming a conservative upper bound of catchment-average depth of the groundwater table at ∼5m, assuming that the lowest groundwater table at each point in the catchment is at the elevation of the nearest stream, a porosity of the silty clay loam soil of 0.4 (Bogena et al., 2018) and field capacity at a relative porewater content of 0.5 suggests an upper limit of SU,max,GW∼1000mm. However, actual roots are very often shallower than these 5 m of the groundwater table. Although sufficient detailed data on root depths are not available in the study catchment, there is no evidence for systematic and widespread roots extending to below 2 m. This is broadly consistent with direct experimental evidence that roots of temperate forests in general (Schenk and Jackson, 2002) and Picea species in particular mostly remain rather shallow (<1m; e.g. Schmid and Kazda, 2001) and with indirect evidence that Picea species rarely tap groundwater and are thus comparatively shallow (e.g. Evaristo and McDonnell, 2017). As a conservative back-of-the-envelope calculation, assuming thus a maximum plausible catchment-average root depth of 2 m, which comes close to the average observed soil depth reported in Graf et al. (2014), rather suggests a physically plausible upper limit of SU,max,RD∼400mm, which is not exceeded by the water balance inferred catchment-scale estimates of SU, max=258±125mm.

Note that the above also suggests the presence of an unsaturated transition zone between the root zone and the groundwater table, i.e. SU,max,TZ=SU,max,GW-SU,max,RD600mm. In the absence of root water uptake and likely negligible soil evaporation in that zone the water content will remain close to field capacity for much of the year, except for days when a wetting front infiltrates towards the groundwater. This transition zone can therefore be considered as hydrologically largely passive so that at timescales of more than a few days dS/dt0. However, this zone also provides a mixing volume that affects tracer circulation and thus water ages (Hrachowitz et al., 2015). Given its hydrologically passive nature and following the idea of a parsimonious model to limit uncertainty, we here, in a simplification, implicitly added the mixing volume SU,max,TZ to the passive groundwater mixing volume SS,p.

For a meaningful interpretation, two specific observations resulting from our analysis warrant special scrutiny. First, model calibration-based estimations of hillslope SU,max,H (Fig. 7a) suggest post-deforestation median SU,max,H reductions of ∼25% as a consequence of clear cutting only ∼10% of the hillslope part of the catchment (Fig. 1). While this may be surprising at first, it can be plausibly explained by considerable further thinning of the remaining forest on the hillslopes in 2015, 2 years after deforestation, and thus by reduced catchment-scale transpiration demand. However, no detailed and systematic data on the degree of forest thinning are available to meaningfully test this hypothesis.

Second, our results suggest that a passive mixing volume SS,p of at least ∼8000mm is necessary for the model to attenuate the amplitudes of the precipitation δ18O signals to those in the stream water. Although SS,p is rather well constrained (Fig. 8h), there has in the past been no hydrogeological evidence for the presence of such a surprisingly large groundwater volume nor for its hydrological relevance in the study catchment. Indeed, the authors are not aware of any catchment-scale study that reported similarly high values for SS,p or functionally equivalent parameters (e.g. Birkel et al., 2011a, b; Hrachowitz et al., 2013, 2015; Benettin et al., 2013, 2015a; Harman, 2015; van der Velde et al., 2015). However, to achieve the degree of damping observed in the stream water, such a volume is necessary if the current understanding of conservative tracer dynamics holds (e.g. Maloszewski and Zuber, 1982; McGuire and McDonnell, 2006). Reflecting our insufficient knowledge of the depth to which exchange with surface water occurs (e.g. Condon et al., 2020), a potential explanation for this observation is that the frequently layered and fractured structure of the Devonian shale bedrock may provide relatively high-permeability pathways for the circulation of and exchange with water at depth. Another, yet, given the current understanding of the Wüstebach (e.g. Graf et al., 2014), less likely hypothesis is the presence of significant lateral groundwater exchange (e.g. Bouaziz et al., 2018; Hulsman et al., 2021b). In other words, the possibility exists that the sub-surface catchment does not match with the surface catchment (Fig. 1) and that older groundwater is imported from “outside” the surface catchment, while an equivalent volume of younger groundwater is exported, maintaining the mass balance. These are hypotheses to be tested in future studies, as the currently available data do not allow a conclusive answer to this question.

7 Conclusions

The small Wüstebach catchment experienced significant deforestation in 2013. Analysing the effects of this deforestation on the hydrology and stable isotope circulation dynamics in the study catchment, our main findings are the following.

Water balance data suggest that deforestation led to a significant increase in streamflow, accompanied by corresponding reductions of evaporative fluxes. This is reflected by an increase in the runoff ratio from CR=0.55 to 0.68 in the post-deforestation period despite similar climatic conditions, supporting previous results based on eddy-covariance measurements (Wiekenkamp et al., 2016).

Based on water balance data, this reduction of evaporative fluxes, as a consequence of reduced vegetation water uptake, could at least partly be linked to a reduction of the catchment-scale water storage volume in the unsaturated soil that is within the reach of active roots and thus accessible for vegetation transpiration from SU, max∼258mm in the pre-deforestation period to SU, max∼101mm in the post-deforestation period.

Estimating SU, max as the calibration parameter of a process-based hydrological model led to similar conclusions. The catchment-average calibrated model parameters representing SU, max for both the pre-deforestation and post-deforestation periods, respectively, correspond to ∼240 and ∼199mm broadly, with SU, max directly estimated from water balance data. Other model parameters, assumed to have a less direct link to vegetation, exhibited much lower levels of systematic change following deforestation.

Using the model to track the age composition of stream water suggested that, in general, water reaching the stream in the pre-deforestation period was rather old, with a mean young water fraction Fyw0.12. In spite of the overall low Fyw, clear shifts in the shape of travel time distributions towards younger water can be seen under wet conditions, with young water fractions increasing up to Fyw0.34.

Deforestation and the associated reduction of SU, max led to shifts in travel time distributions towards younger water. Under wet conditions, this resulted in increases in young water fractions to up to Fyw∼0.37 for individual storms. In contrast, dry-period travel time distributions exhibited only minor changes. Overall the mean fraction of young water in the stream increased to Fyw0.13.

Deforestation resulted in a considerable increase in the sensitivity of young water fractions to discharge under wet conditions from dFyw/dQ=0.25 to 0.36. This implies faster switches towards younger water and thus faster routing of solutes during and shortly after storm events and thus faster routing of solutes with increasing wetness.

The above results suggest that deforestation has not only the potential to affect the partitioning between drainage and evaporation, and thus the fundamental hydrological response characteristics of catchments, but also catchment-scale tracer circulation dynamics. In particular for wet and wet-up conditions, sometimes also referred to as “autumn flush”, deforestation in the Wüstebach caused higher proportions of younger water to reach the stream, implying faster routing of water and plausibly also solutes through the sub-surface, thereby also increasing the risk of faster propagation of contaminants into stream water and groundwater.

Overall, this study demonstrates that post-deforestation changes in both the hydrological response and travel times can to a large extent be traced back and attributed to changes in SU, max, a readily quantifiable catchment-scale sub-surface property (and model parameter) representing the maximum water volume that can be stored within the reach of roots. As such, SU, max and changes therein provide a quantitative, mechanistic hypothesis that can explain why deforestation in the Wüstebach decreased evaporative fluxes, increased streamflow – particularly generated by preferential flows – and reduced travel times. The catchment-scale quantification of SU, max based on water balance data therefore provides a potentially valuable way towards meaningful and data-based catchment-scale representation of vegetation-accessible water where soil and root observations are not available at sufficient spatial and temporal detail to meaningfully represent their respective natural heterogeneities. In addition, and perhaps more importantly, the method may also hold considerable potential for the formulation of temporally adaptive root-zone parameterizations in catchment-scale hydrological models for more reliable predictions in a changing environment.

Code availability

The model code used can be made available by the first author upon request. The equations used in the model are described in the paper.

Data availability

The meteorological and hydrological data of the Wüstebach TERENO site used in this study can be made available by the co-author HB upon request. The model results, including states, fluxes, hydrological signatures, parameter sets, and performance metrics underlying this paper, are available online in the 4TU data repository at (Hrachowitz, 2021).

Author contributions

MH and MS designed the experiment. MH did the analysis and wrote the first draft. All the authors discussed the design, results and the first draft and contributed to writing the final manuscript.

Competing interests

The authors declare that they have no conflict of interest.


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


We would like to thank the editor and four anonymous reviewers for providing a list of critical and very valuable comments that helped to considerably improve the manuscript.

Review statement

This paper was edited by Markus Weiler and reviewed by four anonymous referees.


Ala-aho, P., Tetzlaff, D., McNamara, J. P., Laudon, H., and Soulsby, C.: Using isotopes to constrain water flux and age estimates in snow-influenced catchments using the STARR (Spatially distributed Tracer-Aided RainfallRunoff) model, Hydrol. Earth Syst. Sci., 21, 5089–5110,, 2017. 

Angermann, L., Jackisch, C., Allroggen, N., Sprenger, M., Zehe, E., Tronicke, J., Weiler, M., and Blume, T.: Form and function in hillslope hydrology: characterization of subsurface flow based on response observations, Hydrol. Earth Syst. Sci., 21, 3727–3748,, 2017. 

Arrouays, D., Lagacherie, P., and Hartemink, A. E.: Digital soil mapping across the globe, Geoderma Regional, 9, 1–4,, 2017. 

Arsenault, R., Poissant, D., and Brissette, F.: Parameter dimensionality reduction of a conceptual model for streamflow prediction in Canadian, snowmelt dominated ungauged basins, Adv. Water Resour., 85, 27–44, 2015. 

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

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, 2015a. 

Benettin, P., Rinaldo, A., and Botter, G.: Tracking residence times in hydrological systems: Forward and backward formulations, Hydrol. Process., 29, 5203–5213, 2015b. 

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

Berkowitz, B. and Zehe, E.: Surface water and groundwater: unifying conceptualization and quantification of the two “water worlds”, Hydrol. Earth Syst. Sci., 24, 1831–1858,, 2020. 

Bertuzzo, E., Thomet, M., Botter, G., and Rinaldo, A.: Catchment-scale herbicides transport: Theory and application, Adv. Water Resour., 52, 232–242, 2013. 

Beven, K.: A manifesto for the equifinality thesis. Journal of hydrology, 320, 18–36, 2006. 

Beven, K. J.: Preferential flows and travel time distributions: defining adequate hypothesis tests for hydrological process models, Hydrol. Process., 24, 1537–1547, 2010. 

Beven, K. and Germann, P.: Macropores and water flow in soils revisited, Water Resour. Res., 49, 3071–3092, 2013. 

Birkel, C. and Soulsby, C.: Linking tracers, water age and conceptual models to identify dominant runoff processes in a sparsely monitored humid tropical catchment, Hydrol. Process., 30, 4477–4493, 2016. 

Birkel, C., Tetzlaff, D., Dunn, S. M., and Soulsby, C.: Using lumped conceptual rainfall–runoff models to simulate daily isotope variability with fractionation in a nested mesoscale catchment, Adv. Water Resour., 34, 383–394, 2011a. 

Birkel, C., Soulsby, C., and Tetzlaff, D.: Modelling catchment-scale water storage dynamics: Reconciling dynamic storage with tracer-inferred passive storage, Hydrol. Process., 25, 3924–3936, 2011b. 

Birkel, C., Soulsby, C., and Tetzlaff, D.: Conceptual modelling to assess how the interplay of hydrological connectivity, catchment storage and tracer dynamics controls nonstationary water age estimates, Hydrol. Process., 29, 2956–2969, 2015. 

Blöschl, G., Bierkens, M. F., Chambel, A., et al.: Twenty-three unsolved problems in hydrology (UPH) – a community perspective, Hydrolog. Sci. J., 64, 1141–1158, 2019. 

Bogena, H. R., Bol, R., Borchard, N., Brüggemann, N., Diekkrüger, B., Drüe, C., Groh, J., Gottselig, N., Huisman, J. A., Lücke, A., and Missong, A.: A terrestrial observatory approach to the integrated investigation of the effects of deforestation on water, energy, and matter fluxes, Sci. China Earth Sci., 58, 61–75, 2015. 

Bogena, H.R., Montzka, C., Huisman, J.A., Graf, A., Schmidt, M., Stockinger, M., Von Hebel, C., Hendricks-Franssen, H.J., Van der Kruk, J., Tappe, W. and Lücke, A.: The TERENO-Rur Hydrological Observatory: A multiscale multi-compartment research platform for the advancement of hydrological science, Vadose Zone J., 17, 1–22,, 2018. 

Bogena, H. R., Stockinger, M., and Luecke, A.: Long-term stable water isotope data for the investigation of deforestation effects on the hydrological system of the Wüstebach catchment, Germany, Hydrol. Process., 35, e14006,, 2020. 

Borchardt, H.: Einfluss periglazialer Deckschichten auf Abflusssteuerung am Beispiel des anthropogen überprägten Wüstebaches (Nationalpark Eifel), PhD diss., Lehrstuhl für Physische Geogr. und Geoökologie, Fak. für Georessourcen und Materialtechnik, Rheinisch-Westfälische Tech. Hochsch. Aachen, Aaachen, Germany, 2012. 

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

Bouaziz, L., Weerts, A., Schellekens, J., Sprokkereef, E., Stam, J., Savenije, H., and Hrachowitz, M.: Redressing the balance: quantifying net intercatchment groundwater flows, Hydrol. Earth Syst. Sci., 22, 6415–6434,, 2018. 

Bouaziz, L., Steele-Dunne, S. C., Schellekens, J., Weerts, A., Stam, J., Sprokkereef, E., Winsemius, H. C., Savenije, H. H. G., and Hrachowitz, M.: Improved understanding of the link between catchment-scale vegetation accessible storage and satellite-derived Soil Water Index, Water Resour. Res., 56, e2019WR026365,, 2020. 

Bouaziz, L. J. E., Fenicia, F., Thirel, G., de Boer-Euser, T., Buitink, J., Brauer, C. C., De Niel, J., Dewals, B. J., Drogue, G., Grelier, B., Melsen, L. A., Moustakas, S., Nossent, J., Pereira, F., Sprokkereef, E., Stam, J., Weerts, A. H., Willems, P., Savenije, H. H. G., and Hrachowitz, M.: Behind the scenes of streamflow model performance, Hydrol. Earth Syst. Sci., 25, 1069–1095,, 2021. 

Brand, W. A., Coplen, T. B., Vogl, J., Rosner, M., and Prohaska, T.: Assessment of international reference materials for isotope-ratio analysis (IUPAC Technical Report), Pure Appl. Chem., 86, 425–467, 2014. 

Brooks, J. R., Barnard, H. R., Coulombe, R., and McDonnell, J. J.: Ecohydrologic separation of water between trees and streams in a Mediterranean climate, Nat. Geosci., 3, 100–104, 2010. 

Brooks, P. D., Chorover, J., Fan, Y., Godsey, S. E., Maxwell, R. M., McNamara, J. P., and Tague, C.: Hydrological partitioning in the critical zone: Recent advances and opportunities for developing transferable understanding of water cycle dynamics, Water Resour. Res., 51, 6973–6987, 2015. 

Brunner, I., Herzog, C., Dawes, M. A., Arend, M., and Sperisen, C.: How tree roots respond to drought, Front. Plant Sci., 6, 547, 2015. 

Brutsaert, W.: Daily evaporation from drying soil: Universal parameterization with similarity, Water Resour. Res., 50, 3206–3215, 2014. 

Cain, M. R., Ward, A. S., and Hrachowitz, M.: Ecohydrologic separation alters interpreted hydrologic stores and fluxes in a headwater mountain catchment, Hydrol. Process., 33, 2658–2675, 2019. 

Coenders-Gerrits, A. M. J., Van der Ent, R. J., Bogaard, T. A., Wang-Erlandsson, L., Hrachowitz, M., and Savenije, H. H. G.: Uncertainties in transpiration estimates. Nature, 506, E1-E2, 2014. 

Collins, D. B. G. and Bras, R. L.: Plant rooting strategies in water-limited ecosystems, Water Resour. Res., 43, W06407,, 2007. 

Condon, L. E., Markovich, K. H., Kelleher, C. A., McDonnell, J. J., Ferguson, G., and McIntosh, J. C.: Where is the bottom of a watershed?, Water Resour. Res., 56, e2019WR026010,, 2020. 

Cornelissen, T., Diekkrüger, B., and Bogena, H. R.: Significance of scale and lower boundary condition in the 3D simulation of hydrological processes and soil moisture variability in a forested headwater catchment, J. Hydrol., 516, 140–153, 2014. 

Coxon, G., Freer, J., Wagener, T., Odoni, N. A., and Clark, M.: Diagnostic evaluation of multiple hypotheses of hydrological behaviour in a limits-of-acceptability framework for 24 UK catchments, Hydrol. Process., 28, 6135–6150, 2014. 

Creed, I. F., Spargo, A. T., Jones, J. A., Buttle, J. M., Adams, M. B., Beall, F. D., Booth, E. G., Campbell, J. L., Clow, D., Elder, K., and Green, M. B.: Changing forest water yields in response to climate warming: Results from long-term experimental watershed sites across North America, Glob. Change Biol., 20, 3191–3208, 2014. 

Criss, R. E. and Winston, W. E.: Do Nash values have a value? Discussion and alternate proposals, Hydrol. Process., 22, 2723–2725, 2008. 

Dawson, J. J., Tetzlaff, D., Speed, M., Hrachowitz, M., and Soulsby, C.: Seasonal controls on DOC dynamics in nested upland catchments in NE Scotland, Hydrol. Process., 25, 1647–1658, 2011. 

deBoer-Euser, T., McMillan, H. K., Hrachowitz, M., Winsemius, H. C., and Savenije, H. H. G.: Influence of soil and climate on root zone storage capacity, Water Resour. Res., 52, 2009–2024, 2016. 

Donohue, R. J., Roderick, M. L., and McVicar, T. R.: On the importance of including vegetation dynamics in Budyko's hydrological model, Hydrol. Earth Syst. Sci., 11, 983–995,, 2007. 

Donohue, R. J., Roderick, M. L., and McVicar, T. R.: Roots, storms and soil pores: Incorporating key ecohydrological processes into Budyko's hydrological model, J. Hydrol., 436, 35–50, 2012. 

Douinot, A., Tetzlaff, D., Maneta, M., Kuppel, S., Schulte-Bisping, H., and Soulsby, C.: Ecohydrological modelling with EcH2O-iso to quantify forest and grassland effects on water partitioning and flux ages, Hydrol. Process., 33, 2174–2191, 2019. 

Eagleson, P. S. (Ed.): Land surface processes in atmospheric general circulation models, Cambridge University Press, New York, 1982. 

Ellison, D., Morris, C. E., Locatelli, B., Sheil, D., Cohen, J., Murdiyarso, D., Gutierrez, V., Van Noordwijk, M., Creed, I. F., Pokorny, J. and Gaveau, D.: Trees, forests and water: Cool insights for a hot world, Global Environ. Chang., 43, 51–61, 2017. 

Etmann, M.: Dendrologische Aufnahmen im Wassereinzugsgebiet Oberer Wüstebach anhand verschiedener Mess- und Schätzverfahren, MSc thesis, Institut für Landschaftsökologie, Univ. of Münster, Münster, Germany, 2009. 

Euser, T., Winsemius, H. C., Hrachowitz, M., Fenicia, F., Uhlenbrook, S., and Savenije, H. H. G.: A framework to assess the realism of model structures using hydrological signatures, Hydrol. Earth Syst. Sci., 17, 1893–1912,, 2013. 

Euser, T., Hrachowitz, M., Winsemius, H. C., and Savenije, H. H.: The effect of forcing and landscape distribution on performance and consistency of model structures, Hydrol. Process., 29, 3727–3743, 2015. 

Evaristo, J. and McDonnell, J. J.: Prevalence and magnitude of groundwater use by vegetation: a global stable isotope meta-analysis. Sci. Rep.-UK, 7, 1–12, 2017. 

Evaristo, J., Kim, M., van Haren, J., Pangle, L. A., Harman, C. J., Troch, P. A., and McDonnell, J. J.: Characterizing the fluxes and age distribution of soil water, plant water, and deep percolation in a model tropical ecosystem, Water Resour. Res., 55, 3307–3327, 2019. 

Fan, Y., Miguez-Macho, G., Jobbágy, E. G., Jackson, R. B., and Otero-Casal, C.: Hydrologic regulation of plant rooting depth, P. Natl. Acad. Sci. USA, 114, 10572–10577, 2017. 

Fenicia, F., Savenije, H. H. G., Matgen, P., and Pfister, L.: Understanding catchment behavior through stepwise model concept improvement, Water Resour. Res., 44, W01402,, 2008. 

Fenicia, F., Wrede, S., Kavetski, D., Pfister, L., Hoffmann, L., Savenije, H. H., and McDonnell, J. J.: Assessing the impact of mixing assumptions on the estimation of streamwater mean residence time, Hydrol. Process., 24, 1730–1741, 2010. 

Fovet, O., Ruiz, L., Hrachowitz, M., Faucheux, M., and Gascuel-Odoux, C.: Hydrological hysteresis and its value for assessing process consistency in catchment conceptual models, Hydrol. Earth Syst. Sci., 19, 105–123,, 2015. 

Freer, J., Beven, K., and Ambroise, B.: Bayesian estimation of uncertainty in runoff prediction and the value of data: An application of the GLUE approach, Water Resour. Res., 32, 2161–2173, 1996. 

Gallart, F., von Freyberg, J., Valiente, M., Kirchner, J. W., Llorens, P., and Latron, J.: Technical note: An improved discharge sensitivity metric for young water fractions, Hydrol. Earth Syst. Sci., 24, 1101–1107,, 2020a. 

Gallart, F., Valiente, M., Llorens, P., Cayuela, C., Sprenger, M., and Latron, J.: Investigating young water fractions in a small Mediterranean mountain catchment: both precipitation forcing and sampling frequency matter, Hydrol. Process., 34, 3618–3634,, 2020b. 

Gao, H., Hrachowitz, M., Schymanski, S. J., Fenicia, F., Sriwongsitanon, N., and Savenije, H. H. G.: Climate controls how ecosystems size the root zone storage capacity at catchment scale, Geophys. Res. Lett., 41, 7916–7923, 2014. 

Gao, H., Ding, Y., Zhao, Q., Hrachowitz, M., and Savenije, H. H.: The importance of aspect for modelling the hydrological response in a glacier catchment in Central Asia, Hydrol. Process., 31, 2842–2859, 2017. 

Gentine, P., D'Odorico, P., Lintner, B. R., Sivandran, G., and Salvucci, G.: Interdependence of climate, soil, and vegetation as constrained by the Budyko curve, Geophys. Res. Lett., 39, L19404,, 2012. 

Gerrits, A. M. J., Savenije, H. H. G., Veling, E. J. M., and Pfister, L.: Analytical derivation of the Budyko curve based on rainfall characteristics and a simple evaporation model, Water Resour. Res., 45, W04403,, 2009. 

Graf, A., Bogena, H. R., Drüe, C., Hardelauf, H., Pütz, T., Heinemann, G., and Vereecken, H.: Spatiotemporal relations between water budget components and soil water content in a forested tributary catchment, Water Resour. Res., 50, 4837–4857, 2014. 

Gumbel, E. J.: The return period of flood flows, Ann. Math. Stat., 12, 163–190, 1941. 

Guswa, A. J.: The influence of climate on root depth: A carbon cost-benefit analysis, Water Resour. Res., 44, W02427,, 2008. 

Harman, C. and Troch, P. A.: What makes Darwinian hydrology “Darwinian”? Asking a different kind of question about landscapes, Hydrol. Earth Syst. Sci., 18, 417–433,, 2014. 

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

Hartmann, A., Jasechko, S., Gleeson, T., Wada, Y., Andreo, B., Barberá, J. A., Brielmann, H., Bouchaou, L., Charlier, J.B., Darling, W.G., Filippini, M., Garvelmann, J., Goldscheider, N., Kralik, M., Kunstmann, H., Ladouche, B., Lange, J., Lucianetti, G., Martin, J. F., Mudarra, M., Sanchez, D., Stumpp, C., Zagana, E., and Wagener, T.: Risk of groundwater contamination widely underestimated because of fast flow into aquifers, P. Natl. Acad. Sci. USA, 118, e2024492118,, 2021. 

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

Hengl, T., Mendes de Jesus, J., Heuvelink, G. B., Ruiperez Gonzalez, M., Kilibarda, M., Blagotić, A., Shangguan, W., Wright, M. N., Geng, X., Bauer-Marschallinger, B., Guevara, M. A., Vargas, R., MacMillan, R., Batjes, N., Leenaars, J., Ribeiro, E., Wheeler, I., Mantel, S., and Kempen, B.: SoilGrids250m: Global gridded soil information based on machine learning, PLoS one, 12, e0169748,, 2017. 

Hrachowitz, M.: Model results underlying the paper “Reduction of vegetation-accessible water storage capacity after deforestation affects catchment travel time distributions and increases young water fractions in a headwater catchment”, 4TU.ResearchData [data set],, 2021. 

Hrachowitz, M. and Clark, M. P.: HESS Opinions: The complementary merits of competing modelling philosophies in hydrology, Hydrol. Earth Syst. Sci., 21, 3953–3973,, 2017. 

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, 2009. 

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

Hrachowitz, M., Fovet, O., Ruiz, L., Euser, T., Gharari, S., Nijzink, R., Freer, J., Savenije, H., and Gascuel-Odoux, C.: Process consistency in models: The importance of system signatures, expert knowledge, and process complexity, Water Resour. Res., 50, 7445–7469, 2014. 

Hrachowitz, M., Fovet, O., Ruiz, L., and Savenije, H. H. G.: Transit time distributions, legacy contamination and variability in biogeochemical 1/fα scaling: how are hydrological response dynamics linked to water quality at the catchment scale?, Hydrol. Process., 29, 5241–5256, 2015. 

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

Hulsman, P., Winsemius, H. C., Michailovsky, C. I., Savenije, H. H. G., and Hrachowitz, M.: Using altimetry observations combined with GRACE to select parameter sets of a hydrological model in a data-scarce region, Hydrol. Earth Syst. Sci., 24, 3331–3359,, 2020. 

Hulsman, P., Savenije, H. H. G., and Hrachowitz, M.: Learning from satellite observations: increased understanding of catchment processes through stepwise model improvement, Hydrol. Earth Syst. Sci., 25, 957–982,, 2021a. 

Hulsman, P., Hrachowitz, M., and Savenije, H. H.: Improving the representation of long-term storage variations with conceptual hydrological models in data-scarce regions, Water Resour. Res. 57, e2020WR028837, 2021b. 

Hwang, T., Band, L., and Hales, T. C.: Ecosystem processes at the watershed scale: Extending optimality theory from plot to catchment, Water Resour. Res., 45, W11425,, 2009. 

Jaramillo, F. and Destouni, G.: Developing water change spectra and distinguishing change drivers worldwide, Geophys. Res. Lett., 41, 8377–8386, 2014. 

Jaramillo, F., Cory, N., Arheimer, B., Laudon, H., van der Velde, Y., Hasper, T. B., Teutschbein, C., and Uddling, J.: Dominant effect of increasing forest biomass on evapotranspiration: interpretations of movement in Budyko space, Hydrol. Earth Syst. Sci., 22, 567–580,, 2018. 

Jasechko, S.: Plants turn on the tap, Nat. Clim. Change, 8, 562–563, 2018. 

Jothityangkoon, C., M. Sivapalan, and D. L. Farmer: Process controls of water balance variability in a large semi-arid catchment: Downward approach to hydrological model development, J. Hydrol., 254, 174–198, 2001. 

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

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

Klaus, J., Zehe, E., Elsner, M., Külls, C., and McDonnell, J. J.: Macropore flow of old water revisited: experimental insights from a tile-drained hillslope, Hydrol. Earth Syst. Sci., 17, 103–118,, 2013. 

Kleidon, A.: Global datasets of rooting zone depth inferred from inverse methods, J. Climate, 17, 2714–2722, 2004. 

Knighton, J., Souter-Kline, V., Volkman, T., Troch, P. A., Kim, M., Harman, C., Morris, C., Buchanan, B., and Walter, M. T.: Seasonal and topographic variations in ecohydrological separation within a small, temperate, snow-influenced catchment, Water Resour. Res., 55, 6417–6435, 2019. 

Knighton, J., Singh, K., and Evaristo, J.: Understanding Catchment-Scale Forest Root Water Uptake Strategies Across the Continental United States Through Inverse Ecohydrological Modeling, Geophys. Res. Lett., 47, e2019GL085937,, 2020. 

Kuppel, S., Tetzlaff, D., Maneta, M. P., and Soulsby, C.: EcH2O-iso 1.0: water isotopes and age tracking in a process-based, distributed ecohydrological model, Geosci. Model Dev., 11, 3045–3069,, 2018a. 

Kuppel, S., Tetzlaff, D., Maneta, M. P., and Soulsby, C.: What can we learn from multi-data calibration of a process-based ecohydrological model?, Environ. Modell. Softw., 101, 301–316, 2018b. 

Kuppel, S., Tetzlaff, D., Maneta, M. P., and Soulsby, C.: Critical zone storage controls on the water ages of ecohydrological outputs, Geophys. Res. Lett., 47, e2020GL088897,, 2020. 

Laio, F., D'Odorico, P., and Ridolfi, L.: An analytical model to relate the vertical root distribution to climate and soil properties, Geophys. Res. Lett., 33, L18401,, 2006. 

Loritz, R., Hassler, S. K., Jackisch, C., Allroggen, N., van Schaik, L., Wienhöfer, J., and Zehe, E.: Picturing and modeling catchments by representative hillslopes, Hydrol. Earth Syst. Sci., 21, 1225–1249,, 2017. 

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, 1982. 

Maxwell, R. M., Condon, L. E., Kollet, S. J., Maher, K., Haggerty, R., and Forrester, M. M.: The imprint of climate and geology on the residence times of groundwater, Geophys. Res. Lett., 43, 701–708, 2016. 

Maxwell, R. M., Condon, L. E., Danesh-Yazdi, M., and Bearup, L. A.: Exploring source water mixing and transient residence time distributions of outflow and evapotranspiration with an integrated hydrologic model and Lagrangian particle tracking approach. Ecohydrology, 12, e2042,, 2019. 

Mezentsev, V.: Back to the computation of total evaporation. Meteorologia i Gidologia, 5. 24–26, 1995. 

McDonnell, J. J., McGuire, K., Aggarwal, P., Beven, K. J., Biondi, D., Destouni, G., Dunn, S., James, A., Kirchner, J., Kraft, P. J. H. P., and Lyon, S.: How old is streamwater? Open questions in catchment transit time conceptualization, modelling and analysis, Hydrol. Process., 24, 1745–1754, 2010. 

McGuire, K. J. and McDonnell, J. J.: A review and evaluation of catchment transit time modeling, J. Hydrol., 330, 543–563, 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,, 2012. 

McMurtrie, R. E., Iversen, C. M., Dewar, R. C., Medlyn, B. E., Näsholm, T., Pepper, D. A., and Norby, R. J.: Plant root distributions and nitrogen uptake predicted by a hypothesis of optimal root foraging, Ecol. Evol., 2, 1235–1250, 2012. 

Milly, P. C. and Dunne, K. A.: Sensitivity of the global water cycle to the water-holding capacity of land, J. Climate, 7, 506–526, 1994. 

Mianabadi, A., Coenders-Gerrits, M., Shirazi, P., Ghahraman, B., and Alizadeh, A.: A global Budyko model to partition evaporation into interception and transpiration, Hydrol. Earth Syst. Sci., 23, 4983–5000,, 2019. 

Montanari, A. and Toth, E.: Calibration of hydrological models in the spectral domain: An opportunity for scarcely gauged basins?, Water Resour. Res., 43, W05434,, 2007. 

Morán-Tejeda, E., Zabalza, J., Rahman, K., Gago-Silva, A., López-Moreno, J. I., Vicente-Serrano, S., Lehmann, A., Tague, C., and Beniston, M.: Hydrological impacts of climate and land-use changes in a mountain watershed: Uncertainty estimation based on model comparison, Ecohydrology, 8, 1396–1416, 2015. 

Nash, J. E. and Sutcliffe, J. V.: River flow forecasting through conceptual models: part 1 – A discussion of principles, J. Hydrol., 10, 282–290, 1970. 

Nijzink, R., Hutton, C., Pechlivanidis, I., Capell, R., Arheimer, B., Freer, J., Han, D., Wagener, T., McGuire, K., Savenije, H., and Hrachowitz, M.: The evolution of root-zone moisture capacities after deforestation: a step towards hydrological predictions under change?, Hydrol. Earth Syst. Sci., 20, 4775–4799,, 2016a. 

Nijzink, R. C., Samaniego, L., Mai, J., Kumar, R., Thober, S., Zink, M., Schäfer, D., Savenije, H. H. G., and Hrachowitz, M.: The importance of topography-controlled sub-grid process heterogeneity and semi-quantitative prior constraints in distributed hydrological models, Hydrol. Earth Syst. Sci., 20, 1151–1176,, 2016b. 

Penna, D., Hopp, L., Scandellari, F., Allen, S. T., Benettin, P., Beyer, M., Geris, J., Klaus, J., Marshall, J. D., Schwendenmann, L., Volkmann, T. H. M., von Freyberg, J., Amin, A., Ceperley, N., Engel, M., Frentress, J., Giambastiani, Y., McDonnell, J. J., Zuecco, G., Llorens, P., Siegwolf, R. T. W., Dawson, T. E., and Kirchner, J. W.: Ideas and perspectives: Tracing terrestrial ecosystem water fluxes using hydrogen and oxygen stable isotopes challenges and opportunities from an interdisciplinary perspective, Biogeosciences, 15, 6399–6415,, 2018. 

Penna, D., Geris, J., Hopp, L., and Scandellari, F.: Water sources for root water uptake: Using stable isotopes of hydrogen and oxygen as a research tool in agricultural and agroforestry systems. Agriculture, Ecosystems and Environment, 291, 106790, 2020. 

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, 2018. 

Renner, M., Brust, K., Schwärzel, K., Volk, M., and Bernhofer, C.: Separating the effects of changes in land cover and climate: a hydro-meteorological analysis of the past 60 yr in Saxony, Germany, Hydrol. Earth Syst. Sci., 18, 389–405,, 2014. 

Richter, D.: Ergebnisse methodischer Untersuchungen zur Korrektur des systematische Meßfehlers des Hellmann-Niederschlagmessers, in: Berichte des Deutschen Wetterdienstes, vol. 194, 93 pp, LRG (Lysimeter Research Group), Selbstverlag des Deutschen Wetterdienstes, Offenbach am Main, Germany, 1995. 

Richter, F.: Bodenkarte zur Standorterkundung. Verfahren Quellgebiet Wüstebachtal (Forst), Geologischer Dienst Nordrhein-Westfalen, Krefeld, Germany, 2008. 

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

Rodriguez, N. B. and Klaus, J.: Catchment travel times from composite StorAge Selection functions representing the superposition of streamflow generation processes, Water Resour. Res., 55, 9292–9314, 2019. 

Rodriguez, N. B., McGuire, K. J., and Klaus, J.: Time-varying storage–water age relationships in a catchment with a mediterranean climate, Water Resour. Res., 54, 3988–4008, 2018. 

Rodriguez-Iturbe, I., D'Odorico, P., Laio, F., Ridolfi, L., and Tamea, S.: Challenges in humid land ecohydrology: Interactions of water table and unsaturated zone with climate, soil, and vegetation, Water Resour. Res., 43, W09301,, 2007. 

Roebroek, C. T. J., Melsen, L. A., Hoek van Dijke, A. J., Fan, Y., and Teuling, A. J.: Global distribution of hydrologic controls on forest growth, Hydrol. Earth Syst. Sci., 24, 4625–4639,, 2020. 

Savenije, H. H. G. and Hrachowitz, M.: HESS Opinions “Catchments as meta-organisms a new blueprint for hydrological modelling”, Hydrol. Earth Syst. Sci., 21, 1107–1116,, 2017. 

Sawicz, K., Wagener, T., Sivapalan, M., Troch, P. A., and Carrillo, G.: Catchment classification: empirical analysis of hydrologic similarity based on catchment function in the eastern USA, Hydrol. Earth Syst. Sci., 15, 2895–2911,, 2011. 

Seneviratne, S. I., Wilhelm, M., Stanelle, T., van den Hurk, B., Hagemann, S., Berg, A., Cheruy, F., Higgins, M. E., Meier, A., Brovkin, V., Claussen, M., Ducharne, A., Dufresne, J. M., Findell, K., Ghattas, J., Lawrence, D., Malyshev, S., Rummukainen, M., and Smith, B.: Impact of soil moisture-climate feedbacks on CMIP5 projections: First results from the GLACE-CMIP5 experiment, Geophys. Res. Lett., 40, 5212–5217, 2013. 

Schenk, H. J.: Vertical vegetation structure below ground: scaling from root to globe. In Progress in botany (pp. 341–373). Springer, Berlin, Heidelberg, 2005. 

Schenk, H. J. and Jackson, R. B.: Rooting depths, lateral root spreads and below-ground/above-ground allometries of plants in water-limited ecosystems, J. Ecol., 90, 480–494, 2002. 

Schenk, H. J. and Jackson, R. B.: Global distribution of root profiles in terrestrial ecosystems, ORNL DAAC, Oakridge, Tennessee, USA, 2003. 

Schmid, I. and Kazda, M.: Vertical distribution and radial growth of coarse roots in pure and mixed stands of Fagus sylvatica and Picea abies, Can. J. Forest Res., 31, 539–548, 2001. 

Schymanski, S. J., Sivapalan, M., Roderick, M. L., Beringer, J., and Hutley, L. B.: An optimality-based model of the coupled soil moisture and root dynamics, Hydrol. Earth Syst. Sci., 12, 913–932,, 2008. 

Shamir, E., Imam, B., Gupta, H. V., and Sorooshian, S.: Application of temporal streamflow descriptors in hydrologic model parameter estimation, Water Resour. Res., 41, W06021,, 2005. 

Sivandran, G. and Bras, R. L.: Identifying the optimal spatially and temporally invariant root distribution for a semiarid environment, Water Resour. Res., 48, W12525,, 2012. 

Speich, M. J. R., Lischke, H., and Zappa, M.: Testing an optimality-based model of rooting zone water storage capacity in temperate forests, Hydrol. Earth Syst. Sci., 22, 4097–4124,, 2018. 

Speich, M. J. R., Zappa, M., Scherstjanoi, M., and Lischke, H.: FORests and HYdrology under Climate Change in Switzerland v1.0: a spatially distributed model combining hydrology and forest dynamics, Geosci. Model Dev., 13, 537–564,, 2020. 

Soulsby, C., Tetzlaff, D., and Hrachowitz, M.: Are transit times useful process-based tools for flow prediction and classification in ungauged basins in montane regions?, Hydrol. Process., 24, 1685–1696, 2010. 

Speed, M., Tetzlaff, D., Soulsby, C., Hrachowitz, M., and Waldron, S.: Isotopic and geochemical tracers reveal similarities in transit times in contrasting mesoscale catchments, Hydrol. Process., 24, 1211–1224, 2010. 

Sprenger, M., Leistert, H., Gimbel, K., and Weiler, M.: Illuminating hydrological processes at the soil-vegetation-atmosphere interface with water stable isotopes, Rev. Geophys., 54, 674–704, 2016. 

Sprenger, M., Tetzlaff, D., Buttle, J., Laudon, H., and Soulsby, C.: Water ages in the critical zone of long-term experimental sites in northern latitudes, Hydrol. Earth Syst. Sci., 22, 3965–3981,, 2018. 

Sprenger, M., Llorens, P., Cayuela, C., Gallart, F., and Latron, J.: Mechanisms of consistently disjunct soil water pools over (pore) space and time, Hydrol. Earth Syst. Sci., 23, 2751–2762,, 2019a. 

Sprenger, M., Stumpp, C., Weiler, M., Aeschbach, W., Allen, S. T., Benettin, P., Dubbert, M., Hartmann, A., Hrachowitz, M., Kirchner, J.W., McDonnell, J. J., Orlowski, N., Penna, D., Pfahl, S., Rinderer, M., Rodriguez, N., Schmidt, M., and Werner, C.: The demographics of water: a review of water ages in the critical zone, Rev. Geophys., 57, 800–834, 2019b. 

Stockinger, M. P., Bogena, H. R., Lücke, A., Diekkrüger, 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, 2014. 

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

Stockinger, M. P., Bogena, H. R., Lücke, A., Diekkrüger, B., Cornelissen, T., and Vereecken, H.: Tracer sampling frequency influences estimates of young water fraction and streamwater transit time distribution, J. Hydrol., 541, 952–964, 2016. 

Stockinger, M. P., Lücke, A., Vereecken, H., and Bogena, H. R.: Accounting for seasonal isotopic patterns of forest canopy intercepted precipitation in streamflow modeling, J. Hydrol., 555, 31–40, 2017. 

Stockinger, M. P., Bogena, H. R., Lücke, A., Stumpp, C., and Vereecken, H.: Time variability and uncertainty in the fraction of young water in a small headwater catchment, Hydrol. Earth Syst. Sci., 23, 4333–4347,, 2019. 

Teuling, A. J.: A hot future for European droughts, Nat. Clim. Change, 8, 364–365, 2018. 

Teuling, A. J. and van Dijke, A. J. H.: Forest age and water yield, Nature, 578, E16–E18, 2020. 

Tron, S., Perona, P., Gorla, L., Schwarz, M., Laio, F., and Ridolfi, L.: The signature of randomness in riparian plant root distributions, Geophys. Res. Lett., 42, 7098–7106, 2015. 

Turc, L.: Le bilan d'eau des sols. Relations entre les precipitations, l'evaporation et l'ecoulement, Ann. Agron., 5, 491–596, 1954. 

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

Van der Velde, Y., Vercauteren, N., Jaramillo, F., Dekker, S. C., Destouni, G., and Lyon, S. W.: Exploring hydroclimatic change disparity via the Budyko framework, Hydrol. Process., 28, 4110–4118, 2014. 

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, 2015. 

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

van Oorschot, F., van der Ent, R. J., Hrachowitz, M., and Alessandri, A.: Climate-controlled root zone parameters show potential to improve water flux simulations by land surface models, Earth Syst. Dynam., 12, 725–743,, 2021. 

Wang-Erlandsson, L., Bastiaanssen, W. G. M., Gao, H., Jägermeyr, J., Senay, G. B., van Dijk, A. I. J. M., Guerschman, J. P., Keys, P. W., Gordon, L. J., and Savenije, H. H. G.: Global root zone storage capacity from satellite-based evaporation, Hydrol. Earth Syst. Sci., 20, 1459–1481,, 2016. 

Wang-Erlandsson, L., Fetzer, I., Keys, P. W., van der Ent, R. J., Savenije, H. H. G., and Gordon, L. J.: Remote land use impacts on river flows through atmospheric teleconnections, Hydrol. Earth Syst. Sci., 22, 4311–4328,, 2018. 

Wehrli, K., Guillod, B. P., Hauser, M., Leclair, M., and Seneviratne, S. I.: Identifying Key Driving Processes of Major Recent Heat Waves, J. Geophys. Res.-Atmos., 124, 11746–11765, 2019. 

Weiler, M. and McDonnell, J. J.: Conceptualizing lateral preferential flow and flow networks and simulating the effects on gauged and ungauged hillslopes, Water Resour. Res., 43, W03403,, 2007. 

Weiler, M. and Naef, F.: An experimental tracer study of the role of macropores in infiltration in grassland soils, Hydrol. Process., 17, 477–493, 2003. 

Wiekenkamp, I., Huisman, J. A., Bogena, H. R., Graf, A., Lin, H. S., Drüe, C., and Vereecken, H.: Changes in measured spatiotemporal patterns of hydrological response after partial deforestation in a headwater catchment, J. Hydrol., 542, 648–661, 2016. 

Wiekenkamp, I., Huisman, J. A., Bogena, H. R., and Vereecken, H.: Effects of deforestation on water flow in the vadose zone, Water, 12, 35,, 2020. 

Wilusz, D. C., Harman, C. J., and Ball, W. P.: Sensitivity of catchment transit times to rainfall variability under present and future climates, Water Resour. Res., 53, 10231–10256, 2017. 

Yadav, M., Wagener, T., and Gupta, H. V.: Regionalization of constraints on expected watershed response behavior for improved predictions in ungauged basins, Adv. Water Resour., 30, 1756–1774, 2007. 

Yang, Y., Donohue, R. J., and McVicar, T. R.: Global estimation of effective plant rooting depth: Implications for hydrological modeling, Water Resour. Res., 52, 8260–8276, 2016. 

Yilmaz, K. K., Gupta, H. V., and Wagener, T.: A process-based diagnostic approach to model evaluation: Application to the NWS distributed hydrologic model, Water Resour. Res., 44, W09417,, 2008. 

Zacharias, S., Bogena, H., Samaniego, L., Mauder, M., Fuß, R., Pütz, T., Frenzel, M., Schwank, M., Baessler, C., Butterbach-Bahl, K., Bens, O., Borg, E., Brauer, A., Dietrich, P., Hajnsek, I., Helle, G., Kiese, R., Kunstmann, H., Klotz, S., Munch, J., Papen, H., Priesack, E., Schmid, H., Steinbrecher, R., Rosenbaum, U., Teutsch, G., and Vereecken, H.: A network of terrestrial environmental observatories in Germany, Vadose Zone J., 10, 955–973, 2011. 

Zehe, E., Lee, H., and Sivapalan, M.: Dynamical process upscaling for deriving catchment scale state variables and constitutive relations for meso-scale process models, Hydrol. Earth Syst. Sci., 10, 981–996,, 2006. 

Zehe, E., Elsenbeer, H., Lindenmaier, F., Schulz, K., and Blöschl, G.: Patterns of predictability in hydrological threshold systems, Water Resour. Res., 43, W07434,, 2007. 

Zhang, M., Liu, N., Harper, R., Li, Q., Liu, K., Wei, X., and Liu, S.: A global review on hydrological responses to forest change across multiple spatial scales: Importance of scale, climate, forest type and hydrological regime, J. Hydrol., 546, 44–59, 2017. 

Zuber, A.: On the interpretation of tracer data in variable flow systems, J. Hydrol., 86, 45–57, 1986. 

Short summary
Deforestation affects how catchments store and release water. Here we found that deforestation in the study catchment led to a 20 % increase in mean runoff, while reducing the vegetation-accessible water storage from about 258 to 101 mm. As a consequence, fractions of young water in the stream increased by up to 25 % during wet periods. This implies that water and solutes are more rapidly routed to the stream, which can, after contamination, lead to increased contaminant peak concentrations.