Climate change effects on river droughts in Bavaria using a hydrological large ensemble

Europe and Germany have been hit by severe meteorological droughts in recent years, which have resulted in extreme low flow conditions in streams. Climate projections expect an intensification and increase in the frequency of the associated meteorological drivers. In this study, we investigate the impact of climate change on rare and extreme river droughts in a catchment with a pluvial regime, a northern tributary to the Danube (Wörnitz river) and in a river catchment with nivo-pluvial regime in the Pre-Alps (Ammer river). We employ a unique physically-based modelling chain, where the hydrological model WaSiM (Water balance Simulation Model) is driven by 50 members of the single model initial-condition large ensemble (SMILE) of the Canadian Regional Climate Model version 5 under the high-emission scenario RCP8.5. This results in a hydrological SMILE yielding 1500 simulated years for each of the investigated 30 year periods of a reference (1980–2009; REF), a current (2010–2039; CUR), a near future (2040–2069; NF), and a far future climate (2070–2099; FF). We investigate the seasonality, univariate and bivariate return periods of peak low flow and duration, and explore the climatic drivers causing the alterations.

The Wörnitz catchment shows a summer low flow regime, with climate change affecting the seasonality so that the river droughts are projected to extend further into the autumn. The typical bivariate 100 year event (REF: 7 d peak low flow = 1.96 m³ s⁻¹; event duration = 171 d) shifts to a 30 year (CUR), 17 year (NF), and 6 year (FF) event. In the reference climate of the Ammer catchment, the intensity of winter low flows dominates over summer low flows. However, the low-flow regime is projected to transition from winter-dominated to summer-dominated conditions during the current climate with more intense summer river droughts in the near and far future. While there is a tendency for cold-dry winters to cause low flow conditions in February during the reference climate, future winter low flows shift towards November/December and are triggered by the hot and dry pre-conditions of an antecedent summerly river drought. The most probable bivariate 100 year summer low flow event (REF: 7 d peak low flow = 4.9 m³ s⁻¹, event duration = 60 d) is drastically altered to occur every 34 years (CUR), 8 years (NF), and 2.5 years (FF).

In both catchments, there is an increase in the autocorrelation of peak low flows from one summer to the next, which emphasises the causality and increasing importance of lagged effects and preconditions in the course of climate change. We identify hotter and drier summer seasons as the main driver, with the positive interdependency between heat and drought in climate change intensifying, further exacerbating extremes.

Thereby, the study highlights the opportunities of a hydrological SMILE for the investigation of river droughts. Due to the large sample size of the driving climate simulations, we can robustly assess very rare events in the two catchments and generate bivariate design values narrowing down uncertainties of extreme value statistics in the light of a well characterized internal climate variability. In turn, all results are subject to scenario and model uncertainties, as the simulations are carried out with one hydrological model driven by one climate model under one emission scenario.

1 Introduction

Droughts over Europe trigger a variety of ecological, economic, and societal impacts (AghaKouchak et al., 2023; Stahl et al., 2016). Their frequency and intensity are projected to increase in a warmer climate (Naumann et al., 2021). Generally, drought can be classified into three types (van Loon et al., 2016): meteorological drought, agricultural drought, and hydrological drought referring to a precipitation deficit, soil moisture deficit, and streamflow deficit, respectively. Often, they occur in temporal order (Stephan et al., 2021), where a meteorological drought leads to a soil moisture deficit, which after some delay may induce riverine low flows. The recent European droughts in 2018–2023 (Bakke et al., 2020; Blauhut et al., 2022; Toreti et al., 2022; Toreti et al., 2023) also caused river low flow.

Riverine low flow directly impacts navigability of waterways (Meißner et al., 2017; Conradt et al., 2023), public water supply (Dahlmann et al., 2022; Stahl et al., 2016), agricultural water demand for irrigation (Naumann et al., 2021; Stahl et al., 2016; Stephan et al., 2021), industrial water demand for cooling and production (Stahl et al., 2016), hydropower (Naumann et al., 2021; Stahl et al., 2016), and it may negatively affect the water quality (Otop et al., 2023; Wolff and van Vliet, 2021). Impacts are related to the extremity of the hydrological drought event, which can be diagnosed by intensity (streamflow), duration (time below a certain streamflow threshold), and deficit (cumulative deficit during the event; Brunner et al., 2022; Brunner, 2023). The quantitative measures, how hydrological droughts are assessed cover baseflow, low flow indices, flow percentiles, volume deficit, event duration, annual flow minima, or exceedance probabilities (Brunner et al., 2022; Piniewski et al., 2022).

While trends in hydrological droughts are already detectable in the observations for some European catchments (Brunner et al., 2023; Peña-Angulo et al., 2022), future changes are to be assessed in simulation studies. Decreases in river low flow and increases in drought duration are projected for Southern and Western Europe, while the opposite development is projected for Northern, Eastern Europe and the Alps (Cammalleri et al., 2020; Forzieri et al., 2014; Marx et al., 2018; Muelchi et al., 2021; Roudier et al., 2016). In Central Europe, there is less coherence between studies (Piniewski et al., 2022).

The generation mechanisms of hydrological droughts depend strongly on catchment properties and the local climate (Brunner et al., 2022). Van Loon and van Lanen (2012) and van Loon et al. (2015) qualitatively identify eight different types of river drought generation mechanisms according to the meteorological drivers and seasonality. Brunner et al. (2022) translate this typology into a quantitative classification based on precipitation deficits, air temperature, and anomalies of snow-water-equivalent.

Most large catchments in Europe are heavily managed rather than natural (Brown et al., 2018; Hogeboom et al., 2018), and adjustment of human regulation may serve as adaptation measures to changing climate conditions (Arheimer et al., 2017; Bhaduri et al., 2016, Brunner et al., 2019a). Hence, information about future hydrological drought conditions can contribute to the dimensioning of reservoirs, dams, and channels, where management strategies rely on design values (Ehsani et al., 2017).

For such future low flow projections, hydrological models driven by climate model simulations are applied. In a common model chain, either one or an ensemble of global earth system models are employed, which may be statistically or dynamically downscaled via regional climate models. These climate simulations are run under one or more emission scenarios. Due to model-inherent biases, a bias adjustment is recommended for hydrological impact modelling (Willkofer et al., 2018). Then, one or more hydrological models are set up to simulate streamflow.

Climate simulations are governed by model uncertainty, scenario uncertainty, and internal climate variability (ICV; Bosshard et al., 2013; Lehner et al., 2020), which also propagate non-linearily to the simulated streamflow (Brunner et al., 2021). Furthermore, bias adjustment is to be applied with care (Maraun, 2016), where different methods lead to different streamflow results (Teutschbein and Seibert, 2012; Tootoonchi et al., 2022). Lastly, hydrological models differ, which results in differing streamflow simulations, when driven by the same climate (Melsen et al., 2018; Mendoza et al., 2016). Commonly, hydrological studies address one or more of these sources of uncertainty. Model uncertainties can be addressed by using ensembles of several climate and hydrological models. The limited knowledge about future emissions can be accounted for by applying a range of emission scenarios or warming levels. Single model initial-condition large ensembles (SMILE) are employed to address ICV (Deser et al., 2020). In hydrology, SMILEs are still fairly new, yet they show great potential for exploring variability in mean flows (Poschlod et al., 2020) and floods (Brunner et al., 2021; Faghih and Brissette, 2023; Willkofer et al., 2024). Hydroclimatic extremes are subject to major uncertainties due to ICV, where the SMILE provides a large sample size, which allows for the investigation of very rare events with moderate extreme value statistical uncertainties (Lang and Poschlod, 2024). Hence, in this study, we want to explore the potential of a unique hydrological 50-member SMILE for the assessment of extreme hydrological droughts under climate change. We focus on two example catchments in southern Germany reflecting the dichotomy of low flow drivers, where one catchment in the Pre-Alps is influenced by snow and rain dynamics in the reference climate (Poschlod and Daloz, 2024), while the runoff generation in the other catchment is mainly rainfall-driven. This study aims to investigate the behaviour of very rare hydrological droughts up to 1000 year return periods and the importance of ICV for low flow events. We aim to identify shifts of the dominant hydroclimatic drivers and the climate change effects on seasonality, intensity and frequency of drought events. Further, we illustrate the value of a hydrological SMILE for the detection of changes and the robust assessment of extremes. Thereby, we construct bivariate return periods of hydrological droughts, considering peak and duration of low flow events, illustrating drastic effects of climate change on such bivariate design events.

