Amplified potential for vegetation stress under climate-change-induced intensifying compound extreme events in the Greater Mediterranean Region

The Mediterranean Basin is one of the regions most affected by climate change, which poses significant challenges to agricultural efficiency and food security. While rising temperatures and decreasing precipitation levels already impose great risks, the effects of compound extreme events (CEEs) can be significantly more severe and amplify the risk. It is therefore of high importance to assess these risks under climate change on a regional level to implement efficient adaption strategies. This study focuses on false-spring events (FSEs), which impose a high risk of crop losses during the beginning of the vegetation growing period, as well as heat–drought compound events (HDCEs) in summer, for a high-impact future scenario (Representative Concentration Pathway (RCP) 8.5). The results for 2070–2099 are compared to 1970–1999. In addition, deviations of the near-surface atmospheric state under FSEs and HDCEs are investigated to improve the predictability of these events. We apply a multivariate, trend-conserving bias correction method (MBCn) accounting for temporal coherency between the inspected variables derived from the European branch of the Coordinated Regional Climate Downscaling Experiment (EURO-CORDEX). This method proves to be a suitable choice for the assessment of percentile-threshold-based CEEs. The results show a potential increase in frequency of FSEs for large portions of the study domain, especially impacting later stages of the warming period, caused by disproportionate changes in the behavior of warm phases and frost events. Frost events causing FSEs predominantly occur under high-pressure conditions and northerly to easterly wind flow. HDCEs are projected to significantly increase in frequency, intensity, and duration, mostly driven by dry, continental air masses. This intensification is several times higher than that of the univariate components. This study improves our understanding of the unfolding of climate change in the Mediterranean and shows the need for further, locally refined investigations and adaptation strategies.

1 Introduction

The latest assessment report by the Intergovernmental Panel on Climate Change (IPCC) on the current climate status has yet again made strikingly clear that human actions have significantly impacted weather and climate on a global scale and continue to do so (Eyring et al., 2021). Among the most significant impacts of ongoing climate change is the intensification of extreme events, especially in terms of temperature and precipitation (Seneviratne et al., 2021). These changes are projected to include not only a significant increase in the frequency and intensity of heat-related extremes, as well as intensified heavy-precipitation events, but also an amplified risk of agricultural drought. Additional consequences of human-induced global warming include negative impacts on agriculture. As stated by Bezner Kerr et al. (2022), there is high confidence that climate change is imposing stress on, among others things, agriculture and forestry. Furthermore, the authors point out that extreme events significantly impair food security and that heat and drought are among the main drivers exacerbating the risk of production losses.

While our general understanding of these changes is high on a global level, the level of confidence in these projections varies on a regional to local level. Giorgi (2006) conducted an investigation of regions with the highest exposure to the effects of climate change (climate change “Hot-Spots”) and demonstrated that the Mediterranean, which includes southern Europe as well as parts of northern Africa and western Asia, is one of the two most vulnerable regions worldwide. As this region is characterized by a high population density and high socioeconomic importance, the effects of climate change must be monitored closely to enable suitable protective and mitigative measures.

To assess the potential effects of climate change, general circulation models (GCMs) have been extensively used in the past. These models, driven on a global basis, are capable of providing large-scale information on changes in atmospheric predictors, although their spatial resolution is generally low. While GCMs are able to provide reasonable results under specific circumstances, the quality of GCM output proves to be too low to be applied on a sub-global scale (Di Virgilio et al., 2022; Yang and Villarini, 2021; Wang et al., 2021; Hardiman et al., 2008). To improve the quality of these projections, regional climate models (RCMs) driven by the GCM output can be consulted by conducting simulations on a regional to local level and in a higher temporal and spatial resolution (Giorgi, 2019). While this procedure, also known as dynamical downscaling, can improve the quality of climate projections, systematic model bias may still be inherited, originating from both the GCM and the RCM (Eden et al., 2012). To overcome this limitation, the inherited bias from climate models can be statistically corrected according to a reference data set of a higher quality. This method is referred to as bias correction and can, for example, be performed using station observations, satellite data, or reanalysis. As Cannon (2018), Vrac and Friederichs (2015), and Gudmundsson et al. (2012) have proven, performing statistical bias correction can significantly reduce deviations of one or multiple predictors from a reference data set and therefore improve the reliability of climate projections.

Compound extreme events (CEEs) have moved more and more into scientific focus, as the joint effects of multiple hazards may surpass those of univariate extremes (Zscheischler et al., 2020). In the context of extreme events potentially exacerbating risks to agriculture, water availability, and food security (Bezner Kerr et al., 2022), high demand for robust projections of CEEs in regions with a high vulnerability regarding these aspects becomes evident. Therefore, this study aims at investigating two types of compound event with detrimental effects on agriculture and vegetation in the Greater Mediterranean Region (GMR), a global climate change hot spot. The types of CEEs analyzed in this study include false-spring events (FSEs), which are defined as a freezing event occurring after the start of the crop-related growing season (SGS) (Ault et al., 2013; Gu et al., 2008). Freezing during this highly vulnerable period in the early stages of plant development may cause significant damage, resulting in yield loss or failure. Recently, FSEs have been investigated mostly in moderate mid-latitudinal climate zones. For example, disproportionate changes of the last day of frost (LDF) during spring and the SGS have been shown for central Europe (Vitasse et al., 2018; Zohner et al., 2016) and the United States (Peterson and Abatzoglou, 2014). Chen et al. (2021) demonstrated a negative correlation between the SGS date and mean air temperature over temperate China, indicating an earlier SGS under warmer conditions. Therefore, if the risk of spring freezing events does not proportionally decline, the risk of FSEs increases (Labe et al., 2017; Inouye, 2008; Gu et al., 2008). Less focus has been put on subtropical climate zones, where the risk of frost events is considerably lower than in moderate climate zones. However, as the GMR is characterized by complex topography, with high mountain ridges adjoining wide lowlands that represent a transition zone between the GMR and cool–temperate central Europe, the effects of freezing conditions may potentially reach out to the subtropical parts of the domain. A variety of methods to estimate the LDF and SGS exist. One example is to collect and assess crop-specific empirical data, where specific thresholds for freezing and growing conditions can be taken into account (Chamberlain et al., 2019). This is often applied in regionalized studies or for specific plant species. In terms of heat–drought compound events (HDCEs), many studies have investigated historical periods and demonstrated an increasing trend regarding the frequency and intensity of HDCEs (Ionita et al., 2021; Vogel et al., 2021). Fewer studies have investigated whether these changes will be persistent under future climates. For example, Ruffault et al. (2020) demonstrated that heat–drought-related weather conditions favoring the ignition of wildfires are likely increased under climate change.

In this study, we aim at presenting a large-scale overview of the potential of FSEs occurring in the GMR. Therefore, we apply a simplified method derived from Peterson and Abatzoglou (2014) and Leeper et al. (2021) that uses generalized thermal thresholds to determine the LDF and SGS. In terms of HDCEs, we adopt an approach by Ionita et al. (2021), using the 3-monthly Standardized Precipitation Index (SPI-3) for drought indexing and percentile-based thresholds for daily maximum temperature. By conducting this study, we seek to increase our knowledge of how the effects of global warming will unfold in terms of FSEs and HDCEs in the GMR, regarding both the frequency and the duration of these events.

Figure 1 Study domain covering the Greater Mediterranean Region (GMR) outlined in red and classification of the local climate according to an updated version of the Köppen–Geiger climate classification. Data and visualization obtained from Kottek et al. (2006).

Besides the projected changes in FSEs and HDCEs, we also investigated the deviations of crucial near-surface atmospheric predictors from the mean state during these events in order to improve our level of knowledge on what conditions act favorably towards the occurrence of FSEs and HDCEs. Similar works have been done by Ionita et al. (2021) and Mastrantonas et al. (2021), who investigated connections between extreme and compound events and large-scale atmospheric patterns. The inclusion of weather patterns, however, is aggravated in the context of multivariate bias correction, as the conservation of spatial and temporal coherence of multiple variables requires a significantly higher computational effort (Cannon, 2018). By considering sea level pressure as well as zonal (u) and meridional (v) wind speeds in an only temporally coherent bias correction setup, we seek to reduce computational cost while preserving the ability to assess the origin of air masses and estimate pressure deviations under FSEs and HDCEs.

