Multi-decadal fluctuations in root zone storage capacity through vegetation adaptation to hydro-climatic variability have minor effects on the hydrological response in the Neckar River basin, Germany

Climatic variability can considerably affect catchment-scale root zone storage capacity (S_umax), which is a critical factor regulating latent heat fluxes and thus the moisture exchange between land and atmosphere as well as the hydrological response and biogeochemical processes in terrestrial hydrological systems. However, direct quantification of changes in S_umax over long time periods and the mechanistic drivers thereof at the catchment scale are missing so far. As a consequence, it remains unclear how climatic variability, such as precipitation regime or canopy water demand, affects S_umax and how fluctuations in S_umax may influence the partitioning of water fluxes and therefore also affect the hydrological response at the catchment scale. Based on long-term daily hydrological records (1953–2022) in the upper Neckar River basin in Germany, we found that variability in hydro-climatic conditions, with an aridity index I_A (i.e. $E_{\mathrm{P}}/P$) ranging between ∼ 0.9 and 1.1 over multiple consecutive 20-year periods, was accompanied by deviations ΔI_E between −0.02 and 0.01 from the expected I_E inferred from the long-term parametric Budyko curve. Similarly, fluctuations in S_umax, ranging between ∼ 95 and 115 mm or ∼ 20 %, were observed over the same time period. While uncorrelated with long-term mean precipitation and potential evaporation, it was shown that the magnitude of S_umax is controlled by the ratio of winter precipitation to summer precipitation (p < 0.05). In other words, S_umax in the study region does not depend on the overall wetness condition as for example expressed by I_A, but rather on how water supply by precipitation is distributed over the year. However, fluctuations in S_umax were found to be uncorrelated with observed changes in ΔI_E. Consequently, replacing a long-term average, time-invariant estimate of S_umax with a time-variable, dynamically changing formulation of that parameter in a hydrological model did not result in an improved representation of the long-term partitioning of water fluxes, as expressed by I_E (and fluctuations ΔI_E thereof), or in an improved representation of the shorter-term response dynamics.

Overall, this study provides quantitative mechanistic evidence that S_umax changes significantly over multiple decades, reflecting vegetation adaptation to climatic variability. However, this temporal evolution of S_umax cannot explain long-term fluctuations in the partitioning of water (and thus latent heat) fluxes as expressed by deviations ΔI_E from the parametric Budyko curve over multiple time periods with different climatic conditions. Similarly, it does not have any significant effects on shorter-term hydrological response characteristics of the upper Neckar catchment. This further suggests that accounting for the temporal evolution of S_umax with a time-variable formulation of that parameter in a hydrological model does not improve its ability to reproduce the hydrological response and may therefore be of minor importance for predicting the effects of a changing climate on the hydrological response in the study region over the next decades to come.

1 Introduction

Vegetation is a key component of the terrestrial hydrological cycle as it shapes the hydrological functioning of catchments by regulating the long-term average partitioning of water into drainage and evaporative fluxes (i.e. latent heat), frequently expressed as the runoff ratio C_r = $Q/P$ [–] and evaporative index I_E = $\mathrm{1}-Q/P$ = $E_{\mathrm{A}}/P$ [–], respectively. More specifically, vegetation transpiration, which despite uncertainties (Coenders-Gerrits et al., 2014) globally constitutes the largest fraction of all evaporative fluxes (Jasechko, 2018), is systematically controlled by the interplay between canopy water demand and water supply from the sub-surface (Donohue et al., 2007; Yang et al., 2016; Jaramillo et al., 2018; Mianabadi et al., 2019). To survive, vegetation needs continuous access to water stored in the sub-surface and access to roots to satisfy its canopy water demand. As a consequence, the vegetation present at any moment, and in particular its active root system, reflects its successful adaptation to the prevalent climatic conditions in a region (Laio et al., 2001; Schenk and Jackson, 2002; Rodriguez-Iturbe et al., 2007; Donohue et al., 2007; Gentine et al., 2012; Liancourt et al., 2012). Irrespective of the geometry, distribution, or structure of root systems, the maximum vegetation-accessible water storage volume in the unsaturated root zone of the sub-surface, hereafter referred to as the root zone storage capacity S_umax (mm), represents the hydrologically relevant information of root systems (Rodriguez-Iturbe et al., 2007; Nijzink et al., 2016a; Savenije and Hrachowitz, 2017; Gao et al., 2024). Therefore, S_umax directly reflects the hydrologically relevant information of root systems at the catchment scales. In response to a changing environment, these root systems of vegetation continuously adapt to allow the most efficient use of available energy and resources for survival. The driving factors of changes in root systems are thus also the driving factors of changes in S_umax, as S_umax inherently represents adaptations of the root system.

As a central part of hydrological systems, S_umax is also a critical parameter in hydrological and land surface models. As such, it can, in principle, be estimated as a function of root depth and the sub-surface pore volume between field capacity and permanent wilting point (Scrivner and Ruppert, 1970; Sivandran and Bras, 2012, 2013). However, these data are typically not available at sufficient levels of detail. Alternatively, catchment-scale S_umax can be estimated by three broad approaches. Firstly, it can be obtained by calibration as a parameter of a hydrological model (Nijzink et al., 2018; Bouaziz et al., 2020; Wang et al., 2023; Sriwongsitanon et al., 2023; Roberts et al., 2021; Bahremand and Hosseinalizadeh, 2022; Sadayappan et al., 2023; Tong et al., 2022). Secondly, based on optimality principles, there are some variables like transpiration, nitrogen uptake, or carbon gain that can be maximized to quantify S_umax (Guswa, 2008; McMurtrie et al., 2012; Sivandran and Bras, 2012; Yang et al., 2016; Speich et al., 2018). Thirdly, S_umax can be robustly estimated at the catchment scale directly from annual water deficits based on observed hydro-climatic data, i.e. precipitation and transpiration (e.g. Donohue et al., 2012; Gentine et al., 2012; Gao et al., 2014b; de Boer-Euser et al., 2016; Nijzink et al., 2016a; Dralle et al., 2021; McCormick et al., 2021; Hrachowitz et al., 2021; Stocker et al., 2023; van Oorschot et al., 2021, 2024). For applications of hydrological and land surface models, S_umax (or equivalent parameters) has, except for very few exceptions (Wagener et al., 2003; Merz et al., 2011; Bouaziz et al., 2022; Tempel et al., 2024), been assumed constant over time. As a major knowledge gap, it still remains unknown whether S_umax follows climatic variability and evolves over time, thereby reflecting vegetation adaptation to changing conditions.

In contrast, it is well understood that, due to the importance of vegetation for the hydrological functioning of terrestrial systems, anthropogenic land use management practices, such as deforestation and afforestation (Brown et al., 2005; Brath et al., 2006; Fenicia et al., 2009; Alila et al., 2009; Jaramillo et al., 2018; Teuling et al., 2019; Stephens et al., 2021; Hoek van Dijke et al., 2022; Ellison et al., 2024) or irrigation (e.g. AghaKouchak et al., 2015; Van Loon et al., 2016; Roodari et al., 2021), can induce major shifts in the partitioning between the major components of the terrestrial water and energy cycles and thus between I_E and C_r. Two detailed recent studies with well-documented information on deforestation in several experimental catchments could establish explicit mechanistic links between the reduction in S_umax by > 50 % following deforestation and decreases in I_E (and thus increases in C_r) from ∼ 0.4–0.5 to ∼ 0.1–0.3, depending on the catchment and the scale of deforestation (Nijzink et al., 2016a; Hrachowitz et al., 2021).