2 Study area and data

2.1 Study catchments

Rivers in southern Germany show a variety of streamflow regimes (Poschlod et al., 2020). We select the catchment of the river Wörnitz at gauge Harburg (Fig. 1), a northern tributary to the Danube. It shows a pluvial regime with the highest (lowest) monthly discharge in March (August). The majority of the 1570 km² is elevated between 400 and 550 m, with middle mountain ranges in the west extending until 700 m. The land cover consists of 69 % agricultural areas, 26 % forest, and 5 % artificial surfaces (Loritz et al., 2024). For a reference period of 1980–2009, the mean temperature amounts to 8.3 °C and the annual precipitation is 820 mm.

Figure 1 Topography and location of the two study catchments Wörnitz (Harburg) and Ammer (Weilheim) in southern Germany. Elevation is taken from the European Digital Elevation Model (EU-DEM; EEA, 2018).

Further, we choose the catchment of the river Ammer at gauge Weilheim (Fig. 1) representative for a nivo-pluvial regime with the highest (lowest) monthly discharge in April (February). It covers 608 km² with an elevation range between 550 and 2150 m above sea level. Due to the complex topography and its strong north-south gradient, the landscape can be categorized into the pre-alpine hill country and moorland in the north and the Swabian-Upper Bavarian foothills of the Alps in the southern part of the catchment (Kunstmann et al., 2006). Also, the land cover shows a gradient with agricultural areas (mainly pasture) in the north and forest in the south, resulting in a half-half distribution averaged over the entire catchment area (Ludwig et al., 2003; Marx et al., 2006). For the reference climate, the mean annual temperature is 7.0 °C and the annual precipitation amounts to 1390 mm. The Ammer is the main tributary to Lake Ammersee (Bueche et al., 2020).

The choice of the two river catchments is motivated by the fact that, despite their spatial proximity, they represent two fundamentally different climates and hydrological regimes resulting in deviating precipitation patterns, importance of snow dynamics, flow seasonality, and flood dynamics (Poschlod et al., 2020; Poschlod and Daloz, 2024; Willkofer et al., 2024). Further, both catchments are headwater catchments with little human influence, and their low flow conditions are well reproduced in the hydrological modelling (see next section).

2.2 Hydrological simulations

The hydroclimatic modelling chain was designed within the scope of the ClimEx project (Climate Change and Hydrological Extreme Events) aiming to explore hydrological extremes (Willkofer et al., 2024). For the exploration of extreme events a large homogeneous sample size is beneficial. This is achieved by applying a SMILE as driving climate. The Canadian Earth System Model version 2 large ensemble (CanESM2-LE) consists of 50 independent ensemble members simulating the global climate for 1950–2099 at 2.8° resolution under the high-emission scenario RCP8.5 (Fyfe et al., 2017). As this spatial resolution is too coarse for hydrological impact assessments, these 50 ensemble members are dynamically downscaled with the Canadian Regional Climate Model version 5 (CRCM5) at 0.11° resolution (equalling 12 km; Leduc et al., 2019). The CanESM2-LE setup stems from the CMIP5 phase, while the RCP8.5 was used as the highest-emission scenario from the Representative Concentration Pathways in the IPCC 5th assessment report. Compared to the multi-model ensemble of the coordinated downscaling experiment EURO-CORDEX under RCP8.5, the driving climate of the CRCM5-LE shows stronger signals of drying and increasing temperatures during summer (von Trentini et al., 2019). While the use of the high-emission scenario and the rather hot and dry regional climate model must be considered when interpreting the results, this climate scenario allows for the full range of warming and the exploration of low-probability/high-impact futures as well as the potential warming in the 22nd century. The choice of the CRCM5-LE as driving climate for the hydrological impact simulations is further motivated by the availability of high-resolution large ensembles over Europe, with the CRCM5-LE combining the advantages of high resolution (0.11°) and large ensemble size (50 members) compared to the other available high-resolution large ensembles described by Aalbers et al. (2018; 0.11°, 16 members) and by Brönnimann et al. (2018; 0.44°, 21 members).

For a bias adjustment and further statistical downscaling of the CRCM5 simulations, a sub-daily climate reference (SDCLIREF) data set at 3-hourly temporal resolution and 500 m spatial resolution is created based on interpolated meteorological point observations (Willkofer et al., 2020). It comprises precipitation, air temperature, relative humidity, incoming shortwave radiation, and wind speed. As the temporal and spatial coverage of sub-daily meteorological observations is not sufficient for a robust interpolation, daily observations are disaggregated to 3-hourly resolution applying the Method of Fragments (Poschlod et al., 2018; Westra et al., 2012). The native sub-daily measurements and the disaggregated data are spatially interpolated following the methodology by Rauthe et al. (2013) achieving a final resolution of 500 m. In a next step, the bias of the five meteorological variables in the CRCM5-LE is adjusted via quantile mapping (Mpelasoka and Chiew, 2009) preserving the spread of the 50 ensemble members (Willkofer et al., 2020). Furthermore, the bias-adjusted CRCM5-LE is statistically downscaled to 500 m resolution preserving mass and energy (Marke, 2008).

The hydrological modelling applies the Water balance simulation Model (WaSiM; Schulla, 2021), a spatially distributed hydrological model including modules to simulate evapotranspiration, soil water transfer, groundwater, snow accumulation and snow melt, and glaciers. While most of the modules follow physically-based process descriptions, also conceptual relationships are applied. Based on the water balance, the model computes the generation of discharge and the routing of the streams run at 3-hourly time steps (Schulla, 2021). In our setup, evapotranspiration is modelled with the Penman-Monteith approach (Monteith, 1965). The static parameters governing groundwater, soil, and land cover remain constant over the whole simulation period.

The hydrological years of 2004–2010 and 1995–2002 are used for the calibration and validation of WaSiM, respectively. There, the simulated discharge driven by SDCLIREF is compared to observed discharge, provided by the state water authority (LfU, 2025). Due to the focus of the ClimEx project on high flows, the calibration is carried out using a weighted combination of performance metrics recommended by Moriasi et al. (2007): the Nash and Sutcliffe efficiency (NSE; Nash and Sutcliffe, 1970; 50 % weight) and its logarithmic version (logNSE; 15 % weight), the Kling-Gupta efficiency (KGE; Gupta et al., 2009; 25 % weight), and the ratio of root mean squared error to standard deviation (RSR; 10 % weight; Willkofer et al., 2020). This combined objective function is optimized during the calibration process, which involves manual alternation of WaSiM parameters governing spatially distributed processes as well as automatic parameter estimation for catchment-specific parameters (Willkofer et al., 2020). The validation evaluates the performance metrics NSE, logNSE, KGE, and RSR as well as the percentage bias of mean flow. After the successful validation, the bias-adjusted and statistically downscaled CRCM5-LE simulations are used to drive the hydrological model resulting in the unique hydrological large ensemble WaSiM-LE.

For a detailed description of the setup, we refer the reader to Leduc et al. (2019) for the climate modelling, to Willkofer et al. (2018) for the bias adjustment, to Willkofer et al. (2020) for the WaSiM calibration, and to Poschlod et al. (2020) as well as to Willkofer et al. (2024) for the hydrological modelling. The choice of the objective function and the selection of the time periods for calibration and validation reflect the initial project focus of ClimEx on floods, whereas there is a second follow-up project (“ClimEx II”) with focus on droughts and low flows. This study is conducted at the transition phase from ClimEx I to ClimEx II as a preparation for the second phase in order to explore uncertainties due to climate change and internal variability to provide insights for the hydrological model setup, calibration, bias adjustment, and required ensemble size. Furthermore, the study explores the benefit from the hydrological SMILE for low flow investigation. This enables to adjust the modelling setup accordingly in the next project phase, ClimEx II. Further, we consider the results of this study to be of interest to the scientific community when designing impact studies on droughts and low flows, raising awareness of uncertainties due internal climate variability and low sample size.

