Long-term drought effects on landscape water storage and recovery under contrasting landuses

Recent extreme droughts in Europe have highlighted the urgent need to quantify their effects on ecohydrological fluxes (evapotranspiration, groundwater recharge) and water storage (mainly soil moisture) in the landscape. In response, we combined process-based (EcH 2 O-iso) and machine learning (NARX) models to estimate the enduring effects of long-term drought on water fluxes and storage and to project future short-term groundwater levels and recovery potential under various precipitation scenarios. The work was undertaken at the Demnitz Mill Creek (DMC), a 70 km 2 mixed land use (arable crops and forestry) catchment in northern Germany. Our simulations indicated that the extreme drought years of 2018 and 2022 had the most marked impacts, leading to substantial declines in groundwater recharge ( > 40 %), evapotranspiration (up to 16 %) and soil moisture (up to 6 %). Simulations indicated that groundwater levels may not recover in the next 15 years if recent precipitation anomalies persist. These findings underscore the urgent requirement for enhancing resilience and promoting integrated strategies in managing land and water resources to optimise water retention in landscapes and to better respond to


Introduction
In recent years, the European lowlands have experienced significant increases in the frequency and severity of droughts, with a series of consecutive drought years starting in 2018, which is unprecedented in the past 2110 years (Büntgen et al., 2021).Most recently, there was a 1in-500-year drought in 2022 (Toreti et al., 2022).These extreme drought conditions have made the study of their impacts on landscape water storage and ecosystem resilience both urgent and important.The extensive lowlands of Central Europe, characterized by their flat topography and contrasting landuse patterns, play a crucial role in the provision of water supplies, agricultural production, and the socioeconomic stability of the continent (Buschmann et al., 2020;A. Smith et al., 2021).Specifically, over a third of Central Europe is devoted to extensive agricultural use (Hari et al., 2020), while conifer forests constitute 46 % of all European forests (Forest Europe, 2020).Numerous studies have confirmed that past drought events (e.g. 2003 and 2018) resulted in considerable deficits in water storage, affecting both soil moisture availability and groundwater levels (Hellwig et al., 2020;Samaniego et al., 2018).These events subsequently led to substantial crop yield losses and forest damage across Europe (Ciais et al., 2005;Naumann et al., 2021).According to a cross-sectoral study in Germany, the 2018 drought resulted in 20-40 % reduction in staple crops and forced logging of 25.1 million tons of damaged forests with possibility of other socio-economic instability in future drought events (Conradt et al., 2023).
Given the importance of soil moisture and groundwater as the essential sources of water storage to support vegetation growth and stream flow generation, tracking and resolving ecohydrologicallymediated fluxes and storage such as evapotranspiration, soil moisture and groundwater recharge is of paramount importance to understand water availability and ecosystem functioning at a range of spatial and temporal scales (Smith et al., 2022a).Groundwater, in particular, serving as a principal water resource for drinking and irrigation, requires extensive and continuous monitoring and predictions for sustainable, fair and economically sound use (Chang et al., 2016).Quantitative understanding of the spatio-temporal variations of ecohydrological fluxes and water storage dynamics (e.g.soil moisture and groundwater level) is key to managing water resources particularly during major environmental perturbations such as droughts (Dari et al., 2019;Guse et al., 2019).Furthermore, as climate change further exacerbates the risk of drought, understanding the complex dynamics of water availability is becoming increasingly important for effective mitigation strategies (Walker and Van Loon, 2023).However, the effects of long-term drought on landscape water storage in the European lowlands are poorly understood, despite their importance for maintaining ecosystem services, supporting agricultural production, and ensuring water security for human populations.
Due to the advantages of machine learning models, such as high efficiency in processing large datasets, automation in capturing complex non-linear relationships and their skill in forecasting and prediction based on historical data, they have been increasingly applied to solving ecohydrological problems in recent years (Bergen et al., 2019).However, data-driven machine learning techniques have some disadvantages when compared to traditional process-based ecohydrological models such as limited process interpretability (black box) and difficulties in discerning underlying causal relationships (Tague and Frew, 2021).While process-based models are built on conceptualising the underlying processes, data-driven approaches do not provide such insights.For example, they do not consider the closure of water and energy balances, which are fundamental to process-based ecohydrological modelling (Fatichi et al., 2016;Koppa et al., 2022).Consequently, a promising research direction, and the approach adopted in this study, is integrating process-based understanding with the adaptability and efficiency of machine learning techniques, capitalizing on the strengths of both approaches to better assess drought impacts.
Tracer-aided ecohydrological models offer integrated tools for analyzing the complex interactions between climate, hydrology, and ecosystems, potentially providing crucial insights for policy makers (Vose et al., 2011).These models account for key vegetation influences, such as canopy and rooting properties, as well as management effects (Guswa et al., 2020); and are applicable across various spatial and temporal scales.Through incorporating environmental tracers, e.g. by using stable isotope mass balances, these models allow for more robust representation of the processes influencing ecohydrological fluxes, water partitioning and landscape water storage dynamics under drought conditions (Smith et al., 2022a).This facilitates better tracking of water storage distributions, constraints on ecohydrological fluxes, disaggregation of evapotranspiration components and water age estimates (Birkel et al., 2011;Hrachowitz et al., 2013;Kuppel et al., 2020).
Predicting groundwater levels is commonly used to aid sustainable long-term water resource management, as it facilitates assessment of the direct impacts of environmental change on groundwater resources (Wunsch et al., 2022).In recent decades, artificial neural networks (ANNs) have been widely applied to capture and predict changing groundwater levels in response to environmental drivers (Chang et al., 2016;Coulibaly et al., 2001;Sun, 2013;Wunsch et al., 2021).Similar to other machine learning techniques, the main advantage of ANNs is their capability to model non-linear relationships between input and output variables.A 1D-ANN can be trained to mimic the correlational dynamics between environmental drivers and resultant groundwater level time series data without priori assumptions, and to accomodate noisy input data and adaptation to changing environmental patterns (Sun, 2013).Coulibaly et al. (2001) used insitu hydroclimate measurements and historic groundwater levels to design ANNs for forecasting groundwater levels.Wunsch et al. (2021) developed ANNs to predict groundwater levels using solely climatic input variables.Sun (2013) employed remote sensing data in the formulation of ANN ensembles and found that it enhanced the prediction performance of groundwater levels, particularly during extreme climatic events such as droughts.
By incorporating the outputs of process-based ecohydrological models into the development of ANNs, one can leverage the strengths of both process-based and data-driven approaches and evaluate the potential impacts of various future precipitation scenarios on the recovery of groundwater after extreme droughts to pre-drought levels.However, it is important to recognize that despite the advantages of rapidity and efficiency offered by such hybrid approaches, the improvement in model performance may vary, underscoring the need for evaluations in such ecohydrological applications (Kraft et al., 2020).
To increase understanding of the broader implications of drought on water resource resilience, this research uses data from a 70 km 2 lowland catchment in Northern Germany.The study combines ecohydrological models and ANN to understand changing water fluxes and groundwater recharge in response to climatic variability.As such, the study aims to contribute to the evidence base needed to inform the development of sustainable and adaptive land and water management strategies tailored to the unique challenges posed by long-term droughts in this part of Germany.This study aims to address the following questions: How does drought influence water partitioning and storage in complex heterogenous lowland landscapes?How much time is needed to recover water levels to pre-drought conditions under contrasting landuses for different potential rainfall inputs?Does the incorporation of process-based model outputs into a machine learning model allows predicting recovery times of groundwater levels after a drought?