Mapping the shifts to lower I_E that followed these land conversions from forest to grassland- and rangeland-type vegetation as a function of the aridity index I_A = $E_{\mathrm{P}}/P$ in the Budyko framework (Schreiber, 1904; Ol'Dekop, 1911; Budyko, 1974) corresponds well to the results of previous studies that suggest that, across the world, catchments dominated by grass exhibit a consistently lower I_E at the same I_A than forest environments (e.g. Zhang et al., 2001, 2004; Oudin et al., 2008). These differences in the long-term average I_E are accounted for by parametric reformulations of the Budyko framework, such as the Tixeront–Fu equation (Tixeront, 1964; Fu, 1981). The lumped parameters (here: ω) of these expressions define long-term average catchment-specific positions in the I_A–I_E space. As such, the parameters are typically interpreted as encapsulating vegetation characteristics and all other hydro-climatic and physiographic properties of individual catchments besides I_A (e.g. Roderick and Farquhar, 2011; Berghuijs and Woods, 2016). A frequent assumption is that, with changes in climatic conditions, represented here by I_A, individual catchments can be expected to move to the associated new positions I_E, following their specific trajectories defined by ω (e.g. Zhou et al., 2015; Bouaziz et al., 2022). However, several studies have demonstrated that catchments in many regions worldwide experience deviations ΔI_E from their expected new I_E following a change in I_A (e.g. Jaramillo and Destouni, 2014; van der Velde et al., 2014; Jaramillo et al., 2018; Reaver et al., 2022; Ibrahim et al., 2024; Tempel et al., 2024).

Figure 1 (a) Elevation of the Neckar catchment with discharge and hydro-meteorological stations as well as the water sampling locations used in this study. (b) The spatial distribution of long-term mean annual precipitation in the upper Neckar catchment and the stratification into three distinct precipitation zones P1–P3 (black outline), together with the red outlines indicating three sub-catchments (C1: Rottweil; C2: Plochingen at the Fils River; C3: Horb) within the upper Neckar River basin. (c) Hydrological response units classified according to their land cover and topographic characteristics.

From the above, the following questions arise. (1) Following the notion that vegetation, i.e. individual plants but also the species composition of plant communities, continuously adapts to climatic conditions, does the catchment-scale root zone storage capacity S_umax change over multi-decadal timescales? (2) Do multi-decadal changes in the vegetation response, expressed by changes in S_umax, explain deviations ΔI_E from the expected I_E? (3) Does a time-variable representation of S_umax as a parameter in a hydrological model improve the models' ability to reproduce the hydrological response?

Building on previous studies, the objectives of this study in the upper Neckar River basin in Germany are therefore to provide an analysis of multi-decadal changes in S_umax as a result of changing climatic conditions over a 70-year period (1953–2022) and how this further affects hydrological dynamics. More specifically, we test the hypotheses that (1) S_umax changes significantly over multiple decades (reflecting vegetation adaptation to climatic variability), (2) changes in S_umax affect the long-term partitioning of drainage and evaporation and thus control deviations ΔI_E from the catchment-specific trajectory in the Budyko space, and (3) a time-dynamic implementation of S_umax improves the representation of streamflow in a hydrological model.

Table 1 Characteristics of the Neckar catchment in Germany.

Table 2 Mean annual precipitation P, potential evaporation E_P, temperature T_M, aridity index I_A, evaporative index I_E, parameter ω for the parametric Budyko framework, root zone storage capacity S_umax,WB, and S_umax,cal based on water balance data and hydrological model calibration for Scenario 1 (entire time period T: 1953–2022) and Scenario 2 (four sub-periods t₁: 1953–1972, t₂: 1973–1992, t₃: 1993–2012, and t₄: 2013–2022).

2 Study area

The upper Neckar River basin in south-western Germany covers an area of ∼ 4000 km², with the Black Forest on the western side and the Swabian Jura on the south-eastern side. The river basin has a varying topography, with its elevation ranging from 250 m at the outlet in the north to about 1019 m in the south (Fig. 1a, Table 1). Following the elevation gradient, the landscape is dominated by terrace-like elements and undulating hills (∼ 50 %) with wide valleys used as grasslands and croplands in the lower regions, in particular in the south-eastern parts of the upper Neckar River basin, and by increasingly steep and narrow forested valleys (∼ 40 %) towards the southern parts and the remaining area that includes flat grassland in valley bottoms (∼ 10 %) (Fig. 1c). Annual mean precipitation (P) over the whole river basin has a considerable spatial heterogeneity ranging from ∼ 700 mm yr⁻¹ in the lower parts of the basin to ∼ 1600 mm yr⁻¹ over the Black Forest, with catchment-average long-term mean precipitation (P) reaching ∼ 880 mm yr⁻¹ (Fig. 1b, Table 2). The catchment is characterized by a temperate, humid climate with warm, wet summers and cold, drier winters. Precipitation exhibits some seasonality with ∼ 500 mm yr⁻¹ for the summer months (from May to October) and ∼ 380 mm yr⁻¹ for the winter months (from November to April), respectively (Fig. 3). Although snow is in general not a major component of precipitation in the study region, snowmelt can have a significant influence during individual storm events. The long-term mean temperature is about 8.2 °C and the potential evaporation (E_P) is about ∼ 860 mm yr⁻¹, with an aridity index I_A = $E_{\mathrm{P}}/P$ ∼ 0.97 (Table 2).

Figure 2 (a) Time series of observed monthly precipitation P. (b) Daily cumulative evaporative fluxes for the entire time period (1953–2022), where the dark-brown line indicates the potential evaporation E_P and the yellow lines and light-orange-shaded areas show the actual evaporation E_A modelled using the best-fit parameter sets and the associated 5th and 95th percentiles of all feasible solutions calibrated based on the entire time period. (c) Monthly maximum values of the snow water equivalent (SWE) for the 1953–2022 time period, where the green line indicates the most balanced solution and the light-green shade indicates the 5th and 95th inter-quantile ranges obtained from all the Pareto optimal solutions calibrated based on the entire time period. (d) Observed (blue line) and modelled daily streamflow Q; the red line indicates the most balanced solution and the shaded area indicates the 5th and 95th percentiles of all feasible solutions calibrated based on the entire time period; and the different green background shades from lighter to darker indicate the sub-time periods from t₁ to t₄. Panels (e) and (f) zoom in to the observed and modelled streamflows for the selected wet year (light-grey shade, 1 January 1988–31 December 1988) and the dry year (grey shade, 1 January 2003–31 December 2003), respectively.

Figure 3 The annual and seasonal variability (i.e. winter and summer) of selected climatic indices, including annual averages of precipitation (P), potential evaporation (E_P), estimated snowmelt water, the number of precipitation days (P_num), and the precipitation intensity (P_intensity) for the four sub-time periods (t₁: 1953–1972, t₂: 1973–1992, t₃: 1993–2012, and t₄: 2013–2022). (a–e) The annual variability of selected climatic indices. (f–j) The seasonal variability of selected climatic indices in the winter periods. (k–o) The seasonal variability of selected climatic indices in the summer periods.

3 Data sets

3.1 Data

Daily hydro-meteorological data were available for the period 1 January 1953–31 December 2022 (Fig. 2). Daily precipitation and daily mean air temperature were obtained from stations operated by the German Weather Service (DWD). Precipitation was recorded at 15 stations and temperature measurements were available at 8 stations (Fig. 1) in or close to the study basin. Daily potential evaporation E_P (mm d⁻¹) was estimated using the Hargreaves equation based on the observed daily maximum and minimum temperatures, has been used in many previous studies, and has been shown to be a suitable method for modelling applications (Oudin et al., 2005). Daily mean discharge data for the period 1 January 1953–31 December 2022 at the outlet of the upper Neckar River basin at Plochingen station were provided by the German Federal Institute of Hydrology (BfG). In addition, data of daily mean discharge for the same time period from three sub-catchments within the upper Neckar River basin (Fig. 1) at the gauges Rottweil (C1; 422 km²), Plochingen at the Fils River (C2; 706 km²), and Horb (C3; 1111 km²) were available from the Environmental Agency of the Baden-Württemberg region (LUBW).

Based on the CORINE Land Cover data set of the upper Neckar River basin during the period 1 January 1953–31 December 2022 (https://land.copernicus.eu/en/products/global-dynamic-land-cover (last access: 27 August 2024), there is only very minor change (< 2 %) for all the defined land cover classes (Fig. 1c). The 90 m × 90 m digital elevation model of the study region (Fig. 1a) was obtained from the Hydrologic Derivatives for Modeling and Applications (HDMA) database of the United States Geological Survey (USGS) (Verdin, 2017; https://doi.org/10.5066/F7S180ZP) and used to derive the local topographic indices, including height above nearest drainage (HAND) and slope.

3.2 Data pre-processing

For the subsequent experiment (Sect. 4.2), the study basin was stratified into three zones P1–P3 that are characterized by a distinct long-term precipitation pattern (hereafter precipitation zones), following the approach described and implemented for the Neckar River basin by Wang et al. (2023). Briefly, Goovaerts (2000) and Lloyd (2005) showed that areal precipitation estimates informed by elevation data were often more accurate than those based on precipitation gauge observations alone (e.g. Hrachowitz and Weiler, 2011). Thus, to interpolate and estimate areal precipitation across the basin, we used co-kriging that considered elevation, as a preliminary analysis suggested lower errors. Finally, the individual precipitation estimates for each grid cell were used with k-means clustering to establish three clusters representing the three precipitation zones P1–P3 (see Fig. 1b).

To explore the fluctuations of S_umax over long timescales, we independently estimate the root zone storage capacity S_umax for four subsequent sub-periods of the available data record (t₁–t₄ in Table 2). To survive, root systems of vegetation and the associated vegetation-accessible water storage capacity S_umax respond to the ever-changing conditions of their environment. However, as these changes occur at the landscape scale and are mostly reflected in changes in the composition of plant species present in a specific spatial domain, fluctuations in S_umax largely occur at timescales that reflect the lifecycles of individual plants. Thus, periods of at least 20 years are required to reflect this and to allow for meaningful estimates of S_umax, as also demonstrated by many other studies (e.g. Gao et al., 2014b; Wang-Erlandsson et al., 2016; Singh et al., 2020; Stocker et al., 2023). We therefore had to strike a balance between the number of independent time periods (here: t₁–t₄) and the robustness of the associated S_umax estimates. We deliberately chose to emphasize fewer but longer time periods and thus rather reliable estimates of S_umax.

4 Methods

To test the hypotheses that the key vegetation parameter, i.e. the root zone storage capacity S_umax, evolves over multi-decadal timescales in response to changing hydro-climatic conditions and controls the deviations from expected trajectories in the Budyko space (thereby reflecting the need for time-variable implementations of S_umax as a parameter in a hydrological model), the following stepwise approach is designed. (1) Estimate the observed deviations ΔI_E from the long-term average expected I_E for four consecutive periods t₁–t₄ in the study period (Table 2). (2) Estimate the root zone storage capacity over the entire study period ($S_{\text{umax},\text{WB},\mathrm{T}}$) and the four individual periods t₁–t₄ ($S_{\text{umax},\text{WB},\mathrm{t}}$) based on observed water balance data. (3) Estimate the root zone storage capacity over the entire study period ($S_{\text{umax},\text{cal},\mathrm{T}}$) and the four individual periods t₁–t₄ ($S_{\text{umax},\text{cal},\mathrm{t}}$) by calibration of a hydrological model over the respective time periods to evaluate whether the changes in the calibrated S_umax,cal reflect changes in S_umax,WB directly estimated from water balance data from step (2). (4) Estimate the modelled deviations ($\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},{\mathrm{O}}^{\prime }}$) from the expected I_E using both a long-term average time-invariant $S_{\text{umax},\text{WB},\mathrm{T}}$ and an individual $S_{\text{umax},\text{WB},\mathrm{t}}$ for the four periods t₁–t₄ as model parameters.

