The importance of ecosystem adaptation on hydrological model predictions in response to climate change

To predict future hydrological behavior in a changing world, often use is made of models calibrated on past observations, disregarding that hydrological systems, hence model parameters, will change as well. Yet, ecosystems likely adjust their root-zone storage capacity, which is the key parameter of any hydrological system, in response to climate change. In addition, other species might become dominant, both under natural and anthropogenic influence. In this study, we propose a top-down approach, which directly uses projected climate data to estimate how vegetation adapts its root-zone storage capacity at the 5 catchment scale in response to changes in magnitude and seasonality of hydro-climatic variables. Additionally, the Budyko characteristics of different dominant ecosystems in sub-catchments are used to simulate the hydrological behavior of potential future land-use change, in a space-for-time exchange. We hypothesize that changes in the predicted hydrological response as a result of 2K global warming are more pronounced when explicitly considering changes in the sub-surface system properties induced by vegetation adaptation to changing environmental conditions. We test our hypothesis in the Meuse basin in four 10 scenarios designed to predict the hydrological response to 2K global warming in comparison to current-day conditions using a process-based hydrological model with (a) a stationary system, i.e. no changes in the root-zone storage capacity of vegetation and historical land use, (b) an adapted root-zone storage capacity in response to a changing climate but with historical land use, and (c,d) an adapted root-zone storage capacity considering two hypothetical changes in land use from coniferous plantations/agriculture towards broadleaved forest and vice-versa. We found that the larger root-zone storage capacities (+34 %) 15 in response to a more pronounced seasonality with drier summers under 2K global warming strongly alter seasonal patterns of the hydrological response, with an overall increase in mean annual evaporation (+4 %), a decrease in recharge (-6 %) and a decrease in streamflow (-7 %), compared to predictions with a stationary system. By integrating a time-dynamic representation of changing vegetation properties in hydrological models, we make a potential step towards more reliable hydrological predictions under change. 20 1 https://doi.org/10.5194/hess-2021-204 Preprint. Discussion started: 20 April 2021 c © Author(s) 2021. CC BY 4.0 License.


Introduction
Hydrological models are required to provide robust short-term hydrological forecasts and long-term predictions of the impact of natural and human-induced change on the hydrological response. Common practice is to predict the future using a hydrological model calibrated to the past (Vaze et al., 2010;Blöschl and Montanari, 2010;Peel and Blöschl, 2011;Coron, 2013;Seibert and van Meerveld, 2016). For the near future, it seems acceptable to assume no fundamental change in the hydrological system, 25 although we know that ecosystems, the manager of the hydrological system, have the capacity to adapt to climatic change (Savenije and Hrachowitz, 2017). For longer term predictions, it is therefore not correct to assume an unchanged system within a changing world. This raises the question on the robustness of hydrological predictions, especially in the context of climate change (Coron et al., 2012;Stephens et al., 2019).
For example, Merz et al. (2011) clearly shows the non-stationarity of hydrological model parameters when calibrating 273 30 Austrian catchments in subsequent 5-years periods between 1976 and 2006. Being the core parameter of any hydrological system, Merz et al. (2011) report almost a doubling of the root-zone storage capacity and this gradual increase is assumed to be related to changing climatic conditions, such as increased evaporation and drier conditions in the more recent years. The temporal variability of model parameters could also be attributed to uncertainties in input and model structure or inadequate calibration strategies. However, the observed trends in model parameters are also likely to reflect transient catchment conditions 35 over the historical period.
Under continued global warming, precipitation and temperature extremes are expected to further increase and the hydrological cycle is likely to further accelerate (Allen et al., 2010;Kovats et al., 2014;Stephens et al., 2021). In addition, natural land cover change and anthropogenic activities of land-cover change and land-use management can substantially alter a catchment's water balance (Brown et al., 2005;Wagener, 2007;Fenicia et al., 2009;Jaramillo and Destouni, 2014;Nijzink et al., 2016a;Géographique et Forestière, 2019). In contrast, only 44 % of the 18 th century Walloon forests of Belgium has remained from the original broadleaved forest, the rest being cleared for agriculture on high fertility soils in the North West (30 %) or converted to coniferous plantations (Scots pine, Norway spruce and Douglas-fir) on the poor soils of the Ardennes (26 %, Kervyn et al., 2018). The status of "old growth" forest does not exclude human disturbances, but assumes a relatively limited impact.
Soils are less disturbed and their structure and biochemical composition have been preserved for several centuries. This fa-130 vors a high degree of biodiversity, which is a key element for the resilience of forest ecosystems to perturbations. In contrast, recent short-rotation plantations lack many of these characteristics. Particularly thick canopy plantations, such as the spruce and Douglas-fir, significantly alter the typical biodiversity of forests. Additionally, relatively higher evaporation water use is expected in these recent, short-rotation exotic plantations in comparison to older, more natural forests (Fenicia et al., 2009;Stephens et al., 2021).