Study area
Demnitzer Millcreek is a meso-scale (~70 km 2 ) catchment in the North European Plain (NEP) and typical for these lowlands, with a gentle topography, ranging from 30 to 90 m.a.s.l (Fig. 1).The soils are primarily based on glacial deposits and the water table is usually 1-3 m below the surface.Landuses reflect the hydrogeological units representative of the NEP with two monocultures dominating the catchment.Non-irrigated croplands account for over 50 % of the land area, while conifer forests cover approximately 30 % (Fig. 1c).Dominant soil types include silty brown earths on poorly-drained glacial till, which support arable crops, and podzols on freely-draining sands that sustain forests, mainly comprising Scots Pine (Pinus sylvestris).Along the stream network, poorly drained peat and gley soils provide suitable conditions for pasture for grazing sheep and cattle.Small areas of broadleaved forests, mainly oak and beech, are also present.
The climate is classified as humid continental, with an average annual precipitation of around 600 mm and a mean annual temperature of approximately 10 • C.Over the past decade, precipitation has declined on average by ~20 mm per year, with negative precipitation anomalies becoming more frequent (Smith et al., 2022a).The region also experiences high potential evapotranspiration (PET), at around 600 mm per year, resulting in a negative water balance even in non-drought years.As such, the runoff coefficient is also low (<10 %).

Datasets
The datasets used for forcing, calibration, and validation purposes are provided in Table S1.The daily climate data, which include precipitation, air temperature, relative humidity, and wind speed, were sourced from five nearby weather stations (<10 km) (German Meteorological Service, 2023).These weather stations were strategically chosen to represent the five climate zones in the catchment defined using Thiessen polygons.Daily short-and long-wave radiation data, leaf area index (LAI), along with evapotranspiration within the catchment were extracted from gridded global datasets as listed in Table S1.
Since 2018, soil moisture monitoring has been conducted at a 15-min interval and then aggregated to daily scale using volumetric soil water content sensors (SMT-100, Umwelt-Geräte-Technik GmbH, Germany) with an accuracy of ±3 %.The monitoring has been undertaken at three depths (20, 60 and 100 cm) with 2 replicas per depth under different landuses, including forest, grassland and cropland) across the catchment (two locations at Forest A and one location at Alt Madlitz; Fig. 1a) (Kleine et al., 2020;Landgraf et al., 2022b).Sap flow measurement were conducted with sap flow meters (SFM1 instrument, ICT International, Australia) at tree breast height within a mixed forest stand located at Forest A in 2018 (Fig. 1a) (Landgraf et al., 2022a).Since 2001, daily stream water discharge has been derived from the measured water levels using a rating-curve method at Demnitz and Demnitz Mill, while daily groundwater levels have been monitored at four locations (GW3-GW8) (Fig. 1b).The water levels are recorded with dataloggers (AquiLite ATP 10, AquiTronic Umweltmeßtechnik GmbH, Kirchheim/Teck, Germany) (Smith et al., 2020).
Stable water isotopes have been analysed from bulk soil samples, precipitation, groundwater and stream water.Bulk soil samples were monthly at different, dominant landuse-soil units collected in 2018 and analysed for stable water isotopes using a direct equilibrium method (Kleine et al., 2020).Precipitation water was sampled with an ISCO 3700 autosampler (Teledyne Isco, Lincoln, USA) at Hasenfelde (Fig. 1b), with paraffin oil added to the sampling bottles to prevent evaporation (Landgraf et al., 2022b).Grab samples of stream water isotopes were collected on a weekly/daily basis.Monthly sampling of groundwater was conducted manually using a submersible pump (Kleine et al., 2020).All liquid samples were analysed for deuterium (δ 2 H) and oxygen-18 (δ 18 O) in the laboratory using a Picarro L2130-i cavity ring-down water isotope analyser (Picarro, Inc., Santa Clara, CA, USA) (Landgraf et al., 2022b).