To obtain these results, we apply a multivariate bias correction method based on the N-dimensional probability density function transform (MBCn) presented by Cannon (2018). A total of 13 GCM–RCM combinations were obtained from the Coordinated Regional Climate Downscaling Experiment (CORDEX), as well as ERA5 data as the reference for the bias correction procedure. We inspect a late-century future period (2070–2099) under the high-impact Representative Concentration Pathway (RCP) scenario 8.5 and compare the results to 1970–1999. The choice of scenario was made to demonstrate potential changes in FSEs and HDCEs on the higher end of the range of emission scenarios (Riahi et al., 2011). These research questions are addressed in the following sections:

• a.

  To what extent can the quality of climate model output be improved, in the context of reproducing threshold-based metrics for a multivariate and interdependent task?

• b.

  Are FSEs relevant in the GMR, and how is their occurrence projected to change towards the end of the 21st century? What deviations in the near-surface atmospheric state are connected with FSEs?

• c.

  How will the frequency, intensity, and duration of HDCEs change towards the end of the 21st century? What deviations in the near-surface atmospheric state are connected with HDCEs?

• d.

  What implications will the projected changes in FSEs and HDCEs potentially have for vegetation and crop efficiency?

2 Study domain

In this study, we investigate the potential effects of climate change on the Greater Mediterranean Region (GMR). The domain margins are given as 27.75–50° N and 10° W–37° E and are displayed in Fig. 1. It is important to note that the domain definitions were selected to resemble the Mediterranean domain provided within the Med-CORDEX framework. The label “GMR” as given in this study therefore refers to a model domain definition, rather than a climatic classification of the Mediterranean. The GMR within the context of this study includes most parts of the southern European continent, as well as portions of northern Africa and western Asia. As a result of the Alpine orogeny, this region is characterized by a variety of mountainous formations that merge directly into the Mediterranean Sea. Some of the most prominent ridges are the Alps, the Apennines, and the Pyrenees on the European side; the Atlas formation on the African side; and the Taurus Mountains on the western Asian side. However, these mountain ranges are separated by wide lowlands, where large river systems drain precipitated water from the mountains towards the sea. These lowlands, e.g., the valleys of the Ebro, Rhône, and Po, substantially differ from the mountain ranges in terms of climatic characteristics, adding even more to the diversity and complexity of the GMR.

This climatic complexity is also reflected in the Köppen–Geiger climate classification (Kottek et al., 2006). The northern parts of the study domain are under the influence of the interplay between the polar front and subtropical highs within the zone of westerlies (Cf climate; see Fig. 1). This causes a constant alteration between moderate temperatures and a high likelihood of precipitation under maritime influence and dry air masses of continental origin, causing hot temperatures in summer and the opposite in winter. With decreasing latitude, the influence of the westerlies recedes and the subtropical ridge becomes dominant. This persistent high-pressure system causes calm, warm, and dry conditions under descending air masses. Due to the southward shift of the subtropical ridge following the annual cycle of the Intertropical Convergence Zone (ITCZ), the influence of the westerlies becomes stronger in winter, causing precipitation levels to increase. The resulting warm and dry summer climate (Cs climate), which is most common on the west sides of continents, is also often referred to as Mediterranean climate. The southernmost area of the study domain, where the subtropical ridge is persistently dominant throughout the year and precipitation is rare, is classified as dry–arid, or dry–semi-arid (BW and BS climate). With increasing elevation in the northern mountainous regions, the climate is cooler and more continental, with increased temperature spans and persistently high precipitation (Df climate). The highest portions of the Alps, where the mean temperature of the warmest months is below 10 °C, fall under the polar tundra class (ET; Beck et al., 2018).

3 Data and methods

3.1 Data

To perform bias correction, a reference data set must be obtained, e.g., reanalysis or in situ observations, along with the climate model output. Initially, climate model output data for a historical period must be obtained, as well as the corresponding projection data for a future period under a specific scenario. To optimize the quality of the model output, the reference data set must inherit a higher accuracy regarding the “true” atmospheric state. Depending on the research aim, this can be, for example, observational station data, satellite data, or reanalysis. The latter is especially relevant in regions with a low density of observations, for example, regions with no population or very low population density or regions over the ocean. While a variety of observational data sets are available for the land portions of our study domain, we aimed at specifically including areas over the ocean, thus excluding all purely land-based data sets. Therefore, the latest global reanalysis data set published by the European Centre for Medium-Range Weather Forecasts (ECMWF), ERA5 (Hersbach et al., 2023, 2020), was obtained for this study. Out of the wide variety of variables included in ERA5, hourly 2 m air temperature, the precipitation sum, and sea level pressure, as well as 10 m u and v components of wind speed (eastward and northward, respectively) were obtained for this study and aggregated to daily mean values and, for precipitation, daily sums. Hourly 2 m temperature was aggregated to maximum and minimum daily values. ERA5 data were obtained for the period 1970–2020; however, the 30-year period 1970–1999 is used as historical reference data within the bias correction process and denoted as HIST in the following. There exist weaknesses of the ERA5 data set, as discussed by Velikou et al. (2022) and Xu et al. (2022), for example spatial quality disparities and misrepresentation of extreme events. Nevertheless, as also confirmed by these authors, the overall quality is high, making ERA5 a justifiable choice as the reference.

Table 1 Description of all obtained data sets for this study, including 13 GCM–RCM combinations from EURO-CORDEX for historical and future periods, as well as ERA5 reanalysis data for reference purposes.

In terms of climate model output, a total of 13 different combinations of GCMs and RCMs was obtained from CORDEX, consisting of dynamically downscaled realizations driven by GCMs obtained from the fifth phase of the Coupled Model Intercomparison Project, CMIP5, initiated by the World Climate Research Programme (WCRP; Taylor et al., 2012). For this study, EURO-CORDEX (Cinquini et al., 2014), at the highest available spatial resolution of 0.11 km × 0.11 km and at a daily temporal resolution, was obtained. While multiple RCP scenarios are provided, we focus only on the high-impact scenario RCP8.5 (Riahi et al., 2011) and the distant-future period 2070–2099. Before being subjected to the bias correction procedure, the model data were bilinearly interpolated to match the 0.25 km × 0.25 km spatial resolution of the ERA5 data. The 13 model combinations consist of six GCMs, comprising CNRM-CERFACS-CNRM-CM5 (Voldoire et al., 2013), ICHEC-EC-EARTH (Hazeleger et al., 2012), IPSL-CM5A-MR (Dufresne et al., 2013), MOHC-HadGEM2-ES (Collins et al., 2011), MPI-ESM-LR (Jungclaus et al., 2013; Stevens et al., 2013), and NCC-NorESM1-M (Bentsen et al., 2013), as well as four RCMs, comprising DMI-HIRHAM5 (Christensen et al., 2007), GERICS-REMO2015 (Jacob et al., 2012; Jacob, 2001), KNMI-RACMO22E (van Meijgaard et al., 2008), and SMHI-RCA4 (Strandberg et al., 2014; Samuelsson et al., 2011). The data sets were downloaded from the DKRZ node of the ESGF data portal (Cinquini et al., 2014). All data utilized for this study are described in detail in Table 1.

3.2 Methods

To obtain optimized climate model output, we apply a multivariate and dependency-preserving bias correction method (MBCn). This method is presented below, as well as the definitions of FSEs and HDCEs, for which we adopted threshold-based approaches that have been applied in former studies (Peterson and Abatzoglou, 2014; Leeper et al., 2021; Ionita et al., 2021). We additionally aimed at offering insights into the prevailing near-surface atmospheric conditions during these events. Therefore, in addition to minimum and maximum temperatures and precipitation, we included sea level pressure in the bias correction procedure, as well as u and v components of near-atmospheric wind, to calculate the mean wind direction. The additional information on pressure deviations can give indications of the prevailing type of action center (i.e., high-pressure or low-pressure systems), and the characteristics of the influencing air masses can be derived by means of their origin, i.e., the predominant wind direction. To inspect the statistical significance of near-atmospheric deviations under CEEs from the mean state, we applied a Mann–Whitney U test (Mann and Whitney, 1947; Student, 1908). As the 13 included models are independent, this test is applied separately to each of the models. The sum of models indicating statistical significance is obtained for both positive and negative deviations from the mean and is presented in the results. As the inspection of mean wind directions is aggravated, due to a break in the scale between 359 and 0°, we inspect the mode instead of the mean, i.e., the predominant wind direction with the highest frequency within the sample. The statistical significance of the deviations of wind direction is obtained by applying Fisher's exact test (Fisher, 1935) to the occurrence count distributions of the eight wind directions: N, NE, E, SE, S, SW, W, and NW. For each test, we consider statistical significance at a significance level of 5 %.