Observed historical E-OBS climate data
The E-OBS dataset (v20.0e) includes daily precipitation, temperature and radiation fields for the period 1980-2018 at a 25 km 2 resolution (Cornes et al., 2018). The data are based on station data collated by the European Climate Assessment Dataset (ECA&D) initiative. Temperature is downscaled using the digital elevation model and a fixed lapse rate of 0.0065 • C m −1 .

140
Potential evaporation is estimated using the Makkink formula (Hooghart and Lablans, 1988). There is a relatively large underestimation of precipitation (> 20 %) in the E-OBS dataset in the center of the basin when compared to an operational dataset, which is based on local precipitation data provided by the Service Public de Wallonie for the period 2005-2017 (Bouaziz et al., 2020). A monthly bias-correction factor is applied to improve the consistency between both datasets (Sect. S1 of the Supplement).

Simulated historical and 2K climate data
To study the impact of 2K global warming on the hydrological response of the Meuse basin, we use climate simulations of the historical period 1979-2018 and a 2K global warming simulation, provided by the Royal Netherlands Meteorological Institute (KNMI). The simulations are generated with the regional climate model KNMI-RACMO2 (van Meijgaard et al., 2008) at 12 km x 12 km resolution. RACMO2 uses the land surface scheme HTESSEL (Balsamo et al., 2009), which employs four soil layers with a total depth of 2.9 m. Each land-grid cell includes separate tiles for high and low vegetation (16 vegetation types), bare soil, snow and intercepted water, for which the energy and water balances are solved individually.
The historical simulation uses ERA5 reanalysis data (Hersbach et al., 2020) as initial-and lateral boundary conditions. The 2K simulation is a so-called pseudo-global warming (PGW) simulation (e.g. Schär et al., 1996;Attema et al., 2014;Prein et al., 2017;Brogli et al., 2019), which is an alternative method to generate high-resolution climate change information. Instead of 155 downscaling global climate model (GCM) projections, the historical period is re-simulated, but set against a warmer climate background by adding perturbations to the ERA5 initial-and boundary conditions. The perturbations represent the change in the mean climate state in a globally 2K warmer world, derived from a large initial condition GCM ensemble (Aalbers et al., 2018). The method minimizes biases in the mean climate state of the historical simulation, guaranties a realistic atmospheric circulation under both historical and 'future' conditions and increases the signal-to-noise ratio of the climate change response.

160
A full description of the dataset is provided in Aalbers et al. (2021).

Streamflow
Streamflow data is available for 35 catchments nested within the Meuse basin upstream of Borgharen for the period 2005(Fig. 1, Service Public de Wallonie, 2018Banque Hydro, 2018). The streamflow at Borgharen is a constructed time series which sums the observed streamflow of the Meuse at St Pieter and of the Albert Canal at Kanne to represent the total flow from 165 the tributaries before part of it is extracted in the Albert Canal (de Wit et al., 2007).

Methods
To quantify the importance of reflecting ecosystem adaptation in hydrological models in response to climate change, the following stepwise approach is designed: (1) estimate the long-term runoff coefficient in a 2K warmer world from movements in the Budyko space as a result of a shift in aridity index and a potential shift in dominant land-use from broadleaved forests to 170 coniferous plantation and agriculture and vice-versa by trading space-for-time in the Meuse basin; (2) estimate how the rootzone storage capacity adapts in response to a more pronounced seasonality with drier summers and changing dominant land use using the observed historical and estimated long-term runoff coefficient in a 2K warmer world with potential changes in land use; (3) calibrate a hydrological model with observed historical E-OBS climate data to represent current-day hydrological conditions; (4) test if the historical climate data simulated by the regional climate model leads to a plausible representation of 175 current-day hydrological conditions; (5) run the hydrological model with the 2K climate data in four scenarios describing (a) a stationary system with historical root-zone storage capacity and historical land-use, (b) an adapted root-zone storage capacity in response to a changing climate but a historical land-use, (c,d) an adapted root-zone storage capacity and a shift in dominant land-use; and finally (6) compare the change in hydrological response between 2K and historical conditions for these four scenarios. The long-term partitioning of precipitation (P ) into evaporation (E A ) and streamflow (Q) is mainly controlled by the long-term aridity index (ratio of potential evaporation over precipitation, E P /P ), according to the Budyko hypothesis. To account for additional factors that influence the long-term hydrological partitioning, Fu (1981) introduces a parameter ω to encapsulate the 185 combined influences of climate, soils, vegetation and topography (Equation 1).
We solve Equation 1 to determine the value of ω for each of the 35 catchments of the Meuse basin for observed historical conditions for the period 2005 to 2017 (ω obs ), using the meteorological E-OBS data (P obs and E P,obs ) and observed streamflow (Q obs ). Assuming only a change in long-term mean climate conditions, i.e. aridity index, a catchment will move along its 190 ω obs -parameterized curve from its original position (p 1 ) to a new position (p 2 ) due to the horizontal shift in aridity index (∆E P /∆P , Fig. 2a). Here, we use the simulated historical and 2K climate data to determine how the change in potential evaporation (∆E P = E P,2K − E P,hist ) and precipitation (∆P = P 2K − P hist ) lead to a change in aridity index (Equation 2) and therefore in actual evaporation (∆E A = E A,2K − E A,hist ) and streamflow (∆Q = Q 2K − Q obs ), using Equation 1.
However, land cover and vegetation are likely to also change in response to a changing climate, introducing an additional vertical shift (∆ω) toward a position p 3 on a different ω change curve (Fig. 2a). A downward vertical shift from ω obs to ω change indicates less water use for evaporation, as opposed to an upward shift for higher evaporative water use. These vertical shifts in ω values represent changes in drivers other than aridity index, including e.g. land cover, tree species, forest age, biomass growth and water use efficiency (Jaramillo et al., 2018).