EcH 2 O-iso model
We employed the tracer-aided spatially-distributed ecohydrological model EcH 2 O-iso to quantify ecohydrological fluxes (i.e.evapotranspiration and residual groundwater recharge) and storage dynamics for the entire catchment.EcH 2 O-iso incorporates a tracer module for tracking stable water isotopes (δ 2 H and δ 18 O) transformations and their transport in water fluxes across diverse catchments (Kuppel et al., 2018).The model emcompasses four components coupling energy, water, vegetation and isotopes dynamics.The energy balance is determined using a vertical scheme from the canopy to the land surface, driven by incoming radiation, air temperature, relative humidity, and windspeed, with latent heat fluxes calculated between the canopy and atmosphere, associated with the evaporation of canopy-intercepted water and transpiration.Ground heat flux is calculated for two soil thermal layers and is closely associated with soil heat capacity and thermal conductivity, which depends on soil moisture (Liebethal and Foken, 2007).The water balance follows a multi-layered structure with a vegetation canopy, soil surface, and three subsurface layers, estimating infiltration using the Green-Ampt approximation of the Richards equation.Lateral flow can either occur in the deepest subsurface layer if saturated, or as overland flow if the upper soil layer is saturated or rainfall instensity exceeds infiltration capacity.Vertical and lateral flow in the subsurface is estimated by kinematic routing once field capacity is exceeded.Soil evaporation is restricted to the top layer, while root water uptake is parameterized exponentially across the three subsurface layers from the surface to the saturatured zone.A fully integrated vegetation module employs the K root parameter to determine the transpiration fraction derived from each soil layer and accounts for multiple vegetation types and bare soil, using a Jarvis-type multiplicative model for stomatal conductance estimation (Jarvis, 1976).The tracer module assumes full mixing in each storage compartment, with isotopic fractionation occurring only during soil evaporation (Craig et al., 1965;Gat, 1995).The isotopic composition of fluxes out of each compartment is the same as that stored for the same time step.More information on the model and its recent applications can be found in (A.Smith et al., 2021;Smith et al., 2022a).

Model configuration, calibration and validation
Based on previous calibration of EcH 2 O-iso (fully detailed in A. Smith et al., 2021;Smith et al., 2022b), this study extended the study period to cover the 2022 drought.The model was set up to run at 250 m grid size at daily time steps for a simulation period of 16 years (Jan 2007 -Dec 2022).The initial two years of simulation (2007)(2008) served as a warm-up period to initialise soil moisture, groundwater storage, and discharge, including volumes and isotopic compositions, among other factors.Subsequent to this initial period, the model was calibrated using data from 2009-2014, and 2018-2019, which covered the 2018 drought.The years 2015-2017 and 2020-2022 were set aside and used as validation periods to test the robustness of the model and to ensure its reliability in various climate conditions including the 2022 drought.The model was driven by the comprehensive set of climate data outlined in Sec 2.2, including precipitation, air temperature, relative humidity, wind speed, short-and long-radiation, as well as LAI and precipitation isotopes.The application of forcing data within the catchment was based on the distance to the surrounding measurement locations, which were assigned to their corresponding climate zones.Considering the limited observed spatial variation, precipitation isotopes were uniformly distributed across the entire study domain.
Model multi-criteria calibration and valdation have been conducted using a variety of hydroclimatic and isotopic datasets including field measurements such as soil moisture, sapflow-derived transpiration, groundwater levels, isotopic data of streamwater, bulk soil samples and groundwater, as well as remote-sensing data such as MODIS 8-day evapotranspiration and latent heat (Running et al., 2017) (Table S1).Ten top-performing runs were chosen based on the best overall efficiency, using empirical cumulative distribution functions derived from Nash-Sutcliffe efficiency (NSE) (Nash and Sutcliffe, 1970) and root mean square error (RMSE).These runs were selected from a total of 100,000 parameter configurations, which were generated by combining the parameters that passed the Morris sensitivity analysis for the multi-criteria calibration (Morris, 1991;Sohier et al., 2014).The Morris sensitivity analysis is particularly useful in dealing with high-dimension parameters, which involves a systematic and step-wise examination of the parameters.Each parameter is varied one at a time and other parameters remain at their previous values to measure the effect of the change in this parameter on model outputs and to quantify the degree of variability.In this study, each parameter was systematically varied by 50 % of its range using radial sampling to cover the entire range of possible model parameters, thus identifying the parameters with the highest sensitivity.Table S2 presents the identified most sensitive parameters related to soil, vegetation, and channel properties, as determined by averaging the trajectories' outcomes.

