Large increases of multi-year droughts in north-western Europe in a warmer climate

Three consecutive dry summers in western Europe (2018–2019–2020) had widespread negative impacts on society and ecosystems, and started societal debate on (changing) drought vulnerability and adaptation measures. We investigate the occurrence of multi-year droughts in the Rhine basin, with a focus on event probability in the present and in future warmer climates. Additionally, we investigate the temporally compounding physical drivers of multi-year drought events. A combination of multiple reanalysis datasets and multi-model large ensemble climate model simulations was used to provide a robust analysis of the statistics and physical processes of these rare events. We identify two types of multi-year drought events (consecutive meteorological summer droughts and long-duration hydrological droughts), and show that these occur on average about twice in a 30 year period in the present climate, though natural variability is large (zero to five events can occur in a single 30 year period). Projected decreases in summer precipitation and increases in atmospheric evaporative demand, lead to a doubling of event probability at 1 ∘\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$^\circ$$\end{document}C additional global warming relative to present-day and an increase in the average length of events. Consecutive meteorological summer droughts are forced by two, seemingly independent, summers of lower than normal precipitation and higher than normal evaporative demand. The soil moisture response to this temporally compound meteorological forcing has a clear multi-year imprint, resulting in a relatively larger reduction of soil moisture content in the second year of drought, and potentially more severe drought impacts. Long-duration hydrological droughts start with a severe summer drought followed by lingering meteorologically dry conditions. This limits and slows down the hydrological recovery of soil moisture content, leading to long-lasting drought conditions. This initial exploration provides avenues for further investigation of multi-year drought hazard and vulnerability in the region, which is advised given the projected trends and vulnerability of society and ecosystems.