200
To test the sensitivity of the hydrological response to a change in ω in addition to a change in aridity index, we consider two scenarios. The catchments with relatively high percentages of broadleaved forests (25-38% as in the French part of the basin) receive the ω values of catchments with relatively low percentages of broadleaved forests (1-12% as mainly in the Belgian Ardennes) and vice-versa (Fig. 1b). We denote ω broadleaved for the catchments with relatively high percentages of broadleaved forests and ω coniferous for the catchments where broadleaved forests were largely converted to coniferous plantations or agri-205 culture. These scenarios are meant as a sensitivity analysis in the spirit of trading space-for-time (Singh et al., 2011) to evaluate the effect of potential future land-use management on the overall water balance.
When converting broadleaved forest to coniferous plantations, we expect an increase in water use for evaporation and therefore a vertical upward shift in ω values, as opposed to a downward shift when converting coniferous plantations to more natural broadleaved forests. The described vertical and horizontal movements in the Budyko space are used to estimate the projected  The root-zone storage capacity represents the maximum volume of water which can be held in pores of unsaturated soil and which is accessible to roots of vegetation for transpiration. It is a key element controlling the hydrological response of 220 hydrological systems. The long-term partitioning of precipitation into streamflow and evaporation in a changed climate can only match expectations as estimated from movements in the Budyko space (Sect. 4.1.1) if we consider that vegetation has adapted its root-zone storage capacity to offset hydro-climatic seasonality, by creating a buffer large enough to overcome dry spells (Gentine et al., 2012;Donohue et al., 2012;Gao et al., 2014). This is the main assumption underlying the water-balance method to estimate the root-zone storage capacity at the catchment scale (Gao et al., 2014;Nijzink et al., 2016a;de Boer-Euser 225 et al., 2016;Wang-Erlandsson et al., 2016;Bouaziz et al., 2020;Hrachowitz et al., 2020).
The water-balance method requires daily time series of precipitation, potential evaporation and a long-term runoff coefficient to estimate transpiration, as it depletes the root-zone storage during dry spells. Annual water deficits (S R,def ) stored in the rootzone of vegetation to fulfill canopy water demand for transpiration are estimated on a daily time step as the cumulative sum of daily effective precipitation (P E ) minus transpiration (E R ).

230
First, effective precipitation, i.e. the amount of precipitation that reaches the soil after interception evaporation (E I ), is estimated by solving the water balance of a canopy storage (S I ) with maximum interception storage capacity (I max , here taken as 2.0 mm), according to Equation 4.
Next, the long-term transpiration E R is estimated from the long-term water balance, using mean annual streamflow and 235 effective precipitation (Q and P E , all in mm yr −1 , Equation 5), assuming negligible changes in storage and intercatchment groundwater flows (catchments where E A = P − Q < E P ). The long-term transpiration E R is subsequently scaled to daily transpiration estimates E R , using the daily signal of potential evaporation minus interception evaporation, according to Equation 6 (Nijzink et al., 2016a;Bouaziz et al., 2020).
The maximum annual storage deficits can then be derived from the cumulative difference of effective precipitation (P E ) and transpiration (E R ), assuming an "infinite" storage, according to Equation 7 and illustrated in Fig. 2b. For each year, S R,def represents the amount of water accessible to the roots of vegetation for transpiration during a dry period. Storage deficits are assumed to be zero at the end of the wet period (T 0 , here April) and increase when transpiration exceeds effective precipitation 245 during dry periods, until they become zero again (T 1 ) when excess precipitation is assumed to drain away as direct runoff or recharge.
By fitting the annual maximum storage deficits to the extreme value distribution of Gumbel, the root-zone storage capacity at the catchment scale S R,max can be derived for various return periods. Previous studies used a return period of 20 years 250 for forested areas, meaning that forests develop root systems to survive droughts with a return period of ∼20 years (Nijzink et al., 2016a;de Boer-Euser et al., 2016;Hrachowitz et al., 2020). The root-zone storage capacity of cropland and grasslands is assumed to correspond to deficits with lower return periods of ∼2 years (Wang-Erlandsson et al., 2016). It should be noted that the methodology assumes that vegetation taps its water from the unsaturated zone and not from the groundwater.