Artificial neural networks 2.4.1. Nonlinear autoregressive exogenous (NARX) model
The Nonlinear Autoregressive Exogenous (NARX) model is a type of recurrent ANNs commonly used for time-series predictions (Wunsch et al., 2018).The mathematical expression of NARX model is given in Eq. (1).
where y(t) is the dependent output that is regressed on previous values of the output and previous values of independent (exogenous) input variables x; f is the nonlinear function approximated by a feedforward neural network; n y denotes the feedback delays (FD); n x denotes the input delays (ID).
As shown in Fig. 2, the structure of NARX model employed in this study consists of one input layer, one hidden layer and one output layer.NARX relates the current value of a time series to its past values and exogenous (independant) time series.In open loop mode, the output values are predicted based on previous observations such as climate time series and the true past values of the target (Fig. 2a).This mode is often used during training and testing phases.In closed loop mode, the NARX model incorporates its previous predictions as input (Fig. 2b), and is often used when deploying the model in further predictions.However, closed loop mode is also used for training when Bayesian optimization is used to optimize hyperparameters (Wunsch et al., 2021).Both modes allow for a potential delay to both exogenous variables and target variables.These types of delays are referred to as input delay (ID) and feedback delay (FD), respectively (Fig. 2) (Beale et al., 2010).NARX can also retain information for a duration two to three times longer compared to other recurrent neural networks (Wunsch et al., 2018).Compared to other ANNs and state-of-the-art deep learning techniques such as long short-term memory (LSTM), NARX has demonstrated superior performance in predicting groundwater levels in terms of accuracy (Wunsch et al., 2021).

NARX model: Model selection, training, testing, and projection
We employed the NARX model to project the the duration required to restore average groundwater levels in the catchment under various potential rainfall scenarios.The model was set up for four groundwater wells where data were available daily (Fig. 1b and Table S1).Input variables, aside from groundwater levels (GWL), included widely available climate data, such as precipitation (P), temperature (T), and relative humidity (RH), as well as the derived evapotranspiration data (ET) as an output from the EcH 2 O-iso modelling.Precipitation primarily drives groundwater recharge, whilst temperature and relative humidity account for the conditions affecting groundwater loss and concurrently provide the seasonal information in annual cycles (Wunsch et al., 2021).The analysis was conducted at weekly time steps, with the majority of this process carried out in MATLAB R2022b.
Input data for NARX were splitted into different time blocks: Data from before 2021 were used to select and train models using a Bayesian optimization method (Kraft et al., 2020).Data from Jan 2021 to Dec 2022 were fixed as test and evaluation datasets.During the Bayesian optimization, 80 % of the data were used for training, 10 % for early stopping (dropping non-promising runs) and the remaining 10 % for testing purposes (to disdinguish it from the fixed test dataset in 2021-2022, this testing is denoted as "opt-testing") (Fig. 7b).The training datasets covered both wet years (e.g.2010) and drought years (e.g.2018), which could help model training to better capture the groundwater level changing pattern.All data were normalized between − 1 and 1 before training to reduce the impact of different scales.The training was conducted using the Levenberg-Marquardt algorithm to find the optimal model architecture and hyperparameters (Adamowski and Chan, 2011).Different combinations of input variables (i.e.T, RH and ET; P was a fixed input) and hyperparameters (e.g.hidden neurons, ID and FD) were examined during the optimization.We implemented one hidden layer, which uses the hyperbolic tangent sigmoid activation function, with a variable number of 1-20 hidden neurons to first train the model in closed loop mode.ID and FD were optimized between 1 and 52 (equivalent to 1 year).We established an upper limit of 30 epochs for the training process.Furthermore, the maximum validation failures were set to 5 epochs, after which early stopping of the model training would occur.The number of optimization iterations was between 25 and 50, with each iteration comprising 5 runs to mitigate uncertainties.The objective function was formulated based on the mean of the 5 runs, calculated in terms of the sum of NSE and the squared Pearson correlation coefficient (R 2 ).The model architecture with the highest value of the objective function was selected for further testing and evaluation.To account for uncertainty in testing, ten repetitions were conducted for each groundwater well, following Wunsch et al. (2021).After achieving satisfactory results in testing based on three efficiency criteria-NSE, R 2 , and RMSE-the model loop networks were used for future projections.
Before projections, we examined the long-term (2007-2022) precipitation pattern for alterations using the Pettitt test, a common method for identifying change points in hydrological time series (Luo and Chui, 2020).The precipitation patterns were then categorized into two periods, each exhibiting distinct precipitation pattern with statistical significance (p < 0.05).Recent years showed an exceptionally low average precipitation (<495 mm) and were subsequently designated as dry years.Consequently, the scenarios for groundwaer level forecasting emcompassed: "normal" years (2007-2013), dry years (2014-2022), 2010 as a wet year and 2018 as a drought year.These scenarios were iteratively simulated 2-15 times over a 15-year period.It should be noted that not all input variables were used for projections.Following Bayesian optimization, the input variables were selectively included based on the optimal model architecture, which comprised a fixed P component and optional variables such as T, RH or ET.Instead of using a detrended seasonal component of each variable for projection, the historical datasets preserving the actual variability were directly extracted.Finally, these input datasets were used for predictions with the ten tested models.

