Characteristics of agricultural droughts in CMIP6 historical simulations and future projections

This study explores changes in agricultural drought event characteristics in projections of Earth System Models (ESMs) participating in the Coupled Model Intercomparison Project Phase 6 (CMIP6) for different future scenarios based on three Shared Socioeconomic Pathways (SSP). To quantify the intensity of agricultural droughts, the 6-month Standardized Precipitation Evapotranspiration Index (SPEI6) with a 65-year reference period is applied to the simulations of 18 ESMs. In a first step, these ESMs are evaluated based on performance metrics and pattern correlations of drought related variables including precipitation and approximated reference evapotranspiration with reanalysis datasets including ERA5 and CRU. With this we extend the model benchmarking performed in the third chapter of the IPCC AR6 by 15 years and additional variables. In a second step we analyze global and regional projected SPEI6 distributions to estimate and characterize the changes in agricultural drought in the future based on multi-model means of change rates, distributions and relative area covered by specific events. We quantify the change of drought index values for 42 IPCC AR6 WG1 reference regions individually with a focus on those with most harvest area and find negative trends in water budget and SPEI for higher emission scenarios in most of them, particularly in the Mediterranean and other arid regions. This agrees with other recent studies. Increasing reference evapotranspiration emerges as the dominant driver for drier conditions in these regions. What is considered as the driest 2.3 % months during 1950–2014 is projected to be the new normal or moderate condition in arid regions by 2100, following a high emission future scenario (SSP5-8.5). For this scenario, 40 % of the harvest regions surface is considered to be under extreme drought conditions during Northern Hemisphere autumn. Under a low emission scenario (SSP1-2.6) with an expected global warming of 1.8 °C it would be less than 10 %. Our results show a significant difference between future scenarios regarding distribution shifts and spatial extent of extreme drought conditions in harvesting regions and can serves as a foundation for further impact and mitigation studies.

1 Introduction

The availability of water is an essential requirement for human life, and crucial for food production, freshwater supply and industry. If water is sparse, droughts occur which are dangerous and complex natural hazards. Droughts can affect more people than most other natural disasters (Wilhite, 2000). Changes in their characteristics and their impacts in a changing climate are still not fully understood. Considering a drought as an exceptional dry event requires some normal climatic conditions as reference, which is one of the challenges in defining droughts. A definition further depends on the field of application and can be roughly clustered into different types of drought. A decline of precipitation over several days to months is often referred to as a meteorological drought. Low soil moisture and high plant water stress are indicators for agricultural droughts. Ongoing precipitation deficit in combination with runoff and evapotranspiration can impact freshwater reservoirs like lakes, rivers and groundwater, which are indicators for hydrological droughts (Wilhite, 2000).

Recent events like Northern Hemisphere mega drought in 2022 (Montanari et al., 2023; Xu et al., 2023; Schumacher et al., 2024) or European multi year droughts since 2016 (García-Herrera et al., 2019; Rakovec et al., 2022) have been analyzed and linked to changes in climate (Yu et al., 2023). The Sixth Assessment Report (AR6) of the Intergovernmental Panel on Climate Change (IPCC) found an agreement in recent scientific research with at least medium confidence in an increase of agricultural droughts for several regions of the world (Masson-Delmotte et al., 2021; Seneviratne et al., 2021). The recent increase is more profound than detected for meteorological droughts and may implicate serious consequences for humans.

This motivates the investigation of the impact of changing climate conditions on the processes of formation and characteristics of agricultural droughts in the future. In the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016a) future scenarios for green-house gas (GHG) emissions and land-use change based on Shared Socioeconomic Pathways (SSP) are described to be used in model projections. There is ongoing research in analyzing droughts using different methods and metrics based on such pathways (Bakke et al., 2023; Vicente-Serrano et al., 2022; Balting et al., 2021; Zeng et al., 2021; Zhao and Dai, 2022; Xu et al., 2019).

Studies focusing on (semi) global characteristics agree on a drying trend for the Mediterranean, Europe and parts of America and southern Africa in higher emission scenarios according to the SPEI (Bakke et al., 2023; Vicente-Serrano et al., 2022; Zeng et al., 2021) and other drought indices (Zhao and Dai, 2022; Tabari and Willems, 2022).

We build on this research effort and analyze droughts categorized by SPEI in a low emission scenario following SSP1-2.6 (expected global warming of 1.8 K at the end of the century) in addition to SSP2-4.5 and SSP5-8.5. Further, we investigate projected SPEI distributions for individual regions, to show the regional dependency of responses in drought characteristics to changes in climate. By examining the projected changes in drought distribution and the area covered by normal to extreme drought events in agricultural regions under different climate forcings, this study seeks to establish a scientific foundation for subsequent impact assessments and inform decision-making.

We provide a selective evaluation of the multi-model CMIP6 ensemble with different reanalysis data for the variables that are used throughout the study. The methods we use to achieve that follow recommendations from the World Meteorological Organization (WMO), the Food and Agriculture Organization of the United Nations (FAO) and the American Society of Civil Engineers (ASCE). They are described in detail in Sect. 2. In Sect. 3 the selection of models from the Coupled Model Intercomparison Project Phase 6 (CMIP6; Eyring et al., 2016a) and the datasets used for their evaluation are listed and described. Section 4 describes the results and splits into three parts. The first focus is on the evaluation of historical simulations over the period 1980–2014 followed by a temporal and spatial analysis of SPEI in different regions and future scenarios. In the third part we discuss the limitations of the shown results.

2 Methods

2.1 Reference Evapotranspiration

The Reference evapotranspiration (ET₀) also referred to as Potential Evapotranspiration (PET) or Atmospheric Evaporative Demand (AED) is a measure for the amount of water that the atmosphere can take up from a saturated reference surface with unlimited water supply. Following Allen et al. (1998) we use watered crop land with an albedo of 0.23 as reference surface.