255
Using the above described methodology, we determine several sets of S R,max values for each of the 35 catchments of the Meuse basin to represent the historical and adapted root-zone storage capacity in response to a changing climate and changing/historical land-use, using historical climate observations (E-OBS) and the historical and 2K climate simulations (Table 1).
S R,max,A : Historical root-zone storage capacity from historical land use and observed historical climate 260 The first set is S R,max,A , which represents the historical meteorological and land-use conditions, derived from observed historical E-OBS data (P obs , E P,obs for the period 1980-2018) and observed streamflow data (Q obs for the period 2005-2017).
S R,max,A is used as parameter for three model runs, each forced with a different dataset: historical E-OBS observations, simulated historical and 2K climate data (Sect. 4.4).
In this study, we assume that the observed E-OBS historical climate data is the best available estimate of current-day climate 265 conditions and use this data to estimate historical root-zone storage capacities S R,max,A and to calibrate the hydrological model (Sect. 4.3). The simulated historical climate data is also required to enable a fair comparison with the simulated 2K climate data, as they are both generated with the regional climate model. Despite potential biases in the climate model simulations compared to the observed historical data (here, E-OBS), we do not apply a formal bias-correction of the climate data which may alter the relations between variables in climate models (Ehret et al., 2012). Instead, we force the hydrological model with 270 the native simulated historical climate data, but use the previously determined S R,max,A parameter. An alternative approach would have been to estimate the root-zone storage capacities using the simulated historical climate data, to directly correct for potential biases in the climate data in the estimation of the root-zone storage capacity parameter but with the downside of affecting spatial patterns across catchments. For comparison, this analysis is performed in Sect. S2 of the Supplement.
S R,max,B : Adapted root-zone storage capacity from historical land use and 2K climate 275 We then estimate the root-zone storage capacity S R,max,B based on the 2K climate and historical land use to reflect vegetation adaptation to changing climatic conditions such as differences in seasonality, aridity index (Equation 2) and the resulting runoff coefficient (Equation 3), but under the assumption that the vegetation cover remains unchanged. To account for differences in the observed and simulated historical climate data, S R,max,B is determined by imposing the difference in storage deficits derived from the 2K and historical climate simulations (S R,def,2K − S R,def,hist ) on the observed storage deficit derived with 280 E-OBS data S R,def,obs , as shown in Table 1.
S R,max,C : Adapted root-zone storage capacity from land use conversion to broadleaved forest and 2K climate Subsequently, the root-zone storage capacity is estimated for the 2K climate under two land-use change scenarios, considering that if climate changes, a different vegetation cover might become dominant under natural and anthropogenic influence (Table 1). Making use of a space-for-time exchange, we connect the spatially variable ω parameter of the Budyko curve to dif- The limited number of samples is due to the high computational resources required to run the distributed model. However, our aim is not to find the "optimal" parameter set, but rather to retain an ensemble of plausible parameter sets based on a multiobjective calibration strategy (Hulsman et al., 2019). To best reflect different aspects of the hydrograph, including high flows, 320 low flows and medium-term partitioning of precipitation into drainage and evaporation, parameter sets are selected based on their ability to simultaneously and adequately represent four objective functions, including the Nash-Sutcliffe efficiencies of streamflow, the logarithm of streamflow and, monthly runoff coefficients as well as the Kling-Gupta efficiency of streamflow.
Only parameter sets that exceed a performance threshold of 0.9 for each metric are retained as feasible. The root-zone storage capacity parameter S R,max,A is a fixed parameter, which is derived from annual maximum storage deficits with a return period The performance of the calibrated model for the ensemble of retained parameter sets is also evaluated when the model is forced with the simulated historical climate data, using S R,max,A for the root-zone storage capacity parameter. This is the reference historical run against which the relative effect of 2K global warming is evaluated for different scenarios ( Fig. 4 and Sect. 4.4). In addition, in Sect. S2 of the Supplement, we evaluate the performance of the calibrated model forced with the simulated historical climate data but with a root-zone storage capacity parameter derived directly from this data. While this 335 alternative approach enables to correct for potential biases in the simulated historical climate data directly in the estimation of the root-zone storage capacity parameter, it may also affect the spatial patterns of this parameter across catchments.