EcH 2 O-iso model performances of ecohydrological fluxes
Across different observational data under different landuses, overall uncertainties of EcH 2 O-iso were effectively managed (Figs. 3 and 4).The model successfully captured the different fluxes under contrasting landuses.At both discharge gauge stations Demnitz Mill and Demnitz, the simulation showed good performance, with both calibration and validation NSE exceeding 0.7 and RMSE falling below 0.1 at the former (Fig. 3b).The model also maintained reasonable performance at Demnitz with a calibration NSE of 0.66 and RMSE of 0.11 (Fig. 3c).
The model also demonstrated satisfactory performance in simulating soil moisture variations, with the majority of observations aligning within the ensemble simulation range with the RMSE ≤ 0.05 for both sites (Fig. 3d and e).Importantly, as the Forest A model cell contains different landuses, we only show the comparison between the observation collected in the forested area within this cell with the simulation in a nearby cell with 100 % forest cover.
The seasonal variations of evapotranspiration were also been effectively captured at both the forest and cropland sites with RMSE values of ~1 (Fig. 3f and g).
The calibration of stream isotopes at three locations showed low uncertainties with calibration RMSE values not exceeding 3.6 ‰ and validation RMSE not exceeding 5.2 ‰ (Fig. 4a-c).From 2018 to 2019, the median groundwater isotope simulation results also showed a good agreement with the observed data, exhibiting minimal seasonal variations (Fig. 4d).Multiple monthly observations of soil water isotopes also exhibited reasonable alignment with the simulation ranges at different soil layers (Fig. 4e and f).

Long-term ecohydrological fluxes and storage
Fig. 5 presents the percentage changes in key ecohydrological indicators averaged at the catchment scale since 2018, demonstrating the cummulative effects of recent dry years.Apart from precipitation, the other three ecohydrological indicatorsevapotranspiration (ET), soil moisture in the upper soil layer (top 15 cm) and groundwater recharge were all estimated by EcH 2 O-iso.In recent years, negative anomalies in precipitation have had severe adverse impacts: reduced groundwater recharge, a decline in water storage represented by soil moisture, and consequently, a decrease in ET due to limited water storage.Using the long-term average (i.e.2009-2022) as a Baseline, since 2018 there have been five consecutive dry years, among which 2018 and 2022 were particularly severe drought years with a reduction in precipitation of 30 % and 26 %, respectively.In the remaining three years, the decrease was <10 % but still less than the long-term average values.Corresponding to precipitation, the variable showing the most significant reduction was groundwater recharge, with average declines of 39 % and 44 % in 2018 and 2022, respectively, compared to the long-term average.ET decreased by − 1% to − 16 %, and surface soil moisture decreased by − 1% to − 6%.

Landuse effects on water fluxes and storage
Fig. 6 shows the impacts of different landuses on groundwater recharge, soil moisture in the upper soil layer and ET.The dominant landuse types in the whole NEP, namely conifer forests and croplands, were specifically considered in the analysis.The analysis of all cells provides a comprehensive overview of the catchment as a whole and mitigate any potential bias from cells representative of a single landuse.Fig. 6a presents the distribution of the ecohydrological fluxes (i.e.groundwater recharge and evapotranspiration) and water storage (i.e.soil moisture) in the model domain under different landuses; Fig. 6b shows their percentage changes.In general, the ecohydrological variables exhibited most significant declines in 2018 and 2022 for both conifer forests and croplands.During the Baseline period (2009-2022), croplands showed significantly higher long-term average groundwater recharge than conifer forests.This pattern was not evident during the extreme drought years in 2018 and 2022 due to the substantial decline in groundwater recharge under croplands (Fig. 6a and Figure S2).However, the percentage change observed in croplands was slightly lower than that of conifer forests due to the higher original values and exhibited larger differences in extreme drought years (2018 and 2022) and dry years (2019-2021).For conifer forests, all five years showed a decline in recharge from the long-term mean of >25 %, with the highest decrease exceeding 38 % in 2022 (Fig. 6b).Soil moisture was highest under pasture, as revealed by the Baseline panel (Figs. 1 and S2).Similar to groundwater recharge, the modelled percentage decline in soil moisture was most significant in 2018 and 2022, with croplands showing a higher decline (up to 10 %).ET was estimated as being higher in conifer forests than in croplands (largely due to higher interception and transpiration; Figure S2 and S3), as indicated by the Baseline, and conifer forests experienced a generally greater decline (Fig. 6a), resulting in more uniform, supressed evapotranspiration distribution especially in drought years (2018 and 2022) (Figure S2).Although there was a small increase of less than 1 % in croplands' ET in 2019 and 2020 compared to the Baseline, it was not enough to offset the overall decline in ET observed in both landuses during the dry years (Fig. 6b).