3.2.1 Univariate bias correction

When aiming at an optimization of climate model output, multiple approaches of differing complexity exist. For example, if a linear bias towards a reference data set is to be removed, a delta (additive) or a factor (multiplicative) can be added to the model output (Deque, 2007). In general, the bias-corrected time series of a variable x_bc can be obtained when a statistical transformation function f is applied to the raw model output x_m,p, expressed as

$$\begin{array}{}\text{(1)}& x_{\text{bc}}=f\left(x_{m,p}\right)\end{array}$$

by Piani et al. (2010). However, complex climate models often inherit a more complex bias structure, when, for example, trends and instationarities come into play. To specifically account for differing bias within the distribution of a simulated climate variable, the method of empirical quantile mapping (EQM) was introduced (Piani et al., 2010; Boé et al., 2007; Gudmundsson et al., 2012). Within EQM, the cumulative distribution function (CDF) F of the raw model output m within a historical calibration period c is matched to the corresponding CDF of the reference data (e.g., reanalysis) o within c, given as

$$\begin{array}{}\text{(2)}& F_{m,c}\left(x_{m,c}\right)=F_{o,c}\left(x_{o,c}\right)\end{array}$$

by Gudmundsson et al. (2012) and Tong et al. (2021). In the next step, the CDF of the raw model output within the projected period p is matched to the inverse CDF of the reference data in c, $F_{o,c}^{-\mathrm{1}}$, as given in

$$\begin{array}{}\text{(3)}& x_{\text{bc}}=F_{o,c}^{-\mathrm{1}}\left[F_{m,p}\right(x_{m,p}\left)\right],\end{array}$$

in order to receive the bias-corrected model data x_bc (Gudmundsson et al., 2012; Tong et al., 2021). In a more general sense, the additional knowledge of the reference data regarding the distribution structure of the true atmospheric state is applied to the raw model output for a known historical period. By doing so, information on the performance of the correction can be evaluated. If the correction performance is deemed satisfactory and under the assumption that model bias is stationary within the historical and projection periods (Maraun, 2012; Maraun et al., 2010), the same additional knowledge of the reference data is applied to the projection period.

Trends within the raw model output may have a negative effect on correction performance, and, correspondingly, trends infused by the correction method may not represent the true atmospheric changes, as pointed out by Maurer and Pierce (2014) and Maraun (2013). To account for this, the method of EQM was adjusted to initially extract trends from the raw model output, then perform the correction process, and afterward apply the trend back to the corrected model output as introduced by Cannon et al. (2015). Within this quantile delta mapping (QDM) procedure, the linear trend Δ_m regarding each time step t, defined as

$$\begin{array}{}\text{(4)}& {\mathrm{\Delta }}_{\mathrm{m}}\left(t\right)={\displaystyle \frac{x_{m,p}\left(t\right)}{F_{m,c}^{-\mathrm{1}}\left[F_{m,p}^t\left[x_{m,p}\right(t\left)\right]\right]}},\end{array}$$

is returned to the corrected model output by applying

$$\begin{array}{}\text{(5)}& {\widehat{x}}_{\text{bc},\mathrm{\Delta }}\left(t\right)={\widehat{x}}_{\text{bc}}\left(t\right){\mathrm{\Delta }}_{\mathrm{m}}\left(t\right).\end{array}$$

3.2.2 Multivariate bias correction

All previously described methods are applied independently to each variable. While these methods are not capable of specifically adjusting the day-to-day variability to match that of the reference data, the temporal physical coherence is nevertheless upheld, allowing for long-term climatological inspections (Olschewski et al., 2023). However, the independence of the correction process for multiple variables, where each variable is corrected on its own, aggravates the investigation of multivariate matters such as compound events (Zscheischler et al., 2019; Rocheta et al., 2014). Therefore, methods accounting for the dependency of multiple variables become necessary. As with univariate bias correction, a variety of approaches towards multivariate bias correction have been developed, each differing, for example, in the statistical metric that is used to adjust the dependence or the level of restriction due to specific assumptions (Vrac and Thao, 2020; Cannon, 2016; Vrac and Friederichs, 2015; Bürger et al., 2011). This study applies a multivariate bias correction procedure (MBCn) developed by Cannon (2018), based on an image processing technique using the N-dimensional probability density function transform (N-pdft) presented by Pitie et al. (2005) and Pitié et al. (2007). Within MBCn, a random orthogonal rotation is first applied to the input climate data. Subsequently, QDM is applied to the marginal distributions of the rotated input data as described in the previous section. Finally, the rotated and QDM-adjusted data are inversely rotated to receive the corrected multivariate output (Cannon, 2018). As the author describes subsequently, this process is iterated until the distributions of the model data and the reference data match. In this study, we carried out 100 iterations, which has proven sufficient in previous studies (Dieng et al., 2022; Cannon, 2018). Within the framework of this study, this procedure included the six atmospheric variables listed in Table 2 (wind direction from u and v components) and was individually carried out for every grid cell. The suitability and potential of this method in terms of climate data and their application have been proven in multiple studies, including Cannon (2018), Dieng et al. (2022), Lemus-Canovas and Lopez-Bustins (2021), Singh et al. (2021), and Meng et al. (2022). Additionally, besides MBCn, Cannon (2016) explored other ways of correcting the multivariate dependence structure, based on the Pearson and Spearman correlation dependence structure. However, these approaches are accompanied by pronounced restrictions, making MBCn the preferable choice.

3.3 Compound-event definitions

A comprehensive overview of the types and definitions of compound weather and climate events is provided by Zscheischler et al. (2020). In general, the authors discriminate four different types of compound events, based on their temporal and spatial characteristics. FSEs, as defined for this study, are representative of preconditioned compound events, in which the precondition is the warm anomaly of daily minimum temperature, which is beneficial towards the onset of the growing season, and the hazard is a subsequent frost event, potentially causing damage to crops within the early stages of plant development. Due to its subsequent nature, FSEs may also be categorized as temporally compounding events. FSEs consist of anomalies of one variable only, i.e., daily minimum temperature, and are therefore considered univariate compound events. Heat–drought compound events (HDCEs), by definition in this study, consist of multiple hazards occurring simultaneously, which are classified by Zscheischler et al. (2020) as a multivariate compound event. However, if persistent drought is considered a precondition, HDCEs may also be treated as a preconditioned compound event.

Peterson and Abatzoglou (2014)Leeper et al. (2021)McKee et al. (1993)Perkins and Alexander (2013)Ionita et al. (2021)Peterson and Abatzoglou (2014)Peterson and Abatzoglou (2014)Ionita et al. (2021)

Table 2 Variables, extreme events, and compound-event definitions utilized in this study.

Percentile-based thresholds are a crucial component in this study. These thresholds were calculated for the historical period and subsequently applied to the future period. By applying historical thresholds to future data, the results show changes in extremes that have already been experienced and can therefore be better assessed by potential users. All included variables, as well as extreme event and compound-event definitions, are summarized in Table 2 and described in detail below.

3.3.1 False-spring events

To define FSEs, we take three thresholds into account, including the start of the growing season (SGS), the number of days between the SGS and the frost event, and the thermal definition of frost events. Robeson (2002) suggested the period between spring and fall freezes as the growing season and the SGS correspondingly as the day of the year (DOY) from which daily minimum temperatures are persistently above 0 °C. Leeper et al. (2021) built upon this definition and suggested the use of multiple thermal thresholds, i.e., 0, 5, and 10 °C, and these thresholds are also applied in this study. These thresholds can be considered representative of different phases within the period of continuous warming.