Hydrological change evaluation
We then force the calibrated wflow_FLEX-Topo model for the ensemble of retained parameter sets with the 2K climate data in four scenarios each using a different root-zone storage capacity parameter to represent either stationary or adapted condi- In scenario 2K A (Fig. 4), we assume an unchanged land use and that vegetation has not adapted its root-zone storage capacity to the aridity and seasonality of the 2K climate. This scenario implies stationarity of model parameters by using S R,max,A in both the historical and 2K runs, a common assumption of many climate change impact assessment studies (Booij, 2005;de Wit 350 et al., 2007;Prudhomme et al., 2014;Hakala et al., 2019;Brunner et al., 2019;Gao et al., 2020;Rottler et al., 2020). This is the benchmark scenario against which we compare the hydrological response considering non-stationarity of the system, as in the following three scenarios.
Scenario 2K B : Historical land use and 2K climate adapted root-zone storage capacity (S R,max,B ) In scenario 2K B (Fig. 4), we again assume an unchanged land use (ω obs ). However, we assume that vegetation has adapted its 355 root-zone storage capacity to the aridity and seasonality of the 2K climate conditions by selecting S R,max,B as parameter for the 2K model run, while the historical S R,max,A is used as parameter in the historical run.
Scenario 2K C : Land-use conversion to broadleaved forest and 2K climate adapted root-zone storage capacity (S R,max,C ) In scenario 2K C (Fig. 4), we adapt the root-zone storage capacity to the changing aridity index and seasonality of the 2K In scenario 2K D (Fig. 4), the approach of scenario 2K C is repeated. However, now the broadleaved forest in the French catchments are assumed to be converted to coniferous plantations or agriculture as in the Belgian Ardennes. The parameter S R,max,D is used in the model run forced with the 2K climate, while S R,max,A is used in the historical run. The root-zone storage capacity S R,max,A derived with observed historical E-OBS climate data and observed streamflow is estimated at values of 101 ± 17 mm and 169 ± 24 mm across all study catchments for a 2 year and 20 year return period, respectively (Fig. 5c).
If instead the simulated historical climate data is used to derive the root-zone storage capacity, this results in slightly higher values with 110 ± 18 mm and 180 ± 28 mm for 2 and 20 year return periods, respectively. This overestimation of about +7 % 385 is due to the higher precipitation (on average +9 %) in the simulated historical climate data compared to the observed E-OBS historical data, which leads to relatively lower runoff coefficients and therefore larger evaporative indices and storage deficits in the water balance calculation of the root-zone storage capacity (Sect. S2 of the Supplement).

S R,max,B from historical land use and 2K climate
The adapted root-zone storage capacity S R,max,B , in response to changing climate conditions and an unchanged land use, 390 strongly increases with respect to historical conditions (S R,max,A ) with estimated values of 129 ± 18 mm (+28 %) and 227 ± 27 mm (+34 %) for return periods of 2 year and 20 year, respectively (Fig. 5c). This strong increase is explained by larger storage deficits during summer due to an increase of about +10 % in summer potential evaporation in the 2K climate and, therefore, a more pronounced seasonality (Fig. 2b). In contrast, the change in aridity index between the historical and 2K climate simulations is relatively small with a median of +0.01 across all study catchments. This can be explained by a 395 simultaneous increase in mean annual precipitation (+5 %) and potential evaporation (+7 %) on average over the basin area in the 2K climate compared to the simulated historical climate data. This increase in precipitation mostly occurs during the winter half year (Nov-Apr). In contrast, there is a relatively large variability in precipitation change in summer, characterized by years with wetter and drier summers. year, respectively (Fig. 5c). This corresponds to an increase of +9 % and +7 % for both return periods in comparison with S R,max,B , which does not consider additional land-use changes.

410
The difference in root-zone storage capacity between the 2K and historical climate simulations as a result of a changing climate (aridity and seasonality) is larger (+58 mm or +34 % for a return period of 20 years) than the difference between rootzone storage capacity for a changing climate and additional changes in land use (-8 mm or -3 % for S R,max,C and +16 mm or +7 % for S R,max,D ). This indicates that with the assumed land-use change in scenarios 2K C and 2K D , the strong increase in water demand during summer as a result of a more pronounced seasonality has greater impact on the estimation of the 415 root-zone storage capacity than a change in ω values. However, note that land use is changed in only part of the catchment for both land use change scenarios and that it is plausible to assume that more pronounced changes in land use will reinforce the observed effects.

Model evaluation (historical period)
5.2.1 Model forced with observed historical climate data 420 The ensemble of parameter sets retained as feasible after calibration mimics the observed hydrograph at Borgharen relatively well for the evaluation period (Fig. 6a). Also the seasonal streamflow regime is relatively well reproduced by the model, except for a slight underestimation of about -9 % in the first half year (Fig. 6b). The four objective functions show a relatively similar performance during calibration and evaluation with median values of approximately 0.93 and 0.78 at Borgharen and for the ensemble of nested catchments of the Meuse, respectively (Fig. 7a,b).

Model forced with simulated historical climate data
When the calibrated model is instead forced with the simulated historical climate data, peaks are slightly overestimated in comparison to the model run forced with the observed historical E-OBS data (Fig. 6c). This is due to the on average +9 % overestimation of precipitation in the simulated historical climate data compared to the observed historical E-OBS climate data. This precipitation overestimation results in an overestimation of about +12 % of modeled mean monthly streamflow 430 during the wettest months (Fig. 6d). The streamflow model performance at Borgharen slightly decreases when the simulated historical climate data is used instead of E-OBS, but median values across the ensemble of feasible parameter sets are still above 0.77 for each of the objective functions (Fig. 7c). Although a decrease in model performance is found in a few nested catchments, the performance in the ensemble of nested catchments of the Meuse remains relatively high with median values of around 0.67 (Fig. 7d). The results of the model run forced with the simulated historical data and with the root-zone storage 435 capacity parameter derived directly from this data show a relatively similar behavior, as further detailed in Sect. S2 of the Supplement.
The calibrated model forced with the simulated historical climate data shows a plausible behavior with respect to observed streamflow and is also close to the performance achieved with the observed historical E-OBS climate data. This is important because the effect of the 2K climate on the hydrological response is evaluated with respect to the model run forced with the 440 simulated historical climate data, as they are both generated with the regional climate model. Therefore, the relatively high model performance in the evaluation period enable us to use the retained parameter sets from the calibration with E-OBS data for the subsequent analyses with the simulated historical and 2K climate data. In the 2K A scenario, representing a stationary system with identical parameters in the historical and 2K climate, runoff coefficients are projected to increase with a median of +3 %, the evaporative index (E A /P ) decreases with a median of -2 % and mean annual streamflow increases with a median of +7 %. Maximum annual streamflow is also projected to increase with a median of about +5 %, while the median change in annual minimum of 7-days mean streamflow remains close to zero. The median annual deficit volume below the 90 th percentile historical streamflow increases with +10 %, as shown in Fig. 8.