Recovery times for groundwater levels
To estimate the potential groundwater recovery times after the drought, we first analyzed the long-term precipitation patterns to develop different precipitation scenarios that would be used in the projections.A shift in the precipitation pattern was identified in 2013, with the years before or after 2013 categorized as normal or dry years, with an average annual precipitation of 635 and 495 mm, respectively (Fig. 7a).Fig. 7b displays the gradual decline of groundwater levels at GW4 in the past decade, accompanied by projections for future groundwater levels under different potential rainfall scenarios.Similar declining trends in the past decade were also observed at the other three wells (Figure S1).
The optimal NARX model architecture was determined based on the objective function value derived from the Bayesian optimization (Table S3).The performance during the optimization across all four wells indicated robust model capability with the sum of NSE and R 2 values ranging from 1.2 to 1.6.Regarding the three optional input variables (i.e.RH, T, and ET), no single variable was found to be universally critical for all wells.For instance, ET as the output of EcH 2 O-iso, was not selected for GW4 as it did not enhance the model performance.However, for all other wells, ET was selected as the decisive input variable for groundwater level prediction.Fig. 7b displays the NARX model Bayesian optimization, test-and prediction results, including the testing efficiencies and the projections for one groundwater well (GW4) under different climate scenarios.The NSE was 0.89, while the RMSE was 0.12 m, and the R 2 was 0.83.Performance metrics for the NARX models at other groundwater wells in the catchment are given in Figure S1, all showing satisfactory results.All sites exhibited NSE values no less than 0.64, RMSE values not exceeding 0.22 m, and R 2 values above 0.58.Fig. 7b and Figure S1 also show the most recent groundwater level data from 2023, overlaid with our projections.Except for GW3, which was decommissioned in late 2022, the data from the the remaining wells confirmed that our projections are within a reasonable and justifiable range.
Fig. 7c shows the time needed for the annual average groundwater levels to return to the normal-year average under different climate scenarios.Apart from different scenarios, the analysis of each well also incorporates the median, as well as the lower and upper bounds, derived from ten repetition of the machine learning model to account for the model uncertainties.If future climate patterns resemble those of 2018 (drought year), the groundwater levels at all four locations will see a continuous decline (Fig. 7b and Figure S1) and will not return to normalyear groundwater levels (Fig. 7c).If climate patterns are akin to the dryyear average (post-2013), the median projections suggest that − except for GW8 − the three wells might require an average recovery period of 6 years, with GW4 potentially taking up to 13 years.In contrast, all wells are anticipated to recover under upper bound projections, with the recovery duration shortened to an average of 2.5 years.However, under lower bound projections, none of the wells are expected to recover.If future climate patterns are similar to those in 2010 (wet year), median projections suggest that all sites can return to the normal-year average within 1-2 years.Upper bound projections also suggest a recovery for all wells within 1 year.However, under lower bound projections, only GW3   can recover to the normal-year within 1 year, and the remaining wells will not recover.When climate patterns align with the normal-year average (pre-2013), all sites can recover within a timeframe of 1-3 years and 1-2 years for median and upper bound projections, respectively.The average recovery time is 1.8 years for median projections and 1.5 years for upper bound projections.Similarly to 2010 scenario, under lower bound projections, only GW3 will potentially recover within 1 year.

Value of linking process-based models and machine learning for drought assessment
The advantage of integrating process-based models with machine learning approchaes is to preserve the physical consistency and interpretability inherent in process-based models, while also benefiting from the capabilities of machine learning models; namely providing more efficient, data-driven representations of processes that may be inadequately understood (Koppa et al., 2022).Several applications have demonstrated the effectiveness of such hybrid approaches in modeling ecohydrological processes (Kraft et al., 2022;Lin et al., 2021;Rasp et al., 2018).In our study, the potential benefits of incorporating the output of the process-based model into the machine learning model was demonstrated.During the determination of the machine learning model architecture, the EcH 2 O-iso output (i.e.ET) was identified as an important input variable for predicting groundwater levels in three out of four wells.Thus, incorporating ET estimates from model outputs is valuable in assessing groundwater drought.A study analyzing the groundwater drought in the UK also showed that ET plays a key role in forming and propogating groundwater drought, although estimating ET is challenging (Bloomfield et al., 2019).In alignment with this, the hybrid approach adopted in our study facilitated integrating insights from EcH 2 O-iso, calibrated with extensive, long-term datasets including tracer data, and meanwhile allowing projections using smaller historical datasets.This would not have been possible with process-based models or data-driven models alone.

Effects of drought on water flux partitioning and storage
In recent decades, Europe has witnessed major droughts in 2003, 2015, 2018 and 2022.These events with more frequent anomalous high temperatures and rainfall deficits are now considered not exceptional but a "new normal" to come in the next decades (van der Woude et al., 2023).A recent reconstruction analysis revealed that a longer-term drying trend in Europe started as early as 1850 (An et al., 2023).The extreme drought event in 2018, followed closely by another even more severe event within five years (2022), underscores the escalating frequency and intensity of these extremes.The primary driver of these enhanced summer heatwaves in Europe has been attributed to anthropogenic warming, specifically with effects in the arctic affecting large scale weather patterns (Samaniego et al., 2018;Zhang et al., 2020).When combined with enhanced land-atmosphere feedback (i.e.anomalies in large-scale atmospheric circulation), these factors collectively suppress precipitation, thereby triggering drought conditions (Zhou et al., 2019).Both soil moisture and groundwater recharge declines as observed in our study are predominantly driven by these negative precipitation anomalies.Consequently, diminished water availability in the subsurface (e.g.low soil moisture and deeper groundwater) can restrict evapotranspiration, manisfesting in the reduction in the transpiration, photosynthesis and biomass production (Seneviratne et al., 2010).Figure S3 clearly illustrated this process by showing the differences in fluxes and storages simulated by EcH 2 0-iso in the wet (2010) and the drought years (2018).
Adequate incoming precipitation allows for differentiation in wateruse strategies between different vegetation communities.For instance, forests exhibit higher interception and transpiration losses (Figure S3).Consequently, lower infiltration and groundwater recharge were simulated in the forested parts of the catchment.When precipitation almost halved, the initial impact was on the intercepted water and infiltration, subsequently leading to reduced soil moisture and groundwater recharge.The water scarcity during the drought year narrowed the differences in water fluxes between forests and croplands (shown in Figure S2 with an overview of the entire study catchment).During the two main drought years (2018 and 2022), the spatial variations in groundwater recharge and evapotranspiration were small across the contrasting landuses indicating extensive suppression of fluxes during periods of water limitation.