4.1 Estimation of the temporal trajectory in the Budyko framework

Mapping aridity I_A = $E_{\mathrm{P}}/P$, where E_P is potential evaporation (mm d⁻¹) and P is precipitation (mm d⁻¹), against the evaporative index I_E = $E_{\mathrm{A}}/P$ = $\mathrm{1}-Q/P$, where E_A is actual evaporation (mm d⁻¹) and Q is streamflow (mm d⁻¹), the Budyko framework is an expression of the long-term average water balance for a catchment. It is based on the assumption of negligible storage change over the averaging time period, i.e. $\mathrm{d}S/\mathrm{d}t$ ∼ 0. As demonstrated by Han et al. (2020), this assumption holds for averaging periods ≥ 10 years for a large majority of catchments worldwide. Note that hereafter the term “evaporation” is used to refer to all combined evaporative fluxes, including interception and soil evaporation (E_i) as well as transpiration (E_T), following the terminology proposed by Savenije (2004) and Miralles et al. (2020).

The analysis in this paper is based on the parametric Tixeront–Fu formulation of the Budyko framework (Tixeront, 1964; Fu, 1981):

$$\begin{array}{}\text{(1)}& I_{\mathrm{E},\mathrm{T}}=\mathrm{1}+I_{\mathrm{A},\mathrm{T}}-{\left(\mathrm{1}+I_{\mathrm{A},\mathrm{T}}^{{\mathit{\omega }}_{\mathrm{T}}}\right)}^{\mathrm{1}-{\mathit{\omega }}_{\mathrm{T}}},\end{array}$$

where I_E,T is the observed evaporative index over a chosen averaging period T, I_A,T is the observed aridity index over the same period, and ω_T is the associated catchment-specific parameter that represents all combined catchment properties other than I_A.

In a theoretical catchment that only experiences changes in I_A and no changes in any other hydro-climatic and/or physical catchment characteristics, it can be assumed that ω_T remains constant over time, so that ${\mathit{\omega }}_{\mathrm{T}}={\mathit{\omega }}_{t_i}={\mathit{\omega }}_{t_{i+\mathrm{1}}}$. This implies that, following a disturbance ΔI_A in a subsequent time period t_i+1, the catchment stays on its specific curve defined by ω_T to a new $I_{\mathrm{E}{t}_{i+\mathrm{1}}}$. In such a case, ω_T can thus be used to predict future hydrological partitioning I_E. Based on this assumption, here we use the complete available hydro-climatic data record to estimate the long-term average ω_OT as a reference over the entire 1953–2022 study period. The sub-division into the four time periods t₁–t₄ as shown in Table 2 is then allowed to estimate the expected $I_{\mathrm{E},t_i^{\prime }}$ in the individual periods t₁–t₄: depending on the shift in the observed aridity index along the x axis in t_i ($\mathrm{\Delta }{I}_{\mathrm{A},\mathrm{T},t_i}$ = $I_{\mathrm{A},t_i}-I_{\mathrm{A},\mathrm{T}}$), a catchment will move along its parametric Budyko curve defined by ω_OT to a new expected position $I_{\mathrm{E},t_i^{\prime }}$ (Fig. 4).

Figure 4 Representation of the Budyko space, which shows the evaporative index (I_E = $\mathrm{1}-Q/P$) as a function of the aridity index (I_A = $E_{\mathrm{P}}/P$) and the water and energy limits. A catchment with the long-term mean aridity index I_A,T = $E_{\mathrm{P},\mathrm{T}}/P_{\mathrm{T}}$ and evaporative index I_E,T = $\mathrm{1}-Q_{\mathrm{T}}/P_{\mathrm{T}}$, which is derived from observed data of the entire time period, is plotted at location T on the parametric Budyko curve with ω_T (yellow line) as the baseline. Based on the observed sub-time-period data with the aridity index $I_{\mathrm{A},t_i}$ = $E_{\mathrm{P},t_i}/P_{t_i}$ and evaporative index $I_{\mathrm{E},t_i}$ = $\mathrm{1}-Q_{t_i}/P_{t_i}$, the same catchment is plotted at location t_i on the parametric Budyko curve with ${\mathit{\omega }}_{t_i}$ (green line). A movement in the Budyko space towards $t_{\mathrm{i}}^{\prime }$ along the ω_T curve is shown as a result of a change in the aridity index $I_{\mathrm{A},t_i}$ with the assumption that the long-term mean Budyko curve trajectory and the parameter ω are transferable across time for an individual catchment, which results in a significant deviation $\mathrm{\Delta }{I}_{\mathrm{E},t_i,t_i^{\prime }}$ between the observed evaporative index $I_{\mathrm{E},t_i}$ and the predicted evaporative index $I_{\mathrm{E},t_i^{\prime }}$.

Figure 5 Model structure of the distributed conceptual hydrological model, discretized into three parallel hydrological response units (HRUs), i.e. forest, grassland, and wetland in each precipitation zone from P1 to P3. The light-blue boxes indicate the hydrologically active individual storage volumes. The arrow lines indicate water fluxes and the model parameters are shown in red. All the symbols are described in Table S1 in the Supplement.

Based on the available data, we then estimate the individual observed $I_{\mathrm{E},t_i}$ together with the associated ${\mathit{\omega }}_{t_i}$ for each of the four time periods t₁–t₄ (Fig. 5). For each of the four time periods t₁–t₄, the deviation of $I_{\mathrm{E},t_i}$ from the catchment-specific expected $I_{\mathrm{E},t_i^{\prime }}$, corresponding to a shift from ω_T to ${\mathit{\omega }}_{t_i}\ne {\mathit{\omega }}_{\mathrm{T}}$, was then computed as $\mathrm{\Delta }{I}_{\mathrm{E},t_i,t_i^{\prime }}$ = $I_{\mathrm{E},t_i}-I_{\mathrm{E},t_i^{\prime }}$.

4.2 Estimation of the root zone storage capacity derived by the water balance method S_umax,WB

The root zone storage capacity is the maximum volume of water which can be held in soil pores of the unsaturated zone and which is accessible to root systems of vegetation for transpiration. Here the water balance method that is described in detail in previous papers (e.g. Gao et al., 2014b; Nijzink et al., 2016a; de Boer-Euser et al., 2016; Wang-Erlandsson et al., 2016; Bouaziz et al., 2020; Hrachowitz et al., 2021) is used to determine S_umax,WB. Briefly, S_umax,WB is estimated based on daily observations of precipitation (P), potential evaporation (E_P), and streamflow (Q). As a first step, effective precipitation P_e (mm d⁻¹) that enters the sub-surface is computed by accounting for interception evaporation by

$$\begin{array}{}\text{(2)}& P_{\mathrm{e}}\left(t\right)=P\left(t\right)-E_{\mathrm{i}}\left(t\right)-\mathrm{d}{S}_{\mathrm{i}}/\mathrm{d}t,\end{array}$$

where E_i (mm d⁻¹) is the daily interception evaporation and S_i (mm) is the interception storage. For each time step, E_i is determined by