Applying a time delay between the SGS and the FSE-defining frost event is an attempt to consider the most vulnerable time of leaf tissue, in the phase between budburst and full leafout (Chamberlain et al., 2019). In this context, Peterson and Abatzoglou (2014) applied various numbers of days between 0 and 15, of which the authors found no sensitivity of the results to lag times above 7 d. We adapted the time lag of 7 d, as well as the definition of daily minimum temperatures below −2.2 °C for frost events.

Based on Peterson and Abatzoglou (2014), we investigate changes in the false-spring-event index (FSEI), which indicates the portion of years with an occurrence of FSEs within a selected period. A year is counted towards FSEI when the SGS, which is defined as a period of at least 7 d in which the daily minimum temperature does not fall below 0, 5, or 10 °C, happens before the last frost event, defined as a daily minimum temperature below −2.2 °C. The considered variants of FSEI are therefore dependent on the definition of the SGS and include FSEI₀, FSEI₅, and FSEI₁₀. In general, the data from January through June were considered in the calculation of FSE. However, the occurrence of FSEs is dependent on the DOY of the SGS, which lies within this range of months for the predominant part of the region (see Sect. 4.2.1).

3.3.2 Heat–drought compound events

Ionita et al. (2021) investigated historical compound hot and dry events over Europe, and their approach is adopted in this study. Initially, drought conditions are defined using the Standardized Precipitation Index (SPI; McKee et al., 1993). The SPI takes the standardized accumulated monthly rainfall amount into account and fits it to a gamma distribution, resulting in a mean of 0 and a standard deviation of 1. Values below zero, therefore, describe conditions with below-average accumulated rainfall, whereas values above zero depict wetter-than-average conditions. As we focus on agricultural drought, we choose the 3-monthly aggregated version, i.e., the SPI-3. According to Edwards and McKee (1997), drought conditions are present when the SPI-3 is lower than −1 and when the precipitation level of the preceding 3 months is at least 1 standard deviation lower than average. This is also applied to this study. We define drought conditions as ending with the first day for which the SPI-3 is greater than −1.

Regarding the definitions of heatwave events, as the requirements are diverse and depend on the user and the research question, multiple approaches exist with no general specification. Derived from Perkins and Alexander (2013) and Ionita et al. (2021), we apply a percentile-based threshold to define heatwave events (HWEs). We define an HWE as being present when the daily maximum temperature exceeds the monthly 90th, 95th, and 99th percentiles for at least 6 consecutive days. That is, the sixth day of a period for which this condition is true counts as 1 towards the HWE statistic. If the heat period has a duration of 9 consecutive days, day numbers 6, 7, 8, and 9 will count towards the HWE, as the precondition is given for all of these days. Correspondingly, a heat–drought compound event (HDCE) is defined as the joint occurrence of an HWE and an SPI-3 below −1, and every calendar day for which the HWE and SPI-3 conditions are given counts towards the HDCE statistic. As the hazard potential of HDCEs is highest during the summer months, we analyzed HDCEs for June, July, and August.

4 Results

Firstly, the performance of the MBCn method regarding the overall long-term climatology and specifically the estimation of percentile-based thresholds is evaluated. This is a crucial component in estimating the quality of the projections of FSEs and HDCEs. Subsequently, the prerequisites, the historical and projected frequencies of FSEs and HDCEs, and the corresponding near-surface atmospheric deviations linked to these events are presented.

Figure 2 Evaluation of MBCn performance for long-term climatology and daily-based percentiles. (a) Long-term time series of annual mean maximum temperature for ERA5 and 13 CORDEX models in raw (left) and MBCn-corrected output (middle), with 30-year mean values (right). (b) The same as (a) but for annual precipitation sum. (c) Percentile-based distributions of daily maximum temperature for ERA5 and 13 CORDEX models in historical raw data (left) and MBCn-corrected historical and projected output (middle), with mean absolute error (right). (d) The same as (c) but for daily precipitation.

4.1 Bias correction performance

In terms of long-term climatology for the spatial mean, it becomes clear from Fig. 2 that the regionally downscaled output from CORDEX still inherits a substantial amount of bias regarding ERA5. For the annual mean of daily maximum temperature (Fig. 2a), the raw output of the 13 models lies within a range of roughly 2.5 °C, with most of the models underestimating the value indicated by ERA5. After applying MBCn, the 30-year mean value of all models, and therefore also the model mean, aligns with the ERA5 mean at 19.5 °C. As the daily minimum temperature inherits a close statistical relationship to the daily maximum temperature, the performance of bias correction is similar, and the results are given in Appendix A. For the annual precipitation sum (Fig. 2b), the uncertainty range of the raw output spans almost 200 mm, with an equal number of models over- and underestimating the sum given in ERA5. However, the bias of overestimating models is higher than the bias of underestimating models, leading to a slight positive bias in the model mean of around 75 mm. After performing MBCn, the bias is significantly reduced for all models. While some models show a slight underestimation of the annual rainfall amount after performing MBCn, the uncertainty range of all models is reduced to under 10 mm.

As absolute and percentile-based thresholds are crucial components of the CEE definitions used in this study, we also inspected the ability of MBCn to align the distribution of percentiles derived from the daily data with that of the reference data. As can be seen in Fig. 2c for daily maximum temperature, the quality of raw CORDEX varies between models. In extreme cases, the mean absolute error is more than 1.5 °C. After MBCn, percentile distributions of the models are aligned perfectly, which may also be expected from a quantile fitting method. The indicated percentile values for the MBCn-corrected projection data (illustrated by colored circles) demonstrate the significant changes in daily maximum temperature that are projected by the models, with increases of up to 7 °C. For daily precipitation, the error margin of the raw model output spans up to 0.4 mm. After MBCn, the bias is significantly reduced and lies within a range between 0 and 0.01 mm. The projection data for the future period indicate no clear direction of change, but the majority of models show a reduction in daily precipitation for moderate-precipitation events and an increase in precipitation intensity for extreme-precipitation events.

Figure 3 Historical DOY (a–d) and projected changes in the DOY (e–h) of the last day of frost (LDF) and the start of the growing season (SGS). (a, e) LDF, (b, f) SGS₀, (c, g) SGS₅, and (d, h) SGS₁₀.

4.2 False-spring events

4.2.1 Last day of frost (LDF) and start of the growing season (SGS)

Crucial indicators for changes in the possibility of FSEs occurring are changes in the LDF and the SGS, which are displayed in Fig. 3. In general, three gradients can be derived. Firstly, the LDF and SGS tend to happen earlier in the year with decreasing latitude, which is consistent with an increasing temperature level under a decreasing solar zenith angle. In addition, the LDF and SGS tend to happen earlier in the year in maritime regions under a moderate oceanic influence, whereas more continental areas show a later occurrence. Also, the LDF and SGS tend towards later occurrences with increasing altitude due to decreasing temperature levels. Throughout the low-level areas of southern Europe, the LDF mostly occurs in March or early April. Moving southward, the LDF occurs earlier, around early February, and reaches late January in southern Portugal and Spain, northern Africa, and western Asia. In mountainous areas, depending on altitude, the LDF tends to occur between May and early June. All areas have a projected decrease in the DOY of the LDF within the future period in common. This decrease is mostly around 20 to 30 d, or around 1 month, for most parts of southern Europe, corresponding to increasing temperature levels (Fig. 2). With increasing altitude, the DOY decrease reaches its spatial maximum of around 1.5 months. In coastal southern European areas, northern Africa, and western Asia, where frost events are in general rare, the reduction is less prominent.

The SGS with a temperature threshold of 0 °C (SGS₀, Fig. 3b and f) occurs particularly early in the year in most western and southern parts of the domain, except for in mountainous regions. Moving east, the date increases from middle-to-late January to early March. In elevated areas, SGS₀ starts between April and June, depending on latitude. The projected changes in SGS₀ show the largest decrease in mountainous areas, with reductions in the DOY of up to 2 months. An exception is the Alps, for which the highest situated areas still show a remarkable decrease, but this is less dominant than in the Alpine foreland. In terms of low-level regions, especially the continental European and Turkish regions east of 10° E appear striking, with reductions of more than 30 d.

The spatial characteristics of SGS₅ are similar to SGS₀. However, corresponding to the higher thermal threshold, the occurrence is later within the year. For most regions, SGS₅ occurs around 1.5 to 2.5 months after SGS₀. Within high-elevation areas, the time difference between SGS₅ and SGS₀ is smaller. The projected changes in SGS₅, however, substantially differ. Throughout the domain, the decreases in the DOY reach 30 to 50 d, with larger decreases in the western parts. While mountainous areas are projected to experience a decrease as well, this is less dominant than in the lower-elevation regions. This is particularly true for the Alps.