Over the last century a variety of methods have been developed to approximate the PET, including the methods introduced by Thornthwaite (1948) and Hamon (1963). Although used often, both approximations are not reliable under changing climate conditions (Shaw and Riha, 2011). In contrast, Eq. (1) provides a more complex combination of parameters to approximate ET₀ with a closer link to the underlying physical processes. This approximation is more adequate and recommended by FAO and ASCE (Pereira et al., 2015; Allen et al., 1998; Technical Committee on Standardization of Reference Evapotranspiration, 2005). It is derived from Penman (1948) and Monteith (1965) and includes several variables like the soil heat flux density G, mean daily air temperature at 2 m T_d, wind speed at 2 m u₂, the psychometric constant γ, saturation and actual vapor pressure e_s and e_a, such as net radiation at plants surface R_n and heat flux density G. Δ is the slope of the vapor pressure curve (Allen et al., 1998). Instead of 0.12 m short grassland we assume a crop height of 0.5 m following Walter et al. (2001). With adjusted crop coefficients (1600 and 0.38) for monthly time steps this leads to the ASCE equation for standardized reference crop evapotranspiration for tall surfaces ET_RS (Technical Committee on Standardization of Reference Evapotranspiration, 2005):

$$\begin{array}{}\text{(1)}& {\text{ET}}_{\mathrm{0}}={\mathrm{ET}}_{\mathrm{RS}}={\displaystyle \frac{\mathrm{0.408}\mathrm{\Delta }(R_{\mathrm{n}}-G)+\mathit{\gamma }\frac{\mathrm{1600}}{T+\mathrm{273}}{u}_{\mathrm{2}}(e_{\mathrm{s}}-e_{\mathrm{a}})}{\mathrm{\Delta }+\mathit{\gamma }(\mathrm{1}+\mathrm{0.38}{u}_{\mathrm{2}})}}\end{array}$$

Allen et al. (1998) and Zotarelli et al. (2010) provide detailed instructions on the calculation of the psychometric constant and how to approximate variables that may not be available in observational or model data (i.e. u₂ from wind speed at 10 m, R_n from downwards short wave radiation and T_d from monthly mean daily minimum and maximum temperature).

2.2 Standardized Precipitation Evapotranspiration Index

In contrast to its predecessor SPI, which is solely based on precipitation, the SPEI includes additional atmospheric variables through potential or reference evapotranspiration ET₀. While the SPEI, as originally proposed by Vicente-Serrano et al. (2010), uses only temperature data to approximate ET₀ (Thornthwaite, 1948), we make use of a more complex approximation described above. The SPEI is calculated from ET₀ and P by fitting a probability distributions of water budget (ET₀−P) over a chosen reference period and transfer it to a standard distribution. In contrast to precipitation in the SPI the water budget can be negative, therefore a three parameter distribution like the generalized log-logistic is used (Beguería et al., 2014):

$$\begin{array}{}\text{(2)}& F\left(D\right)={\left[\mathrm{1}+{\left({\displaystyle \frac{\mathit{\alpha }}{D-\mathit{\gamma }}}\right)}^{\mathit{\beta }}\right]}^{-\mathrm{1}}\end{array}$$

In Eq. (2) F is the probability distribution of a variable D with α, β and γ being scale, shape and location parameters. The fit is performed for each grid cell and each month of the year individually. It is further possible to apply the index on different time scales by defining the number of accumulated months. Similar to (Yu et al., 2023) we choose 6 months, which is within the range of 3–6 months of accumulation for agricultural droughts recommended by the WMO (2012). To quantify the strength of a drought (or wet spell) the probabilities are mapped to SPEI values and can be assigned to categories shown in Table 1 (Guenang and Kamga, 2014). A more detailed guide on the described calculations can be found in McKee et al. (1993) and Guenang and Kamga (2014).

Table 1 Mapping of SPEI values to their category and occurrence. These probabilities are used to fit the probability function used by the SPEI (Vicente-Serrano et al., 2010). SPEI values between 0 and −1 for the SPEI have originally been labeled as mild droughts, but are considered as normal for simplification (McKee et al., 1993).

Due to changes in precipitation and potential evapotranspiration in many regions over the last decades, the definition of the reference period is important for the resulting SPEI values. It is recommended to use a reference period of 30 years or more, aligning well with 30 year climate normals which are standard for benchmarking in climate change assessments (WMO, 2017). However, we found that calibrating the SPEI on short reference periods (30 values per month of the year) can lead to an artificial increase in the occurrence of extreme events. Prolonging the period into the past on the other hand decreases the accuracy and availability of global observations of atmospheric variables. As a compromise we choose a reference period of 65 years (1950–2014) and show the evaluation for its last 35 years (1980–2014).

For this study we used the open source SPEI R-package to calculate ET₀ and the SPEI (Beguería and Vicente-Serrano, 2023).

2.3 Model evaluation

The metrics used to evaluate the models used in this study are listed below with their corresponding equations. Spatial statistics are generally area weighted, but for temporal statistics the variable months lengths is neglected and monthly accumulated units are avoided when possible (i.e. precipitation in mm d⁻¹ rather than mm per month). All metrics and the SPEI calculation is applied for each model individually and multi-model means are calculated afterwards.

2.3.1 Pattern correlation

The Pearson Correlation Coefficient R_p can be used to measure the linear dependency between two variables X and Y. When R_p is applied along spatial coordinates we refer to it as pattern correlation R_pc. The pattern correlation between the temporal mean of a simulated variable X_sim and an observed reference X_ref combined with area weighting can be written as

$$\begin{array}{}\text{(3)}& R_{\mathrm{pc}}={\displaystyle \frac{{\sum }_i{a}_i(X_i-\overline{X})(Y_i-\overline{Y})}{\sqrt{{\sum }_i{a}_i(X_i-\overline{X}{)}^{\mathrm{2}}}\sqrt{{\sum }_i{a}_i(Y_i-\overline{Y}{)}^{\mathrm{2}}}}}.\end{array}$$

2.3.2 Root mean squared distance

In addition to pattern correlation the root mean squared distance σ is used to measure the similarity between simulations and references. To make models inter-comparable across multiple variables of different units, σ is centered and normalized by the median distance across all simulations for each variable. The normalized centered root mean squared distance

$$\begin{array}{}\text{(4)}& \widehat{\mathit{\sigma }}={\displaystyle \frac{\mathit{\sigma }-M\left(\mathit{\sigma }\right)}{M\left(\mathit{\sigma }\right)}}\text{, with }{\mathit{\sigma }}_{\mathrm{RMSD}}=\sqrt{{\displaystyle \frac{\mathrm{1}}{N}}\sum _i{\left(X_i-Y_i\right)}^{\mathrm{2}}}\end{array}$$