For this study, we focus on the hydrological simulations between 1980–2099 for the Wörnitz and Ammer. We select 1980–2009 as reference period (REF) and divide the simulations into 2010–2039 (current climate: CUR), 2040–2069 (near future: NF), and 2070–2099 (far future: FF). With 50 ensemble members, this results in 1500 simulated years per 30 year period. Both catchments are head catchments, where no major water management structures, such as reservoirs, dams or channels, are implemented. The spatial distribution of land use and land cover is assumed to stay constant. As the simulations stemming from the ClimEx project are designed to represent floods, we have to ensure their applicability for low flow assessments. The validation results for the simulations of these two catchments are given in Table 1 indicating a sufficiently good model performance. The logNSE values above 0.7 diagnose that low flow conditions are well reproduced during the validation (Parra et al., 2019). Flow duration curves (see Fig. S1 in the Supplement) demonstrate that the hydrological model chain is able to well represent the characteristics of low flows for the two catchments in this study, while mean flows are slightly underestimated reflecting the bias reported in Table 1.

Table 1 Performance metrics of the discharge validating WaSiM simulations driven by the climatological reference SDCLIREF against discharge observations.

3 Methods

3.1 Summer and winter low flow

The German hydrological year extends from November to October, whereas for low flow also April–March is a common assessment period (LfU, 2021). Here, we divide the hydrological year into winter half (November–April) and summer half (May–October).

For the quantification of low flow intensity, we apply the lowest discharge per hydrological (half) year, considering 7 d windows (L7Q), where the river flow is averaged over a moving window of 7 d and the minimum value is chosen as L7Q. The mean of the annual L7Q indices over a longer time period (e.g. 30 years) is referred to as ML7Q. These two metrics are used operationally by the low flow monitoring (LfU, 2020). We further define the low flow event duration. Therefore, we count the number of consecutive days, where L7Q is below the ML7Q of the reference period 1980–2009. We treat events as separate if they are separated by at least one day on which L7Q > ML7Q, as the L7Q index is based on 7 d moving windows. We assign an event to the winter (summer) period, if the majority of the event duration is in the winter (summer) half year. We assess summer and winter low flows in the Ammer catchment, while only investigating summer low flows in the Wörnitz catchment due to its pronounced regime with wet winter seasons.

A comparison between discharge observations, WaSiM driven by meteorological observations and the WaSiM-LE driven by the CRCM5-LE reveals that the L7Q intensity is underestimated by 5 % to 10 % in the simulations (Fig. S2). The observed summerly event duration in the Wörnitz catchment is slightly overestimated by up to 30 d for the longest events, while observed event durations in the Ammer catchment are well within the ensemble spread for summer and winter (Fig. S3).

3.2 Extreme value statistical assessment

The flowchart in Fig. S4 provides an overview of the methodological workflow of the extreme value statistical assessment. The theoretical foundation of this framework is explained in the following sections.

3.2.1 Univariate return levels

Extreme Value Theory (EVT) is a statistical framework addressing extreme deviations from the median of probability distributions. Hence, it is applied to study the behaviour of processes, when they take unusually large or small values (Coles, 2001). The aim is to derive occurrence probabilities of rare events and their according intensity, size, duration or any other feature. In Earth system sciences, it is commonly applied to unusually large values (Hamdi et al., 2021), e.g. riverine flood discharge (Basso et al., 2021; Villarini et al., 2011), coastal flood water levels (Caruso and Marani, 2022), extreme rainfall intensities (Poschlod and Koh, 2024), extreme snow depths (Le Roux et al., 2020), or high temperature (Parey et al., 2019) but is also applicable to small values such as low flows (Raynal-Villasenor and Raynal-Gutierrez, 2014) or drought (Naumann et al., 2018).

There are two common strategies in classical EVT: Block maxima or minima sampling associated with the Generalized Extreme Value (GEV) distribution and Peak over (under) threshold sampling associated with the Generalized Pareto (GP) distribution. Note that the practical application of the GEV on minima only slightly differs from the more common application on maxima: the minima are negatively transformed by multiplication with −1 before the distribution parameters are fitted. Applying the negative transformation, the lowest value of the minima series becomes the highest value of the transformed series and vice versa. Then, the common workflow of GEV for block maxima can be applied (Hewa et al., 2007). Further note that data censoring due to zero flows has to be handled, which however does not occur in our study catchments.

For the low flow intensity, we sample the L7Q per hydrological half year and apply a negative transformation, which is per definition a block maximum sampling. However, the L7Q is not necessarily considered to be “extreme” in any half year. Hence, we allow for larger block sizes, such as two or three consecutive winter (summer) half years. Following the Fisher–Tippett–Gnedenko theorem (Fisher and Tippett, 1928; Gnedenko, 1943), the GEV distribution describes the block maxima M for a sufficiently large sample size. It approximates the exceedance probability