The spatial characteristics of SGS₁₀ closely resemble those of SGS₅. The shift between the two is around 1 to 2.5 months and therefore almost linear when comparing SGS₀, SGS₅, and SGS₁₀. The highest-elevation parts of the Alps show no occurrence where nightly temperatures do not remain above 10 °C within the inspected period. In terms of projected changes, the spatial structure of the continental regions resembles that of SGS₅, although the decrease is in general on a lower level, reaching 20 to 30 d. Most noticeable here are most of the coastal regions, where the projected decrease is the highest, with reductions of up to 2 months.

It becomes apparent from Fig. 3 that the reductions in the DOY for the LDF are disproportional to the reductions in the DOY of the SGS. The level of disproportionality is dependent on the thermal definition of the SGS as well as on regional characteristics such as latitude, altitude, and distance from the sea. From a statistical point of view, only SGS₀ tends to occur before the LDF, leading to a potential expectancy that FSEI₀ may be high in the study domain but not FSEI₅ or FSEI₁₀. However, due to disproportional reductions of the SGS compared to the LDF and with many regions projected to experience a larger reduction in the DOY of the SGS than for the LDF, the risk of an increased FSEI rises.

Figure 4 Historical occurrence and projected changes of the false-spring-event index (FSEI): FSEI₀ in the first column, FSEI₅ in the second column, and FSEI₁₀ in the third column. Historical occurrence in the first row (a, e, and i), the difference in projected changes in the DOY of the LDF and SGS in the second row (b, f, and k), projected changes in FSEI in the third row (c, g, and l), and projected changes in FSEI in percent in the fourth row (d, h, and m).

4.2.2 False-spring-event index (FSEI)

The historical and projected changes of FSEI are displayed in Fig. 4. Following the assumed distributions described above, FSEI₀ is high in almost all of southern Europe and western Asia, as well as in some parts of northwestern Africa, and lies mostly in a range between 15 and 25 (Fig. 4a). By definition, this equals an occurrence ratio of 50 %–83 %. Over the ocean, northeast Africa, and mountainous regions, FSEI is low. On the contrary, FSEI₅ (Fig. 4e) and FSEI₁₀ (Fig. 4i) are in general low within the study domain, with the former reaching around 5 in many parts of the domain and the latter mainly not occurring over land areas in the historical period. An exception within FSEI₅ are the coastal areas of western Europe and Türkiye, where the occurrence is highest at around 10 to 15, equaling 33 %–50 % of all years.

Regarding the projected changes in the future period, Fig. 4b, f, and k display the difference between the decrease in the DOY of the SGS and LDF. In areas marked red, the SGS is projected to retract more quickly towards the beginning of the year compared to the LDF, therefore potentially increasing the risk of FSEs. The opposite is projected for areas marked blue, i.e., a faster retraction of the LDF compared to the SGS, thus potentially lowering the risk of FSEs. The third row of Fig. 4 shows the actual projected changes in FSEI₀, FSEI₅, and FSEI₁₀, and a comparison of the two metrics shows a close agreement for most of the inspected regions. In terms of FSEI₀, a decrease of up to 20 events per 30-year period is projected for most of the western and southern parts of the domain. For most of the eastern domain and the mountainous areas except the Atlas, however, up to 10 events per 30-year period are projected. For eastern Europe, a similar picture is apparent for FSEI₅. However, in the western parts of the domain, the signal is reversed and depicts an increase in FSEs of up to 5. Only the most maritime parts of France, as well as parts of Spain and Portugal, show matching signals for FSEI₀ and FSEI₅. For FSEI₁₀, especially the non-mountainous areas of the study region show a potential increase in the number of events, with mostly one additional event per 30-year period. Low-level southwestern France, where FSEI₀ and FSEI₅ are projected to predominantly decrease, becomes most dominant for FSEI₁₀, however, with the highest increase of up to two events. It should be noted, especially for FSEI₅ and FSEI₁₀, that some regions appear blue regarding the disproportionate retraction of the LDF and SGS but already show no occurrence of FSEs in the historical period. Therefore, no further decline in FSEI can be detected.

The percentage change of FSEI, shown in the bottom row of Fig. 4, shows a more pronounced manifestation of changes in FSEI₅, reaching up to a 5-fold increase in the number of FSEs, whereas increases amount to around 50 %–100 % for FSEI₀. The decreases projected for FSEI₀ amount to 50 %–75 % of the historical count. While the occurrence of FSEI₁₀ is extremely rare in the historical period, Fig. 4m shows that this event is projected to occur in almost every portion of the study domain in the future period.

Figure 5 Mean deviation of the near-surface atmosphere from the 30-year mean during frost events after the SGS for daily mean sea level pressure (SLP; a, e, and h) and daily minimum temperature (T_n; b, f, and i). Number of models indicating a statistically significant positive deviation (e, f) and statistically significant negative deviation (h, i). Mode of wind direction within the 30-year historical period (c) and only for frost events after the SGS (d). Number of models indicating significantly different distributions of wind directions (g).

4.2.3 Atmospheric variations related to FSEs

The deviation from the mean atmospheric state for sea level pressure (SLP) and daily minimum temperature (T_n) for frost events after SGS₀ is shown in Fig. 5 (SLP – Fig. 5a, e, and h; T_n – Fig. 5b, f, and i). Note that the deviations in the historical period are shown, whereas the deviations in the future period are given in Appendix B. In general, frost events after the SGS are accompanied by above-average SLP, indicating the potential predominance of high-pressure systems. There is a high level of agreement within the 13 models, as shown by the high number of models indicating this positive SLP deviation (Fig. 5e). In addition, these deviations are statistically significant in all 13 models for the majority of areas. A notable exception is the high-elevation mountainous regions, where SLP deviations show a negative deviation. As can be expected, T_n deviations are significantly below average in most low- and mid-elevation areas, where SGS₀ happens particularly early within the year. In higher-elevation areas, where SGS₀ is closer to summer, the deviations turn towards positive values, as frost events occurring post-SGS tend to be warmer than in winter.

Following the climatology of the study domain (see Sect. 2), the wind most often originates from southwesterly to northwesterly directions in the less elevated areas. This changes to northerly to easterly directions with decreasing latitude, where the influence of the westerlies decreases. Over the mountain ridges, however, local anomalies appear dominant in the mode of the wind direction. Here, wind most often blows from the main ridge towards the lowlands, resulting in, for example, predominant southwesterly winds to the north and east and northeasterly winds to the south and west of the Apennines. The interplay of large-scale atmospheric flow and local anomalies causes the complex structure shown in Fig. 5c. During frost events post-SGS, this structure is significantly changed (Fig. 5g), with predominant easterly and northeasterly winds throughout the lowlands (Fig. 5d). Along mountainous areas, the mode of the wind direction appears similar, with winds blowing from the ridge to the lowlands. However, a high number of models still indicate statistically significant shifts in the wind direction distribution.

Figure 6 The same as Fig. 5 but for FSEI₅.

The general picture is similar for frost events after SGS₅ (Fig. 6). But with a decreasing number of frost events and therefore a reduced sample size, fewer models indicate statistical significance. For SGS₁₀, with even fewer events, no clear picture can be derived, and the corresponding plot is therefore given in Appendix B. The overall picture indicates that frost events after the SGS are most often accompanied by cold, easterly airflow under high-pressure conditions, which all significantly deviate from the 30-year mean. In mountainous regions, descending air masses flowing from the ridge to the lowlands are predominant.

Figure 7 Projected changes in the deviations of daily mean sea level pressure (SLP; a, c, and e) and daily minimum temperature (T_n; b, d, and f) from the 30-year mean during frost events after the SGS. SGS₀ (a, b), SGS₅ (c, d), and SGS₁₀ (e, f).

In the future, there is no indication of a change to the general picture of these characteristics. While cold, easterly flows under high pressure are also predominant in the projections, fewer models indicate a statistically significant deviation from the 30-year mean (Appendix B). Comparing the mean atmospheric deviation in the future and historical periods, it can be seen in Fig. 7 that the models project an increase in the positive deviation of SLP from the mean state for most of eastern Europe and western Asia.