$$\begin{array}{}\text{(3)}& E_{\mathrm{i}}\left(t\right)=\left\{\begin{array}{ll}{E}_{\mathrm{P}}\left(t\right)& \text{ if }{E}_{\mathrm{P}}\mathrm{d}t<S_{\mathrm{i}},\\ \frac{S_{\mathrm{i}}}{\mathrm{d}t}& \text{ if }{E}_{\mathrm{P}}\mathrm{d}t\ge S_{\mathrm{i}}.\end{array}\right.\end{array}$$

Then, the effective precipitation P_e (mm d⁻¹) is further estimated according to

$$\begin{array}{}\text{(4)}& P_{\mathrm{e}}\left(t\right)=\left\{\begin{array}{ll}\mathrm{0}& \text{ if }{S}_{\mathrm{i}}<S_{\text{imax}},\\ \frac{S_{\mathrm{i}}-S_{\text{imax}}}{\mathrm{d}t}& \text{ if }{S}_{\mathrm{i}}\ge S_{\text{imax}},\end{array}\right.\end{array}$$

where S_imax (mm) is the maximum interception storage. As S_umax is not very sensitive to the choice of S_imax, as previously shown by Hrachowitz et al. (2021) and Bouaziz et al. (2022), here we used a value of S_imax = 2 mm, which was previously also used by de Boer-Euser et al. (2016).

Hereafter, the long-term mean transpiration ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{r}}$ (mm d⁻¹) is estimated from the long-term water balance, with the assumption of no additional gains or losses:

$$\begin{array}{}\text{(5)}& {\stackrel{\mathrm{‾}}{E}}_{\mathrm{r}}={\stackrel{\mathrm{‾}}{P}}_{\mathrm{e}}-{\stackrel{\mathrm{‾}}{Q}}_{\mathrm{O}},\end{array}$$

where ${\stackrel{\mathrm{‾}}{P}}_{\mathrm{e}}$ (mm d⁻¹) is the long-term mean effective precipitation and ${\stackrel{\mathrm{‾}}{Q}}_{\mathrm{O}}$ (mm d⁻¹) is the long-term mean observed streamflow. Considering the seasonal fluctuation of energy input, the daily transpiration E_r (mm d⁻¹) is estimated by subsequently scaling the daily potential evaporation E_P (mm d⁻¹) minus the interception evaporation E_i (mm d⁻¹) (see Eqs. 2 and 3) by the long-term mean transpiration ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{r}}$ (mm d⁻¹), according to (Bouaziz et al., 2022; Hrachowitz et al., 2021)

$$\begin{array}{}\text{(6)}& E_{\mathrm{r}}\left(t\right)={\displaystyle \frac{{\stackrel{\mathrm{‾}}{E}}_{\mathrm{r}}}{\left({\stackrel{\mathrm{‾}}{E}}_{\mathrm{P}}-{\stackrel{\mathrm{‾}}{E}}_{\mathrm{i}}\right)}}(E_{\mathrm{P}}-E_{\mathrm{i}}),\end{array}$$

where ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{P}}$ (mm d⁻¹) is the long-term mean potential evaporation and ${\stackrel{\mathrm{‾}}{E}}_{\mathrm{i}}$ (mm d⁻¹) is the long-term mean interception evaporation.

From daily storage deficits S_rd,n(t) (mm) during dry periods, estimated as the cumulative sum of daily effective precipitation P_e (mm d⁻¹) minus transpiration E_r (mm d⁻¹), the maximum storage deficit S_rd,n of a specific year n is then computed as follows:

$$\begin{array}{}\text{(7)}& {\displaystyle }\begin{array}{rl}& S_{\text{rd},n}\left(t\right)\\ & =\left\{\begin{array}{ll}{\int }_{t_{\mathrm{0},w}}^{t_{\mathrm{0},d}}\left(P_{\mathrm{e}}\right(t)-E_{\mathrm{r}}(t\left)\right)\mathrm{d}t,& \text{ if }{\int }_{t_{\mathrm{0},w}}^{t_{\mathrm{0},d}}\left(P_{\mathrm{e}}\right(t)-E_{\mathrm{r}}(t\left)\right)\mathrm{d}t\le \mathrm{0},\\ \mathrm{0},& \text{ if }{\int }_{t_{\mathrm{0},w}}^{t_{\mathrm{0},d}}\left(P_{\mathrm{e}}\right(t)-E_{\mathrm{r}}(t\left)\right)\mathrm{d}t>\mathrm{0},\end{array}\right.\end{array}\text{(8)}& {\displaystyle }{S}_{\text{rd},n}=max\left(\mathrm{|}{S}_{\text{rd},n}\right(t\left)\mathrm{|}\right).\end{array}$$

where t is the time step (d), t_0,w is the day at the end of the wet period when the storage deficits are zero but P_e(t)−E_r(t) < 0, and t_0,d is the day when the storage deficits return to zero again after the beginning of the next wet period when the water supply exceeds the canopy water demand, i.e. (P_e(t)−E_r(t)) > 0. Any cumulative precipitation surplus is assumed to be drained from the root zone and is released from the system either directly as streamflow or via recharge of the groundwater.

Table 3 Flow signatures and the associated performance metrics used for model calibration and evaluation. The performance metrics include the Nash–Sutcliffe efficiency (NSE) and the relative error (RE).

The Gumbel extreme value distribution (Gumbel, 1941) was used previously by several other studies (Gao et al., 2014b; Nijzink et al., 2016a; de Boer-Euser et al., 2016; Bouaziz et al., 2020, 2022; Hrachowitz et al., 2021) to estimate the root zone storage capacity through the water balance approach. Based on fitting the Gumbel distribution to the maximum annual storage deficits for all n years during one of the four time periods t₁–t₄, the root zone storage capacity S_umax,WB can be derived from various return periods of the sequence of n maximum annual storage deficits S_rd. Previous studies suggested that vegetation develops root zone storage capacities large enough to survive in dry spells with return periods of ∼ 20–40 years (Gao et al., 2014b; de Boer-Euser et al., 2016; Wang-Erlandsson et al., 2016; Hrachowitz et al., 2021). Therefore, here we define S_umax,WB as the maximum storage deficit in a 40-year period so that $S_{\text{umax},\text{WB}}=S_{\text{rd},\mathrm{40}\text{years}}$.

Using the above water-balance-based method, we determine S_umax,WB for the entire study period 1953–2022 ($S_{\text{umax},\text{WB},\mathrm{T}}$) as well as individually for the four time periods t₁–t₄ ($S_{\text{umax},\text{WB},\mathrm{t}}$) to quantify potential fluctuations of the root zone storage capacity reflecting the adaptation to changing climatic conditions.

4.3 Hydrological model

4.3.1 Model architecture

Loosely based on the flexible DYNAMITE modular modelling framework (e.g. Hrachowitz et al., 2014), here we used a semi-distributed, process-based model that was previously successfully implemented and tested for the Neckar study basin (Wang et al., 2023) and for many other contrasting environments worldwide (e.g. Prenner et al., 2018; Hulsman et al., 2021a, b; Hanus et al., 2021; Bouaziz et al., 2022). Briefly, this hydrological model consists of three parallel hydrological response units (HRUs), i.e. forest, grassland/cropland, and wetland, which are linked through a common storage component representing the groundwater system (Fig. 5). The classification into the three HRUs was based on the HAND (Gharari et al., 2011) metric and land cover similar to previous studies (e.g. Gao et al., 2014a; Gharari et al., 2014; Nijzink et al., 2016b; Bouaziz et al., 2021). The model was further spatially discretized by stratification into 100 m elevation bands for a more detailed representation of the snow storage (S_snow) and was finally implemented in parallel, i.e. individually for each of the three precipitation zones P1–P3, to balance to a certain degree spatial differences in precipitation with computational requirements. Rain (P_rain) and meltwater (M_snow) from the different elevation zones were aggregated according to their associated spatial weights in each elevation zone as further input to the subsequent layers of the model in each HRU. The outflows from each HRU in each precipitation zone and finally the outflows from each precipitation zone were likewise aggregated according to their respective spatial weights to represent the catchment-aggregated outflows. While the three HRUs are characterized by distinct parameters that reflect their respective functioning, the parameters between the individual zones P1–P3 were, in the spirit of model parsimony, kept the same in what is elsewhere referred to as a distributed moisture-accounting approach (e.g. Ajami et al., 2004; Fenicia et al., 2008; Euser et al., 2015). Overall, the model consists of snow (S_snow), interception (S_i), the unsaturated root zone (S_u), and fast-responding (S_f) and slow-responding storage (S_s) components for each HRU and precipitation zone. The maximum storage volume in the unsaturated root zone component in each HRU is defined by the corresponding calibration parameters S_umax,F, S_umax,G, and S_umax,W, respectively. The catchment average S_umax,cal is then inferred by aggregating these parameters according to their spatial weights. Water can be released from unsaturated root zones as the combined soil evaporation and transpiration flux E_t (mm d⁻¹), which is a frequently applied way of representing vegetation water stress (e.g. Bouaziz et al., 2021; Gharari et al., 2013; Gao et al., 2014a). The equations of the model are provided in Table S1 in the Supplement, and more detailed descriptions of the model are provided by Wang et al. (2023) and the other earlier implementations referred to above.

4.3.2 Model calibration

The model was run with a daily time step and has 18 calibration parameters. Briefly, the model parameters were calibrated by using the Borg_MOEA algorithm (Borg Multi-objective evolutionary algorithm; Hadka and Reed, 2013) and were based on uniform prior distributions (Table S2 in the Supplement). To best reflect different aspects of the hydrograph (including high flows, low flows, and the partitioning of precipitation into runoff and evaporation), the parameters were calibrated using a multi-criteria approach that includes seven objective functions as performance metrics E_Q,n (Table 3). There are multiple ways of dealing with sets of Pareto front solutions, as described in detail by e.g. Efstratiadis and Koutsoyiannis (2010) or Gharari et al. (2013). We chose to use all solutions on the Pareto front to obtain a conservative estimate of uncertainty. The seven performance metrics were subsequently also combined into an overall performance metric based on the Euclidian distance (D_E), where D_E = 1 indicates a perfect fit. To find a somewhat balanced solution in the absence of more detailed information, all individual performance metrics were equally weighted here (e.g. Hrachowitz et al., 2021; Hulsman et al., 2021b; Wang et al., 2023):

$$\begin{array}{}\text{(9)}& D_{\mathrm{E}}=\mathrm{1}-\sqrt{{\displaystyle \frac{{\sum }_{n=\mathrm{1}}^N(\mathrm{1}-E_{\mathrm{Q},n}{)}^{\mathrm{2}}}{N}}},\end{array}$$

where N=7 is the number of performance metrics with respect to streamflow (E_Q,n). Note that the different units and thus different magnitudes of residuals in the individual performance metrics introduce some subjectivity in finding the most balanced overall solution according to D_E (Eq. 9). However, a preliminary sensitivity analysis with varying weights for the individual performance metrics in D_E suggested limited influence on the overall results and is thus not further reported here. In addition, the model was tested for its ability to represent spatial differences in the hydrological response by evaluating it against streamflow observations in three sub-catchments (C1–C3) of the upper Neckar catchment without further recalibration, where each one of the sub-catchments largely represents the hydrological response from one of the precipitation zones (Fig. 1).

The model is calibrated following two distinct calibration scenarios as indicated in Table 2. In the first scenario, the model and thus also S_umax,F, S_umax,G, and S_umax,W are calibrated over the full length of the 70-year study period from 1953 to 2022. This reflects the common assumption of a system that is stable over time. By extension, this also implies that the role of vegetation and thus S_umax does not change and that vegetation does not adapt to climatic variability. In the second scenario, individual calibration to the four time periods t₁–t₄ allowed us to estimate fluctuations in the parameters S_umax,F, S_umax,G, and S_umax,W between the time periods as an indicator of vegetation adaptation to changing climatic conditions.

Figure 6 (a) Green dots from the light dot to the dark dot indicate the observed positions for the four sub-time periods t₁–t₄. The black dot $t_{\mathrm{4}}^{\prime }$ indicates the expected location on the parametric Budyko curve with ω_T derived from the entire observed time period. We select time period t₄ as an example to present the modelled positions in the zoom-in plot (b). The grey dot t_4,mT indicates the modelled position based on Scenario 1, which is with $S_{\text{umax},\text{WB},\mathrm{T}}$, and the orange dot t_4,mt indicates the modelled position based on Scenario 2, which is with $S_{\text{umax},\text{WB},\mathrm{t}\mathrm{4}}$. $\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},{\mathrm{O}}^{\prime }}$ (black arrow) indicates the deviation between the modelled I_E,mT based on Scenario 1 and the expected $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$. $\mathrm{\Delta }{I}_{\mathrm{E},\text{mt},{\mathrm{O}}^{\prime }}$ (orange arrow) indicates the deviation between the modelled I_E,mt based on Scenario 2 and the expected $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$. $\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},\mathrm{O}}$ (grey arrow) indicates the deviation between the modelled I_E,mT based on Scenario 1 and the observed I_E,O. $\mathrm{\Delta }{I}_{\mathrm{E},\text{mt},\mathrm{O}}$ (green arrow) indicates the deviation between the modelled I_E,mt based on Scenario 2 and the observed I_E,O.