describes the relative performance of a model in a multi model ensemble and cannot be interpreted independently. M is the multi-model median.

2.3.3 Seasonal cycle

For a variable X the seasonal cycle of a region r is calculated from the area weighted mean over all grid cells i in that region. The average of spatial means for each month of the year m is calculated for the 35 year evaluation period 1980–2014. With a_i being the cells area weight (area relative to total area of the region) it can be written as:

$$\begin{array}{}\text{(5)}& X_{\mathrm{m},\mathrm{r}}={\displaystyle \frac{\mathrm{1}}{\mathrm{35}}}\sum _{y=\mathrm{1980}}^{\mathrm{2014}}\sum _i{a}_i{X}_i\end{array}$$

2.4 Drought characteristics

To analyze the drought characteristics, we calculate ET₀ and SPEI as described on the same regular 1° × 1° for 18 CMIP6 models, described in Sect. 3. We calculate multi-model means of the following metrics to compare them across different scenarios.

2.4.1 Decadal change rate

As a metric for the non-linear change over time we use a decadal rate for evaluating the model simulations and analyzing projections. The decadal change rate is calculated by subtracting the average of the periods first N years from the last N years and then normalizing the difference to an average change per 10 years. The result is roughly comparable to linear regressed trends.

$$\begin{array}{}\text{(6)}& D_N={\displaystyle \frac{\mathrm{120}}{M-N}}\left(\sum _{t=M-N}^M{x}_t-\sum _{t=\mathrm{0}}^N{x}_t\right)\end{array}$$

2.4.2 Event area fraction

Frequency, severity and duration of events are often used to describe drought characteristics. For large changes in climate, however, it becomes difficult to distinguish individual drought events when they emerge in time or space to large long lasting periods that are considered as a single event. This way a general shift towards drier conditions can result in an artificial decrease of frequency as shown in Xu et al. (2019). Therefore these metrics are not suited to analyze SPEI-based droughts in long-term future projections considering high emission scenarios. As an alternative approach to quantify and visualize the temporal development of global drought characteristics, we calculate the relative land surface area for certain SPEI ranges including near normal, moderate, severe and extreme droughts and wet spells for each month as a metric for spatial extent. The area fraction A_e is calculated from the number of cells in the regular 1° × 1° grid weighted by its area a_x,y for 18 CMIP6 models.