Figure 8 Historical (a–c) and projected changes (d–f) in the 90th (a, d), 95th (b, e), and 99th (c, f) percentile of daily maximum temperature. Projected deviation of the SPI-3 from the historical mean in the future period (right column) annually (g) and in summer (h).

The homogeneous spatial indication of this change for SGS₀ is less uniform for SGS₅, where the majority of cells indicate an increase in positive deviation mixed with cells indicating a decrease. In Spain, Portugal, and northwestern Africa, however, the projected development is the opposite, indicating a general decrease in positive SLP deviation. In terms of T_n, regions with varying changes in the deviation mix closely for SGS₀ with no clear spatial distinction. For SGS₅, which is in general projected to occur earlier within a year, the temperature deviation is projected to decrease in correspondence to a higher potential for lower temperatures closer to the winter period. In general, not accounting for regional variations, the projections indicate an intensification of the high-pressure anomaly during frost events post-SGS, with a higher potential for lower temperature anomalies than in the historical period.

4.3 Compound events of heat and drought

4.3.1 SPI-3 and percentiles of T_x

Prerequisites for the investigation of heat- and drought-related compound events are changes in the percentiles of daily maximum temperature (T_x) and the corresponding change in the exceedance of these thresholds, as well as changes in drought behavior indicated by the SPI-3. These projected changes are shown in Fig. 8.

Following the climatological description in Sect. 2, the thermal gradients regarding latitude, elevation, and distance from the sea become apparent for the historical percentile-based thresholds of T_x including the 90th, 95th, and 99th percentiles in Fig. 8a–c. The given values all refer to the meteorological summer, including June, July, and August (JJA). The corresponding projected changes under RCP8.5 are shown in Fig. 8d–f and depict a minimum increase of 4.5–5 °C for a majority of land areas in the study domain. Additionally, a southward increase in the deviation can be detected, reaching from around +5 °C in central Europe to over +6 °C in northwestern Africa. Thirdly, some mountainous areas, e.g., the Alps and the Pyrenees, show intensified warming compared to the surrounding lowlands. Besides the spatial gradients, an increase in warming can be detected for increasing percentiles, which is true for most parts of the study domain. In general, when applying the percentile-based thresholds derived from the historical period (see Fig. 2c) to the future daily T_x, the historical 90th percentile is projected to shift to approximately the 66th percentile. Therefore, while this threshold was historically exceeded on 1 of 10 d on average, 1 out of 3 d would exceed the historical threshold in the future. This tripling of the occurrence probability is surpassed by the 99th percentile, for which a 15-fold increase is projected.

Figure 9 Historical occurrence and projected changes of heat–drought compound events (HDCEs). HDCE₉₀ in (a–d), HDCE₉₅ in (e–h), and HDCE₉₉ in (i–m). Historical occurrence rate (a, e, and i) and projected future occurrence rate (b, f, and k). Average historical length of consecutive days meeting the HWE criteria (HDCE period) (c, g, and l) and projected changes in the average length of HDCE periods (d), h, and m).

The SPI-3 indicates an increased probability for drought conditions in the study domain, excluding parts of central and eastern Europe, which is reflected both annually (Fig. 8g) and for JJA (Fig. 8h). However, lower values for JJA indicate an aggravated drought situation for the summer months, compared to the annual mean. For most regions indicating a decrease in the SPI-3, the reductions are projected at around 0.5. This equals a reduction in the mean precipitation sum of around 0.5 standard deviations of the historical period. However, in the western parts of the domain, including most parts of Spain, Portugal, and Morocco, the SPI-3 is projected to decrease by more than 1 standard deviation, indicating prevailing drought conditions in the future. In the regions where ΔSPI-3 falls below −1, on average, every month is expected to meet the drought-based requisite for the occurrence of HDCEs.

Figure 10 Mean deviation of the near-surface atmosphere from the 30-year mean under HDCE₉₀ for daily mean sea level pressure (SLP; a, e, and h) and daily maximum temperature (T_x; b, f, and i). Number of models indicating a statistically significant positive deviation (e, f) and statistically significant negative deviation (h, i). Mode of wind direction within the 30-year historical period (c) and only for HDCEs (d). Number of models indicating significantly differing distributions of wind directions (g).

4.3.2 Heat–drought compound events (HDCEs)

The historical and future spatial characteristics of HDCEs, as well as the corresponding length of consecutive days meeting the HWE criteria, are shown in Fig. 9. During the historical period (Fig. 9a, e, and i), HDCE₉₀ most often occurred in the northern half of the domain, including Italy, southern France, and the Balkans, reaching occurrence rates of once every 3 to 5 years, on average. While the spatial differences are similar for HDCE₉₅, the frequency of events is reduced, mostly occurring less than once every 10 years. HDCE₉₉ is virtually not simulated by the models, showing only singular appearances within the 30-year period. Due to the low historical occurrence rates compared to the future period, the absolute occurrences, instead of Δ, are shown in Fig. 9b, f, and k. For every grid cell of the domain, increases in the occurrence of all types of HDCEs are projected. The largest frequency over land areas appears in regions close to, or bordering, the Mediterranean Sea. For HDCE₉₀, up to 20–25 d per summer is projected to occur in these regions, with frequencies dropping to around 5–10 events per year in central Europe and northeastern Africa. As with the historical occurrence, the spatial structure of the projected increases is similar for HDCE₉₅ and HDCE₉₉, with decreasing frequencies for higher thermal thresholds. Within the aforementioned regions with the highest increase, HDCE₉₅ is projected to occur around 10–20 times per summer and 5–10 times per summer for HDCE₉₉.

In accordance with the definition of HDCEs, the number of consecutive days meeting the HWE criteria (HDCE period) is at least 6. In the historical period, the average HDCE period has a length of around 7 to 9 d for HDCE₉₀, with no apparent spatial discrimination (Fig. 9c). The length is reduced to around 7 to 8 d for HDCE₉₅. The projected changes of the HDCE period are highest in the regions with the highest projected increase in HDCEs (Fig. 9d, h, and m). For HDCE₉₀, gains reach up to 15 d, resulting in an overall projected length of around 3 weeks, on average. Increases in the length of the HDCE period reach 10 d for HDCE₉₅, resulting in more than 2 weeks with consecutive days reaching the HDCE₉₅ criteria. Most land areas show an increase of 6 d for HDCE₉₉, which is mostly due to HDCE₉₉ being projected for the first time within the future period.

4.3.3 Atmospheric variations related to HDCEs

For this investigation, all days meeting the HDCE criteria were extracted, and their mean climatology was compared to that of the entire 30-year period. Due to the low number of events (see Fig. 9), the results for near-surface atmospheric variations during the historical period are given in the Appendix. The general picture of the variations detected in the future period similarly applies to the historical period (cf. Figs. 10–12 and Appendix B). As shown by Fig. 10, HDCE₉₀ is accompanied by below-average SLP conditions over most of the low-elevation areas. These negative deviations are often statistically significant, although the number of models varies.

Over mountainous regions, however, the SLP deviations are positive, indicating above-average pressure conditions. Over the Alps, Pyrenees, and Atlas Mountains, these deviations are statistically significant in a majority of models. As expectable, the T_x deviations are significantly above average. The mode of the wind direction is similar for both the entire 30-year period and for HDCEs. As with FSEs, westerly to northerly wind directions are predominant in the eastern half of the study domain. In the low-elevation areas of the northwestern parts of the domain, where westerly winds are predominant for FSEs (see Fig. 5), winds predominantly originate from the north and northeast. Mountainous areas, similarly to FSEs, are determined by regional characteristics depending on the orientation of the main ridge, causing a complex picture in the mode of the wind direction. Contrarily to FSEs, near-surface flows tend to move towards the mountain ridge, as can be seen for the Atlas Mountains, Pyrenees, and Apennines (Fig. 10c). For many land areas in the study domain, the projected distribution of the wind directions under HDCE₉₀ in the future differs significantly from the total distribution (Fig. 10g).

Figure 11 The same as Fig. 9 but for HDCE₉₅.