Drought recovery
Our study revealed that the diverse landuse patterns, ranging from agricultural to forest ecosystems, exhibited varying degrees of vulnerability to drought.As indicated by the more substantial decline during extreme drought years (Fig. 6a and S2), in our study region croplands seem more vulnearable/sensitive to the effects of drought on groundwater recharge but also with the potential to recover quicker after drought.Concurrently, croplands encountered a greater decline in soil moisture, with a slight rise in evapotranspiration during wetter years in the study period (Fig. 6b).On the other hand, conifer forests demonstrated longer-lasting negative effects from drought on groundwater recharge since the extreme drought year 2018 (Fig. 6b).In the following years, simulations showed that rewetting from 2019 to 2021 was probably inadequate for groundwater recharge to recover under the conifer forests.The differences in groundwater recharge recovery between the two landuses can be attributed to the nonlinear relationship between the simulated soil moisture and groundwater recharge (Figure S4).In both forests and croplands, groundwater recharge remains low until soil moisture exceeds a specific threshold, after which it increases sharply.This pattern is consistent with the nonlinear exponential relationships established in Apurv et al. (2017) and Berthelin et al. (2023).Interestingly, both the soil moisture and recharge ranges at our site in Central Europe were of the same magnitude as those in an arid climate (Apurv et al., (2017).In contrast, the soil moisture range was about half of that in a more humid and permeable karst aquifer in Germany and groundwater recharge was also one order of magnitude smaller (Berthelin et al., 2023).During drought years when soil moisture further decreased it would result in greater reduction in groundwater recharge in the cropland than under conifer forest (Fig. 6a).This nonlinear relationship is particularly informative on how groundwater droughts can develop locally under different landuses.
The observed variations in vulnerability across different landuses can also be attributed to their distinct subsurface storage characteristics (in turn strongly dependent on landscape evolution and pedogenesis).Particularly in catchments receiving lower precipitation inputs, storage capacity becomes crucial in retaining water within the catchment soils and aquifers and facilitating its gradual release (Apurv and Cai, 2021;Carey et al., 2010).Both cropland and conifer forest systems are located on soils with relatively low storage capacities; though the brown soils of the croplands are more silty throughout their profile and more retentive than the sandier forest podzols.However, for the croplands, low soil organic carbon contributes to the low water-holding capacity in the top soil, thereby reducing the overall soil moisture (L.C.Smith et al., 2021).This partly explains why evapotranspiration and groundwater recharge of croplands could recover slightly after the 2018 drought, but the decline in soil moisture was greater and lasted longer (Fig. 6b).Unlike croplands, conifer forests typically have denser canopy cover and high interception and transpiration rates, generally consume more water, with the former reducing net rainfall reaching the ground and infiltrating (then percolating) into the subsurface all year round, thus limiting the subsurface water storage capacity in the system (Douinot et al., 2019;Zal et al., 2015).This also explains the slower recovery groundwater of conifer forests if negative precipitation anomalies persist after drought (Fig. 6).
The projected recovery times under different precipitation scenarios raise serious concerns regarding groundwater resources over the coming decade if future precipitation negative anomalies continue.The upper bound projections represent best-case scenarios, while the lower bound projections represent more conservative estimates.In the worst-case scenarios groundwater recovery is nearly impossible for all sites, except for GW3 in wetter years.Here, we discuss median projections as the maximum drought recovery time.According to our models, given the current landuse patterns and recent average precipitation levels (i.e., 495 mm per year), optimistically all sites will have the potential to recover to pre-drought groundwater levels with an average of 2.5 years.Under longer-term average precipitation levels (635 mm), the simulations suggest it will only 1-2 years for all sites to recover.A study encompassing all of Germany revealed that the northern regions tend to experience longer-lasting groundwater droughts as a result of lower precipitation than Germany's south due to lower hydraulic conductivity and flatter topography.These droughts effects have the potential to persist for years, averaging ~2.4 years in the North, which aligns with our findings (Hellwig et al., 2020).Notably, the projected recovery time could be longer.Under the post-2013 scenario, median projections for groundwater recovery, with an average of 6 years, exhibited large variance from 1 year at GW7 to 13 years at GW4.Such variance may be attributed to the local differences in recharge and hydrogeology (Bloomfield et al.,2015).Bloomfield and Marchant (2013) used a Standardized Groundwater level Index (SGI) and found that maximum groundwater droughts in England lasted for 12-73 months.Globally, droughts may endure for years and even longer, particularly in the context of anthropogenic warming ( (Bloomfield et al., 2019;Walker and Van Loon, 2023).

Wider implications
Catchments with low resilience are sensitive to changes in external perturbations such as droughts (Carey et al., 2010).In this context, both monoculture landuses showed limited resilience to external changes, as monocultures lack biodiversity and can reduce the stability of ecosystems (de Groot et al., 2021;Wright et al., 2021).Nevertheless, arable lands and conifer forests stand out as the predominant types of landuse across Central Europe.In the face of frequent reccuring drought events, it is very likely that drought "memory effects" on soil moisture occur, which may in turn influence plant phenology and biomass productivity, eventually decreasing the ecosystem resistance to subsequent droughts (Hoover et al., 2021).The results emphasize the need for landuse management strategies that enhance the resilience of landuses to drought events by increasing water retention.This may include transitioning away from monocultures and exploring alternative land management plans, such as agroforestry and mixed forests.These strategies could enhance the landscape's storage capacity, allowing for the gradual release of input water, thereby better addressing the challenges posed by climate change (Liu et al., 2018).Our findings highlight the vulnerability of groundwater resources and landscape water storage to precipitation changes and vegetation cover.They also underscore the importance of adaptive water management and conservation strategies.Policymakers and stakeholders need such an evidence base to assess the potential long-term impacts of precipitation anomalies on groundwater levels, and develop landuse strategies to ensure food and water security in increaseingly drought prone areas such the extensive lowlands of Central Europe.

Conclusion
We merged the advantages of both process-based (EcH 2 O-iso) and machine learning (NARX) models to evaluate the impact of increasingly frequent drought events in Central Europe on landscape ecohydrological flux partitioning and water storage, and to predict the recovery of groundwater under different precipitation scenarios.This hybrid approach effectively incorporates the established physical principles inherent in the process-based model, while also leveraging the more straightforward and flexible features of the machine learning model.This is especially useful for assessing future scenarios where data for process-based models may be lacking.Our results revealed that recent droughts have had a profoundly negative impact on the ecohydrology of a lowland catchment in Northern Germany.Catchment-scale groundwater recharge, soil moisture, and evapotranspiration have seen substantial reductions, especially during the main drought years of 2018 and 2022, when we observed that water fluxes across the entire catchment exhibited minimal spatial variations (across different land uses) under sustained water scarcity.However, we still observed distinct vulnerabilities in the ecohydrological pathways to drought between contrasting landuses of croplands and conifer forests, primarily due to their different subsurface storage capacities.Assuming the continuation of recent negative precipitation anomalies, our predictions for groundwater recovery is pessimistic, with one well unlikely to recover and the other three wells may require an average of 6 years to regain predrought water levels.These findings underscore the urgent need and importance of exploring alternative land management strategies to augment water storage in the landscape and resilience to drought, given the escalating frequency and severity of droughts.Our research thus provides a valuable science-base to support evaluations of drought impacts on ecohydrological processes and future land management studies, serving the needs of policy-makers and stakeholders.

Fig. 1 .
Fig. 1.Demnitzer Millcreek catchment: (a) soils, (b) topography and (c) landuse types.Soil sampling sites are shown on the soil map; groundwater and stream sampling sites are shown on the topography map.

Fig. 3 .
Fig. 3. Performance evaluation of the EcH 2 O-iso model under contrasting landuses and gauge stations.(a) precipitation; (b) and (c) discharge at two discharge gauge stations (Demnitz Mill and Demnitz); (d) and (e) soil moisture of layer 1 (top 15 cm) at forest (Forest A) and cropland sites (Alt Madlitz), respectively; (f) and (g) evapotranspiration at forest and cropland sites, respectively.The "*" after Forest A represent a nearby cell with 100 % conifer forest composition.The lines denote the median simulation results, and the shaded areas depict the bounds of ensemble simulation uncertainty.The black dots represent the observations in the calibration periods (2009-2014 & 2018-2019), and the pink dots represent the observations in the validation periods (2015-2017 & 2020-2022).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 4 .
Fig. 4. Performance evaluation of isotope simulation of EcH 2 O-iso model across different landuses.(a-c) stream isotopes at Demnitz Mill, Bruch Mill and Peat North, respectively; (d) groundwater isotopes at GW4; (e) and (f) soil water isotopes of layer 1 (top 15 cm) and layer 2 (top 15-35 cm) at forest site.The lines denote the median simulation results, and the shaded areas depict the bounds of ensemble simulation uncertainty.The black dots represent the observations in the calibration periods (2018-2019), and the pink dots represent the observations in the validation periods (2020-2022).(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 6 .
Fig. 6.(a) and percentage change (b) of daily average groundwater recharge (GWr), soil moisture in upper soil layer (SMC-L1), and evapotranspiration (ET) in two dominant landuses (i.e.conifer forests and croplands) of the study catchment and the North European Plain.Baseline represents the long-term (2009-2022) average.The percentage change was computed based on Baseline.

Fig. 7 .
Fig. 7. Forecasts of groundwater levels using a Nonlinear autoregressive exogenous (NARX) model.(a) Annual precipitation pattern divided into two periods (normal and dry years) using the Pettitt test; (b) Left panel (until 2020): Historical groundwater levels at GW4 that were used for training, early stopping and opttesting during the Bayesian optimization.Middle panel (2021-2022): Test period when the model was run ten times to account for uncertainties (more transparent lines).The solid pink line represents the median of the ten repeats.Right Panel: Projections generated by the ten tested models under various climate scenarios.The dashed lines depict trendlines corresponding to the solid lines of matching colors in the time series; (c) Time needed to return to normal-year average groundwater level under different climate scenarios.'Median', 'Lower', and 'Upper' denote the median value and the lower and upper bounds of the projections from ten repititions of the machine learning model, respectively.Empty columns = no bars indicate that the groundwater levels could not recover within the next 15 years of simulation.(For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)