$$\begin{array}{}\text{(7)}& \begin{array}{rl}& A_{\mathrm{e}}\left(t\right)={\displaystyle \frac{\mathrm{1}}{A_{\mathrm{land}}}}\sum _{y=\mathrm{0}}^{\mathrm{180}}\sum _{x=\mathrm{0}}^{\mathrm{360}}{a}_{x,y}\cdot i_{e,x,y}\left(t\right)\\ & \text{with }i\left(t\right)=\left\{\begin{array}{ll}\mathrm{1},& \text{if SPEI in category}\\ \mathrm{0},& \text{otherwise}\end{array}\right.\end{array}\end{array}$$

Through the area fraction A_e we can visualize multidimensional SPEI data (time, space, model dimensions) for different scenarios, while maintaining some temporal and spatial information which would get lost in multi-model temporal or spatial means. This makes the event area fraction a well suited metric to analyze the change of drought characteristics in future projections, especially for comparably rare events like extreme droughts.

2.5 Harvest regions

We make use of the IPCC AR6 WG1 Land reference regions (Iturbide et al., 2020) to apply regionally integrated statistics. Our main focus are regions with a significant amount of crop production and agricultural land use. We identified these regions by their relative amount of harvest area in 2010 according to the GFSAD1KCM dataset (Teluguntla et al., 2016). Figure 1 shows a regridded combination of major and minor rain-fed and irrigated harvest area masks. The reference regions are shown as overlay in panel (a). Panel (b) shows the relative amount of harvest area for each region. Regions above 33.3 % harvest area are considered as highly relevant for agriculture throughout this study. The agricultural most relevant regions are namely Western&Central-Europe (WCE), C. North-America (CNA), S. Asia (SAS), E. Asia (EAS), E. Europe (EEU), Western-Africa (WAF), S. E. South-America (SES), W. Siberia (WSB), Mediterranean (MED), N. E. South-America (NES), S. Australia (SAU). They are highlighted in Fig. 1b.

Figure 1 Selected regions with at least 33 % harvest area. In panel (a) the land subset of the IPCC AR6 WG1 reference regions is shown as overlay (Iturbide et al., 2020). The color represents the regridded relative harvest area derived from a combination of irrigated and rain fed masks in the 1 km crop mask of GFSAD (Teluguntla et al., 2016). Panel (b) shows the relative amount of harvest area per reference region. Selected regions are outlined in black.

2.6 ESMValTool

To load and pre-process the CMIP6 model simulation output and reanalysis data, described in the next section, ESMValTool is used. ESMValTool is an open source software for community-developed diagnostic and performance metrics to evaluate ESM simulations including CMIP contributions (Eyring et al., 2020; Weigel et al., 2021; Lauer et al., 2020). Since its first release in 2016 (Eyring et al., 2016b), the software has been developed continuously and experienced technical improvements in terms of performance and usability (Righi et al., 2020). The new benchmarking and model evaluation capabilities of ESMValTool are shown by Lauer et al. (2025). Existing diagnostics for extreme event analysis are described by Weigel et al. (2021).

The ESMValTool framework contains a dedicated python module, called ESMValCore, for data handling and general processing, and a collection of scientific diagnostics. Recipes can be used to setup experiments as a specification of input data, pre-processing methods and a list of diagnostics to be applied to the data. We used the ESMValCore v2.11 to load the described data, converted them to comparable physical units, masked and regridded them to the same 1° × 1°regular grid (Andela et al., 2025).

The methods described in this section and multiple plotting routines are added to the ESMValTool framework as a set of diagnostics in the scope of this study. This allows to reproduce the results and re-use the methods for further studies i.e. analysis for upcoming CMIP7 simulations.

3 Data

To asses the severity of agricultural droughts in the future we choose 18 models that participated in the CMIP6 ScenarioMIP and provided simulations of the three Shared Socioeconomic Pathways (SSPs; O'Neill et al., 2016; Meinshausen et al., 2020) that are analyzed in this study:

• SSP1-2.6 as a low greenhouse gas (GHG) emission scenario corresponding to a forcing of 2.6 W m⁻² and an expected global warming of 1.8 K at the end of the century with respect to the pre-industrial period,

• SSP2-4.5 as the middle of the road scenario in terms of GHG emissions with an end of the century forcing level of 4.5 W m⁻² and an expected warming of 2.7 K,

• SSP5-8.5 as a high GHG emission scenario corresponding to a forcing of 8.5 W m⁻² in 2100 resulting in an expected warming of 4.4 K.

The minimum requirement for a model to be included in this study was the availability of all variables required to calculate the reference evapotranspiration (wind speed, surface pressure, total cloud cover, daily minimum and maximum temperature) and precipitation. We limited the selection to only one model per institute, to account each ESMs atmosphere and land component only once. The final selection of models is listed with some additional information in Table 2. Relative humidity is not considered as a selection requirement because we found that we get similar results by approximating actual vapor pressure as $e_{\mathrm{a}}=\mathrm{0.611}\exp \frac{\mathrm{17.27}{T}_{\min }}{T_{\min }+\mathrm{237.3}}$ following Allen et al. (1998, Eq. 48).

Throughout this study we also discuss our results with respect to projected soil moisture, which is arguably a drought indicator on its own. Vicente-Serrano et al. (2022) used simulated total column soil moisture as a reference to discuss droughts. It might be a useful indicator for hydrological drought types at long time scales. However, the uncertainties between the predictions of different CMIP6 models increase with depth and the total depth of all layers is not the same between the models. This makes a meaningful comparison of integrated soil moisture between CMIP models difficult, and also leads to even more uncertainties. Similar discrepancies for CMIP5 are described by Berg et al. (2017). Since the maximum root depths of the biosphere is less than 2 m in most regions of the world (Müller Schmied et al., 2021) and more than 50 % of the roots of most crops are in the upper 20 cm of the soil (Fan et al., 2016), we assume the water content in the upper 10 cm of soil (surface soil moisture) is a more reliable indicator for agricultural droughts in the context of model simulations than the total column. Another benefit of surface soil moisture is the availability of this quantity from satellite observations and reanalysis to validate models for the historical period. We assumed a density of 1000 kg m⁻³ for converting the surface soil moisture to volumetric percent. This way surface soil moisture values from the CMIP6 model simulations are comparable to the reanalysis and observations listed below. The 1 km crop mask of the GFSAD dataset Teluguntla et al. (2016) has been used solely for selecting regions of interest.

Dix et al. (2019)Semmler et al. (2019)Xin et al. (2019)Swart et al. (2019)Lovato et al. (2021)Voldoire (2019)EC-Earth (2020)Li (2019)Song et al. (2019)John et al. (2018)NASA/GISS (2020)Volodin et al. (2019)Boucher et al. (2019)Byun et al. (2019)Shiogama et al. (2019)Schupfner et al. (2021)Yukimoto et al. (2019)Sellar et al. (2019)

Table 2 Available CMIP6 models providing variables to calculate ET₀ and mrsos.

To evaluate the model data we choose ECMWF's monthly averaged ERA5 reanalysis as primary reference (Hersbach et al., 2020). The successor of the former commonly used ERA-Interim dataset does not provide true measurements, but a wide range of variables and global coverage. It is based on the Integrated Forecasting System (IFS) Cy41r2 operating on a 31 km spatial and hourly temporal resolution on 137 vertical levels. In the ERA5 reanalysis billions of observations including data from over 200 satellite mounted instruments, in situ observations, balloons and aircraft measurements (Hersbach et al., 2020) are assimilated. Full spatial and temporal coverage enables direct comparison with CMIP6 model output. Most variables are available in the ERA5 monthly averaged data on single levels from 1940 to present dataset provided by the CDS (Hersbach et al., 2023). Daily maximum and minimum temperatures have been derived from ERA5 hourly data on single levels (Copernicus Climate Change Service, Climate Data Store, 2023) using the CDS toolbox¹.

In addition to the ERA5 dataset we use the Climatic Research Unit (CRU) gridded Time Series version 4.07 (TS4.07) alternative reference for temperature and precipitation in the evaluation process (Harris et al., 2020). The CRU dataset also provides approximated ET₀ based on the Penman-Monteith equation, similar to Eq. (1). Using pre-calculated ET₀ values in one reference dataset also serves as a sanity check for the implementation of Eq. (1). For soil moisture the CDS's combined satellite dataset (CDS-SM; Dorigo et al., 2019) is used. The dataset is blended from a set of passive (starting 1978) and active (starting 1991) microwave measurements with a resolution of 0.25 °. A detailed list of used sensors can be found in the datasets user guide (Preimesberger et al., 2025).

4 Results

In the following we evaluate the last 35 years (1980–2014) of historical simulations by 18 CMIP6 models with ERA5 as primary reanalysis dataset. Characteristics of agricultural droughts in future projections (2015–2100) by the same models for three different scenarios are analyzed and discussed in Sect. 4.2.

4.1 Evaluation of historical simulations

We are interested in the capabilities of models to project long-term characteristics like mean and change rates of drought relevant variables. Therefore we compare these variables with observations and reanalysis with time integrated metrics for the period 1980–2014 to estimate the reliability of these simulated variables. The pattern correlations of all 18 analyzed models, CRU TS v. 4.07 (Harris et al., 2020) and CDS-SM data (Dorigo et al., 2019) with ERA5 (Hersbach et al., 2023) are shown in Fig. 2.

Figure 2 Overview of the CMIP6 multi-model performance compared to ERA5 reanalysis over land surface. In panel (a) pattern correlation of historical CMIP6 simulations with ERA5 for precipitation (P), reference evapotranspiration ET₀, daily minimum and maximum temperature (T_min, T_max), wind speed at 10 m (V_10 m), surface down welling shortwave radiation (RS_down), surface pressure (p_surf) and soil moisture (SM_surf) are shown. The data is averaged over the period 1980 to 2014. The additional markers show pattern correlation of ERA5 to CRU and CDS-SM. In panel (b) the relative model performance for the same variables is shown. The root mean squared distance (RMSD) of each model to ERA5 is centered and normalized at the median RMSD. Negative (blue) values indicate that the model predicts results more similar with ERA5 compared to the other models. Red colors indicate worse performance. For some variables the performance compared to CRU is shown as well in the bottom right triangle. The alternative reference dataset for SM_surf is CDS-SM.

The temporal averaged pattern of CMIP6 models over the period 1980–2014 show high similarity (P>0.8 for most models) with the ERA5 reference data for surface downwelling shortwave radiation (RS_down), reference evapotranspiration (ET₀), surface pressure (p_surf) and daily minimum and maximum temperature (T_min, T_max). Lower pattern correlation ($\mathrm{0.6}<P<\mathrm{0.8}$) is found for precipitation (P). The agreement of wind speed at 10 m (V_10 m) and soil moisture (SM_surf) patterns is the lowest with correlation coefficients P<0.6 for most models, also with a higher spread between the models. The averaged soil moisture pattern of the CDS-SM dataset based on satellite observation does not show higher correlation to ERA5 than the models, which makes it difficult to judge model performance based solely on ERA5. However, we see high agreements in the most important input variables for ET₀ (RS_down, p_surf, T_min and T_max). This propagates to a ET₀ pattern correlation above 0.8 for 12 out of 18 analyzed models.

Figure 2b shows the centered median between each models root mean squared distance $\widehat{\mathit{\sigma }}$ (Eq. 4) when compared with ERA5 (upper left rectangles). For the models providing soil moisture the CDS-SM dataset is used as an alternative reference REF2 (bottom right). For the other variables CRU is used as alternative, where available. Due to the normalization, this plot lacks information about the actual distance to the reference, but instead highlights the relative performance between the models for multiple variables. A blue color indicates smaller errors to the reference and therefore better performance compared to the other models. Errors larger than multi-model median are shown in red. Biases in model simulations directly contribute to RMSD, while the normalized SPEI is able to ignore them. Therefore, it is possible that models with high RMSD are still able to project drought conditions reliably. The AWI-CM-1-1-MR and MIROC6 models for example show the highest RMSD to ERA5 for ET₀ and T_min. BCC-CSM2-MR, CanESM5 and FGOALS-g3 simulate ET₀ more similar to ERA5 but with large differences to the alternative reference dataset CRU. Throughout all models ACCESS-CM2, CMCC-ESM2, CNRM-CM6-1, EC-Earth3-Veg-LR, FIO-ESM-2-0, GFDL-ESM4, and MPI-ESM1-2-LR perform better in simulating ET₀ than the ensemble median in respect to both reference datasets. For soil moisture we find low agreement in RMSDs between reference datasets ERA5 and the CDS-SM satellite observations. We are not able to comment on model performance in these cases.

In the IPCC AR6 a similar method has been used to evaluate CMIP5 and CMIP6 model performance for the period 1980–1999 (Eyring et al., 2021, Fig. 3.42). Comparing precipitation, where ERA5 has been used for reference, we find that our centralized RMSD values are generally higher, because the median model in our subset is performing better for precipitation than the median of the 58 model ensemble analyzed in Eyring et al. (2021). Despite a few exceptions, this is in agreement with the relative ranking between the subset of 18 models for the period 1980–2014 used in this study.

Due to the importance of precipitation and reference evapotranspiration for drought formation, they are evaluated further in the following. Similar figures for other variables can be found in Appendix A (Figs. A1, A2, A3).

Figure 3 Global pattern of temporal mean precipitation based on (a) ERA5, (b) CRU and (c) CMIP6 multi-model mean over the period 1980–2014.

Precipitation is the main source of fresh water for most agricultural regions and therefore a crucial variable in drought analysis. In the multi-model mean of CMIP6 models and both reanalysis datasets we find mean patterns that qualitatively agree over the period 1980–2014 in Fig. 3. Average daily precipitation close to 0 mm can be found for dry regions in northern Africa and South-west Asia. Regions of high precipitation caused by topological boundaries like mountain ranges can be seen clearly for the Andes in South America and for the Himalayas, which forms the northern boundary of the Indian monsoon climate in the CMIP6 multi-model mean and both reference datasets. The CMIP6 simulations also agree with ERA5 and CRU reanalysis precipitation over rain forest regions of South America (8–10 mm d⁻¹) and Central Africa (3–6 mm d⁻¹).

For the eleven crop producing regions, selected in Sect. 2.5, seasonal cycles of 18 simulations and two reanalysis are shown in Fig. 4. The regional average for each month of the year is calculated according to Eq. (5) over the period 1980–2014. For each region at least 15 of the 18 models are able to reproduce the main features of the seasonal cycles derived from ERA5 and CRU.

All simulations for regions of Northern Hemisphere summer monsoon – Western-Africa (WAF), Eastern-Asia (EAS) and Southern-Asia (SAS) – agree in less than 1 mm d⁻¹ precipitation in the winter months (DJF) and maximum precipitation during the summer (JJA). For WAF the average August precipitation spreads throughout 3 mm d⁻¹ (FGOALS-g3) and 11.5 mm d⁻¹ (GFDL-ESM4), with most simulations being closer to the reanalysis around 7.5 mm d⁻¹. For SAS the precipitation during JJA is slightly underestimated by many models. The highest discrepancy to ERA5 and CRU is found in CanESM5 with 4 mm d⁻¹ less precipitation in July. EAS seasonal cycle is captured by all simulations within a maximum distance of 2 mm d⁻¹ to ERA5 or CRU.

The seasonal precipitation in Southern Hemisphere summer monsoon regions N. E. South-America (NES) and S. E. South-America is also qualitatively captured by all simulations. Mentionable is that both reanalysis show significant lower precipitation (6–7 mm d⁻¹) for NES over the first three months of the year. Pereira et al. (2017) presents seasonal cycles of Northeast Brazil over the period 1971–2000 for different reanalysis and observations including the NCEP/NCAR dataset, with 8 mm d⁻¹ average precipitation in February.

For the other regions North-American, European and Australian regions the seasonal signal of precipitation is much lower and captured by most simulations. Exceptions are the historical simulations of the FGOALS-g3 and KACE-1-0-G models. Both underestimate summer precipitation in higher latitudes and do not represent the seasonal cycle of Western & Central-Europe (WCE), E. Europe (EEU) and W. Siberia (WSB) as described by ERA5 or CRU. However, overall the regional features of the annual precipitation cycle are captured by most of the models in most of the regions, that are relevant to this study. The seasonal temperature cycles are recreated much closer by all models, with only one exception (MIROC6) overestimating temperatures in a few regions. A figure similar to Fig. 4, but for temperature, is appended as Fig. A4.

Figure 4 Seasonal cycles of precipitation for historical simulations of 18 CMIP6 models and two reference datasets calculated over the time period 1980–2014.

Figure 5 Global pattern of temporal mean of reference evapotranspiration (ET₀) based on (a) ERA5, (b) CRU and (c) CMIP6 multi-model mean over the period 1980–2014. For CRU the variable is taken from the dataset. For ERA5 and CMIP6 ET₀ is calculated from other atmospheric variables (Eq. 1) similar to the method used for CRU.

ET₀ is calculated for ERA5 and the CMIP6 model output using Eq. (1) based on the same set of variables. The ET₀ provided by the CRU dataset is derived using a similar method, but with a crop height of 0.12 m and included relative humidity to calculate the vapor pressure deficit (Harris et al., 2020; Ekström et al., 2007). The differences in the methods might cause some systematical differences, but we can still see a general agreement in the ET₀ pattern over the period 1980–2014 between all datasets in Fig. 5. Further, we can find that models reproduce expected temporal mean features like decreasing ET₀ towards high latitudes due to decreasing temperatures, low ET₀ at the high altitude Tibetan Plateau and highest ET₀ values at the Sahara due to high temperature and low cloud coverage in an arid region.

Figure 6 Annual precipitation, reference evapotranspiration and the derived SPEI. All variables are averaged over global land surface. The multi-model means of 18 CMIP6 models are drawn as solid lines, with shaded areas indicating the multi-model standard deviation. For the period (2015–2100) different future scenarios can be distinguished by color.

4.2 Projections

For this study only non glaciated land area is considered and referred to as global in the following. The global mean of P and ET₀ for the base period and different future scenarios are shown in Fig. 6. For both variables a higher positive trend can be noticed with higher emission scenarios from SSP1-2.6 to SSP5-8.5. The increase of ET₀ is significantly higher compared to precipitation. This results in a decreasing water budget for future projections and a general shift of SPEI towards drier conditions, which we can see in the corresponding annual SPEI time series, the third plot in Fig. 6. All 18 analyzed models agree in a projected decrease in global mean water budget over the period 2015–2100 in SSP2-4.5 and SSP5-8.5 scenarios. The general downward trend of water budget and a shift towards drier climate in higher emission scenarios agrees with several other studies (Balting et al., 2021; Vicente-Serrano et al., 2022).

Figure 7 Decadal change rate of SPEI with reference period 1950–2014 derived from CMIP6 projections for three future scenarios: (a) SSP1-2.6, (b) SSP2-4.5 and (c) SSP5-8.5. The change is defined as the absolute difference between 2015–2034 and 2081–2100 normalized to the average change over 10 years.

However, these changes are not distributed uniformly over the globe. We found significant regional differences in quantity of water budget and therefore SPEI changes in future scenarios as shown in Fig. 7. It shows maps of the multi-model mean decadal change of SPEI over the period 2015–2100 for three different future scenarios. The SPEI is calculated and calibrated for each model individually using 1950–2014 as common baseline for all three scenarios. While the change rates are different between the scenarios, the general pattern and signs of the trends are the same for most regions. Beside the decrease in most parts of North and South America, Africa and Southwest Asia small regions with increasing water budget can be identified, mainly in the Arctic and Subarctic as well as some mountain regions with relatively low crop production. Positive and negative changes are intensified for the SSP5-8.5 scenario compared to SSP1-2.6.

Change rates on their own are not sufficient to describe the development of projected drought characteristics. Due to regional differences we analyze the statistical distribution of the drought index for individual IPCC AR6 WG1 reference regions of similar climatic conditions (Iturbide et al., 2020). We further select eleven harvest regions based on their relative land use for crop production. The selection method and regions are described in Sect. 2. An overview of the distribution for combined and individual harvest regions is given in Fig. 8. The same figure for global distribution and all non glaciated land regions is appended as Fig. A7 and for semi arid regions of Europe and Asia as Fig. A9.

Figure 8 SPEI distributions for the CMIP6 multi-model ensemble in different future scenarios. The SPEI distributions for all grid cells in selected harvest regions are shown in panel (a) and regional distributions in form of boxplots for IPCC AR6 WG1 reference regions (Antarctic excluded) in panel (c). The hexagons in panel (b) show the mean decadal change of mean SPEI for each of the selected region but only for the SSP5-8.5 scenario. The regions are positioned relative to each other based on their real locations.

In Fig. 8a we can see the distribution of index values for the reference period 1950–2014 and the last 30 years of the projected century 2070–2100 for three future scenarios, collected over the most relevant agricultural regions. While the values in the reference period follow a Gaussian distribution, which is centered at zero by definition, the distributions of the future scenarios are shifted towards negative indices (drier climate), which agrees with the observed global downward trend in Fig. 6. We also expect more extremes for higher emission scenarios, according to widening distributions.

The hexagons in panel (b) are shown to provide a spatial overview of the regions and their projected SPEI decadal change in the SSP5-8.5 scenario. Due to the symmetric Gaussian distribution of the SPEI the decadal change of the mean is roughly proportional to the shift of the distribution. Similar to Fig. 7 the lower emission scenarios SSP1-2.6 and SSP2-4.5 show smaller change rates (Fig. A6) but qualitatively agree with Fig. 8b for most regions.

The regional statistics in Fig. 8c show the distributions as box plots for the IPCC AR6 WG1 land reference regions, ordered by median in the SSP5-8.5 scenario. For nine of the 11 regions a clear shift towards lower median SPEI values enhancing with higher emission scenarios can be found. This implies a general shift towards drier conditions and we find more severe and extreme droughts towards the end of the century in all three future scenarios based on SPEI projections of 18 ESMs. For the MED region we see the median of the drought index for 2070–2100 being below −2 following the SSP5-8.5. Which means the 2.3 % driest conditions during the simulated period 1950–2014 can be considered as normal in the Mediterranean during the projected period 2070–2100. The Mediterranean is mentioned as a hotspot for increase of extreme events including heat waves and droughts in other studies (Balting et al., 2021; Paçal et al., 2023). For East and South Asia no general shift can be found, but the number of months classified as extreme droughts is also higher compared to historical reference.

The first five of all global non-glaciated land regions (Fig. A7) SAH, ARP, MED, WCA and ECA are neighbors and part of northern Africa, southern Europe and western Asia. A general shift towards drier conditions with increasing forcing scenarios can be seen for 32 of the 42 regions. Seven regions (GIC, RAR, RFE, NEN, NWN, NEU, CAF) show an opposite shift towards wetter conditions. Six of them are located at high latitudes. While the widening of the index distribution for future scenarios can be found for all of the analyzed regions, it is most prominent in Western Africa (WAF) region.

Figure 9 Timeseries of relative drought event area according to SPEI. The color indicates how much of the agricultural relevant regions area is covered by certain drought conditions for each month. A 20 year interval from the baseline period 1950–2014 is shown on the left, while the plots on the right display relative area by drought category over the last 20 years (2080–2100) of different future scenarios.

Using the ability of the SPEI to assign a category for the severity of drought based on climatic conditions we calculated the area that is covered by drought conditions of a certain category at any time. This allows us to quantify the relative area that might be affected by droughts under changing atmospheric conditions in future projections as given in Eq. (7).

In contrast to spatial mean time series the area fraction plots show the increase in wet and dry conditions at the same time without compensating positive and negative index values.

Furthermore, area fraction plots are not influenced by localized extreme SPEI values, as they focus on the proportion of cells exceeding a certain threshold rather than the magnitude of the values themselves, unlike mean time series which can be skewed by a few cells with exceptionally high trends.

Figure 9 compares the event area fraction for an intervals during the historical period 1970–1990 with the end of century 2070–2100 over the selected harvest area. The complete series 1950–2100 is appended as Fig. A8.

The first interval in Fig. 9 presents the relative surface area for each index category during the calibration period. As each individual grid cell and month of the year is calibrated to match a normal distribution with certain shape, it is expected that we see 2.3 % of the area under extreme drought conditions without seasonal cycle on average. The interval 1970–1990 roughly represents this. However, the 2080–2100 interval selected from the end of the projected future shows significantly more area under moderate to extreme droughts for all three scenarios especially in autumn. Further the temporal resolution reveals an increasing seasonal dependency of the spatial extent of droughts.

According to the SPEI for future projections following a worst case scenario SSP5-8.5 between 20 % and 40 % of the land area in agricultural active regions would experience extreme drought conditions at the end of century. In contrast, for the SSP1-2.6 the extreme drought area does not exceed 10 %. For the drought categories we find peaks in Northern Hemisphere autumn (September–November), resulting from high evapotranspiration and low precipitation accumulated over the previous six months in spring and summer. In the extreme drought category we identify the largest differences between the scenarios. However, the peak area under at least severe (moderate) conditions also ranges from 20 % (35 %) to 50 % (60 %) throughout the scenarios.

For wet spells we find a similar increasing seasonality with opposite phase peaking in early spring. The area of wet spells with SPEI >1.5 does not clearly expand from the expected occurrence rate of 16.9 % (Table 1). It stays below 20 % on yearly multi-model average in all scenarios. Nevertheless we see a significant growth in area of extreme wet spells, which is consistent with the broadening and shifting of the global SPEI distribution for different scenarios (Fig. 8a), which overlap at SPEI between −1 and −2.

4.3 Discussion

The historical evaluation in Sect. 4.1 is meant to be a sanity check for the ET₀ calculation methods and a validation for the model projections in the first place. Considering the spread of ET₀ and P between CMIP6 model projections some sort of performance based weighting in the multi model ensemble could improve the results based on a more sophisticated evaluation taking seasonal distributions and biases into account. The differences in reanalysis and lack of reliable global observations especially for wind and precipitation is an additional challenge (Soci et al., 2024).

We decided to use the SPEI for our analysis of the development of agricultural drought characteristics, due to its ability to be applied globally while being spatially comparable over different regions. However, the statistical approach of standardized indices implies that the intensity or severity of a drought is calculated based on its rarity and not based on the actual physical soil parameters or how it impacts crop growths. Analyzing the actual water balance levels and how they impact yield for different crop types on a regional level would be a valuable addition to this study.

Beside the relative nature of the index we also have to consider that most of the area with highest increase in droughts in high emission scenarios are already arid regions i.e. the Sahara desert. Nevertheless, Singh et al. (2022) found a tenfold higher agricultural exposure to severe compound droughts towards the end of century using model projections from RCP8.5 in CMIP5.

Figure 10 Temporal average (a, b) and decadal change (c, d) in (a, c) global water budget (P−ET₀) and (b, d) soil moisture projected from CMIP6 multi-model ensemble for the SSP2.4-5 scenario (2015–2100).

While root zone soil moisture is a good indicator for agricultural droughts, the availability of reliable global soil moisture data is limited and simulations show high inter model variability (Fig. 2). Further, available water is not the only parameter responsible for crop failure. A high atmospheric evaporative demand can also cause increasing plant water stress. Figure 10 illustrates the difference between a purely hydrological approach to quantify droughts like soil moisture and approaches like SPEI, that consider the influence of ET₀. In the upper panels (a, b), we can see a that regions with an average water budget below −5 mm d⁻¹ correspond to arid regions with low soil moisture. Especially Northern Africa, South West Asia and Australia. However, under the SSP2-4.5 future scenario, soil moisture is expected to decrease over Europe and South America, but shows no or even positive changes in arid regions (Fig. 10d), because there is not much water available to evaporate or runoff (other scenarios can be found in Fig. A5). The atmospheric demand of water, shown in Fig. 10c, can further increase even if the soil does not contain water at all. However, the influence of ET₀ on crop failure and agricultural drought severity might need further research to be quantified exactly.

Due to the above mentioned reasons the SPEI might be suited to investigate the impact of changes in climate to the characteristics of droughts, but this is not directly transferable to any ecological impact regardless of the focuses on regions with lot of harvest area. With this study we hope to provide further motivation for research on regional impact. Balting et al. (2021) analyzed SPEI characteristics projected by CMIP6 models for the same three future scenarios, with a different set of models and focus on the Northern Hemisphere in summer (JJA). We do not show occurrence rates in this work, but based on regional distribution shifts and the increase of extreme dry area we can confirm the highest increase in projected droughts for the MED region, less but still increased drying in WCE EEU and WSB. By investigating the change of extreme drought (SPEI $<-\mathrm{2}$) in addition to severe droughts, we find higher difference between the future scenarios for regions where the SPEI stays below −1.5 for many months. We also complement their study by providing SPEI analysis for additional regions and over all seasons, even though it is most important for agriculture during growing season.

Although, the SPEI is symmetric, which makes it possible to apply it to wet spells in a similar way, it is not well suited to identify extreme events such as floods. Therefore, we focus the analysis on droughts, while still discussing the whole distribution of the index over all categories.

5 Conclusions

To evaluate the performance of 18 CMIP6 models for simulation of drought related variables, we compare their simulations of precipitation P, surface wind speed V_10 m, downwelling solar radiation RS_down, surface pressure p_surf and daily minimum and maximum temperature T_min, T_max monthly averaged for the last 65 years of the historical experiment (1950–2014). As reference for the evaluation we use the monthly and hourly ERA5 reanalysis datasets (Hersbach et al., 2023; Copernicus Climate Change Service, Climate Data Store, 2023). We find generally high agreement between reanalysis and all CMIP6 models for T_min, T_max p_surf and RS_down. They show similar temporal mean pattern correlations of 0.8 and higher. In contrast precipitation shows low agreement between different models and reanalysis while the trends in several regions are significant and directly impact SPEI. Since we found that CRU and ERA5 show very different precipitation patterns for that time period we are not rating the performance of individual models and rather incorporate the knowledge that precipitation comes with huge uncertainties in reanalysis and simulated variables in our analysis. Similar to the study by Vicente-Serrano et al. (2022), we find that uncertainty in precipitation data limits the quality of global drought assessment. The approximated ET₀ values based on CMIP6 model predictions, are considered to be reliable even though, their significance depend on the operational definition of a drought. For the 6-month SPEI we identified ET₀ as the main driver of decrease in water budget for many regions.

In agreement with the results of Balting et al. (2021) we find negative trends of water budget in the Mediterranean and several other regions. The trends are significantly stronger towards dry climate in scenarios of higher warming levels. The drought conditions that have been considered as the driest 2.3 % of each month 1950–2014 become the new normal in several arid regions at the end of this century following the SSP5-8.5. In these scenarios last two decades the projected SPEI shows extreme drought conditions with peak spatial extent of 40 % of the land surface in agricultural relevant regions. With an estimated temperature increase of 2.7 °C at the end of the century following the SSP2-4.5 pathway we calculated area fraction of maximum 20 %, and maximum 10 % in the SSP1-2.6 projection experiment for an estimated temperature increase of 2 °C. The increasing ET₀ is the dominant driver for drier conditions in arid and semi arid regions. Even with high uncertainties in precipitation projections it can be shown that the number of extreme droughts globally increases in SSP2-4.5 and SSP5-8.5, with seasonal peaks in autumn.

Relative indices such as the SPEI require a baseline period as reference. The length and the starting point of this period has an impact on the resulting Index values. Therefore, the results always need to be interpreted with respect to the reference period and are not directly comparable between analysis with different reference periods. The SPEI maps changes in the water budget to strengths of droughts only by taking their distribution in the past into account. For regions with low variance extreme conditions can be found for small water budget changes, even if crops would be capable of adapting to such small changes.

We found significant higher global increase in moderate and extreme atmospheric drying conditions on a six months timescale, in higher emission future scenarios. This also holds for most agricultural relevant IPCC WG1 reference regions.

Appendix A

Figure A1 Global pattern of temporal mean (a,b,c) and average change (d, e, f) of temperature at 2 m height based on ERA5 (a, d), CRU (b, e) and CMIP6 multi-model mean (c, f) over the period 1950–2014.

Figure A2 Global pattern of temporal mean (a, b, e, f) and average decadal change (c, d, g, h) of surface pressure p_sfc (a–d) and wind speed at 10 m (e–g) based on ERA5-Land (a, c, e, g) and CMIP6 multi-model mean (b, d, f) over the period 1950–2014.

Figure A3 Global pattern of temporal mean (a, b) and average change (c, d) of downward surface solar radiation based on ERA5 and CMIP6 multi-model mean over the period 1979–2014.

Figure A4 Seasonal cycles of average temperature at 2 m height for historical simulations of 18 CMIP6 models and two reference datasets calculated over the time period 1980–2014. The shown temperature is the average of monthly mean values of daily minimum and maximum temperature.

Figure A5 Change rate and average of soil moisture (a–f) and water budget (g–l) in different future scenarios (a, d, g, j: SSP1-2.6, b, e, h, k: SSP2-4.5, c, f, i, l: SSP5-8.5) from 2015–2100.

Figure A6 Average decadal change between 2015–2034 and 2081–2100 of SPEI in different future scenarios. Similar plot for SSP5-8.5 is shown in Fig. 8b.

Figure A7 Global overview of SPEI values for the CMIP6 multi-model ensemble for different future scenarios. The global SPEI distributions for all non glaciated land grid cells for the historical period (1950–2014) and end of century (2070–2100) are shown in panel (a) and regional distributions for the same periods in form of boxplots for IPCC AR6 WG1 reference regions in panel (c). The hexagons in panel (b) show the regional mean average decadal change of SPEI for each region only for SSP5-8.5. The regions are positioned relative to each other based on their real locations.

Figure A8 Complete time series of relative drought event area in harvest regions according to SPEI. The color indicates how much of the total non glaciated land surface area is covered by certain drought conditions from near normal.

Figure A9 Analysis analogous to Figa. 8 and 9 for 5 selected semi arid european and asian regions: WCE, EEU, MED, WCA, ECA.