450
Streamflow is projected to increase from December until August with +8 % and decrease between September and November with -7 %. In the months where evaporation demand exceeds precipitation, the root-zone soil moisture decreases, with a maximum of -22 % in September. Actual evaporation increases throughout the year with +3 % except in July and August (-4 %) when the availability of water in the root-zone of vegetation is not sufficient to supply canopy water demand. Recharge to the groundwater storage increases with approximately +5 % in all months except November, as shown in Fig. 9. Changes are substantially different in the 2K B scenario which considers that the root-zone storage capacity of vegetation has adapted to the change in aridity and seasonality of the 2K climate. Runoff coefficients are instead projected to decrease with a median of -2 %, while the evaporative index increases with a median of +2 % and the median change of mean annual streamflow is close to zero (Fig. 8). Also the median change of, both, annual maximum streamflow and minimum 7-days mean 460 streamflow remain close to zero. However, there is a substantial increase of +38 % in median annual deficit volume below the 90 th percentile historical streamflow. This result suggests that while the minimum streamflow remains relatively similar, the length of the low flow period strongly increases if we consider that the root-zone storage capacity has adapted to the 2K climate ( Fig. 8).
Seasonal changes indicate a decrease in streamflow of -19% between September and November, which is longer and more 465 pronounced than in the 2K A scenario (Fig. 9). The root-zone soil moisture increases throughout the year with an average of +34 % due to the larger root-zone storage capacities. Actual evaporation is no longer reduced as a result of moisture stress in the root-zone and strongly increases with approximately +7 % from May to October to supply canopy water demand. The increase in evaporation during summer strongly reduces the groundwater recharge with -5 % from October to February (Fig. 9).
5.3.3 Scenario 2K C : Land-use conversion to broadleaved forest and adapted root-zone storage capacity (S R,max,C )

470
The predicted hydrological response in the 2K C scenario is very similar to the response of the 2K B scenario, despite considering additional changes in the root-zone storage capacity as a result of a land-use conversion from coniferous plantations and agriculture to broadleaved forest (Figs. 8 and 9). This is in line with the limited differences in root-zone storage capacities of approximately +3 % between both scenarios (Sect. 5.1).

capacity (S R,max,D )
In contrast, the change in hydrological response is most pronounced for the scenario S R,max,D , which considers land use conversion of the broadleaved forests in the French part of the basin to coniferous plantations and agriculture (Figs. 8 and 9).
Runoff coefficients decrease with a median of -4 %, while the evaporative index increases with a median value of +4 % and mean annual streamflow decreases with a median of -2 %. If the median change in streamflow extremes remains relatively 480 close to zero, there is a strong increase of +54 % in the median annual deficit volume, suggesting a strong increase in the length of the low flow period (Fig. 8).

Streamflow decreases from August to January with an average of -23 % and evaporation strongly increases from May to
October with an average of +9 %. This increased evaporation during summer further reduces recharge from October to February with -7 % (Fig. 9). In comparison with the hydrological response of scenario 2K B , the additional land-use conversion in 485 scenario 2K D results in relatively similar patterns of change but with an additional +2 % increase in evaporation, -2 % decrease in streamflow and -2 % decrease in recharge on average throughout the year.

Stationary versus adaptive ecosystems
There is a difference of -7 % in the change of mean annual streamflow between the scenarios 2K B , 2K C , 2K D with adaptive ecosystems and the stationary 2K A scenario. Additionally, the scenarios with adaptive ecosystems show a more pronounced 490 decrease in streamflow from September to January and a delay in the occurrence of the lowest streamflow from September to October. The change in mean annual actual evaporation is approximately +4 % higher in the scenarios with adaptive ecosystems and the increase mainly occurs between May and October. Instead of a year-round increase in recharge in the 2K A scenario, there is a decrease in winter recharge in the three other scenarios, resulting in a mean annual difference of -6 % between the scenarios with ecosystem adaption and the stationary scenario 2K A . Hence, the hydrological response in the 2K climate of the 495 stationary scenario 2K A is substantially different from the responses of the three scenarios 2K B , 2K C , 2K D , which consider a change in the root-zone storage capacity to reflect ecosystem adaptation in response to climate change.