The overall picture of HDCE₉₅ and HDCE₉₉ is similar to HDCE₉₀, with below-average sea level pressure conditions in the low-lying areas and above-average conditions over the western mountain ridges (HDCE₉₅ in Fig. 11, HDCE₉₉ in Appendix C). Winds predominantly originate from northern to eastern directions, excluding regional phenomena in the vicinity of mountain ranges. While a majority of the 13 inspected models prove that the atmospheric variations during HDCEs differ from the overall mean with statistical significance, the number of indicating models decreases with decreasing sample size and higher thermal thresholds.

Figure 12 Projected changes in the deviations of daily mean sea level pressure (SLP; a, c, and e) and daily maximum temperature (T_x; b, d, and f) from the 30-year mean during HDCEs. HDCE₉₀ (a, b), HDCE₉₅ (c, d), and HDCE₉₉ (e, f).

The projected changes in the mean near-atmospheric state during HDCEs are shown in Fig. 12. The results described above for the future period are therefore the result of a decrease in sea level pressure for most regions. Still, some regions show an increase in sea level pressure, adding to the complexity of the projected changes. The deviation in mean T_x during HDCEs, as indicated by the general projected warming, is shown to increase by around 2.5–5 °C, and often the projected warming is more intense for HDCE₉₅ than for HDCE₉₀. Due to the low occurrence rate in the historical period, comparisons for HDCE₉₉ could not be made.

5 Discussion

The use of dynamically downscaled climate model output offers several advantages over the output of general circulation models, foremost a higher spatial and temporal resolution, as well as a higher data accuracy. However, as the results of this study prove for the GMR, a high amount of bias remains in the output of RCMs when compared to state-of-the-art reanalysis data. This bias varies for different variables and is dependent on the underlying GCM and RCM, but it is not distinguishable whether the choice of GCM or RCM contributes more bias. With one exception, temperatures are underestimated in the raw model output, which is in line with former experiences with CORDEX output (Dosio et al., 2022; Top et al., 2021). Precipitation sums, on the other hand, show an almost equal number of over- and underestimating models. However, the magnitude of the bias is higher for models with an overestimation of the precipitation sum. Reasons for the over- and underestimation of precipitation, including the “drizzle effect”, have been discussed before, e.g., by Chen et al. (2021), Demory et al. (2020), and DeMott et al. (2007). A bias is present in both the climatology of daily data, i.e., percentiles, and the long-term climatology. After applying QDM, both aforementioned forms of bias are removed to a significant degree, resulting in nearly zero bias in terms of percentiles and minimal deviations in long-term climatology, depending on the selected GCM–RCM combination and variable. Therefore, and by preserving the dependence structure of multiple corrected variables, MBCn proves to be a suitable choice for assessing percentile-threshold-based compound events. The performance of several different multivariate bias correction methods could be evaluated in follow-up studies, similarly to what has been done before for univariate correction methods applied to CORDEX (Olschewski et al., 2023; Laux et al., 2021). An additional limitation of this study lies in the choice of reference data. While ERA5 is proven to be a suitable state-of-the-art reference for temperature and precipitation (Lavers et al., 2022; Mistry et al., 2022; Velikou et al., 2022; Hassler and Lauer, 2021), some level of uncertainty regarding extreme events remains, especially for precipitation. When conducting local refinement of the presented results in potential follow-up studies, additional reference data sets, including station-based observations, could be compared to further improve the robustness of the projections.

The disproportionate behavior of the LDF and the SGS shown in this study has also been demonstrated, for example, for the United States (Allstadt et al., 2015; Peterson and Abatzoglou, 2014; Marino et al., 2011). These authors conclude that FSE-related risk varies significantly for different regions. In this context, Chamberlain et al. (2021) discuss potential reasons for the large regional variability, including elevation (Vitasse et al., 2018) and distance from the sea (Ma et al., 2019). The findings of this study also demonstrate these two drivers as decisive when determining the potential future risk of FSEs. However, the nature of the projected change is also strongly dependent on the selected SGS threshold. Regarding the SGS, the presented thresholds are strictly based on atmospheric temperatures and do not include data on, for example, soil temperatures and observations of vegetation. Apart from highly elevated regions, where the onset of the growing season is generally later in the year, the three selected thermal thresholds for the SGS are proven in this study to represent three different sections within the warming period between January and June. This offers a benefit for a future refinement of the presented results, for example for specific crops with specific thermal characteristics referring to one of the sections. However the broad nature of the thermal definitions still imposes a limitation on the results of this study. A more sophisticated approach, e.g., including crop-related modeling or crop-specific thresholds of budburst and freezing damage, for example, discussed by Allstadt et al. (2015) and Chamberlain et al. (2019), may lead to increased preciseness and applicability for practitioners. On the other hand, this study aims to provide a general overview of the potential risk of FSEs within a large domain that has not been previously closely investigated in this regard. This was done using simplified and straightforward thermal definitions. This study proves the spatial resolution to be sufficient to distinguish regional characteristics of FSEs. In terms of agricultural applications on a local level, where crop-specific metrics and local climate characteristics gain in importance, the presented results may be the starting point of local refinement.

The results regarding the atmospheric deviations under FSEs are in line with our general understanding of climatology in the study domain, and this especially applies to the low-lying areas. During spring, air masses originating from westerly and southerly directions are more and more likely to cause warm conditions, as the level of warming closely follows the northward shift of the subtropical ridge (see Sect. 2). However, especially when high-pressure conditions cause clear skies and the nights are still long, a high level of outgoing longwave radiation over land masses causes temperature levels to decrease. Under easterly flows, which this study proves to be predominant under FSEs, these cooled air masses are forced into the study domain and increase the risk of frost events. In a recent study for Europe, Quesada et al. (2023) determined the advection of northerly to easterly air masses to be crucial drivers of cold waves over Europe. These results can be directly linked to the findings of this study, increasing the predictability of FSEs. Over mountainous areas, the predominant winds differ, indicating flows from the mountain ridges towards the lowland. In general, descending air masses are linked to dissolving cloud coverage, resulting in clear skies. As described above, the longwave outgoing radiation can therefore not be retained in the lower atmosphere, contributing to low near-surface temperatures and the occurrence of FSEs. On the other hand, for example, under foehn conditions, descending air masses on the leeward side of a mountain ridge may be linked to warm conditions (Jansing et al., 2022) and therefore contradict the results of this study. To gain further insights into the local conditions caused by the large-scale atmospheric steer flow, additional vertical layers of atmospheric flow that are less sensitive to local orography, e.g., the 500 hPa layer, could be assessed in further studies. For this study, we specifically focused on the near-surface atmospheric conditions, as these are most relevant in the context of crops. A similar gain may also be achieved by considering circulation patterns, for example, based on Jenkinson and Collison (1977) or Hess and Brezowsky (James, 2007). However, the application of large-scale circulation patterns is aggravated in the context of multivariate bias correction, as the spatial coherence, alongside the temporal coherence, of all variables must be preserved for the entire study domain. As Cannon (2018) pointed out, this is, in general, possible for MBCn. However, considering the size of the study domain, the additional consideration of spatial coherence required a significantly higher computational effort that was not available within the means of this study.

The nature of compound extreme events including heat and drought over Europe and the GMR has been extensively studied for historical periods (Vogel et al., 2021; Ionita et al., 2021; De Luca et al., 2020; Russo et al., 2019). The results presented in this study embody the continuation of the historical trend of increasing HDCEs that all former studies share. Inspecting the relevant variables independently already shows the tendency of each variable to favor the occurrence of HDCEs in the future. As described above, the projected increase in daily maximum temperature for the spatial mean causes extremes exceeding the 99th percentile to be 15 times more likely. However, as a majority of regions show the potential for persistent drought during summer (JJA), the projected change in the occurrence rate of HDCEs is several times higher. To further increase the robustness of these projections, additional metrics and thresholds to define heat and drought could be considered in a follow-up study.

The aggravated risk of compounding extremes has, for example, been indicated by Zscheischler et al. (2018). The potential effects of the projected changes in heat and drought behavior are diverse. An increased risk of crop failure, for example, has been shown by He et al. (2022) on a global level and by Ribeiro et al. (2020) for a Mediterranean domain. While it must be noted that crop failure risk is highly dependent on the crop type and local climate factors, accumulated drought and heat nevertheless have negative impacts on the amount of water available for irrigation. Apart from agricultural aspects and natural ecosystems, the projected level of warming will significantly impact human health and mortality (Raymond et al., 2020; Gasparrini et al., 2015). The impact of increased temperatures is also shown to be aggravated in highly urbanized areas (Merkenschlager et al., 2023) and further amplifies socioeconomic disparities (Hsu et al., 2021).