5 Results

5.1 Observed multi-decadal hydro-climatic variability

Based on the initial analysis of water balance data for the four sub-time periods, significant differences were observed in the variability of different hydro-climatic indicators over the 1953–2022 study period (Fig. 3). While periods t₁ and t₄ were characterized by rather low mean annual precipitation of ∼ 870 and 811 mm yr⁻¹, respectively, periods t₂ and t₃ were subject to, on average, higher precipitation of ∼ 911 mm yr⁻¹. While summer precipitation remained rather stable over the study period (Fig. 3f), the above was mostly caused by fluctuations in winter precipitation (Fig. 3k). In contrast, potential evaporation E_P has gradually increased by 7 % from 836 to 906 mm yr⁻¹ (Fig. 3b). Similarly reflecting increases in temperature (Fig. 3b), the annual snowpack and associated snowmelt have continuously decreased from around 98 mm yr⁻¹ to around 50 mm yr⁻¹ between t₁ and t₄ (Fig. 3c). A slight decrease in the number of days with precipitation from ∼ 264 to 251 (Fig. 3d), on average, mostly due to changes in the summer months (Fig. 3n), was accompanied by some rather limited variability in precipitation intensities (Fig. 3e), mostly during winter (Fig. 3j). Overall, the comparatively humid periods t₁–t₃ that were characterized by I_A fluctuating between 0.93 and 0.97 were followed by a markedly more arid period t₄ with I_A = 1.12 (Table 2, Fig. 6). In response to the multi-decadal variability in I_A, expressed as movement along the x axis in the Budyko framework, the catchment experienced I_E varying between 0.56 and 0.59 (Table 2, Fig. 6). However, this observed variability was somewhat lower than the variability I_E,ωT = 0.55–0.61 that would have been expected based on ω_T. This illustrates that the hydrological response did not consistently follow its long-term trajectory defined by ω_T. Instead, deviations $\mathrm{\Delta }{I}_{\mathrm{E},t_i}$ from the expected positions, and thus values of ${\mathit{\omega }}_{t_i}$ that are different to ω_T, were observed for the individual periods. More specifically, the deviations gradually decreased from $\mathrm{\Delta }{I}_{\mathrm{E},t_{\mathrm{1}}}$ = 0.01 in t₁ to $\mathrm{\Delta }{I}_{\mathrm{E},t_{\mathrm{4}}}$ = −0.02 in t₄ (Fig. 6). This systematic shift towards lower (more negative) $\mathrm{\Delta }{I}_{\mathrm{E},t_i}$ and thus also lower ${\mathit{\omega }}_{t_i}$ indicates that, at the same I_A, a smaller fraction of precipitation is released as evaporation, i.e. I_E, now than at the start of the 70-year study period. Although the magnitude of deviations of $\mathrm{\Delta }{I}_{\mathrm{E},t_i}\le$ ± 0.02 remains rather minor, similar to what was recently reported elsewhere (Ibrahim et al., 2024; Tempel et al., 2024), their systematic shift in particular to one direction implies that changes in the system other than I_A have a visible effect on the hydrological response pattern.

Figure 7 The time series of storage deficits as calculated by Eq. (7) for the entire time period T (1953–2022) and for the four sub-time periods (green shades from light to dark for the time period from t₁ to t₄) (a, c, e, g, and i). The maximum annual deficits are indicated by the dots. Estimation of S_umax as the storage deficit associated with a 40-year return period using the Gumble extreme value distribution for different time periods (b, d, f, h, and j). The orange crosses indicate the values of S_umax for the different time periods.

5.2 Root zone storage capacity S_umax,WB estimated from water balance data

As the baseline of our study, the annual maximum storage deficits fluctuate between 97 mm in 2022 and 16 mm in 1970 (Fig. 7a). Assuming an adaptation to dry spells with 40-year return periods, the root zone storage capacity over the entire 1953–2022 study period (Scenario 1) was estimated to be $S_{\text{umax},\text{WB},\mathrm{T}}$ = 105 mm (Table 2, Fig. 7b). In the next step, the storage deficits and the associated root zone storage capacity for each period from t₁ to t₄ were estimated (Scenario 2). $S_{\text{umax},\text{WB},\mathrm{t}\mathrm{1}}$ and $S_{\text{umax},\text{WB},\mathrm{t}\mathrm{3}}$ for periods t₁ and t₃, respectively, are estimated at the same value of 95 mm. In contrast, and somewhat counterintuitively, the highest value over the study period is found in the wettest period (t₂) and reaches $S_{\text{umax},\text{WB},\mathrm{t}\mathrm{2}}$ = 115 mm, while the driest period (t₄) is characterized by $S_{\text{umax},\text{WB},\mathrm{t}\mathrm{4}}$ = 100 mm (Table 2, Fig. 7c–j). These patterns suggest that S_umax,WB did vary by ∼ 20 mm, equivalent to ∼ 20 % throughout the 1953–2022 period. In contrast to $\mathrm{\Delta }{I}_{\mathrm{E},t_i}$, which was characterized by a systematic shift towards more negative deviations over time, no evidence was found of a systematic, one-directional shift in S_umax,WB. Instead, S_umax,WB evolved following a somewhat cyclic pattern.