Implications
The hydrological response under 2K global warming with respect to historical conditions shows distinct patterns of change if 500 we explicitly consider the non-stationarity of climate-vegetation interactions in a process-based hydrological model. We implement a dynamic root-zone storage capacity parameter, which is directly inferred from long-term and seasonal water balances of historical observations in combination with simulated historical and future climate conditions. A time-dynamic parameterization of the root-zone storage capacity was previously introduced by Nijzink et al. (2016a) in the context of deforestation, while it was implemented by Speich et al. (2020) in the context of climate change. In the latter study, forest growth in response to 505 climate change leads to a six times higher reduction of streamflow if a dynamic representation of, both, the Leaf Area Index and the root-zone storage capacity is implemented as opposed to a study in which only the Leaf Area Index varies (Schattan et al., 2013). These results are more pronounced than our findings but point towards the same direction of change. While Speich et al.
(2020) combine a forest landscape model with a hydrological model to simultaneously represent the spatio-temporal forest and water balance dynamics, we rely on a simpler approach of movements in the Budyko framework to include potential land-use 510 change.
The concept of trading space-for-time, which uses space as a proxy for time (Singh et al., 2011) could be further explored by selecting a region outside the Meuse basin with a current climate similar to the projected climate. This approach is also commonly referred to as climate-analogue mapping, i.e. statistical techniques to quantify the similarity between the future climate of a given location and the current climate of another location (Rohat et al., 2018;Bastin et al., 2019;Fitzpatrick and 515 Dunn, 2019). Finding a climate analogue for future projections in present conditions, may allow us to estimate future ω or rootzone storage capacity values in a region where the future climate may resemble today's climate elsewhere. These methods are intuitive but not straightforward, as they rely on the selection and combination of relevant climate variables and their relation with vegetation, despite non-linear vegetation responses to climate change (Reu et al., 2014).
In comparing several scenarios for the root-zone storage capacity parameter, we include some form of system representation 520 uncertainty, which improves our understanding in the modeled changes by placing them in a broader context (Blöschl and Montanari, 2010). Actual evaporation in the study catchments is projected to decrease if the historical root-zone storage capacity is used as a result of moisture stress in the root-zone of vegetation. In contrast, the increased water demand during summer is met when we assume that vegetation has adapted its root-zone storage capacity. This is an indication of disagreements among model representations on processes that become relevant in the future, in line with findings of Magand et al. (2015); Melsen 525 and Guse (2020).
The impact of climate change on low flows in the Meuse basin has been previously studied by de Wit et al. (2007). Using simulations from regional climate models which project wetter winters and drier summers, they question if the increase in winter precipitation reduces the occurrence of summer low flows due to an increase in groundwater recharge. However, they were unable to address this question with their model due to its poor low-flow performance. Our results indicate an increase in 530 groundwater recharge during winter if the historical root-zone storage capacity parameter is used, as opposed to a decrease for the models with a time dynamic root-zone storage capacity, as a result of an increased water demand for evaporation during summer. This further emphasizes the major impact of vegetation in regulating the water cycle (Luo et al., 2020;Wang et al., 2020;Stephens et al., 2021).
The land surface scheme HTESSEL, that is used in the regional climate model RACMO2 to generate the historical and 2K 535 climate simulations, assumes, as most land surface models, a fixed root-zone storage capacity. Ideally, this discrepancy between the land surface model and the hydrological model could be reduced by updating the adapted root-zone storage capacity from one model to the other in several iteration steps, thereby including soil moisture -atmosphere feedback mechanisms.

Limitations and knowledge gaps
Our study relies on the assumption that vegetation has had the time to adapt its root-zone storage capacity in a changing climate.
Yet, considering the unprecedented scale and rate of change (Gleeson et al., 2020), it is unclear how ecosystems will cope with climate change, also considering the impact of storms, heatwaves, fires and biotic infestations as a result of water stress on forest ecosystems (Lebourgeois and Mérian, 2011;Allen et al., 2010;Latte et al., 2017;Stephens et al., 2021). Additionally, when exposed to different environmental conditions, ecosystems may adapt their behavior by reducing or increasing their 545 water use to the water availability . Similarly, direct human interventions, such as the conversion of natural forests to fast-growing monoculture plantations in many parts of the world has significantly altered forests, making them more susceptible and vulnerable to disturbances (Schelhaas et al., 2003;Levia et al., 2020). However, humans also have the ability to positively influence the water cycle through vegetation, by promoting sustainable agricultural practices and integrated forest management with a simultaneous focus on biodiversity, recreation and timber production. Additionally, the conversion of exotic 550 to native species may also increase biodiversity and the resilience of ecosystems to climate change (Klingen, 2017). However, these processes are slow, implying that current management practices shape the forests of decades and centuries to come in an uncertain future. Increasing our understanding on how to include these changes in hydrological models to reliably quantify their impact is a way forward in the development of strategies to mitigate the adverse effects of climate change.
We quantify the changes in the hydrological response as a result of a changing climate in combination with several land use 555 scenarios (historical, conversion of broadleaved forests to coniferous plantations and agriculture and vice-versa). The relatively limited additional effects of land-use change on the hydrological response should be understood in the context of the relatively limited areal fraction under potential land-use conversion. The land-use changes are integrated in the root-zone storage capacity as single parameter. However, climate and land use changes likely affect other aspects of catchment functioning (Seibert and van Meerveld, 2016). For example, changes in the maximum interception storage capacity (Calder et al., 2003) are not explicitly 560 considered in the estimation of the adapted root-zone storage capacities, as the impact was shown to be relatively minor in Bouaziz et al. (2020). Additional effects of soil compaction and artificial drainage on peak flows as a result of future land conversion (Buytaert and Beven, 2009;Seibert and van Meerveld, 2016) are difficult to quantify but may partly be captured in the changed ω values.
In a first step towards temporally adaptive models, and trading space-for-time for different land-use scenarios, we did not 565 consider any additional vertical movements in the Budyko space due to the effects of increasing CO 2 levels in terms of increased productivity through fertilization, on the one hand, or water use efficiency on the other hand (Keenan et al., 2013;van der Velde et al., 2014;van Der Sleen et al., 2015;Ukkola et al., 2016;Jaramillo et al., 2018;Yang et al., 2019;Stephens et al., 2020), as these effects remain problematic to quantify in a meaningful way. Neither did we, for the same reason, consider how the relatively high ω values may be related to intercatchment groundwater losses (Bouaziz et al., 2018). Note that as our analyses 570 should be understood in the context of a sensitivity analyses of the impact of potential additional vertical shifts in the Budyko space as a result of a changing land use (Fig. 2), the potential effects on groundwater losses on the results are likely to be minor.
In addition, we performed a limited calibration of the hydrological model to retain an ensemble of plausible solutions and only used a single climate simulation despite the uncertainty in initial and boundary conditions of regional climate models. Our analyses are intended as a proof-of-concept to introduce a top-down methodology to quantify non-stationarity of the root-zone 575 storage capacity parameter through optimal use of projected climate data, rather than a comprehensive climate change impact assessment of the Meuse basin.