The results of this study depict that below-average sea level pressure conditions and predominantly northeasterly to easterly flows are accompanied by HDCEs. A reason for below-average pressure could be near-surface thermal lows induced by intense surface warming, as previously described by Lavaysse et al. (2016) and Hoinka and Castro (2003). This, in combination with the advection of dry, continental air masses from eastern Europe, may be crucial for the occurrence of HDCEs. The limitations in terms of drawing conclusions on the large-scale atmospheric circulation, as discussed above, similarly apply to HDCEs. While the results clearly depict significantly differing conditions during HDCEs and offer benefits regarding the predictability of these events, a follow-up study including circulation patterns could improve the quality of these predictions even more. This also applies to the mountainous regions, which show different characteristics during HDCEs than the surrounding lowlands. For both FSEs and HDCEs, the climate projections show an intensification of the pressure anomalies. For FSEs, the already above-average conditions show a projected increase in deviation from the current mean atmospheric state. For HDCEs, there is a less clear picture, although more regions show below-average conditions for historical HDCEs. This below-average deviation is amplified in the future projections for a majority of regions in the study domain. Whether centers of action (Osman et al., 2021) and related blocking systems will intensify in the future is subject to the scientific discourse, as discussed by Kautz et al. (2022).

6 Summary and conclusions

This study aims at demonstrating characteristics and potential changes in false-spring events (FSEs) and heat–drought compound events (HDCEs) in a high-impact future scenario (RCP8.5) for the end of the 21st century. The inspected periods are 1970–1999 and 2070–2099. We applied a multivariate, i.e., a variable dependence-preserving, bias correction method based on the N-dimensional probability density function transform (MBCn) to regional climate model output obtained from CORDEX. The results prove that MBCn can successfully remove biases in the model output and therefore increase the reliability of projections. Due to MBCn representing a percentile-adjusting method, this is specifically true for the percentile distributions of jointly corrected climate variables, making it a valuable method for the derivation and inspection of percentile-threshold-based compound extreme events.

The presented results indicate that, while only FSEs of the lowest thermal threshold (SGS₀) are historically relevant on a widespread basis in the Greater Mediterranean Region, FSEs are projected to gain in frequency and relevance in the future period. This is mainly due to a disproportionate reduction in the day of the year (DOY) of the last day of frost (LDF) and the start of the growing season (SGS). Therefore, warm periods have an increased probability of occurring particularly early in the year, whereas the risk of frost events remains. In low-lying areas, FSEs are mostly coupled with high-pressure anomalies and northerly to easterly winds. This is in line with cold, continental air masses flowing into the GMR and causing an increased risk of frost events. For mountainous areas, local winds appear predominant, mostly flowing from the mountain ridge towards the lowlands.

All investigated metrics of HDCEs indicate significant changes in the future. These include an increase in intensity, with more severe phases of agricultural drought and increased daily maximum temperatures of up to 6 °C. Also, the frequency of HDCEs and the length of consecutive days meeting the HDCE criteria are projected to increase. In general, the changes in HDCEs are several times more intense than those of the underlying univariate extremes. Local thermal lows are shown as the predominant near-surface condition linked to the occurrence of HDCEs, together with the inflow of dry, continental air masses from northern to eastern directions.

The projected changes will potentially have severe negative effects on vegetation and crop efficiency if no additional adaptive measures are applied. An increased risk of frost exposure due to an earlier SGS can have adverse effects on crop yield by causing damage in the early stages of development. This is particularly true for plant species with only rare historical exposure and therefore a low degree of adaption to frost. The projected level of heat- and drought-induced stress could potentially aggravate the agricultural and socioeconomic stress even more. This could, for example, manifest itself in the form of severe water shortages for humans, livestock, and vegetation; stress on food security due to the increased risk of crop losses or failures; increased tree and forest mortality endangering natural ecosystems; and deteriorated human health and increased mortality.

While the public focus is often on the effects of increased temperature levels due to global warming, this study demonstrates the additional need for adaptation on the other end of the thermal spectrum. To prevent vegetation and crop damage or loss, which would further amplify the already-aggravated ecological and nutritional stress, efforts to protect vegetation and crops should equally focus on heat, drought, and frost. The results of this study can offer benefits in this regard to allow for better-informed adaptation strategies in agriculture. An additional benefit is given by increasing the level of predictability of FSEs and HDCEs. This will allow local actors to react to forecasts that indicate a potential occurrence in a timely manner. The data sets generated in this study can additionally be used for further climate change impact analyses that, for example, focus on a local refinement or specific crop species in the GMR. By applying more complex definitions of compound extreme events, for example, when additional predictors or higher-resolution data sets are available, potential follow-up studies will be able to further improve the assessment capability of compound extreme events.

Appendix A: List of abbreviations used in this study

Abbreviation	Name
CDF	Cumulative distribution function
CEE	Compound extreme event
CMIP5	Fifth phase of the Coupled Model Intercomparison Project
CORDEX	Coordinated Regional Climate Downscaling Experiment
DIR	Daily mean wind direction
DOY	Day of year
ECMWF	European Centre for Medium-Range Weather Forecasts
EQM	Empirical quantile mapping
ERA5	Reanalysis product provided by ECMWF
FSE	False-spring event
FSEI	False-spring-event index
GCM	General circulation model
GMR	Greater Mediterranean Region
HDCE	Heat–drought compound event
HIST	Historical time period (1970–1999)
HWE	Heatwave event
IPCC	Intergovernmental Panel on Climate Change
LDF	Last day of frost
MBCn	Multivariate bias correction using the N-dimensional probability density function transform
N-pdft	N-dimensional probability density function transform
PR	Daily precipitation sum
QDM	Quantile delta mapping
RCM	Regional climate model
RCP	Representative Concentration Pathway
SGS	Start of the growing season
SLP	Daily mean sea level pressure
SPI-3	3-monthly Standardized Precipitation Index
Tn	Daily minimum temperature
Tx	Daily maximum temperature
WCRP	World Climate Research Programme

Appendix B: Additional material for bias correction performance

Figure B1 Evaluation of MBCn performance for long-term climatology and daily-based percentiles. (a) Long-term time series of annual mean minimum temperature for ERA5 and 13 CORDEX models in historical raw data (left) and historical MBCn-corrected output (middle), with 30-year mean values (right). (b) Percentile-based distributions of daily minimum temperature for ERA5 and 13 CORDEX models in historical raw data (left) and MBCn-corrected historical and projected output (middle), with mean absolute error (right).

Appendix C: Additional material for FSEs

Figure C1 Mean deviation of the near-surface atmosphere from the 30-year mean during frost events after SGS₀ for daily mean sea level pressure (SLP; a, e, and h) and daily minimum temperature (T_n; b, f, and i). Future deviations under RCP8.5. Number of models indicating a statistically significant positive deviation (e, f) and statistically significant negative deviation (h, i). Mode of wind direction within the 30-year future period (c) and only for frost events after SGS₀ (d). Number of models indicating significantly differing distributions of wind directions (g).

Figure C2 The same as Fig. C1 but for SGS₅.

Figure C3 The same as Fig. C1 but for SGS₁₀ in the historical period.

Figure C4 The same as Fig. C1 but for SGS₁₀ in the future period.

Appendix D: Additional material for HDCEs

Figure D1 Mean deviation of the near-surface atmosphere from the 30-year mean under HDCE₉₀ for daily mean sea level pressure (SLP; a, e, and h) and daily maximum temperature (T_x; b, f, and i). Future deviations under RCP8.5. Number of models indicating a statistically significant positive deviation (e, f) and statistically significant negative deviation (h, i). Mode of wind direction within the 30-year future period (c) and only for HDCEs (d). Number of models indicating significantly differing distributions of wind directions (g).

Figure D2 The same as Fig. D1 but for HDCE₉₅.

Figure D3 The same as Fig. D1 but for HDCE₉₉ in the historical period.

Figure D4 The same as Fig. D1 but for HDCE₉₉ in the future period.