Figure 8 Selected model performance metrics in the entire time period 1 January 1953–31 December 2022 of the upper Neckar River basin against the model performance in uncalibrated sub-catchment-based (C1: Rottweil; C2: Plochingen at the Fils River; C3: Horb) parameter sets derived from the calibration for the entire time period. The dots indicate all the Pareto optimal solutions in the multi-objective model performance space. The shades from dark to light indicate the overall model performance based on the Euclidean distance D_E, with the black solutions representing the overall better solutions (i.e. larger D_E).

5.3 Root zone storage capacity S_umax,cal estimated as a calibration parameter

5.3.1 Model calibration for 1953–2022 (Scenario 1)

The model parameter sets determined as feasible after calibration over the entire 1953–2022 study period in Scenario 1 reproduce the main features of the hydrological response (Fig. 2d). More specifically, the modelled hydrographs in particular describe the timing of high flows well while somewhat underestimating the flow peaks for the best-performing model in terms of D_E (Eq. 9). The low flows and the shapes of the recessions are generally captured well (NSE_log Q = 0.67). Crucially, the model also reproduces the other observed streamflow signatures well, e.g. the flow duration curves (NSE_log FDC = 0.96), the autocorrelation function (NSE_AC = 0.99), and the long-term and seasonal runoff coefficients (${\mathrm{NSE}}_{C_{\mathrm{r}}}$ = 0.90, ${\mathrm{RE}}_{C_{\mathrm{r}},\text{summer}}$ = 0.83, and ${\mathrm{RE}}_{C_{\mathrm{r}},\text{winter}}$ = 0.91). The latter further implies that the modelled long-term actual evaporative fluxes E_A (Fig. 2b) and thus I_E,ωT are, on average, consistent with the observed ones, which can be seen in Fig. 6. The model, calibrated on the overall response of the upper Neckar River basin, also exhibited considerable skill in representing spatial differences in the hydrological response by reproducing observed streamflow in the three sub-catchments (C1–C3) similarly well (Fig. 8) without any further recalibration. The overall model skill to mimic the hydrological response corresponds well to a similar implementation of the model in the greater study area by Wang et al. (2023). A detailed list of the performance metrics is provided in Table S4 in the Supplement.

Figure 9 (a) S_umax values derived from the water balance method and the hydrological model for different time periods. The yellow boxes from light to dark indicate the range of S_umax,cal for the sub-time period from t₁ to t₄ and the entire time period T based on the corresponding parameter sets derived from the model. Yellow dots indicate the corresponding $S_{\text{umax},\text{cal},\text{bal}}$ based on the most balanced solution, and green dots indicate the corresponding S_umax,WB derived from the water balance method. (b) Values of $S_{\text{umax},\text{cal},\text{bal}}$ against S_umax,WB. The yellow–green dots from light to dark indicate the values of S_umax for the sub-time period from t₁ to t₄ and the entire time period T.

The model calibration resulted in pronounced differences in the root zone storage capacity parameters for three individual landscape classes. While for forest-dominated land it was estimated at S_umax,F = 158 mm for the best-performing model (5th and 95th percentiles of all feasible solutions: 138–168 mm), it reached S_umax,G = 95 mm (5th and 95th percentiles: 71–123 mm) for grassland/cropland and S_umax,W = 61 mm for wetland (5th and 95th percentiles: 49–68 mm), which reflects differences in the vegetation type and position in the landscape (Fan et al., 2017). Remarkably, the catchment root zone storage capacity, estimated by aggregating the individual values according to their areal fractions, at S_umax,cal = 116 mm (5th and 95th percentiles: 99–130 mm; Fig. 9a) came very close to the estimate of S_umax,WB = 105 mm that is directly derived from the water balance method without any calibration, as described in Sect. 5.2.

5.3.2 Model calibration for the individual periods t₁–t₄ (Scenario 2)

The model parameter sets obtained from the individual calibration for each period from t₁ to t₄ reproduce the hydrographs of the corresponding periods as well or slightly better than when using the long-term average parameters from Scenario 1 (see the detailed performance metrics in Table S4). In particular, the runoff coefficients of ${\mathrm{NSE}}_{C_{\mathrm{r}}}$ ∼ 0.86–0.91, ${\mathrm{RE}}_{C_{\mathrm{r}},\text{summer}}$ ∼ 0.84–0.90, and ${\mathrm{RE}}_{C_{\mathrm{r}},\text{winter}}$ ∼ 0.88–0.92 could be mimicked rather well. Similarly, the daily dynamics with NSE_log Q ∼ 0.63–0.72 for the best-performing model of each period and most other hydrological signatures could be reproduced marginally better.

The individual calibration over each period t₁–t₄ resulted in associated differences in the catchment-scale root zone storage capacity of each period. Based on the best-performing models, the calibrated values varied between low values for t₁ and t₂, with $S_{\text{umax},\text{cal},t_{\mathrm{1}}}$ = 98 mm and $S_{\text{umax},\text{cal},t_{\mathrm{3}}}$ = 99 mm, and higher values for the other two periods, with $S_{\text{umax},\text{cal},t_{\mathrm{2}}}$ = 122 mm and $S_{\text{umax},\text{cal},t_{\mathrm{4}}}$ = 107 mm (Table 2, Fig. 9a). The magnitudes of $S_{\text{umax},\text{cal},t_i}$ obtained by calibration in the individual time periods t₁–t₄ are, with a difference of 5 mm (∼ 5 %), on average very close to the $S_{\text{umax},\text{WB},t_i}$ estimated on the basis of the water balance for the same periods. Perhaps even more notably, the temporal evolution of $S_{\text{umax},\text{cal},t_i}$ and $S_{\text{umax},\text{WB},t_i}$ follows the same sequence over time (R² = 0.95, p = 0.05; Fig. 9b).

5.4 Effect of S_umax on temporal fluctuation in the trajectories of the Budyko curve

The deviations $\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ between the expected evaporative index $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$ and the observed evaporative index I_E,O for all the periods t₁–t₄ become gradually more negative from t₁ ($\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ = 0.013) to t₄ ($\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ = −0.020), which is consistent with decreases in ${\mathit{\omega }}_{t_i}$ and downward shifts of the associated parametric Budyko curves over time as described in Sect. 5.1. These systematic reductions in $\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ over the 70-year study period are not reflected in the fluctuations of root zone storage capacities, irrespective of how they were estimated, i.e. $S_{\text{umax},\text{WB},t_i}$ (p = 0.85) or $S_{\text{umax},\text{cal},t_i}$ (p = 0.96), as illustrated in Fig. 10.

Figure 10 Relationships between the deviations $\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ and the values of S_umax,WB and S_umax,cal for the four sub-time periods (t₁–t₄) which are derived from the water balance method (green circles) and the hydrological model calibration (yellow circles), respectively.

Figure 11 (a) The grey boxes (I_E,mT) indicate the modelled evaporative index based on all the Pareto front solutions for the four sub-time periods based on Scenario 1 with a stationary $S_{\text{umax},\text{WB},\mathrm{T}}$ and grey dots based on the most balanced solution based on Scenario 1. The green circles from light to dark in panels (a) and (b) indicate the observed evaporative index for the four sub-time periods from t₁ to t₄. (b) The orange circles (I_E,mt) indicate the modelled evaporative index based on all the Pareto front solutions for the four sub-time periods (from lighter to darker shades) based on Scenario 2 with time-variant $S_{\text{umax},\text{WB},t_i}$ and the orange circles based on the most balanced solution for Scenario 2. The black boxes in panel (c) indicate the deviations $\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},{\mathrm{O}}^{\prime }}$ = I_E,mT − $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$ (Δ i) based on all the Pareto front solutions for the four sub-time periods and the dark-grey circles based on the most balanced solution based on Scenario 1. The orange boxes in panel (c) indicate the deviations $\mathrm{\Delta }{I}_{\mathrm{E},\text{mt},{\mathrm{O}}^{\prime }}$ = I_E,mt − $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$ (Δ iii) based on all the Pareto front solutions for the four sub-time periods and the orange circles based on the most balanced solution for Scenario 2. The grey dots indicate the deviation $\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ (Δ ii) between the observed I_E,O for each sub-time period and the corresponding expected $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$. The light-grey boxes in panel (d) indicate the deviations $\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},\mathrm{O}}$ = I_E,mT − I_E,O (Δ j) based on all the Pareto front solutions for the four sub-time periods and the grey circles based on the most balanced solution. The green boxes in panel (d) indicate the deviations $\mathrm{\Delta }{I}_{\mathrm{E},\text{mt},\mathrm{O}}$ = I_E,mt − I_E,O (Δ jj) based on all the Pareto front solutions for the four sub-time periods and the green circles based on the most balanced solution.