Conclusions
Understanding non-stationarity of hydrological systems under climate and environmental changes has been recognized as a major challenge in hydrology (Blöschl et al., 2019). Despite our strong awareness of non-stationarity of hydrological param-580 eters, we often lack knowledge to implement system changes in hydrological models. In this proof-of-concept study in the Meuse basin, we propose a top-down approach to introduce a time-dynamic representation of the root-zone storage capacity parameter within process-based hydrological models, using regional climate model simulations. Our approach relies, on the one hand, on a space-for-time exchange of Budyko characteristics of dominant land-use types to estimate the hydrological behavior of potential land-use changes and, on the other hand, on the interplay between the long-term and seasonal water budgets 585 to represent climate-vegetation interactions under climate and land-use change. Despite knowledge gaps on future ecosystem water use, we implement potential system changes in a hydrological model based on our current understanding of hydrological systems. The predicted hydrological response to 2K warming is strongly altered if we consider that vegetation has adapted its root-zone storage capacity to offset the more pronounced hydro-climatic seasonality under 2K global warming compared to a stationary system. The increased vegetation water demand under global warming results on average annually in -7 % less 590 streamflow, +4 % more evaporation and -6 % less recharge for the scenarios assuming non-stationary conditions compared to a stationary system. These differences even lead to a distinct change of sign of median annual streamflow. Our study contributes to the quest for more plausible representations of catchment properties under change and, therefore, more reliable long-term hydrological predictions.    Figure 2. (a) Representation of the Budyko space, which shows the evaporative index (EA/P ) as a function of the aridity index (EP/P ) and the water and energy limit. A catchment with aridity index (EP/P ) obs and evaporative index (EA/P ) obs , derived from observed historical data, plots at location p1 on the parametric Budyko curve with ω obs . A movement in the Budyko space towards p2 along the ω obs curve is shown as a result of a change in aridity index (EP/P )∆ towards a projected (EA/P ) 2K,ωobs associated with aridity (EP/P )2K.
An additional vertical shift ∆ω towards a location p3 on a ω change curve is shown if additional factors (e.g. land use) are projected to change besides aridity index. Here, the represented downward shift in ω reduces the change in evaporative index to (EA/P ) 2K,ωchange .
(b) Cumulative storage deficits (S R,def ) derived from effective precipitation (PE) and transpiration (ER) using the simulated historical and 2K climate data. Estimates of transpiration (ER) are derived from long-term water balance projections as a result of movements within the Budyko framework in response to climate and potential land use changes.  . Model scenarios using the observed historical and the simulated historical and 2K climate data. The model is calibrated using observed E-OBS data and the historical root-zone storage capacity SR,max,A. The model is then forced with the simulated historical climate data using SR,max,A as root-zone storage capacity parameter. We then define four scenarios to compare the change in hydrological response to 2K global warming in comparison to historical conditions for the ensemble of feasible parameter sets. In scenario 2KA, we assume an unchanged system (no changes in land use, nor root-zone storage capacity). In scenario 2KB, we assume that vegetation has adapted its root-zone storage capacity to the 2K climate, but no changes in land use. In scenario 2KC, we test the combination of an adapted root-zone storage capacity in response to the changed climate and a hypothetical conversion of coniferous plantations and agriculture to broadleaved forests in part of the catchment. A similar but reversed approach in land-use changes is assumed in scenario 2KD.