$$\begin{array}{}\text{(1)}& \begin{array}{rl}\mathrm{Pr}& (M\le x)\approx G(x;\mathit{\xi })\\ & =\left\{\begin{array}{ll}\exp \left(-{\left[\mathrm{1}+\mathit{\xi }\left(\frac{x-\mathit{\mu }}{\mathit{\sigma }}\right)\right]}^{-\mathrm{1}/\mathit{\xi }}\right),& \mathit{\xi }\ne \mathrm{0}\\ \exp \left(-\exp \left(-\frac{x-\mathit{\mu }}{\mathit{\sigma }}\right)\right),& \mathit{\xi }=\mathrm{0},\end{array}\right.x\in \mathbb{R},\mathit{\sigma }>\mathrm{0},\end{array}\end{array}$$

with the parameters μ (location), σ (scale) and ξ (shape). The shape parameter governs the tail behaviour of the distribution. The return level z can then be calculated by inverting Eq. (1) (Coles, 2001) and transforming the discharges back to positive values.

For the event durations, we apply a peak over threshold sampling. There, we sample all values s above a threshold u considering them as extreme values (Balkema and de Haan, 1974; Pickands, 1975). For a random variable S, the GP distribution (Davison and Smith, 1990; Martins and Stedinger, 2001) describes the exceedances $y=s-u$ with the corresponding cumulative distribution function (CDF)

$$\begin{array}{}\text{(2)}& \begin{array}{rl}\mathrm{Pr}& (S>s\mathrm{|}S>u)\approx H(y;\mathit{\xi })\\ & =\left\{\begin{array}{ll}\mathrm{1}-{\left(\mathrm{1}+\frac{\mathit{\xi }y}{\mathit{\beta }}\right)}^{-\mathrm{1}/\mathit{\xi }},& \mathit{\xi }\ne \mathrm{0},\mathit{\beta }>\mathrm{0},y>\mathrm{0}\\ \mathrm{1}-\exp \left(\frac{-y}{\mathit{\beta }}\right),& \mathit{\xi }=\mathrm{0},\mathit{\beta }>\mathrm{0},y>\mathrm{0}\end{array}\right..\end{array}\end{array}$$

Thereby, y is defined as the event duration excess over the threshold u. The event duration threshold u must be chosen suitably: selecting u too low may violate the asymptotic assumptions of the GP distribution, while choosing u too high may result in too few exceedances, leading to high estimation uncertainty and reduced statistical power. The scale parameter β describes the spread of the distribution, and the shape parameter ξ governs the tail behaviour of the GP distribution (Coles, 2001). The calculation of return levels requires to account for the exceedance rate, which is the number of events per year n_Y. Then, the N-year return level z_N is derived by

$$\begin{array}{}\text{(3)}& z_N=\left\{\begin{array}{ll}u+\frac{\mathit{\beta }}{\mathit{\xi }}\left[(N{n}_Y{\mathit{\zeta }}_u{)}^{\mathit{\xi }}-\mathrm{1}\right],& \mathit{\xi }\ne \mathrm{0},\mathit{\beta }>\mathrm{0}\\ u+\mathit{\beta }\log \left(N{n}_Y{\mathit{\zeta }}_u\right),& \mathit{\xi }=\mathrm{0},\mathit{\beta }>\mathrm{0},y>\mathrm{0}\end{array}\right.,\end{array}$$

with ${\mathit{\zeta }}_u=\mathrm{Pr}(S>u)$ (Coles, 2001).

We estimate the GEV and GP parameters in a Bayesian framework applying a Markov chain Monte Carlo algorithm (Bocharov, 2022; Foreman-Mackey et al., 2013). We check the goodness of fit via the Anderson–Darling test at the significance level of α=0.05. Within each 30 year period, we assume stationarity as simplification.

The two theoretical frameworks of the GEV and GP do not contradict each other, but rather correspond as one framework can be expressed by the other (Goda, 2011; Serinaldi and Kilsby, 2014; Wang and Holmes, 2020). Here, the choice of the respective framework is motivated by the data and the definition of intensity (L7Q) and event duration (consecutive days, where L7Q > ML7Q).

3.2.2 Bivariate return levels

For the combination of the extremal behaviour of low flow intensity and event duration, we calculate bivariate return periods. Gräler et al. (2013) provide an overview of different approaches to derive multivariate return periods in hydrology with an application to flood events. Tu et al. (2018) propose a multivariate design combination of drought properties to better capture the impact on regional water supply. Following Tu et al. (2018), we use copulas to describe the bivariate dependence structure. With I as the unite interval [0,1], the bivariate copula C is the joint distribution on $I^{\mathrm{2}}=[\mathrm{0},\mathrm{1}{]}^{\mathrm{2}}$ with uniform margins (Zscheischler and Fischer, 2020). Hence, for two random variables X and Y with the CDFs $F_X\left(x\right)=\mathrm{Pr}(X\le x)$ and $F_Y\left(y\right)=\mathrm{Pr}(Y\le y)$, the copula C expresses the joint distribution function as (Sklar, 1959; Zscheischler and Seneviratne, 2017)

$$\begin{array}{}\text{(4)}& F(x,y)=\mathrm{Pr}(X\le x,Y\le y)=C\left(F_X\right(x),F_Y(y\left)\right),\end{array}$$

with the marginal CDFs $F_X\left(x\right)=\mathrm{Pr}(X\le x)$ and $F_Y\left(y\right)=\mathrm{Pr}(Y\le y)$ Based on the copula fitted to the data, one can directly define four “hazard scenarios”, where Salvadori et al. (2016) provide an intuitive graphical representation (Fig. 1 therein). Here, we explore bivariate exceedance probabilities with the joint survival distribution

$$\begin{array}{}\text{(5)}& \stackrel{\mathrm{‾}}{F}(x,y)=\widehat{C}\left({\stackrel{\mathrm{‾}}{F}}_X\right(x),{\stackrel{\mathrm{‾}}{F}}_Y(y\left)\right),\end{array}$$

following Salvadori et al. (2013) with the marginal survival functions ${\stackrel{\mathrm{‾}}{F}}_X=\mathrm{1}-F_X$ and ${\stackrel{\mathrm{‾}}{F}}_Y=\mathrm{1}-F_Y$ and the survival copula $\widehat{C}$ (AghaKouchak et al., 2014; Zscheischler et al., 2017). For any given $(x,y)\in R^{\mathrm{2}}$, there is a unique survival critical layer (or isoline), on which the set of realizations of X and Y have the same probability $t\in (\mathrm{0},\mathrm{1})$: $L_t^{\stackrel{\mathrm{‾}}{F}}=\mathit{\left\{}\right(x,y)\in R^{\mathrm{2}}:\stackrel{\mathrm{‾}}{F}(x,y)=t\mathit{\}}$, where $L_t^{\stackrel{\mathrm{‾}}{F}}$ is the survival critical layer associated with the probability t (AghaKouchak et al., 2014; Salvadori et al., 2011).

$\stackrel{\mathrm{‾}}{K}$ is the Kendall's survival function associated with $\stackrel{\mathrm{‾}}{F}$ and the probability t defined as

$$\begin{array}{}\text{(6)}& \stackrel{\mathrm{‾}}{K}\left(t\right)=\mathrm{Pr}\left(\stackrel{\mathrm{‾}}{F}\right(X,Y)\ge t)=\mathrm{Pr}\left(\widehat{C}\right({\stackrel{\mathrm{‾}}{F}}_X\left(x\right),{\stackrel{\mathrm{‾}}{F}}_Y\left(y\right))\ge t).\end{array}$$

With an average interarrival time of 1 for X and Y, the survival Kendall's return period ${\stackrel{\mathrm{‾}}{\mathit{\kappa }}}_{XY}$ is calculated as (AghaKouchak et al., 2014)

$$\begin{array}{}\text{(7)}& {\stackrel{\mathrm{‾}}{\mathit{\kappa }}}_{XY}={\displaystyle \frac{\mathrm{1}}{\mathrm{1}-\stackrel{\mathrm{‾}}{K}\left(t\right)}}.\end{array}$$

The bivariate survival Kendall's return period has been applied by AghaKouchak et al. (2014) in order to analyse the 2014 drought in California and by Zscheischler et al. (2017) in order to explain crop yield variability in Europe based on temperature and precipitation.

In this study, we use the Multivariate Copula Analysis Toolbox (MvCAT; Sadegh et al., 2017) to fit the copula as it considers 26 different bivariate copula families. The copula parameters are estimated via a gradient-based local optimization (Sadegh et al., 2017), and the best copula is chosen according to the Bayesian Information Criterion (BIC, Schwarz, 1978), which measures the trade-off between model fit and model complexity.

3.2.3 Bivariate design events

For a given probability t, there are numerous bivariate combinations on survival the critical layer $L_t^{\stackrel{\mathrm{‾}}{F}}$, which are associated with the same bivariate survival Kendall's return period. For the selection of a bivariate design event, a common approach is to choose the “most likely design scenario” (Gräler et al., 2013; Sadegh et al., 2018). Therefore, weighted samples are drawn from the critical layer $L_t^{\stackrel{\mathrm{‾}}{F}}$, with the copula density function as weight (Sadegh et al., 2018). The bivariate combination with the highest density is then chosen as the most likely design scenario. Note that this design pair of peak low flow and low flow duration is not necessarily the most impactful realization. However, it can be interpreted as the “typical” bivariate realization, where the water management should be designed accordingly to withstand this design scenario for a given return period (Salvadori et al., 2014).

3.3 The Budyko framework

Budyko (1951) found that climate aridity governs the long-term partitioning of precipitation into river discharge and evaporation formalizing the ratio of mean annual evaporation E and mean annual precipitation P as:

$$\begin{array}{}\text{(8)}& {\displaystyle \frac{E}{P}}={\left[{\displaystyle \frac{E_{\mathrm{p}}}{P}}\tanh \left({\displaystyle \frac{P}{E_{\mathrm{p}}}}\right)\left(\mathrm{1}-e^{-\frac{E_{\mathrm{p}}}{P}}\right)\right]}^{\mathrm{0.5}},\end{array}$$

with the mean annual potential evaporation E_p. The framework is based on an empirical foundation, where the resulting Budyko curve is found to capture around 90 % of the inter-catchment variability in evaporative fractions across various catchments in Europe (Budyko, 1951). We refer the reader to Berghuijs et al. (2020) for an insightful commentary on the Budyko framework.

4 Results

4.1 Climatic and hydrological drivers

River discharge is driven by catchment-specific features, human management and the climatic conditions. As the first two factors are kept constant in this simulation experiment, the changes in river discharge and low flow conditions can be attributed to climate change. Figure 2 presents the seasonal course of air temperature, precipitation, evaporation, soil moisture in the rootzone, snow, and snow storage for the reference climate and future projections. As WaSiM is a spatially distributed model run at 500 m resolution, the climatic and hydrological drivers are averaged over the two catchment areas upstream of the river gauges. Compared to the reference climate, the annual temperature in the Wörnitz catchment is projected to increase by +1.0 °C (CUR), +2.3 °C (NF), and +4 °C (FF), respectively. In the Ammer catchment, slightly higher increases are simulated with +1.1 °C (CUR), +2.5 °C (NF), and +4.3 °C (FF). This leads to an increase of evaporation, especially in the summer season, which in turn induces lower soil moisture conditions. The annual precipitation sums change by less than 6 % in the future periods for both catchments. However, the seasonality shifts, where precipitation is projected to increase in the winter and decrease in the summer.

Figure 2 Climatic and hydrological drivers in the (a) Wörnitz and (b) Ammer catchments. The centerline of the boxplots shows the median over all 50 members and 30 years per period. The boxes represent the inter-quartile range of the 1500 years, where the whiskers extend to the 5th and 95th percentiles, respectively. An alternative figure focussing on the changes compared to the reference period is given in the Supplement (Fig. S5).

Due to the temperature increase and the summerly precipitation decrease, soil moisture conditions during summer are projected to be drier in the future. This depletion of soil moisture is even more pronounced for drought years (see the lower whiskers for the 5th percentile). In the Ammer catchment, the median soil moisture conditions are close to saturation over the whole year in the reference climate. For the near and far future, the soils are projected to dry out stronger during summer. In the Wörnitz catchment, the soil moisture is projected to decrease in summer and shift the seasonality towards drier soils later in the year. Snow dynamics play a minor role in the Wörnitz catchment, with a projected decease of snowfall. In the Ammer catchment, snow storage reaches over 50 mm snow water equivalent during February and March of the reference climate. Snowfall and snow storage are projected to strongly decrease.

4.2 Drying and Seasonality, duration

The alterations of the hydroclimatic conditions affect the riverine flow regime and low flow intensities. For the annual L7Q with no separation of summer and winter low flow, Fig. 3 shows the median drying tendencies with an earlier onset for Wörnitz than for Ammer. The annual median over 50 members results in a relatively smooth temporal progression despite using an extreme metric (L7Q).

Figure 3 “Drying stripes” (inspired by the warming stripes by Hawkins, 2018): With the ML7Q of 1980–2009 as reference, the temporal course of the annual L7Q visualizes the projections of low flow intensity. Each annual stripe pools all 50 members and represents the median of these 50 years.

Summer events dominate the annual L7Q in the Wörnitz catchment, with constant shifts towards later in the year with future warming (Fig. 4). In the Ammer catchment, the low flow events in the reference period are mostly in winter with most events in February. In the current climate, the simulations show a broad distribution of events throughout the year, while summer events dominate in the two future periods.

Figure 4 Seasonality of the annual L7Q events. The violin plots provide kernel density estimates, where the grey line shows the median of all events in the respective 30 year period. Per time period, all 30 years times 50 members are pooled, allowing the kernel density estimate to be based on 1500 years.

While Figs. 3 and 4 focus on the peak low flow intensity and its date, the event duration and its seasonality are also affected by a changing climate. For eight categories of event duration (see Table 2), Fig. 5 shows the event seasonality for summer low flows in the Wörnitz catchment as well as summer and winter low flows in the Ammer catchment.

Table 2 Categories of event duration visualized in Fig. 5.

Figure 5 Seasonality and event duration for summer low flow events in the (a) Wörnitz catchment during summer and Ammer catchment during (b) summer and (c) winter. Per time period, all events from the 30 years times 50 members are pooled. The density colour code depicts the occurrence of low flow events within the simulated 1500 years. Inset Numbers (1–8) indicate the event duration for 8 categories ranging from category 1 (<30 d) to category 8 (>211 d). See Table 2 for details.

The summer low flow events are projected to extend towards longer event durations and to shift later in the year (Fig. 5a and b). For the Wörnitz, the typical low flow events (highest density in Fig. 5a) in the reference period cover 31–120 d centred around August and September. In the far future period, prolonged events of 91–180 d centred around August to October are most common. In the Ammer catchment during the reference period, typical summer low flows cover 1–90 d at any time between May and November (Fig. 5b). Until the far future, they are projected to occur mainly from July to November and extend to durations ranging from 31–120 d. Winter low flows are shorter and typically last 1–60 d at any time between November and March.

4.3 Return levels of low flow intensity and duration

The L7Q indices are well represented by the GEV distribution (Fig. 6) passing the goodness-of-fit test for all time periods and catchments. Due to the large sample size of 50 members times 30 annual indices, the extreme value statistical uncertainties are small. The L7Q is projected to intensify during summer for both catchments and across all return periods. For the winter low flow in the Ammer catchment, the intensification is less pronounced. The current climate only slightly differs from the reference period. The near future shows a lower L7Q for 5 year and rarer events. In the far future, a strong decline of the low flow discharge is projected across all event probabilities.

Figure 6 Univariate return levels of 7 d low-flow discharge peaks L7Q (a, c, e) and event duration (b, d, f) for both study catchments. Per time period, all 50 members are pooled to a total of 1500 years. The coloured dots are plotted via Weibull plotting positions.

For the event duration, the GP distribution is applicable above a duration threshold, which varies over the catchments, seasons and time periods. The threshold is in the range of the respective 2 to 15 year return periods (see Fig. 6b, d, and f; where the line starts). Due to the sample size of the WaSiM-LE, the extreme value statistical evaluation shows relatively small uncertainties. The width of the confidence intervals depends on the sample size above the chosen threshold.

During summer, event durations are increasing in a warmer climate for both catchments. It is notable, that the GP distributions show a light tail (ξ<0) in the Wörnitz catchment, indicating that there is an upper limit for the duration of very rare events caused by the interruption of low flow events in the wet months of February and March (Fig. 5a). Three periods (CUR, NF, FF) show similar durations for very rare events above the 100 year return period (Fig. 6b). However, the moderate extremes up to 5 year events show a strong increase with warming.

In the Ammer catchment, the event durations increase across the different return periods (Fig. 6d), where three periods (REF, CUR, NF) show a slightly heavy tailed behaviour (ξ>0). With the strong warming until FF, the durations increase, but the distribution shifts towards a light tail. During winter, no clear changes are detected, only a tendency for longer events in the far future period.

The bivariate distribution of L7Q and event duration reveals a high dependence between both event characteristics during the summer season (Fig. 7). All possible realizations of bivariate 10 year and 100 year events are given by the grey isolines, where the black star marker represents the most likely combination. Hence, a typical 10 year summer low flow event in the Wörnitz catchment in the reference period (L7Q=2.44 m³ s⁻¹; event duration = 95 d) is projected to occur every 5.2 a (CUR), 3.0 a (NF), and 1.7 a (FF). The typical 100 year event (L7Q=1.96 m³ s⁻¹; event duration = 171 d) shifts to a 30 year (CUR), 17 year (NF), and 6 year (FF) event.

Figure 7 Bivariate survival Kendall's return levels of event duration and L7Q. The grey isolines show the 10 year and 100 year return levels, where the black star marker represents the most likely combination of L7Q and event duration for the respective bivariate return periods. For a direct visual comparison of the changes in the isolines between the time periods, Fig. S6 is provided.

In the Ammer catchment, the projected changes during summer are even more drastic. The most probable 10 year event is characterized by a 7 d peak low flow of 6.2 m³ s⁻¹ and a duration of 24 d. Those characteristics are reached every 6 years (CUR), 2.5 years (NF), and 1.3 years (FF). This indicates that a moderate extreme event of the reference period is projected to become the “new normal” until 2070–2099 under high greenhouse gas emissions. The 100 year event (L7Q=4.9 m³ s⁻¹, event duration = 60 d) is altered to occur every 34 years (CUR), 8 years (NF), and 2.5 years (FF).

In the winter, the dependence between event intensity and duration is less pronounced. The most likely bivariate characteristics of a 10 year event during the reference period show similar values as the summerly 10 year event with 6.4 m³ s⁻¹ 7 d peak low flow and 18 d. The 100 year event is only slightly more intense (5.8 m³ s⁻¹) extending to 43 d. Due to the differing dependence structure, the bivariate return level values deviate stronger from their univariate counterparts of the marginal distributions. Under warmer climate conditions, the 10 year (100 year) events of the reference period occur every 15 (92) years in the current climate, 12 (62) years in the near future, and 5 (20) in the far future, respectively.

4.4 Summer low flow

In a warmer climate, the results show more intense peak low flows, longer event durations and an extension of events towards later in the year for both study catchments. We refer those changes to the alterations of the hydroclimatic drivers (Fig. 2). In the Ammer catchment during the reference period, the rank correlation between the bivariate event rarity and hydroclimatic drivers (averaged over the summer half year) shows a weak correlation for temperature, moderate correlations for evaporation, precipitation, and the climatic water balance, and strong correlation for the soil moisture (Fig. 8b). The strength of the rank correlation is projected to increase in a warmer climate, with strong correlation for all drivers in the far future. We reason that events prolong in a warmer climate covering larger parts of the summer half year (Fig. 5b), which is why the dependence with the governing hydroclimatic drivers increases. Secondly, we find an increase of the autocorrelation of L7Q comparing the previous and following summer season (Fig. 9b).

Figure 8 Rank correlation between climatic and hydrological drivers and the bivariate event rarity in the summer half year for the (a) Wörnitz catchment and (b) Ammer catchment. The colours refer to the respective time periods as in the previous Figs. 2, 4, and 6. Per time period, all 50 members are pooled to a total of 1500 years. Drivers are averaged over the summer half year. The climatic water balance is defined as precipitation minus potential evaporation. The transparent ranges provide 95 % confidence intervals for Kendall's τ.

Figure 9 Lag-1-autocorrelation for L7Q and low flow duration in the summer half year for the (a) Wörnitz catchment and (b) Ammer catchment. The colours refer to the respective time periods as in the previous Figs. 2, 4, and 6. Per time period, all 50 members are pooled to a total of 1500 years. The lag relates to the previous summer half year. The coloured transparent ranges provide 95 % confidence intervals for the autocorrelation, while the grey transparent range represents no significant autocorrelation (α=0.05).

This indicates that the pre-conditions of the previous year gain influence on the conditions in the following year with the potential for multi-year river droughts highlighting the importance of lagged effects under climate change for the assessment of drought risks (Bevacqua et al., 2024). Furthermore, we find an increase of inter-driver dependence, notably for temperature and precipitation in both catchments. Figure 10 illustrates that hotter (colder) summers are more likely to be drier (wetter) in the future periods. This stronger dependence has also been found for compound hot and dry summers during June, July, and August over Europe (Felsche et al., 2024) propagating to soil moisture droughts (Boehnisch et al., 2025).

Figure 10 Bivariate scatter plots of summer (May to October) temperature and precipitation. The scatter plot shows one dot per year, adding up to a total of 1500 dots pooling the 50 members. The color scale connects the climatic conditions to the respective summer L7Q.

In the Wörnitz catchment, a positive autocorrelation of L7Q is already significant for the reference climate and is projected to further increase in the near and far future (Fig. 9a). We also find increases of the dependence between the bivariate event rarity and hydroclimatic drivers (Fig. 8a), where the changes are however less pronounced than in the Ammer catchment. The rank correlation of bivariate event rarity and evaporation drops from a moderate level in the reference and current climate to almost zero in the far future (Fig. 8a). We argue that under the high emission scenario, the water availability for evaporation is limited especially during dry summers, where the potential evaporation largely exceeds the precipitation. This shift is illustrated by Fig. 11 within the Budyko framework. Here, we find that both catchments move along the Budyko curve, towards higher evaporative fractions. In the far future, the climate aridity of the Wörnitz catchment shows a spread around 1.0, indicating that it is projected to shift to arid conditions until the far future.

Figure 11 Evaporative fractions in the Budyko framework of both study catchments for four climate periods. Each dot refers to the average conditions of one member per 30 year climate period, where the spread of the dots represents the range of internal climate variability between the 50 members.

4.5 Winter low flow

River droughts in cold climates can be caused by different drivers (van Loon et al., 2015). Van Loon and van Lanen (2012) have identified eight categories, where three subtypes of cold snow season droughts are governed by cold temperature anomalies. Cold temperatures can either induce an earlier onset of the snow season, or extend the season with later snowmelt. In mild winter climates, a cold anomaly can affect the water availability during the whole season, as precipitation falls as snow preventing from infiltration to the groundwater and runoff to the river (van Loon et al., 2015). Warm snow season droughts are governed by temperature and precipitation. In harsh climates, an earlier than normal snowmelt peak can induce a drought afterwards. In milder winter climates, joint negative snowmelt anomalies and rainfall deficit can cause low flows. Snowmelt droughts and glaciermelt droughts are induced by negative anomalies of meltflow during the respective seasons. Lastly, a rain-to-snow season drought develops as rainfall deficit during the summer and autumn, which is then acerbated by a rather sudden temperature drop during winter, which does not allow the system to recover (van Loon et al., 2015). Compared to the investigated catchments in cold climates by van Loon et al. (2015), the reference climate of the Ammer would be a mild winter climate with a mean annual air temperature of 7 °C without glacier influence. Even though the catchment is considered a mild winter climate compared to cold climates, the seasonality of low flows and their intensity shows that winter low flows dominate in the reference climate (Figs. 4 and 6). However, there is a shift towards dominating summer events in the future climates. Such shifts are also projected for nivo-pluvial catchments in Austria under RCP8.5 until end of the century (Hanus et al., 2021).

Hence, during the reference climate, the Ammer catchment is mostly prone to rain-to-snow season droughts, cold snow season droughts with negative temperature anomalies during the winter season, and warm snow season droughts with joint negative snowmelt anomalies and rainfall deficit. Due to the mix of potential drivers for these different river drought types, a correlation analysis between hydroclimatic drivers averaged over the whole winter season and the drought metric will be less conclusive than for the summer events (Haslinger et al., 2014). As the dependence structure between L7Q and event duration is weaker in winter than in summer (see Fig. 7b and c), a correlation analysis is conducted separately for L7Q and event duration instead of a single analysis for the bivariate event rarity.

There, a weak positive rank correlation between temperature and L7Q as well as precipitation and L7Q, and a negative correlation between temperature and event duration as well as precipitation and event duration in the reference climate indicate that there is a tendency for colder and drier winters intensifying and prolonging river droughts (Fig. 12). The sign of correlation is projected to switch in the future periods, where warmer winters weakly correlate with more extreme drought events. We interpret this as a projected decline of cold snow season droughts in a warmer climate. Furthermore, the same analysis for the hydroclimatic drivers of the antecedent summer (Fig. 13) reveals that there is a strong projected increase in the correlation of all antecedent drivers and the following winter L7Q. We refer this behaviour to a higher number of L7Q events, which are preconditioned by intense summer river droughts. Either these river droughts start in summer and extend over the winter, so that the longer part is in the winter season and the L7Q event is counted as a “winter L7Q event”. Or an intense summer drought extends into the winter, with the longer part in the summer, which is why it would not be counted as a “winter event”. However, this event can be interrupted by a wet spell, which increases the L7Q above the ML7Q, followed by another dry spell, causing the discharge to decrease below the ML7Q as the system has not recovered yet. Then, the second part of the event counts as “winter event”. Both described types of the river drought progressions explain the increasing importance of antecedent summer conditions for the winter L7Q. However, for the winter event duration the antecedent driver conditions are weakly correlated (Fig. 13b). There, the moisture conditions during winter govern the event duration (Fig. 12b). Generally, winter event durations stay shorter than their summerly counterparts. In the reference and current climate, snowmelt and rainfall lead to a recovery until April for almost all years (see Figs. 2 and 5c). This also applies for the future climates, where the combination of higher rainfall and low evaporation ends eventual drought events in April (see Figs. 2 and 5c).

Figure 12 Kendall's rank correlation between winter climatic drivers and (a) winter low-flow intensity L7Q, as well as between winter climatic drivers and (b) winter low-flow event duration. The colours refer to the respective time periods as in the previous Figs. 2, 4, and 6. Per time period, all 50 members are pooled to a total of 1500 years. The transparent ranges provide 95 % confidence intervals for Kendall's τ.

Figure 13 Kendall's rank correlation between climatic drivers of the antecedent summer and (a) winter low-flow intensity L7Q, as well as between climatic drivers of the antecedent summer and (b) winter low-flow event duration. The colours refer to the respective time periods as in the previous Figs. 2, 4, and 6. Per time period, all 50 members are pooled to a total of 1500 years. The transparent ranges provide 95 % confidence intervals for Kendall's τ.

5 Discussion

5.1 Comparison to existing scientific literature

Mauser et al. (2008) use a statistical climate simulator to generate climate projections under the medium-emission scenario A1B driving the physically-based hydrological model PROMET in the catchments of the Upper Danube. Comparing 1971–2003 to 2011–2060, they find a strong intensification of river low flows, where 50 year peak low flow is less than half of the reference until 2060. Huang et al. (2015) apply a multi-model RCM ensemble featuring 16 simulations (12 simulations from the project ENSEMBLES and 4 additional simulations) under the medium-emission scenario A1B to drive the hydrological model SWIM over central Europe. For the catchment of the Upper Danube, they find that the 16 simulations do not agree on changes in 50 year river droughts until 2021–2060 and 2061–2100. However, more than 60 % of the simulations suggest an intensification of river droughts in the non-alpine catchments in the far future. Building on the 12-member multi-model ensemble from the project ENSEMBLES under the A1B emission scenario, Forzieri et al. (2014) apply the hydrological model LISFLOOD to investigate climate change effects on the annual minimum flow of 446 river gauges in Europe. Across Europe, they find a pronounced spatial gradient for 20 year low flow, with intensification of river droughts in the south and west and less extreme peak low flow in the north to north-east. For the catchments of the Upper Danube, increases and decreases of river droughts are projected, with the tendencies for intensification in the southern parts.

Stanzel and Kling (2018) compare the effect of the driving climate of the ENSEMBLES project (21 simulations) and the EURO-CORDEX ensemble (16 simulations) under medium (RCP4.5) and high (RCP8.5) emission scenarios in the Upper Danube catchment. They apply the conceptual hydrological model COSERO and find no significant changes in the 10th percentile of streamflow until 2021–2050. In the far future of 2071–2100, they report a decrease under RCP8.5 and based on the ENSEMBLES simulations, while no changes are simulated for RCP4.5.

Probst and Mauser (2022) select three simulations from EURO-CORDEX under RCP2.6 and RCP8.5 to drive the hydrological model PROMET. They find less extreme low flows in the future for the Upper Danube (gauge Achleiten) for all scenarios, where the gauge in Achleiten drains the whole Upper Danube catchment and is also strongly governed by the Inn. Hence, the comparability to the small catchments in this study is limited.

Marx et al. (2018) devise a multi-model study over Europe featuring five global climate models under three emission scenarios and three hydrological models equalling to 45 hydrological simulations. They analyse the 10th percentile of streamflow as low flow for the future warming levels of +1.5, +2, and +3 K. Compared to the baseline of 1971–2000, they find no significant changes in low flow in the Upper Danube for the catchments north of the Danube, while the projections suggest less extreme river droughts in alpine and pre-alpine river catchments south of the Danube.

Kiesel et al. (2019) apply the hydrological model Soil and Water Assessment Tool (SWAT) driven by 16 climate projections from EURO-CORDEX and explore different bias adjustment approaches. They find a 32 % decrease of the L7Q in the Ammer catchment until the end of the century under RCP8.5 compared to 2010. They also detect the shift of the dominating low flow regime from winter low flow to summer low flow, however occurring between the 2050s and the 2090s.

In sum, the existing scientific literature on river droughts in the Upper Danube area is ambiguous as the review of Piniewski et al. (2022) follows for the whole of central Europe. Our study tends to show a rather strong signal of low flow intensification compared to other studies in the Upper Danube. However, the selection of river catchments, low flow metrics, the differentiation between winter and summer low flow, driving climate models, emissions scenarios, time periods, and hydrological models limits the comparability of the studies.

5.2 Limitations and benefits of the hydrological SMILE for river drought assessments

Uncertainties in future climate projections stem from three specific sources: (1) climate model uncertainty, (2) emission scenario uncertainty, and (3) internal climate variability (ICV; Lehner et al., 2020). For projections of streamflow, further sources of uncertainty from add to the overall uncertainty, such as the bias adjustment (Willkofer et al., 2018), hydrological model calibration (Herrera et al., 2022), and the model uncertainty of hydrological models (Thébault et al., 2024). In our study, we apply only one climate model under one emission scenario, one bias adjustment methodology, one hydrological model calibrated by one algorithm on one calibration period. However, the SMILE as driving climate allows us to quantify ICV (Deser et al., 2020; Maher et al., 2021) and its large sample size enables a robust assessment of extreme low flow events.

The RCP 8.5 emission scenario represents high-emission scenario (Riahi et al., 2011). Following Hausfather and Peters (2020) we would recommend to see this scenario as a worst-case estimate, rather than the most likely scenario. While the RCP 8.5 tracks the global cumulative emissions well until 2020 (Schwalm et al., 2020), the IPCC AR6 report stated that the likelihood of such a high-emission scenario is considered low due to recent developments in the energy sector (Chen et al., 2021). Sarofim et al. (2024) still stress the utility of high-emission scenarios to assess low-probability/high-impact futures and to represent the climate of the 22nd century. In consequence, it is the only scenario to allow for an analysis of extreme event dynamics for the full range of global warming levels.

Compared to the multi-model ensemble of EURO-CORDEX under RCP8.5, the CRCM5-LE shows stronger signals of drying and increasing temperatures during summer. During winter, the CRCM5 projects a stronger increase of precipitation, where the temperature is close to the mean of the 22 EURO-CORDEX models (von Trentini et al., 2019). Hence, the driving climate of the CRCM5-LE induces a stronger intensification of summer low flow events and a less intense increase of winter low flows in comparison to the EURO-CORDEX. In turn, EURO-CORDEX simulations are found to underestimate the historical and projected warming compared to their driving global climate models of the CMIP5 phase (Schumacher et al., 2024). Furthermore, CMIP6 models have a wider range of equilibrium climate sensitivity than CMIP5 also resulting in a wider range of warming responses (Meehl et al., 2020).

With their multi-model study, Marx et al. (2018) can assess the contribution of the driving climate and the hydrological modelling to the overall uncertainty within low flows. They find a dominating contribution of the driving climate in large parts of central Europe, whereas hydrological modelling governs uncertainty in the (pre-)alpine catchments. Parajka et al. (2016) differentiate between summer low flow and winter low flow regimes for Austrian catchments. They find that both low flow regimes show a varying sensitivity to different sources of uncertainty. Investigating different climate scenarios, calibration periods, and calibration objective functions, they find that the driving climate is the dominating source of uncertainty for summer low flows contributing 76 % of overall uncertainty. For winter low flows, the driving climate is also the largest contributor (26 %), accompanied by the selection of the calibration period (13 %) and the objective function (13 %). Willkofer et al. (2018) show that the choice of the bias adjustment method has a moderate effect on projections of low flow, which ranges up to 10 % for a 2 year low flow.

Troin et al. (2022) assess how the hydrological model structure in the representation of potential evapotranspiration, snow dynamics, vertical and lateral flow can affect low flow simulations. They find that the lateral and vertical flow schemes are most important sources of uncertainty for low flows. This is because aquifers contribute to the streamflow via the baseflow (Beck et al., 2013), which is particularly relevant during periods of low precipitation and low surface runoff. As the upward movement of groundwater can influence soil-moisture-climate interactions, depending on the soil type (Vogelbacher et al., 2024; Shokri and Salvucci, 2011), integrating information on soil characteristics, groundwater level, as well as capillary rise and its contribution to soil moisture within the hydro-climatological model chain is particularly valuable for future climate projections. Shallow groundwater tables can enhance soil moisture storage, potentially mitigating the intensity and frequency of droughts (Van Loon and Van Lanen, 2012; Van Loon et al., 2024). In the setup of the study, the hydrological model WaSiM is able to represent the aforementioned processes, as it captures infiltration, surface runoff and water flow in the unsaturated zone using the Richards equation as well as deep percolation and capillary rise of water from the saturated to the unsaturated zone (Schulla, 2021). Lateral groundwater flow is represented by calculating the baseflow conceptually based on the saturated hydraulic conductivity, elevation of the water table and the soil surface (Neukum and Azzam, 2012).

In addition, changes in land cover and land use may further influence these dynamical interactions between groundwater, soil moisture and climate, whereby our model setup assumes static land use and land cover. Teuling et al. (2010) analyse the differing responses of grasslands and forests, emphasizing the mitigating effects of forests on the persistence of increased temperatures for long-lasting events. Zipper et al. (2019) suggest that land use and land cover changes after the Soviet Union collapse and European Union formation might have induced drying in the Mediterranean and Eastern Europe regions. For a detailed analysis of groundwater, soil moisture, and climate interactions, a fully coupled model chain (e.g. Fersch et al., 2020; Maxwell and Kollet, 2008; Naz et al., 2023; Rummler et al., 2022) as well as dynamic scenarios of land cover and land use (Huber García et al., 2018) would be beneficial. However, in this study, the climate simulations by the CRCM5-LE drive the hydrological model WaSiM without feedback from the hydrological conditions into the climate modelling. Given the scope of this analysis on large ensembles, a fully coupled model chain featuring dynamic land use and land cover was not incorporated. Nevertheless, we emphasize potential changes in land use and land cover as well as the potential contribution of groundwater to soil moisture and the interaction with the climate in the interpretation of the results.

In addition to the sources of uncertainty of the scenario uncertainty and various components of the model uncertainty of the climate models and hydrological models discussed above, ICV affects future climate and streamflow projections. Lehner et al. (2020) compare the contribution of internal variability to overall uncertainty for mean temperature and mean precipitation globally and for different regions. They demonstrate that ICV influences projections of precipitation more than of temperature and that it can be the dominating source of uncertainty for small regions and local assessments. For regional to local extreme events (Poschlod, 2022) and compound events (Bevacqua et al., 2023; Lehner, 2024), the contribution of ICV is even higher than for mean values. The propagation of this uncertainty source to hydrology is less investigated. Using the WaSiM-LE, Willkofer et al. (2024) report a range of −67 % (−51 %) to +171 % (+91 %), when estimating 100 year flood peaks in 98 catchments in central Europe based on 30 (100) years. Also based on a hydrological SMILE, Faghih and Brissette (2023) confirm the strong influence of ICV on flood projections under climate change for 133 catchments across the eastern and north-eastern United States. Further, they emphasize that floods in small catchments show a higher sensitivity to ICV than floods in large catchments. For future projections of low flow peaks, they also find a considerable influence of ICV, however with a more even sensitivity across catchment sizes.

In our study, we focus less on quantifying the ICV influence, but rather pool the 50 members of the hydro-SMILE to generate a robust database for analysing very rare and extreme low flows in terms of their univariate (Fig. 6) and bivariate (Fig. 7) characteristics. If we were to use only a single 30 year simulation per time period instead of the large ensemble, we would not be able to detect the climate change-induced changes as robustly. When resampling 30 years per time period with replacement for 100 times, the 100 year low flow intensities and duration are systematically underestimated (Fig. S7). The inner 90-percent-range of estimated 100 year L7Q spans from −21 % to +32 %, whereas uncertainty is larger with −62 % to +81 % for 100 year durations. While the “true” sign of change is detected, with the pooled ensemble representing the “truth”, only a small partition of the bootstrapped samples can reveal the significance of future changes. Also, for the bivariate return levels (Fig. 7) or trivariate assessments of temperature, precipitation and low flow (Fig. 10) the large sample size of the ensemble is paramount.

5.3 Implications of the results for the local hydrology and water management

The projected changes in river drought characteristics for the Wörnitz and Ammer rivers are expected to have substantial implications for the water availability. While both rivers are not contributing to public water supply, major shipping infrastructure or cooling of thermal powerplants, the projected river droughts would affect hydropower generation, agricultural irrigation, recreational activities, and riverine ecosystems, particularly in terms of water temperature and quality (FGG Donau, 2021; LfU, 2016; WWA Donauwörth, 2013). The drastic scale of the projected changes would put water management and adaptation planning under intense pressure to act. However, due to the reliance on a single combination of global climate model, regional climate model, and hydrological model, as well as a single emission scenario, the results are subject to considerable model and scenario uncertainty. Given that RCP8.5 represents the highest-emission scenario and that the CRCM5 is positioned within the hotter and drier range relative to the multi-model EURO-CORDEX ensemble, the projected drastic changes for the period 2070–2099 can be interpreted as a worst-case scenario by the end of the century. Nevertheless, significant alterations in river drought frequency and severity are projected for the period 2040–2069, corresponding to a global warming of slightly above 3 °C above preindustrial levels, which aligns with the middle-of-the-road SSP4.5 scenario for 2100. Hence, we argue that the low flow events that are simulated for the period 2040–2069 are in the range of possible events occurring until the end of the century. We recommend that water management and adaptation planning consider the projected intensity and duration of these rare events for a well-informed and anticipatory drought management.

6 Conclusion

Based on the above discussion, we do not recommend a hydro-SMILE as replacement of multi-model assessments, but for the focus on very rare events of multivariate extreme event properties. Laaha et al. (2016) suggest a three-pillar approach for the projections of low flow. They recommend extrapolation of observational trends, hydrological modelling driven by climate simulations under emission scenarios, and extrapolation of changing stochastic rainfall characteristics into the future combined with rainfall–runoff modelling. They consider the consistency between the three approaches as an indicator of the uncertainty of the future projections. Within their third pillar, they apply a stochastic rainfall generator to create an ensemble of 100 runs, which are intended to capture climate variability. We argue, that a process-based hydrological SMILE could either enrich the third pillar or form a separate fourth pillar, if there is a focus on very rare and extreme events or multivariate event characteristics.

We conclude that hydro-SMILEs are an important tool for the investigation of rare river droughts and complex event characteristics. As such events are most impactful, a hydro-SMILE, especially when based on a spatially explicit, process-based impact model, can also be used to design physically-plausible storylines of drought events (Shepherd, 2019; Sillmann et al., 2021). Many studies only analyse moderate extreme levels of peak low flow, which might hide climate change induced alterations of the low flow. Hence, we recommend to distinguish between summer low flow and winter low flow, as rainfall-driven and melt-driven catchments show distinct reactions to climate change (Brunner et al., 2019b). A low flow typology could be applied to even better disentangle the effects of the climatic and hydrological drivers on different types of low flow genesis (Brunner et al., 2022). Further we suggest to analyse the seasonality, low flow duration, or bivariate return levels of peak flow and duration. They show a highly varying degree of sensitivity to climate change alterations, but are relevant for the water management. Climate change projections of these characteristics can inform decision-makers and adaptation planners for the dimensioning of reservoirs, renaturation of rivers, and regulations of groundwater and surface water withdrawal.

Given the high-emission scenario and strong warming in the CRCM5, the projected drastic changes in river droughts may be interpreted as a worst-case scenario. However, using global warming levels as reference, the projected local changes can be linked to global emissions. The projected climate of 2040–2069 corresponds to a global warming of slightly more than 3 °C above pre-industrial levels. Hence, the projected river drought events for the period 2040–2069 are within the range of possibilities that water managers must consider for adaptation planning until the end of the century.