The above is further corroborated by comparing the modelled I_E from Scenarios 1 and 2 for each period t₁–t₄. More specifically, in Fig. 11c it can be seen that Scenario 1, based on a time-invariant, long-term average $S_{\text{umax},\text{WB},\mathrm{T}}$ obtained over the entire 1953–2022 period generates deviations $\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},{\mathrm{O}}^{\prime }}$ from the expected long-term average $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$ for each period t₁–t₄. In this case, the modelled I_E does not follow the expected $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$. However, nor does it follow the sequence of increasingly negative deviations from 0.013 in t₁ to −0.020 in t₄ as observed in reality ($\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$). Instead, $\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},{\mathrm{O}}^{\prime }}$ remains negative for all the time periods and fluctuates between $\mathrm{\Delta }{I}_{\mathrm{E},\text{mT},{\mathrm{O}}^{\prime }}$ = −0.005 and −0.029 (white boxplots in Fig. 11c). Replacing the time-invariant $S_{\text{umax},\text{WB},\mathrm{T}}$ with an individual $S_{\text{umax},\text{WB},t_i}$ for each period t₁–t₄ in Scenario 2 accounts for the different effects of vegetation in the individual periods. If S_umax controlled the observed deviations from the expected $I_{\mathrm{E},{\mathrm{O}}^{\prime }}$, Scenario 2 would generate estimates of $\mathrm{\Delta }{I}_{\mathrm{E},{\text{mt}}_i,{\mathrm{O}}^{\prime }}$ that are closer to the observed ones than those of Scenario 1. However, no evidence was found of that: the deviations $\mathrm{\Delta }{I}_{\mathrm{E},\text{mt},{\mathrm{O}}^{\prime }}$ obtained by Scenario 2 with time-variable S_umax for each period t₁–t₄ are largely indistinguishable (orange boxplots in Fig. 11c) from those generated by Scenario 1 with time-invariant S_umax. As a consequence, the evaporative index I_E modelled with time-variable $S_{\text{umax},\text{WB},t_i}$ is not found to be closer to the observed I_E for Scenario 2 than for Scenario 1. On the contrary, the deviations $\mathrm{\Delta }{I}_{\mathrm{E},m,\mathrm{O}}$ from the observed I_E obtained from the time-invariant Scenario 1 are in most time periods, albeit only slightly, less pronounced (Fig. 11d).

Figure 12 The relative errors of observed and modelled high (Q^5th), median (Q^50th), and low (Q^95th) flow quantiles and the mean Q for different time periods based on two scenarios. The grey shades and the corresponding dots indicate the relative errors based on Scenario 1 with all the Pareto front solutions, and the grey circles indicate the most balanced solution. The corresponding values for Scenario 2 are shown in green.

5.5 Effect of S_umax on streamflow

Corresponding to the above findings, there is no significant difference in the modelled average streamflow between Scenario 1, using the long-term average S_umax, and Scenario 2, using individual S_umax values for each time period (Fig. 12d). While the model for both scenarios consistently and similarly underestimates high flows (Q^5th, Fig. 12a) by ∼ 10 %, it overestimates median flows by ∼ 15 % with both time-invariant and time-variable S_umax for all time periods (Fig. 12b). Interestingly, the low flows are overpredicted by ∼ 10 %–20 % in the first two periods, while they are underpredicted by up to ∼ 20 % in the later periods in both scenarios (Fig. 12c). In addition, it was also found that using a time-variable S_umax in Scenario 2 did not have any discernible effect on the seasonal flow pattern (not shown). The fact that both scenarios generate similar estimates over different flow percentiles and, in particular, that the time-variable Scenario 2 reflects the same systematic shift in the ability of the model to reproduce low flows as Scenario 1 suggests, together with the very minor effects of the time-variable S_umax in Scenario 2 on the model performance metrics, that the adaptation of S_umax to changing climatic conditions does not significantly affect the average hydrological response pattern in the Neckar River basin.

6 Discussion

6.1 Multi-decadal changes in root zone storage capacity S_umax

As Gao et al. (2024) suggested, considering the terrestrial ecosystem structure can improve our understanding of hydrological processes and how the ecosystem can survive and be developed. It is valuable to explore how ecosystems adapt to climatic variability (which is reflected in fluctuations in S_umax) and how this affects the long-term partitioning of drainage, evaporation, and hydrological response. This study is the first to systematically and explicitly quantify how the root zone storage capacity S_umax changes with changing climatic conditions over time. The values of root zone storage capacity, estimated from both water balance data and as a model calibration parameter, indeed show significant and corresponding fluctuations over multiple decades, varying by up to ±20 %. At S_umax ∼ 95–115 mm, the overall estimated magnitudes fall well within the range of long-term average values reported previously for the greater region (e.g. Bouaziz et al., 2021; Hrachowitz et al., 2021; Tempel et al., 2024) and other temperate, humid environments (e.g. Kleidon, 2004; Gao et al., 2014b; de Boer et al., 2016; Wang-Erlandsson et al., 2016; Stocker et al., 2023; van Oorschot et al., 2024).

The values of S_umax obtained from both methods are very similar and within an error margin of only ∼ 5 %. In addition, they both follow a comparable change over time. Together, this lends support to the underlying assumption that this temporal evolution of S_umax may indeed be a fingerprint of vegetation adaptation to changing climatic conditions. More specifically, as S_umax,WB is explicitly based on the estimates of transpiration E_r (Eq. 9), it could be plausibly argued that during specific years only more water is used for E_r but that the size of the water storage volume accessible for roots may not necessarily change. In that case, changes in S_umax,WB would not reflect actual changes in the active root system but only in how much water was used by them. In contrast, S_umax,cal inferred as a calibration parameter of a hydrological model regulates not only transpiration, but, critically, also the generation of streamflow. If therefore the active root system in reality did not change and fluctuations in S_umax,WB were a mere artifact of changes in water uptake from a fixed-size volume instead of an actual change in maximum vegetation-accessible sub-surface water volumes, fluctuations in S_umax,cal would not mirror those of S_umax,WB, and the use of S_umax,WB in the hydrological model would, due to the non-linear character of the flow generation function in the model (Eq. S20 in Table S1), lead to misrepresentations of streamflow dynamics. However, no deteriorations of the model performance with a changing S_umax were found here. Even more, the fact that S_umax,WB and S_umax,cal are characterized by very similar magnitudes and fluctuations adds further evidence that their evolution over time is a manifestation of vegetation adapting its active root system to changing climatic conditions.

Figure 13 Relationships between the temporal fluctuation of the deviations ($\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$, black dots in Fig. 10a), S_umax,WB, and climate indices, including precipitation (P), potential evaporation (E_P), the aridity index (I_A), the estimated snowmelt water (Snowmelt), the number of precipitation days (P_num), and the precipitation intensity (P_intensity) for the four sub-time periods. Relationships between the temporal fluctuation of the deviations ($\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$, black dots in Fig. 10a) and the climate indices for the four sub-time periods are shown in the first two rows in panels (a–l) (i.e. (a–f) $\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ vs. annual climatic indices and (g–l) $\mathrm{\Delta }{I}_{\mathrm{E},\mathrm{O},{\mathrm{O}}^{\prime }}$ vs. seasonal climatic indices). Relationships between the temporal fluctuation of S_umax,WB derived from the water balance method and climate indices for the four sub-time periods are shown in the last two rows in panels (m–x) (i.e. (m–r) S_umax,WB vs. annual climatic indices and (s–x) S_umax,WB vs. seasonal climatic indices).

Several previous studies in similar environments found that the root zone storage capacity S_umax can decrease by 50 % or more after deforestation and that these changes not only cause reductions in I_E by −0.2 or more, which reflect changes in ω and thus the overall functioning of the system, but also influence hydrological dynamics on short timescales, such as the magnitudes of flow peaks (Nijzink et al., 2016a; Hrachowitz et al., 2021). In contrast to the above studies, the ±20 % fluctuations of S_umax here did not lead to similarly marked shifts in I_E or ω. This is further corroborated by an analysis of different variables as potential controls on S_umax and ΔI_E as shown in Fig. 13. Fluctuations of S_umax can to a large extent be attributed to the variability in the ratio of winter precipitation to summer precipitation (Fig. 13s) as a simplified metric for precipitation seasonality. This comes as no surprise, as the computation of S_umax,WB is explicitly based on the seasonal water deficit (Eq. 7). It merely shows that, when more precipitation falls in summer, at a time when evaporative demand is highest, the lower S_umax needs to provide vegetation with sufficient and continuous access to water for continuous vegetation transpiration. All the other tested variables do not exert any major influence on S_umax in the study region. Conversely, it was found that the deviations ΔI_E are largely independent of the seasonality of precipitation (Fig. 13g). Instead, increases in summer E_P are correlated with decreases in ΔI_E (Fig. 13h) and thus with a reduction in E_T. The observed systematic shift towards a more negative ΔI_E, which indicates proportionally less evaporation, thus coincides with the gradually increasing summer E_P over time. This points towards different controls on ΔI_E than on S_umax and the potential role of increased vegetation water stress in summer as the main driver of ΔI_E. Thus, while there is compelling evidence of fluctuations in S_umax, the above illustrates that these changes cannot explain the observed deviations from the expected long-term Budyko trajectory in the study region.