Introduction
Three consecutive dry springs and/or summers in the Netherlands and large parts of north-western Europe (2018,2019,2020), without full drought recovery in the winters in between, led to widespread problems regarding water availability and water quality. In the aftermath, water managers and policy makers raised questions on the role of climate change in creating such events, especially regarding the likelihood of multi-year drought (MYD) events versus single year or 'normal' drought events. Here we present an exploratory study into MYD occurrence in the Rhine basin, based on state-of-the-art large ensemble climate model simulations. The focus of our analysis is on the likelihood of MYD events in the present climate and in a warmer future, relating projected changes in probability to anthropogenically forced changes in mean climate and in climate variability. Additionally we describe the temporally compound physical processes in a typical MYD event.
The summer of 2018 was one of exceptional drought in western Europe ( Fig. 1a; e.g. Philip et al. 2020;Zscheischler and Fischer 2020;. Low precipitation in the summer half year combined with high 1 3 atmospheric evaporative demand due to high temperatures and high incoming solar radiation (Fig. 2a,d,g), led to large scale drought conditions decreasing soil moisture availability, lowering ground water tables and reducing river discharge ( Fig. 1c-f). In many regions the recovery of ground water tables was not completed in the subsequent winter (Brakkee et al. 2021), making ecosystems and different societal sectors more vulnerable to dry conditions in the following summer (pre-conditioning, ). The summer of 2019 was again drier than normal, which magnified existing drought impacts in vulnerable ecosystems and societal sectors. Dry conditions in the spring of 2020 completed this sequence of consecutive meteorological summer droughts (Figs. 1, 2).
The severity of drought conditions and impacts varied widely between regions, for example due to differences in soil characteristics and topography, which lead to differences in the temporal disconnect between precipitation and groundwater levels (Van Loon et al. 2014;Brakkee et al. 2021). Therefore, in some locations ground water tables and ecosystems recovered fully every winter, meaning that those droughts could maybe be analysed as individual droughts. However, in locations without full recovery in the winter season, failing to recognise the multi-year nature of the events can lead to an underestimation of impacts and  , c preciptation [mm/month], d Makkink reference evapotranspiration [mm/month], e volumetric soil moisture 10 3 m 3 ∕m 3 , f runoff [mm/ month]. All data based on based on ERA5 reanalysis; a, b SPEI indi-ces computed using calibration period 1961-2020, c-f anomalies relative to a 1991-2020 climatology, grey shading shows the 10-25-75-90 percentiles of this climatology. Orange shading and bars in a, b show identified MYD events risks, as for example ecosystems are more vulnerable to repeated droughts (e.g. Anderegg et al. 2020;Peters et al. 2020). In the Netherlands, the full hydrological recovery of the 2018-2019-2020 drought event was announced in the spring of 2021 (LCW 2020; NOS 2021). Note that system stress, either in natural or socio-economic systems, may continue to influence the system also after 'full' meteorological/ hydrological drought recovery.
From a climatological perspective, summers in western Europe are characterised by a precipitation deficit. End of summer soil moisture and ground water deficits are recharged during winter seasons which have a (larger) precipitation surplus. Natural variability in the climate system causes variability in the timing and severity of the summer precipitation deficit, when it is larger than normal we speak of 'meteorological drought' conditions (Van Loon 2015). i 2020. All data based on ERA5 reanalysis regridded to a 2 • × 2 • grid; anomalies relative to a 1991-2020 climatology. Black rectangle denotes the Rhine basin box used for time series analysis in the study Meteorological drought can lead to lower than normal soil moisture initiating a 'soil moisture drought', which in turn may develop in a 'hydrological drought' characterised by lower than normal ground water levels or river discharge. This propagation of drought through the hydrological cycle is strongly affected by the local soil characteristics (Van Loon et al. 2014). The duration and severity of a drought are strongly impacted by both the meteorological forcing as well as the soil and groundwater characteristics (Van Lanen et al. 2013).
The 2018-2020 drought had large societal and ecological impacts, putting drought management prominent on the agenda of policy makers. Consequently, questions on the effects of climate change on drought, the likelihood of a sequence of droughts, and how to manage (future) droughts were asked. In the Netherlands, where many of the current water management systems are designed to prevent floods and limit the impact of wet periods, these drought events started high-level discussions on how to prepare for, manage and adapt to future droughts (Witte et al. 2020).
The aim of this paper is twofold: we provide a historical context for the 2018-2020 MYD event, and investigate changes in the likelihood of similar events due to anthropogenic climate change. We use multiple reanalysis products and large ensemble climate model simulations to robustly assess (changes in) MYD frequency and duration. Large ensemble data are uniquely suited for the analysis of rare extreme events like MYDs (Deser et al. 2020). Their main advantage lies in the simulation of the full distribution of climate variability and extreme events, which means we can analyse extremes without relying on fitted (sometimes extrapolated) statistical models of the tail of the distribution of variability (Van der Wiel et al. 2019). We limit our analysis to the large-scale meteorological processes leading to MYD, and because of that refer to this study as 'exploratory'. It is however important to note that local meteorological and hydrological processes influence drought development (Van Loon et al. 2014;Brakkee et al. 2021), further study of local hydrological MYD processes is therefore warranted.

Data
We used the ERA5 reanalysis dataset to investigate MYD events in the recent history. ERA5, including its backward extension, covers the time period 1950-2020 (Hersbach et al. 2020;Bell et al. 2021). From this archive monthly precipitation, near-surface air temperatures, incoming solar radiation, soil moisture content and runoff were used to characterise MYD events and quantify MYD hazard risk. We limit the analysis to landpoints in a box containing the Rhine river basin (46-54 • N, 2-12 • E, Fig. 2). Other reanalysis datasets (JRA-55, ERA-20C, 20CRv3, detailed in Table 1) were used to investigate uncertainty in reanalysis data. ERA5 and JRA-55 are both modern reanalyses using up-to-date meteorological models and 3-D observations for data-assimilation. In contrast, ERA-20C and the 20CRv3 only assimilate surface fields, and are therefore expected to have inferior performance. However, these century-long reanalysis products provide more reliable statistics with regards to rare events, highly relevant for analysing rare MYD events. Finally, we use the gridded observations dataset E-OBS v23.1e (Cornes et al. 2018) to quantify biases in reanalysis data.
We used the climate model experiments available in the Multi-Model Large Ensemble Archive (MMLEA, Deser et al. 2020). The available models and the number of ensemble members per model experiment are shown in Table 1. All models use historical forcing for the period up to 2005, and the strong forcing scenario RCP8.5 from 2006 onwards (Van Vuuren et al. 2011;Schwalm et al. 2020a;Hausfather and Peters 2020;Schwalm et al. 2020b). These are quite coarse resolution models, with horizontal resolution in the atmospheric model component ranging from about 1 • to 2.8 • . To enable fair comparison between reanalyses and models, all data were regridded to a common 2 • × 2 • grid. Finally, we limit most of the analysis to the years 1950-2100, this period is covered by the climate models and the ERA5 reanalysis (up to 2020).

Model evaluation
We first investigate the ability of the different climate models employed in this study to simulate the observed regional hydroclimate. Since we study MYD events in the Rhine basin, we quantify potential model biases in the simulated time series over this region. We assess the annual mean bias, and the root mean square error (RMSE) and temporal correlation of the simulated mean annual cycle, taking the ERA5 reanalysis as observational truth. To assess observational uncertainty we include the same quantities for the other reanalysis datasets. Then we assess if the models reproduce the linear trends in the recent historic period, and finally we compare persistence in simulated time series. Additional figures and tables for model evaluation are provided in Section S2 of the Supplementary Information (SI).
All reanalysis products and all climate models simulate too high annual mean precipitation values. The modern reanalysis products (ERA5, JRA-55) reproduce the mean annual cycle quite well, the models have larger biases. In terms of the simulated mean annual cycle, CESM1-CAM5 most closely follows the cycle of E-OBS (Table 2). Individual model biases are between one and two times larger than the biases in the modern reanalyses. For reference evapotranspiration (ET 0 ) differences between ERA5 and JRA-55 are substantial, caused by differences in the simulation of clouds and aerosols and consequently differences in incoming solar radiation. ERA5 reproduces the observed values from E-OBS well, JRA-55 is positively biased (too high First three columns show annual mean bias, root mean square error (RMSE) of the annual cycle and the temporal correlation coefficient (Corr.) of the annual cycle, all relative to the E-OBS data set. Evaluation period is 1991-2020, variables shown precipitation (PR), reference evapotranspiration (ET 0 ) and their difference, the precipitation deficit (ET 0 −PR). Next two columns show the linear trend over the period 1958-2010, where precipitation and reference evapotranspiration are regressed onto a smoothed global mean surface temperature. Final column shows the 1-month autocorrelation for precipitation and reference evapotranspiration anomalies over the period 1960-2020. Similar statistics for additional variables are included in the SI For monthly precipitation deficit, defined to be the difference between reference evapotranspiration and precipitation, CESM1-CAM5 stands out as the model with smallest annual bias and RMSE, CANESM2 and GFDL-CM3 additionally stand out as better performing models, though have somewhat compensating errors. CSIRO-Mk3-6-0 and GFDL-ESM2M (biased precipitation), and EC-Earth 2.3 (biased reference evapotranspiration) have lowest realism in terms of the mean annual cycle of precipitation deficit. All climate models reproduce the observed historic trend of precipitation (Table 2), the JRA-55 reanalyses is an outlier as the only dataset showing a negative trend of precipitation. The range of trends within a single model ensemble is due to internal variability of the climate system (each member is in a different phase of e.g. decadal variability), but we do not expect a range of trends between the observations and reanalysis products, given that these all represent the historic period which is a single, but importantly the same, realisation of internal variability. The trend of reference evaporation is overestimated by ERA5 and JRA-55 compared to E-OBS, caused by too high trend in incoming solar radiation. Six of the models are in agreement with this observed trend, CSIRO-Mk3-6-0 overestimates the observed trend also due to too high trend in radiation. There are quite large model differences, caused by differences in local temperature trends (too low for GFDL-ESM2M) and incoming radiation (too high for CSIRO-Mk3-6-0 and GFDL-CM3).
Finally we compare 1-month auto-correlation values of simulated time series. There is no auto-correlation in monthly precipitation, the climate models have slightly higher values than E-OBS and the reanalyses. Auto-correlation is larger for temperature, incoming radiation and reference evapotranspiration. Also here the climate models show slightly larger values, with the exception of GFDL-ESM2M. Because SPEI-6 and SPEI-12 time series are based on running mean type calculations, auto-correlation in those time series are very high (0.85 and 0.95 respectively). The climate models reproduce these values (SI Section S2).

Multi-year drought definitions
We used the Standardized Precipitation Evapotranspiration Index (SPEI, Vicente-Serrano et al. 2010;Beguería et al. 2014) to characterize drought. The SPEI computes a water balance using monthly values of precipitation (PR) and reference evapotranspiration (ET 0 ) similar to the precipitation deficit discussed in Sect. 2.1. By including reference evapotranspiration, SPEI takes into account the influence of atmospheric evaporative demand on drought development. We used the Makkink formula to compute reference evapotranspiration (SI Section S1), this choice was made with consideration of the availability of different variables in the MMLEA, local standards and sensitivity of calculations of evapotranspiration for trends in temperature. The SPEI was calculated using the R package SPEI (Beguería et al. 2014), assuming a loglogistic probability distribution, and calibrated based on the period 1961-2020 (for ERA-20C and 20CRv3 using all available years in this period).
The SPEI can be computed on different timescales, enabling the identification of different types of drought. We used the 6-month SPEI (SPEI-6) to quantify meteorological drought (Fig. 1a), and the 12-month SPEI (SPEI-12) as a measure of hydrological drought (Fig. 1b). The propagation from meteorological drought to hydrological drought takes time, hence the longer SPEI time scale, in ERA5 SPEI-12 has a high Pearson correlation with runoff anomalies of 0.66, indicating a that SPEI-12 time series are good proxy for hydrological droughts. Due to the variety of societal and natural impacts, caused by different vulnerability to drought conditions, it is impossible to categorise the 2018-2020 drought event in western Europe as a certain type of drought. We therefore quantify the changing hazard risk of MYD using two metrics: the first based on sequential or consecutive meteorological summer droughts, the second based on long-duration hydrological droughts.

Consecutive meteorological summer droughts
We focused our investigation of meteorological drought to the summer months, because the vulnerability to lower than normal water availability is substantially larger in these months than in the winter season due to on average higher water demands from nature and society and lower water availability. We selected years in which the SPEI-6 time series drops below minus one, often classified as 'moderate drought' (lowest 15% of data, dashed blue line in Fig. 1a), for at least two months in the period June to November. Note that this period is shifted from the traditional summer drought season (April-September) to take into account the lagged effect of SPEI-6. If a minimum of two such years with meteorological droughts occur in a sequence, we refer to this event as a Consecutive Meteorological Summer Drought (CMSD, orange shading in Fig. 1a). We studied both the frequency and the duration of such events, with a sequence of meteorological droughts longer than two years counted as a single event, and the length of such a sequence (in years) as the measure of duration (i.e. 2018-2020 is one event of 3 year duration, Fig. 1a).

Long-duration hydrological droughts
We based the identification of Long-Duration Hydrological Droughts (LDHD) on run theory (Yevjevich 1967). We used a SPEI-12 threshold of minus one, 'moderate drought', to select drought events (dashed blue line in Fig. 1b), if these events last a minimum of 12 months we classify the event as a MYD event (orange shading in Fig. 1b). Comparable methods of determining drought length, i.e. periods under a (SPI/ SPEI) threshold, are frequently used in drought research (e.g. Serinaldi et al. 2009;Barker et al. 2019). Here we took the number of selected events for the frequency count, and the total length of the period (in months) in which SPEI-12 stays below the threshold as a measure of duration of a single event.

Drought 2018-2020 in historical perspective
The dry and warm northwestern European spring and summer of 2018 was caused by a persistent positive phase of the North Atlantic Oscillation, regional climate change, and, for part of the summer, a stationary Rossby wave (Drouard et al. 2019;Kornhuber et al. 2019;Philip et al. 2019). From April to November precipitation was much lower than normal precipitation (Figs. 1c, 2a) and higher than normal reference evapotranspiration (Figs. 1d, 2d, g) resulted in a SPEI-6 value lower than −2 (exceedance probability < 2.5 %, Fig. 1a). Soil moisture and runoff responded to the meteorological drought with a slight delay (Fig. 1e, f), it is this hydrological response that caused the widespread negative impacts of the drought on ecosystems and society (Ecorys et al. 2018;Peters et al. 2020;Brakkee et al. 2021). The summer of 2019 and spring of 2020 were also drier and warmer than normal (SPEI-6 < −1 for both years in ERA5 and JRA55, Fig. 1a and SI Fig. S3; Fig. 2). Based on this meteorological perspective and looking at the events in isolation, these droughts were thus not as severe as the 2018 drought. However, when looked at from a compound event perspective (Zscheischler et al. 2018), and if we consider 2018 as a pre-conditioning event to the droughts of 2019 and 2020 or alternatively if we consider the three droughts together as a temporally compounding event , indeed a MYD event becomes apparent as was also found by Brakkee et al. (2021). The effect of this temporal compounding can, for example, be seen in the time series of runoff, which remains below normal also in the period between the 2018 and 2019 meteorological droughts (Fig. 1f). The SPEI-12 time series captures this compounding effect from 2018 to 2019 (Fig. 1b), though, on the regional scale investigated here, the hydrological event as measured by SPEI-12 or runoff, does not continue into 2020.
The 2018-2020 event is not the first MYD event in the region. In the period covered by ERA5 we find two periods of consecutive meteorological summer droughts and five periods of long-duration hydrological drought (Table 3, Fig. 3). The identified droughts are notorious dry periods which have each been topic of many studies (e.g. Namias 1964; Zaidman et al. 2002), and have been identified in many other objective classifications of large droughts in general (e.g. Sheffield et al. 2009;van der Schrier et al. 2013;Spinoni et al. 2015;Zink et al. 2016;Barker et al. 2019). Note that for the start dates of long-duration hydrological droughts in Table 3 show the first month in which SPEI-12 is smaller than −1 , but that this takes into account accumulated precipitation and reference evapotranspiration values over the 12 months leading up to that month.
Any estimation of return periods for MYD events in the Rhine basin based on the historical record is very uncertain. This is due to the small number of events in a 30 year period, but also due to sensitivity to small differences as for example shown by differences between reanalyses products. A similar analysis using JRA-55 data for example additionally identifies 1995-1996-1997 as a MYD event, both from meteorological and hydrological point of view, but does not identify 2018-2020 because it is slightly wetter (e.g. only one month under SPEI-6 threshold in 2019). A full overview of historic MYD events based on the four reanalysis products is provided in SI Section S3.

Multi-year drought hazard risk in present-day climate
To overcome the limitation of the relatively short observational record (for this purpose), we use the large ensemble climate model simulations (Table 1) to quantify the hazard risk of MYD events in the Rhine basin in the current climate. Robust risk quantification is possible using large ensemble experiments, because the many realisations of an experiment allow us to compute a detailed empirical probability distribution from which the probability of extreme  Pascale et al. (2020Pascale et al. ( , 2021 and Van der Wiel et al. (2019. For each climate model and each ensemble member, we compute the SPEI-6 and SPEI-12 time series (calibrated for the period 1961-2020) and identify all consecutive meteorological summer droughts and long-duration hydrological droughts that start in the period 1991-2020 (but do not necessarily end in this period). The models indicate that the climatological (statistical mean) occurrence of consecutive meteorological summer droughts is two events in that 30 year period (Table 4), the estimated probabilities of individual models vary between 1.06 and 2.30 (EC-Earth 2.3 and CSIRO-Mk3-6-0, respectively). The actual occurrence of consecutive meteorological summer droughts in a single 30-year record (such as the observational record) can vary quite substantially from this mean climatological value because of internal variability in the climate system (i.e. random processes, chance). Individual model members, representing different but equally likely alternatives to the observed historical record, have between zero and five events in the 30-year period. In plain language, based on these simulations, one would expect about two events in the period 1991-2020, but wouldn't be surprised to have zero or five events. The event count in the historical record, here most reliably analysed using the two modern reanalyses (ERA5 and JRA-55, top rows in Table 4), is in agreement with these model results.
Similarly, for the long-duration hydrological droughts we find an average climatological mean occurrence of two events in 30 years in the present day climate (Table 4). The intermodel variability is comparable, between 1.12 and 2.50 events (EC-Earth 2.3 and CSIRO-Mk3-6-0). Models with fewer consecutive meteorological summer droughts also have fewer long-duration hydrological droughts (Pearson correlation r = 0.90 ), suggesting that there is a link between our two definitions of MYDs and that similar events are identified in both metrics (indeed the 2018-2020 event is identified as both a consecutive meteorological summer drought and a long-duration hydrological drought in ERA5, Table 3).
The average duration of MYD events depends strongly on our selection method and event definition. Duration statistics can therefore not be compared between consecutive meteorological summer droughts and long-duration hydrological  Fig. 1a, b, but here for the full ERA5 reanalysis record. The SPEIintensity threshold that was used for selection is shown as a dashed blue line, selected a consecutive meteorological summer droughts and b long-duration hydrological drought events are highlighted with orange shading and an orange bar, background shading shows different classifications of dry and wet periods based on SPEI levels droughts. Consecutive meteorological summer droughts last on average 2.5 years (by definition a minimum of 2 years), i.e. two or three dry summers in a row. Long-duration hydrological droughts last on average 16 months (by definition minimum 12 months), though this takes only the slowresponding hydrological system into account (as measured by SPEI-12). Combining the climatological values of event frequency and event duration, we find that the Rhine basin is in a state of MYD for 17 or 9% of the time respectively.

Consecutive meteorological summer droughts: chance or multi-year physical processes?
When investigating MYD events, it is natural to ask whether or not there are multi-year physical processes within the climate system that have a driving role in such long-lasting events. To this end we have developed an experiment to test whether or not the occurrence of multi-year events is the result of random processes causing a summer drought in two consecutive years or whether there are causal links in the climate system linking two consecutive summer seasons with drought. Such causal links might for example be found in long-lasting soil moisture anomalies or slow oceanic processes that influence the likelihood of atmospheric circulations conducive to drought. In our method of identifying consecutive meteorological summer droughts (Sect. 3) we do not take into account any relation between the first year and any subsequent year of a MYD event, dry summers are found in isolation, after which we look for consecutive events. To investigate the existence of multi-year physical processes driving consecutive meteorological summer droughts we artificially remove any multi-year process. In each model ensemble, we shuffle individual years between ensemble members. This removes any pre-conditioning effects from December into January of the following year (e.g. soil moisture anomalies) or any drivers that influence two consecutive years (e.g. slow oceanic processes), allowing us to estimate the mean frequency of consecutive meteorological summer droughts without multi-year physical processes. Similar shuffling methods were applied in e.g. Van den Hurk et al. (2015) and Goulart et al. (2021).
The shuffling experiment provides no evidence of multiyear physical processes that drive consecutive meteorological summer droughts in the Rhine basin. We performed 100 repetitions of the described shuffling experiment for each model ensemble, generating a distribution of mean frequencies in a world without multi-year physical processes (Fig. 4). The mean frequency of consecutive meteorological summer droughts in the original climate model time series is not significantly different from the range of mean frequencies derived from the shuffling experiment. We can therefore not disregard the hypothesis that two consecutive dry summers are the result of random chance. Comparatively, the EC-Earth model lies furthest in the tail of the shuffled distribution, however this is in the opposite tail as expected and might imply the model includes multi-year processes opposing MYD events. Based on this experiment, three summers of meteorological drought in a row, like 2018-2020, seems to be mostly an unfortunate coincidence, though it could of course be that the climate models systematically underestimate multi-year physical processes. Note that the impact of independent consecutive meteorological summer droughts on the hydrological system can be of multi-year nature, e.g. through incomplete recovery of wintertime ground water tables, this effect is investigated in this work through the long-duration hydrological drought events.

Increasing probability of multi-year droughts in a warmer climate
Anthropogenic climate change has already influenced temperatures, precipitation and incoming solar radiation in north-western Europe (e.g. KNMI 2021; Philip et al. 2020;Krauskopf and Huth 2020;Twardosz et al. 2021;Pfeifroth et al. 2018), which has led to changes in the probability of drought (Wanders et al. 2015;Samaniego et al. 2018;Marx et al. 2018;Hari et al. 2020). In this section we investigate what the effect of projected future trends in these quantities has been on longer duration MYD events.
To determine the statistics of MYD under climate change, we use the same procedure as was used for the present day climate (Sect. 5, Table 4), but now evaluating moving windows of each 30 years wide. Under the high-emission scenario investigated here, all models project a substantial increase of the frequency of both consecutive meteorological summer droughts and long-duration hydrological droughts. On average this increase is from two events in the period 1991-2020, to four events (consecutive meteorological summer droughts or long-duration hydrological droughts) per 30 year period in 2021-2050. There are however substantial model differences in these values, which can be explained, at least partly, by differences in climate sensitivity (i.e. the amount of global temperature rise at a given external forcing, SI Fig. S7) and differences in event frequency in the present-day climate ( Table 4). The relative increase in frequency, captured by the Probability Ratio, at different warming levels removes these differences. From that perspective, the models project an average doubling of both types of MYD events (Probability Ratio = 2 ) at 1 • C additional global warming (Fig. 5a, b, Table 5).
For long-duration hydrological droughts, model differences can be explained by annual mean trends in precipitation and reference evapotranspiration. Models with larger (non-compensating) trends in these variables, i.e. with larger trends in the precipitation deficit (ET 0 −PR), have a relatively larger trend in SPEI-12 (Fig. 6, SI Section S3). Changes in the characteristics of the SPEI-12 time series directly impact the probability and length of long-duration hydrological droughts by design of the method (Sect. 3). The seven models investigated show a continuum of changes, with MPI-ESM-LR projecting largest changes (due to a relatively strong decrease of summer precipitation and a relatively weak increase in winter precipitation) and CanESM2 the smallest changes (due to a relatively weak trend in summer precipitation and a strong trend in winter precipitation).
Consecutive meteorological summer droughts show very little model dependence in the change in MYD event frequency, model differences show in differences in the change of event duration (Fig. 5c). Models with larger trends in summer precipitation, reference evapotranspiration, and SPEI-6 show larger increases in event duration. The EC-Earth 2.3 model is an outlier in this group of climate models, showing much larger trend in event frequency but a lower trend in event duration.
For large values of global warming (e.g. exceeding + 1 • C global temperature rise in many mdels), or in the second half of the 21 st century (SI Fig. S8), MYD frequency continues to rise but our analysis breaks down. The computed trends from that point onwards are incorrect, noted by the use of dotted lines in Fig. 5. The breakdown of the analysis is related to the limit of the total number of droughts of a certain length that can occur in a 30 year period. For example, one cannot fit four droughts of ten year duration in a 30 year period. In the case of consecutive meteorological summer droughts this limit is very clearly approached in five models (CanESM2, CESM1-CAM5, CSIRO-MK3.6, GFDL-CM3 and MPI-ESL-LR, Fig. 5e). For long-duration hydrological droughts it is less obvious, but also here this limit corrupts the results slightly. We have masked all data in the figure from the point in time that we find the system is in a state of MYD for more than 75% of the time, note that this threshold is a somewhat arbitrary choice.
A physical reason for the break down of the analysis is the strong trend in SPEI (Fig. 6d, SI Fig. S10). There is such a strong general drying trend, that identifying events below the SPEI ≤ −1 threshold does not make sense anymore, invalidating the MYD selection approaches. A new hydrological regime has formed as a result of climate change (Wanders et al. 2015). The model projections used here are based on the strong forcing scenario (RCP8.5), more moderate scenarios of future greenhouse gas emission would likely give weaker trends and consequently smaller increases of MYD occurrence and duration throughout the 21st century. Assuming there is no path dependence, i.e. the relation between the speed of climate change and projected changes in MYD at warming levels, the described changes at warming levels discussed here are not impacted by the use of the strong forcing scenario.
For this analysis of the influence of projected climate change on MYD events, we have compared future event frequencies to  Fig. S7). Probability Ratio is the relative increase of frequency, calculated as the ratio of frequencies at two warming levels, e.g. PR +1C = F +1C ∕F PD . Duration change shown as percentage of event duration in the present-day period, with 100% indicating no change. a, c, e For consecutive meteorological summer droughts (CMSD), b, d, f for long-duration hydrological droughts (LDHD). Coloured lines show the individual models, black dashed line shows the multi-model mean. Lines are dotted when total time in drought exceeds 75% (see main text for details) present-day frequencies . Some of the large ensemble simulations used in this study also include data for the period longbefore the present-day period (Table 1). These models show that in the period 1850-1950 the number of MYD events was approximately stable, at about half the frequency as compared to the present-day (probability ratio of approximately 0.5) (SI Fig. S8). From approximately the 1950s the number of events starts increasing in response to the radiative forcing (historical + RCP8.5).

Contribution of forced changes in climate variability to the total increase of multi-year droughts
In the previous section we have discussed the total effect of climate change on MYD occurrence, this includes both changes in mean climate and changes in climate variability. Trends in climate variability can be different to the trends in mean climate, and the effect of those trends on extreme event occurrence isn't necessarily similar (Van der Wiel and Bintanja 2021). The models in this study project a clear shift in mean climate towards wetter winters and drier summers (Fig. 6). However, any additional effects from forced changes in climate variability on MYD occurrence remain unclear. We therefore repeat the analysis using detrended time series (i.e. the mean climate change signal removed), following the framework of Van der Wiel and Bintanja (2021), to separate the effects attributable to a changing mean climate from the effects attributable to changes in variability. We detrend the Rhine basin mean time series of simulated precipitation, temperature and incoming solar radiation using a 4 th order polynomial, recalculated the time series of reference evapotranspiration and SPEI, and identified MYD events in these time series. Any trend in MYD frequency or duration in these time series can then be attributed to forced changes in climate variability. The computation of the 4 th order polynomial was done over all ensemble members taken together but for each month individually to get the best estimate of the forced change; the results are similar if a 3 rd order polynomial is taken.
There is a small contribution of changes in climate variability to the increase of MYD events in the Rhine basin (Fig. 7), though it is much smaller than the total trend. The multi-model mean gives an increase of 0.6 consecutive meteorological summer droughts per 30 year period over a century ( +0.6 #/30 years /100 years) and an increase of 0.7 long-duration meteorological droughts per 30 year period per century, both relative to the present-day climate (Table 4). All individual models have a positive trend, which are statistically significant ( p = 0.05 ) with the exception of CanESM2 for long-duration hydrological droughts. CanESM2 has the smallest trend due to changing climate variability, CESM1-CAM5 and MPI-ESM-LR have the largest trends due to variability changes. Quantitative estimates of the trend in individual models are comparibly more uncertain than the multimodel mean due to large fluctuations, which are related partly to the number of members in the model's ensemble (Table 1). The multi-model mean estimate is robustly positive, indicating trends in climate variability add to the changes in MYD hazard risk due to changes in mean climate. Investigating the dynamical processes behind this trends is beyond the scope of this study, possibly changes in land-atmosphere feedbacks at lower mean soil moisture levels or changes in atmospheric persistence play a role (Berg and Sheffield 2018;Pfleiderer et al. 2019).

Physical processes during a typical multi-year drought
In the previous sections we have discussed the statistics of MYD events in the Rhine basin, and shown that anthropogenic climate change strongly impacts event probability and event duration. In this section we move away from the  discussion of statistics, and switch focus to local physical processes. Through composite analysis on identified simulated events in the period 1991-2020 we describe a 'typical' consecutive meteorological summer drought and a typical long-duration hydrological drought.
To investigate the local physical processes contributing to or being impacted by MYD events, we have aligned time series of all identified MYD events in the model ensembles that start in the period 1991-2020 (as in Sect. 5). We then took the average over these time series, to calculate the composite mean MYD time series, describing the behaviour of a typical MYD (Fig. 8). Consecutive meteorological summer droughts are somewhat aligned in time and have the largest SPEI-6 deviation in the summer months (by experimental design, Sect. 3), we therefore take the average over all identified events. The number of events in a composite mean varies from 17 to 227, depending on a model's ensemble size (Table 1) and mean event frequency (Table 4). The  Fig. 5a, b but here for trends in MYD attributable to changes in climate variability (see main text for details). X-axis shows time [years] rather than smoothed GMST change long-duration hydrological droughts do not have this seasonal-lock, we therefore limited the composited events to those events that start in the period June-November (number of events in a composite varies from 11 to 95). In both cases, calendar months are aligned in the resulting composite calculation; the start of the long-duration hydrological drought according to our selection criteria (Sect. 3) can vary from June to November.
The composite mean time series of SPEI-6 and SPEI-12 ( Fig. 8a, b) show the direct result of the selection process as defined in Sect. 3. For consecutive meteorological summer droughts two summers with SPEI-6 values below the threshold for moderate drought have been identified. There are no model differences here, and no obvious differences between the first and second year of drought. In the year before the consecutive droughts, SPEI-6 is above zero which can be explained by the fact that these summers (by design of the method) cannot have a meteorological summer drought (SPEI-6 < −1 for two months), leading to a small positive shift of the mean. Interestingly, in the year after the consecutive droughts (i.e. year 3) the mean SPEI-6 is still slightly negative, which may indicate lingering drought effects due to incomplete recovery. For long-duration hydrological droughts SPEI-12 lies below the − 1 threshold for twelve months as per the event definition, recovery to near-zero SPEI-12 values takes almost another twelve months.
Consecutive meteorological summer droughts are driven by two summers of lower than normal precipitation, and higher than normal temperature and incoming solar radiation (Fig. 8c, g, i). This leads to higher than normal reference evapotranspiration and a larger than normal precipitation deficit. From a meteorological perspective the two summers are very similar, though there is a hint of a slightly warmer and sunnier conditions in the second summer. The soil moisture content decreases in the first summer of drought, and doesn't fully recover to normal climatological levels in the winter following this summer (Fig. 8k). Consequently, the soil moisture response in the second summer of drought is larger (pre-conditioned compound event, . The larger soil moisture anomaly in the second summer may have impacted the temperature and incoming radiation through land-atmosphere feedbacks, but such an effect cannot be attributed in the present experimental setup. The drivers of long-duration hydrological droughts are similar, though the temporal sequence is different. In the first summer of drought, the departures from normal are larger in terms of precipitation, temperature and incoming radiation compared to the consecutive meteorological summer droughts (Fig. 8d, h, j), and precipitation remains low throughout the winter season. This results in a larger decrease of soil moisture content, and a much reduced recovery of soil moisture in the winter season compared to the consecutive meteorological summer droughts (Fig. 8l). The second summer also shows a minor precipitation deficit but much less pronounced than the first summer for this longduration hydrological event or compared to the summers in the consecutive meteorological summer drought case, but the higher than normal precipitation deficit does hinder drought recovery and hence increases the duration of the event.
The described typical MYD event is comparable to the observed 2018-2019-2020 event in many ways (Fig. 1), for example in the soil moisture response. Note that we expect less temporal variability in the composite mean event (Fig. 8) due to the averaging procedure which removes noise from internal variability and thus removes between-event variability. This has advantages because, for example, driving mechanisms can more easily be identified, but event specifics can be very important in determining impacts ) and should therefore also be taken into account in further studies.

Discussion
The presented results depend on the assumption that the climate models adequately capture the processes leading to MYDs in the region. In Sect. 2.1 we assessed potential biases in simulated annual mean cycles, historic trends and auto-correlation, but a thorough assessment of simulated MYD processes is not possible due to sparse observational evidence (too little observed events). Here we discuss some potential model and experimental deficiencies that may influence our results, and provide some suggestions for future work.
The seven climate models used in this study (Table 1) are members of a family of coupled atmosphere-ocean models, and contributed to IPCC's CMIP5 project. Interactions between atmosphere, ocean, land and sea ice are simulated, which is essential as drought is a coupled land-atmosphere phenomenon (Van Loon 2015). Reduced soil moisture influences air temperatures, evapotranspiration, and moisture recycling, and through this pathway a positive feedback exists towards further drought development (e.g. Zampieri   . a, c, e, g, u, k for consecutive meteorological summer droughts (CMSD), b, d, f, h, j, l for longduration hydrological droughts (LDHD). To align droughts in time we only include long-duration hydrological droughts events that start in the months June-November. Coloured lines show the individual models, bold dashed line shows the multi-model mean, not all models have data for soil moisture content available. Yellow shading shows the multi-model mean 25th-75th percentile range. Grey/white vertical shading show different years ◂ et al. 2009; Berg and Sheffield 2018;Herrera-Estrada et al. 2019;Benedict et al. 2021). Though the climate models do simulate soil moisture, the model formulations used are much simplified compared to state-of-the-art hydrological understanding. E.g. processes like ground water feedbacks, surface runoff, and vegetation feedbacks, all shown to be of influence to drought development (e.g. Dirmeyer 1994;Evans et al. 2017), are not yet considered in these models. More detailed assessment of the contribution of local hydrological processes in the context of multi-year events is therefore recommended.
We have identified MYD events using SPEI-6 and SPEI-12 time series, and though these time series correlate with soil moisture, ground water levels and river discharge, it is a top-down method of drought identification. Besides (largescale) meteorological forcing, local hydrological processes have a large influence on local drought impacts (Van Loon et al. 2014;Brakkee et al. 2021). The above recommendation for more detailed consideration of local processes can therefore also be applied to the event selection procedure, for example by choosing (hydrological) variables or indices closer to actual societal and/or natural impacts ( Van der Wiel et al. 2020). This would also circumvent the use of reference evapotranspiration, a measure of atmospheric evaporative demand, which during times of drought is not causally related to actual evapotranspiration (Berg and Sheffield 2018).

Conclusions
This study aimed to investigate the historical context of the 2018-2020 multi-year drought in western Europe, and investigate whether anthropogenic climate change is expected to change the likelihood of similar events. We combined multiple reanalyses datasets and large ensemble climate model simulations to explore present-day and future MYD frequency, duration and characteristics. This large data approach allowed us to robustly investigate (changing) extreme event statistics in the presence of natural variability Deser et al. 2020), and furthermore describe the multi-variate and temporally compounding drivers of MYD events in the region (Zscheischler et al. 2018. We investigated two types of MYD: consecutive meteorological summer droughts and long-duration hydrological droughts, identified through time series analysis of the 6-month and 12-month SPEI respectively. In the present-day climate around two of these events are expected in a 30 year period, though because of natural variability (internal variability in the climate system) the real number of events in a specific 30 year period can vary between zero and five events. Drought events like 2018-2020 are thus part of the realm of possibilities in the present-day climate (as is obvious by it occurring), and a comparable event could have been expected to occur based on the average frequency or return period.
The models project a strong increase of MYD event frequency and duration in response to further global climate change. The probability of MYD occurrence will double if the global mean surface temperature rises with 1 • C relative to the present-day value (approximately 2 • C warming relative to pre-industrial levels). Changes in MYD occurrence are caused by large changes in mean climate, for example due to reductions in mean summer precipitation and increases in summer reference evapotranspiration. On top of these changes in mean regional climate, forced changes in climate variability add slightly to the increase of MYD event occurrence (approximately 10-20% of the total change, Van der Wiel and Bintanja 2021). An increase of climate variability results in larger deviations from mean climate, either in amplitude or in duration. This can have influences on extreme events (as shown here for MYD frequency), providing additional challenges for societal and natural adaptation.
MYD events are the result of 'classic' meteorological forcing for drought: lower than normal precipitation and higher than normal atmospheric evaporative demand (reference evapotranspiration, due to higher than normal temperatures and incoming solar radiation). For consecutive meteorological summer droughts, at least two summers with such conditions occur in succession. We have not found evidence for multi-year processes in atmospheric forcing, which would indicate that a second dry meteorological summer is an unfortunate coincidence from a meteorological point of view (in terms of precipitation and atmospheric evaporative demand). In contrast, the hydrological response to such consecutive meteorological summer droughts can not be seen as separate droughts: soil moisture does not recover fully in the winter following the first drought. This makes the system more vulnerable to a second drought, due to reduced initial spring soil moisture before the second dry summer (Fig. 9a). Long-duration hydrological droughts are kick-started by a summer of drought, followed by a drier than normal autumn and winter season, and another relatively dry summer season. Soil moisture recovery after the first summer drought is limited by the lingering meteorological dry conditions, resulting in a longer than normal period of drought (Fig. 9b). Both types of MYD event investigated here can be viewed as compound events, they match some of the characteristics of both pre-conditioned events and temporally compounding events .
We have investigated MYD events in the Rhine basin using large ensemble climate model simulations, a similar approach can be applied in other regions where MYD events have occurred in the past (e.g. U.S. Dust Bowl, Australian Millennium Drought, South African day-zero drought ;Worster 2004;Schubert et al. 2004;van Dijk et al. 2013;Otto et al. 2018;Pascale et al. 2020). Large ensemble climate modelling is uniquely suitable as a tool for analysis of complex extreme weather or climate events, as it removes the need to rely on hypothetical statistical (extrapolation) models (Van der Wiel et al. 2019), as also shown by other application of such data for the study of drought events (e.g. Diffenbaugh et al. 2015;Pascale et al. 2020Pascale et al. , 2021. The next step following from this exploratory study is to better focus on the hydrological response during MYD events, as for example regional topographic-geological differences can lead to large differences in societal or natural consequences (Brakkee et al. 2021). For such an activity the use of large ensembles of regional climate models, with more regional details in meteorological forcing, can also be of use. Furthermore, human actions (ground water extraction, drought management) and vegetation feedbacks potentially impact MYD development (Cowan et al. 2020). Through dedicated modelling of these effects, a more bottom-up identification method for MYD events can be developed (Van der Wiel et al. 2020), leading to higher quality, more relevant and more useable climate change information (Doblas-Reyes et al. 2021).
Author contributions KvdW: conceptualisation, methodology, formal analysis, data curation, visualisation, writing-original draft, writing-review and editing. TJB: formal analysis, data curation, writingreview and editing. NW: methodology, writing-review and editing.  Higher risk due to climate change Lower risk due to climate change + Fig. 9 Physical elements of compound multi-year droughts in the Rhine basin, for a consecutive meteorological summer droughts, b long-duration hydrological droughts. Yellow boxes show meteorolog-ical forcing, blue boxes the hydrological response, red boxes societal and natural impacts; arrows indicate causal links Funding This work was supported by funding from H2020 project EUCP, grant agreement number 776613. NW acknowledges funding from NWO 016.Veni.181.049.

Data availability
The large ensemble climate model simulation data was obtained from the NCAR Climate Data Gateway (http:// www. cesm. ucar. edu/ proje cts/ commu nity-proje cts/ MMLEA/) and the MPI grand ensemble archive (https:// mpimet. mpg. de/ en/ grand-ensem ble). Reanalysis data sets and E-OBS data are available online at the source.

Conflict of interest
The authors declare that they have no conflict of interest.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.