It is also important to note that the temporal fluctuations of both S_umax and ΔI_E can be subject to uncertainties. In spite of the findings reported by Han et al. (2020) that, for most river catchments worldwide, $\mathrm{d}S/\mathrm{d}t$ ∼ 0 holds over averaging periods similar to the ones used here (t₁–t₄), this assumption may not completely hold in the study region. In relation to this, we did not consider potential effects of unobserved groundwater import or export on the long-term water balance either (Bouaziz et al., 2018).

As only < 2 % of the study area experienced documented land use change over the 1953–2022 period and no major reservoirs are present upstream of the study basin outlet, here we interpret fluctuations in S_umax as being a reflection of adaptations of root systems to changing hydro-climatic conditions. However, some of the fluctuations may be a consequence of land management practices not quantified by available gridded land cover products such as CORINE, including forest thinning (Hrachowitz et al., 2021) or rejuvenation (Teuling and Hoek van Dijke, 2020). In addition, although here we mainly attribute changes in S_umax to changes in root systems, these may be complemented by additional effects of changes in vegetation water use due to feedbacks with increases in atmospheric CO₂ (e.g. Berghuijs et al., 2017; Jaramillo et al., 2018).

6.2 Effect of changing S_umax on the representation of streamflow in a model

Reflecting their lack of explanatory power for the changes in ΔI_E, our results correspondingly indicate that signatures of both annual flow, such as the average Q^5th, Q^50th, or Q^95th, and seasonal flow are not reproduced better by the hydrological model when replacing a time-invariant, long-term average S_umax with a temporally dynamic S_umax. Overall, these results are in contrast to previous studies that quantified the effect of a time-variable S_umax parameter following deforestation. For example, Nijzink et al. (2016a) reported that adjusting the parameter S_umax to a lower value improves a model's ability to reproduce streamflow after deforestation. These findings were strongly supported by Hrachowitz et al. (2021), who found that post-deforestation model recalibration resulted in a lower S_umax and a significantly better performance compared to using parameters from pre-deforestation calibration. However, our results are also different from those reported by Duethmann et al. (2020), who found that accounting for vegetation dynamics in a model in the form of changing surface resistances to vegetation improved the long-term performance of the model. Similarly, Bouaziz et al. (2022) estimated the future S_umax based on projected future hydro-climatic conditions. In a somewhat more humid environment, they found that an estimated ∼ 25 % future increase in S_umax from ∼ 170 to 226 mm may lead to reductions in the mean and maximum annual Q of ∼ 5 %. More pronounced effects were reported on the intra-annual timescale, with reductions in autumn and winter Q of up to ∼ 15 %. This was accompanied by increases of up to ∼ 15 % in summer evaporation and decreases of 10 % in winter groundwater levels. Irrespective of the additional uncertainties in their study introduced by future projections, the much less pronounced effects we found in our analysis are most likely a consequence of the lower absolute magnitude of S_umax that remains below 115 mm in the study region. These lower S_umax values reflect lower storage requirements in summer, due to a precipitation pattern in the Neckar River basin that is more evenly spread throughout the year. In other words, the fact that here ∼ 55 %–60 % of the annual precipitation falls in summer (Fig. 3f and k) when it is needed most by vegetation due to high E_P removes the need for a larger S_umax as a water storage buffer to allow vegetation to survive. However, the lower the magnitude of S_umax, the more frequently storage deficits can be overcome by even rather small rainstorms and the less water is (or needs to be) stored. Thus, even if the relative changes are similar between Bouaziz et al. (2022) and Tempel et al. (2024) in a somewhat more humid environment and our study, abundant summer precipitation causes absolute S_umax fluctuations of less than ±20 mm over time in the Neckar catchment. This in turn limits the influence of the changes on the hydrological response, which has wider implications for the use of models in the Neckar River basin and potentially in other temperate regions with similar hydro-climatic characteristics. More specifically, it has been argued that a changing climate will affect the properties of terrestrial hydrological systems (e.g. Stevens et al., 2020). As these properties are represented by typically time-invariant parameters in hydrological or land surface models, accounting for changing system properties with time-variable formulations of parameters may facilitate more reliable predictions. For many model parameters such a time-variable formulation to estimate their future values is not trivial due to frequently insufficient data and a general lack of mechanistic understanding of the underlying processes. The estimation of S_umax and its temporal evolution based on observed historical or projected future water balance data opens up an opportunity to estimate its magnitude under future conditions for use in models. However, in contrast to the findings in other regions (e.g. Merz et al., 2011) and as discussed above, adapting S_umax to changing conditions in the Neckar River basin does not lead to improved modelled representations of the hydrological response. It is therefore plausible to assume that the use of a time-variable parameter S_umax does not substantially improve future predictions and is thus not necessarily required for at least the next few decades to come and that the use of a long-term average S_umax, obtained either by calibration or based on the water balance, is sufficient in the Neckar River basin and in hydro-climatically similar regions. This in itself is already an interesting finding as it gives modellers process-based evidence that the use of a time-invariant S_umax as a model parameter will also be sufficient for meaningful hydrological predictions in the near future in such an environment. However, it can also be expected that, in more arid regions with less summer precipitation and generally higher S_umax (see e.g. Gao et al., 2014b; Stocker et al., 2023), changes in S_umax will play a more prominent role.

7 Conclusions

The catchment-scale root zone storage capacity (S_umax) is a critical factor reflecting the moisture exchange between land and atmosphere as well as the hydrological response in terrestrial hydrological systems. However, as a major knowledge gap, it is unclear whether S_umax at the catchment scale evolves over time, reflecting vegetation adaptation to changing climatic conditions. As a consequence, it also remains unclear how potential changes in S_umax may affect the partitioning of water fluxes and, as a consequence, the catchment-scale hydrological response. In this study, for the upper Neckar catchment, based on long-term daily hydrological data (1953–2022), we quantify and analyse how S_umax dynamically evolves over multiple decades, reflecting vegetation adaptation to climate variability and the effects on the hydrological system.

The main findings of our analysis are the following:

1. S_umax has fluctuated by ±20 % between 95 and 115 mm, in response to climatic variability over the 70-year study period.

2. Estimates of S_umax obtained from both methods, i.e. based on water balance data and model calibration parameters, respectively, were, with differences of ∼ 5 %, highly consistent with each other and correlated in time (R² = 0.95, p = 0.05). Findings (1) and (2) support the hypothesis that S_umax, even in temperate, humid climates such as in the Neckar River basin, significantly changes over multiple decades, reflecting vegetation adaptation to climatic variability.

3. The estimated fluctuations of S_umax were inconsistent with the temporal sequence of observed deviations ΔI_E ∼ ± 0.02 from the expected I_E over the study period (R² = 0.02, p = 0.85).

4. As a consequence, replacing a long-term average, time-invariant parameter S_umax in a hydrological model with a time-variable formulation of S_umax does not lead to a better representation of the observed ΔI_E. Based on (3) and (4), the hypothesis that S_umax affects the long-term partitioning of drainage and evaporation and thus controls deviations ΔI_E from the catchment-specific trajectory in the Budyko space needs to be rejected for the Neckar River basin.

5. Replacing the time-invariant S_umax with a time-variable S_umax in the hydrological model leads to only very minor improvements of the model to reproduce streamflow dynamics. The hypothesis that a time-dynamic implementation of S_umax improves the representation of streamflow in the hydrological model therefore also needs to be rejected for the Neckar River basin.

Overall, our study is the first to systematically document the temporal evolution of S_umax, and although limited to the Neckar River basin, it provides clear quantitative evidence that S_umax can significantly change over multiple decades, reflecting vegetation adaptation to climate variability. However, these changes do not cause deviations from the long-term average Budyko curve under changing climatic conditions. This implies that the temporal evolution of S_umax does not control variation in the partitioning of water fluxes and has no significant effects on fundamental hydrological response characteristics of the upper Neckar River basin. As the use of time-variable S_umax over the 70-year study period does not improve the performance of the hydrological model, it can plausibly be assumed that in the study region the use of time-invariant S_umax as a model parameter will be sufficient for meaningful predictions over at least the next few decades.