Near-future CORDEX projections of the SPEI6 drought index over the Mediterranean region

5.2.8. Near-future CORDEX projections of the SPEI6 drought index over the Mediterranean region#

Production date: 25-03-2026

Produced by: CMCC foundation - Euro-Mediterranean Center on Climate Change. Albert Martinez Boti and Lorenzo Sangelantoni.

🌍 Use case: Calculating drought indices from climate projections data to inform water management#

❓ Quality assessment question#

What are the projected future changes and associated uncertainties in droughts over the Mediterranean?

The Standardized Precipitation-Evapotranspiration Index (SPEI) [1] is a widely used drought indicator that integrates precipitation and potential evapotranspiration, providing a more complete measure of water deficits than precipitation alone. Its multiscalar formulation allows the evaluation of different types of drought [2], from short-term agricultural or soil moisture deficits to long-term hydrological or ecological events. Being standardised, SPEI is dimensionless, facilitating comparisons across regions and climates. By capturing variations in frequency, duration, and severity, it effectively reflects the complex nature of drought, making it a robust and physically meaningful proxy for drought assessment in the context of climate change [3].

Understanding SPEI is particularly relevant for the water management sector. By capturing multiple types of drought SPEI serves as a versatile indicator for various users. Climate projections using SPEI enable planners to anticipate future water deficits, guide allocation, and design effective strategies for reservoir management, irrigation, and drought preparedness, helping to reduce socio-economic and ecological impacts [4][5].

This assessment uses data from a subset of models from CORDEX Regional Climate Models (RCMs) to explore the signal and the uncertainty in future projections of the SPEI6 index. The analysis focuses on the Mediterranean region, evaluating changes in SPEI6 for the near-future period (2026–2055; under the RCP8.5 scenario) relative to the 1971–2000 baseline, which serves as the reference period for standardisation. The projected uncertainty is assessed by considering the ensemble inter-model spread of projected changes.

This notebook complements the evaluation of the biases of the same subset of CORDEX models performed through the assessment named “CORDEX biases in the SPEI6 drought index over the Mediterranean region”

📢 Quality assessment statement#

These are the key outcomes of this assessment

The CORDEX ensemble median for the subset of models considered in this assessment indicates an overall increase in projected drought conditions for the near-term future (2026–2055) under the RCP8.5 scenario.
There is substantial inter-model spread, with differences in both the magnitude and, in some cases, the sign of projected SPEI6 changes, particularly at regional and seasonal scales, indicating that projections are highly model-dependent and should be interpreted in an ensemble-based framework.
Despite this spread, SPEI6 anomalies remain predominantly negative at the Mediterranean scale, indicating a consistent basin-wide tendency towards drier conditions.
The projected intensification of drought may be conservative, as indicated by a complementary historical assessment showing that the analysed CORDEX models tend to underestimate drought conditions relative to ERA5.

Visual abstract — Fig. 5.2.8.1 Spatial patterns of SPEI6 change over the Mediterranean region. Panels show (a) CORDEX ensemble median (median of the selected subset of models at each grid cell) and (b) ensemble spread (standard deviation across the selected models). SPEI6 is standardised over the reference period 1971–2000. SPEI6 change is calculated at each grid point as the difference between the monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000), and then averaged to obtain the ANNUAL mean. Future projections used for the target period follow the RCP8.5 scenario.#

📋 Methodology#

The Standardized Precipitation-Evapotranspiration Index accumulated to six months (SPEI6) has been calculated for this assessment — through the use of xclim — as a proxy to evaluate projected drought conditions for a subset of CORDEX models. The study covers the Mediterranean domain, following IPCC-AR6 regional definitions [6]. Results include spatial maps of SPEI6 changes for the near-future period (2026–2055) relative to the 1971–2000 baseline, which serves as the reference period for standardisation. Additionally, time series are computed to show the temporal evolution of future SPEI6 changes. Uncertainty is assessed through the inter-model ensemble spread, defined as the standard deviation across the selected subset of models, of projected SPEI6 changes.

Potential evapotranspiration (PET) is calculated with the Hargreaves method [7], which requires only maximum and minimum temperature (and optionally mean temperature), providing a computationally efficient yet robust estimate (e.g., [8]). The accumulated series of precipitation minus PET is then standardised over the 1971–2000 reference period.

Five key concepts are defined:

Reference period: The period used to standardise SPEI values, here 1971–2000.
Baseline period: The period against which changes are calculated. In this notebook the baseline period is selected to be the same as the reference period.
Target period: The future period for which changes are evaluated, here 2026–2055. The Representative Concentration Pathway is set to the RCP8.5 scenario for this assessment.
SPEI6 change: Differences in SPEI6 between the target period and the baseline, representing anomalies in drought conditions. Since, in this notebook, the baseline and the reference period are the same, these changes directly represent anomalies relative to the standardisation climatology used to compute SPEI6.

The analysis and results follow the next outline:

1. Parameters, requests and functions definition

1.1. Import packages
1.2. Define Parameters
1.3. Define models
1.4. Define land-sea mask and ERA5 precipitation requests
1.5. Define model requests
1.6. Functions to cache

2. Downloading and processing

2.1. Download and transform the regridding model
2.2. Download and transform models
2.3. Apply land-sea mask and change some attributes

3. Plot and describe results

3.1. Define plotting functions
3.2. Plot ensemble maps
3.3. Plot ensemble maps - Seasonal aggregations
3.4. Plot model maps
3.5. Timeseries
3.6. Results summary
3.7. Implications for the users

NOTE ON REFERENCE PERIOD SELECTION:
This notebook follows the same methodology as the assessment ‘CORDEX biases in the SPEI6 drought index over the Mediterranean region’. However, the reference period (1971–2000) is shorter than the 1971–2010 period used in the historical assessment. This choice avoids including scenario data in the standardisation (as CORDEX historical simulations extend only until 2005), ensures computational efficiency due to the shorter period, and maintains consistency with the 30-year length of the future target period.

📈 Analysis and results#

1. Parameters, requests and functions definition#

1.1. Import packages#

1.2. Define Parameters#

In the “Define Parameters” section, various customisable options for the notebook are specified:

The reference and baseline periods are set to be the same in this notebook and can be adjusted by modifying reference_slice; the default range is 1971-2000.
The target period can be adjusted by modifying target_slice; the default range is 2026-2055.
A check is performed on the target_slice and reference_slice to ensure consistency with CDS CORDEX requests and to verify that the reference period falls within the historical experiment, i.e. the last year of reference_slice does not exceed 2005.
collection_id is not customisable for this assessment and is set to ‘CORDEX’.
The area parameter specifies the geographical domain of interest. By default, the Mediterranean domain is used, following IPCC-AR6 regional definitions [6]
The time_agg allows selecting the time aggregation for the analysis. Options include: annual, DJF, MAM, JJA, SON, Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec.
SPEI accumulation window (window): Defines the accumulation period in months for the SPEI calculation. Common options include 3, 6, or 12 months, depending on whether short-term or long-term droughts are of interest. In this assessment, a six-month accumulation is used.
rcm_model_regrid and gcm_driven_model_regrid specify the EURO-CORDEX regional climate model and the driving global climate model, respectively, whose grid will be used as the target for remapping. Note that both models must be included in the matrix of RCM-GCM combinations defined by the models_cordex dictionary in Section 1.3.
The chunks selection allows the user to define if dividing into chunks when downloading the data on their local machine. Although it does not significantly affect the analysis, it is recommended to keep the default value for optimal performance.

1.3. Define models#

The following climate analyses are performed using a subset of CORDEX models that provide both the historical and RCP8.5 experiments, as well as maximum, minimum and mean temperature and precipitation, within the Climate Data Store (CDS). All Global Climate Models (GCMs) driving the selected CORDEX simulations were included, and the selection ensures a representative coverage of the available Regional Climate Models (RCMs). The selected models are the same as those used for the historical evaluation in the assessment titled “CORDEX biases in the SPEI6 drought index over the Mediterranean region”.

1.4. Define land-sea mask and ERA5 precipitation requests#

Within this notebook, ERA5 will be used to download the land-sea mask when plotting. In this section, we set the required parameters for the cds-api data-request of ERA5 land-sea mask. A request for monthly ERA5 precipitation data is also defined to compute a bare-earth (barrean) mask, which needs to be applied to ensure meaningful results in arid regions.

1.5. Define model requests#

In this section we set the required parameters for the cds-api data-request.

Show code cell source

Hide code cell source

request_cordex = {
    "format": "zip",
    "domain": "europe",
    "horizontal_resolution": "0_11_degree_x_0_11_degree",
    "temporal_resolution": "daily_mean",
    "variable": [
        "2m_air_temperature",
        "maximum_2m_temperature_in_the_last_24_hours",
        "minimum_2m_temperature_in_the_last_24_hours",
        "mean_precipitation_flux"
    ],
    "ensemble_member": "r1i1p1",
}


def get_cordex_years(
    year_start,
    year_stop,
    start_years=list(range(1951, 2097, 5)),
    end_years=list(range(1955, 2101, 5)),
):
    start_year = []
    end_year = []
    years = set(range(year_start, year_stop + 1))
    for start, end in zip(start_years, end_years):
        if years & set(range(start, end + 1)):
            start_year.append(start)
            end_year.append(end)
    return start_year, end_year



model_requests = {}
if collection_id == "CORDEX":
    for gcm_model in gcm_models:
        for model in models_cordex[gcm_model]:
                # Historical part (until 2005)
                start_years_hist = list(range(baseline_slice.start, 2006, 5))
                end_years_hist = list(range(baseline_slice.start + 4, 2010, 5))
                requests = [
                    {
                        "rcm_model": model,
                        "gcm_model": gcm_model,
                        **request_cordex,
                        "experiment": "historical",
                        "start_year": [str(start_year)],
                        "end_year": [str(end_year)],
                    }
                    for start_year, end_year in zip(
                        *get_cordex_years(
                            baseline_slice.start,
                            baseline_slice.stop,
                            start_years=start_years_hist,
                            end_years=end_years_hist,
                        )
                    )
                ]

                # Future part (RCP8.5 from 2006 onwards)
                start_years_fut = list(range(target_slice.start, target_slice.stop + 1, 5))
                end_years_fut = list(range(target_slice.start + 4, target_slice.stop + 5, 5))
                requests += [
                    {
                        "rcm_model": model,
                        "gcm_model": gcm_model,
                        **request_cordex,
                        "experiment": "rcp_8_5",
                        "start_year": [str(start_year)],
                        "end_year": [str(end_year)],
                    }
                    for start_year, end_year in zip(
                        *get_cordex_years(
                            target_slice.start,
                            target_slice.stop,
                            start_years=start_years_fut,
                            end_years=end_years_fut,
                        )
                    )
                ]
                model_requests[f"{gcm_model}_{model}"] = (
                    "projections-cordex-domains-single-levels",
                    requests,
                )
else:
    raise NotImplementedError

request_grid_out = model_requests[model_regrid]

1.6. Functions to cache#

In this section, functions that will be executed in the caching phase are defined. Caching is the process of storing copies of files in a temporary storage location, so that they can be accessed more quickly. This process also checks if the user has already downloaded a file, avoiding redundant downloads.

Description of the main functions:

compute_spei_cordex_fut: Uses the xclim package to calculate the Standardized Precipitation-Evapotranspiration Index (SPEI). The reference_slice argument specifies the reference period used for standardization, while the window argument defines the number of months over which precipitation and evapotranspiration are accumulated (6 months in this assessment).
get_mask_pr: Generates a mask identifying grid points where precipitation is below 0.3 mm/day, allowing filtering of extremely dry regions.

Show code cell source

Hide code cell source

# Ensure CF-compliant attributes for model datasets
def ensure_cf_attrs(ds):
    ds["longitude"].attrs.update({
        "standard_name": "longitude",
        "units": "degrees_east"
    })
    ds["latitude"].attrs.update({
        "standard_name": "latitude",
        "units": "degrees_north"
    })
    return ds


def add_bounds(ds):
    for coord in {"latitude", "longitude"} - set(ds.cf.bounds):
        ds = ds.cf.add_bounds(coord)
    return ds


def get_grid_out(request_grid_out, method):
    ds_regrid = download.download_and_transform(*request_grid_out)
    coords = ["latitude", "longitude"]
    if method == "conservative":
        ds_regrid = add_bounds(ds_regrid)
        for coord in list(coords):
            coords.extend(ds_regrid.cf.bounds[coord])
    grid_out = ds_regrid[coords]
    coords_to_drop = set(grid_out.coords) - set(coords) - set(grid_out.dims)
    grid_out = ds_regrid[coords].reset_coords(coords_to_drop, drop=True)
    grid_out.attrs = {}
    return grid_out


def compute_spei_cordex_fut(ds,
                            reference_slice,
                            future_slice,
                            window=6,
                            request_grid_out = None,
                            area=None,
                            **regrid_kwargs
                           ):


    ds=ds.convert_calendar('proleptic_gregorian', align_on='date', use_cftime=False)
    # Cutout
    if area:
        regionalise_kwargs = {
            "lon_slice": slice(area[1], area[3]),
            "lat_slice": slice(area[0], area[2]),
        }
        ds = utils.regionalise(ds, **regionalise_kwargs)

    # Original bounds for conservative interpolation
    if regrid_kwargs.get("method") == "conservative":
        regrid_kwargs.setdefault("ignore_degenerate", True)
        ds = add_bounds(ds)
        bounds = [
                ds.cf.get_bounds(coord).reset_coords(drop=True)
                for coord in ("latitude", "longitude")
            ]
    else:
        bounds = []

    if request_grid_out:
        ds = ds.drop_vars([
            var for var, da in ds.data_vars.items()
            if len(da.dims) != 3 or da.dtype == "object"
        ])  
        ds = ds[list(ds.data_vars)]
        vars_to_drop = [
            "rotated_latitude_longitude",
            "rotated_pole",
        ]
        for var in vars_to_drop:
            if var in ds.data_vars:
                ds = ds.drop_vars(var)
        ds = diagnostics.regrid(
                ds,
                grid_out=get_grid_out(request_grid_out, regrid_kwargs["method"]),
                **regrid_kwargs,
        )

    pr = ds["pr"] * 86400
    pr.attrs["units"] = "mm/day"


    
    da_pet=xclim.indices.potential_evapotranspiration(tas = ds["tas"],
                                                      tasmin = ds["tasmin"],
                                                      tasmax = ds["tasmax"],
                                                      method='HG85')*86400
    da_pet.attrs["units"] = "mm/day"

    wb = pr - da_pet
    wb.attrs["units"] = "mm/day"
    wb = wb.convert_calendar("proleptic_gregorian", align_on="date", use_cftime=False)    
    wb = wb.sortby("time")
    

    
    da_spei = xclim.indices.standardized_precipitation_evapotranspiration_index(wb,
                                                                                freq='MS',
                                                                                window=window,
                                                                                cal_start=f"{reference_slice.start}-01-01",
                                                                                cal_end=f"{reference_slice.stop}-12-31")
    # Wrap into dataset
    da_spei_hist = da_spei.sel(time=slice(str(reference_slice.start),str(reference_slice.stop)))
    da_spei_fut = da_spei.sel(time=slice(str(future_slice.start),str(future_slice.stop)))
    da_spei = xr.concat([da_spei_hist, da_spei_fut], dim="time")

    ds_out = xr.Dataset({"spei": da_spei})
    ds_out["spei"].attrs = da_spei.attrs

    return ds_out

    
def get_mask_pr(ds,
                threshold):
    # Select pr         
    pr = ds["tp"] * 1000   
        
    pr.attrs["units"] = "mm/day"
    pr_mean =pr.mean(dim = 'time')
        
    pr_mask = 1 - xr.where(pr_mean<threshold, 1, 0)
    ds_out = xr.Dataset({"pr_mask": pr_mask})
    return ds_out

2. Downloading and processing#

2.1. Download and transform the regridding model#

In this section, the download.download_and_transform function from the c3s_eqc_automatic_quality_control package is employed to download daily data from the selected EURO-CORDEX regridding model, select the subregion of interest, compute daily potential evapotranspiration, calculate the daily water balance, and derive the SPEI (SPEI6 for this assessment). Results are cached to avoid redundant downloads and repeated processing.

The regridding model here refers to the model whose grid is used as the target grid for remapping the other models. This ensures all datasets share a common grid, facilitating direct comparison at each grid cell. Within this notebook, the regridding model is set to "ichec_ec_earth_dmi_hirham5" but it can be changed by modifying the gcm_driven_model_regrid and rcm_model_regrid parameters in Section 1.2. Note that both models must be included in the matrix of RCM-GCM combinations defined by the models_cordex dictionary in Section 1.3. It is important to highlight that the choice of the target grid can impact the analysis depending on the specific application.

2.2. Download and transform models#

In this section, the download.download_and_transform function from the ‘c3s_eqc_automatic_quality_control’ package is used to download daily data from CORDEX models, compute daily potential evapotranspiration, calculate the daily water balance, and derive the SPEI (SPEI6 for this assessment). SPEI6 is computed on each model’s native grid and also on the "ichec_ec_earth_dmi_hirham5" grid to enable bias assessment. When regridding is required, it is performed for each essential climate variable before the index calculation to preserve standardisation. Results are cached to avoid redundant downloads and repeated processing.

model='cccma_canesm2_clmcom_clm_cclm4_8_17'
model='cccma_canesm2_gerics_remo2015'
model='cnrm_cerfacs_cm5_ipsl_wrf381p'
model='cnrm_cerfacs_cm5_smhi_rca4'
model='ichec_ec_earth_dmi_hirham5'
model='ichec_ec_earth_knmi_racmo22e'
model='ipsl_cm5a_mr_ipsl_wrf381p'
model='ipsl_cm5a_mr_knmi_racmo22e'
model='miroc_miroc5_clmcom_clm_cclm4_8_17'
model='miroc_miroc5_gerics_remo2015'
model='mohc_hadgem2_es_cnrm_aladin63'
model='mohc_hadgem2_es_mohc_hadrem3_ga7_05'
model='mpi_m_mpi_esm_lr_mpi_csc_remo2009'
model='mpi_m_mpi_esm_lr_cnrm_aladin63'
model='ncc_noresm1_m_dmi_hirham5'
model='ncc_noresm1_m_mohc_hadrem3_ga7_05'
model='ncc_noresm1_m_smhi_rca4'

2.3. Apply land-sea mask and change some attributes#

This section changes some attributes and applies a land–sea mask to all models. A bare-earth (barrean) mask is also applied to ensure meaningful results in arid regions, as in the C3S atlas. In particular, we follow the methodology available in their Github repository. While their implementation of the mask also considers factors such as snow depth and vegetation cover, here we only apply the filter based on annual mean precipitation for 1971–2000 being less than 0.3 mm day⁻¹. This simplification is justified because, in the region under consideration, the bare-earth mask obtained with the full set of filters closely matches the one derived from the precipitation threshold alone.

Note: ds_interpolated contains data from the models regridded to the "ichec_ec_earth_dmi_hirham5" grid. model_datasets contain the same data on the original grid of each model. Regridding is performed for every essential climate variable prior to the index calculation in order to avoid compromising standardisation.

100%|██████████| 1/1 [00:00<00:00, 25.66it/s]

3. Plot and describe results#

This section will display the following results:

SPEI6 change maps. The layout includes the ensemble median (defined as the median of the change values across the selected subset of models at each grid cell) and the ensemble spread (calculated as the standard deviation of the change across the selected subset of models).
SPEI6 change maps for each individual model.
Time series of SPEI6 change (averaged over the Mediterranean region).

Note: The projected SPEI6 change is calculated at each grid point as the arithmetic difference between climatologies in the target period (2026–2055) and the baseline period (1971–2000).

3.1. Define plotting functions#

The functions presented here are used to calculate SPEI6 changes and generate corresponding layout plots. Three types of layout can be displayed, depending on the plotting function:

Reference and ensemble summary: Includes the ensemble median and the ensemble spread. This is generated using plot_ensemble().
Individual models: Displays all models individually using plot_models().

Calculation of SPEI6 change:

compute_spei_change() calculates SPEI6 change at each grid point as the arithmetic difference between monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000). It then aggregates the change according to the user-specified time_agg in Section 1.2.
compute_change_4timeseries() computes SPEI6 anomalies relative to the baseline climatology for each month of the timeseries and aggregates them according to the time_agg parameter. This provides monthly, seasonal, or annual-level information depending on the user’s selection.

Show code cell source

Hide code cell source

def compute_spei_change(ds, baseline_slice, target_slice, time_agg="annual"):
    # Mapping months to names
    month_map = {
        "Jan": [1], "Feb": [2], "Mar": [3], "Apr": [4],"May": [5], "Jun": [6],
        "Jul": [7], "Aug": [8],"Sep": [9], "Oct": [10], "Nov": [11], "Dec": [12]}

    # Mapping seasons to months
    season_map = {"DJF": [12, 1, 2],"MAM": [3, 4, 5],
                  "JJA": [6, 7, 8],"SON": [9, 10, 11]}
    
    valid_opts = ["annual"] + list(season_map.keys()) + list(month_map.keys())
    assert time_agg in valid_opts, f"Invalid option: {time_agg}"
    # Select baseline and target periods
    baseline = ds.sel(time=slice(f"{baseline_slice.start}-01", f"{baseline_slice.stop}-12"))["spei"]
    target   = ds.sel(time=slice(f"{target_slice.start}-01", f"{target_slice.stop}-12"))["spei"]
    # Monthly climatologies 
    clim_baseline = baseline.groupby("time.month").mean("time", keep_attrs=True)
    clim_target   = target.groupby("time.month").mean("time", keep_attrs=True)
    #get delta
    delta_time = clim_target - clim_baseline
    delta_time.attrs = ds["spei"].attrs  
    # Aggregate according to user choice
    if time_agg == "annual":
        result = delta_time.mean("month", keep_attrs=True)
    elif time_agg in season_map:
        result = delta_time.sel(month=season_map[time_agg]).mean("month", keep_attrs=True)
    elif time_agg in month_map:
        result = delta_time.sel(month=month_map[time_agg]).squeeze("month")
        result.attrs = ds["spei"].attrs  

    return result


def compute_change_4timeseries(ds, baseline_slice, time_agg="annual"):
    # Mapping months to names
    month_map = {
        "Jan": [1], "Feb": [2], "Mar": [3], "Apr": [4],"May": [5], "Jun": [6],
        "Jul": [7], "Aug": [8],"Sep": [9], "Oct": [10], "Nov": [11], "Dec": [12]}

    # Mapping seasons to months
    season_map = {"DJF": [12, 1, 2],"MAM": [3, 4, 5],
                  "JJA": [6, 7, 8],"SON": [9, 10, 11]}
    valid_opts = ["annual"] + list(season_map.keys()) + list(month_map.keys())
    assert time_agg in valid_opts, f"Invalid option: {time_agg}"

    # Select baseline period
    baseline = ds.sel(time=slice(f"{baseline_slice.start}-01", f"{baseline_slice.stop}-12"))["spei"]

    # Monthly climatology for the baseline
    clim_baseline = baseline.groupby("time.month").mean("time", keep_attrs=True)

    # Compute monthly anomalies relative to baseline climatology
    anomalies = ds["spei"].groupby("time.month") - clim_baseline

    # Preserve attributes
    anomalies.attrs = ds["spei"].attrs

    # Aggregate anomalies
    if time_agg == "annual":
        # Group by year
        result = anomalies.groupby("time.year").mean("time", keep_attrs=True)
    elif time_agg in season_map:
        # Select months for the season, then group by year
        sel_months = season_map[time_agg]
        season_anom = anomalies.sel(time=anomalies["time.month"].isin(sel_months))
        result = season_anom.groupby("time.year").mean("time", keep_attrs=True)
    elif time_agg in month_map:
        # Select the month, keep time as is
        result = anomalies.sel(time=anomalies["time.month"].isin(month_map[time_agg]))

    return result
    
def plot_trend(time_series, var_values, label, color):
    time_numeric = pd.to_numeric(time_series.to_series())
    slope, intercept, _, p_value, _ = linregress(time_numeric, var_values)
    trend = intercept + slope * time_numeric
    slope_per_decade = slope * 10 * (60 * 60 * 24 * 365.25 * 10**9)  # nanoseconds to decade
    plt.plot(time_series, trend, linestyle='--', color=color, label=label)
    return slope_per_decade, p_value

def set_extent(da, axs, area):
    extent = [area[i] for i in (1, 3, 2, 0)]
    extent[3] += -3 
    for ax in axs:
        # Grid
        subplotspec = ax.get_subplotspec()
        gridliners = [a for a in ax.artists if isinstance(a, Gridliner)]
        for gl in gridliners:
            gl.remove()
            gl = ax.gridlines(draw_labels=True, linestyle="--", linewidth=1, color="black", alpha=0.15)
            gl.top_labels = gl.right_labels = False
            gl.left_labels = subplotspec.is_first_col()
            gl.bottom_labels = subplotspec.is_last_row()
        # Extent
        ax.set_extent(extent)

def plot_models(
    data,
    da_for_kwargs=None,
    col_wrap=4,
    subplot_kw={"projection": ccrs.PlateCarree()},
    figsize=None,
    layout="constrained",
    area=area,
    flip_cmaps=False,
    cmap_params=None,
    gcm_models=None,         
    split_model_title=True,   
    **kwargs,
):
    if isinstance(data, dict):
        assert da_for_kwargs is not None
        model_dataarrays = data
    else:
        da_for_kwargs = da_for_kwargs or data
        model_dataarrays = dict(data.groupby("model"))

    # --- Determine colour parameters ---
    default_kwargs = {"robust": True, "extend": "both"}
    if cmap_params is None:
        kwargs = default_kwargs | kwargs
        cmap_params = xr.plot.utils._determine_cmap_params(
            da_for_kwargs.values, **kwargs
        )

    fig, axs = plt.subplots(
        *(col_wrap, math.ceil(len(model_dataarrays) / col_wrap)),
        subplot_kw=subplot_kw,
        figsize=figsize,
        layout=layout,
    )
    axs = axs.flatten()

    # --- Helper function to split GCM and RCM ---
    def format_model_title(model_name):
        if gcm_models is None:
            return model_name  # fallback
        for gcm in gcm_models:
            if model_name.startswith(gcm):
                rcm = model_name[len(gcm)+1:]  # +1 to remove underscore
                return f"{gcm}\n{rcm}"
        return model_name  # fallback if no match

    # --- Plot each model ---
    for (model, da), ax in zip(model_dataarrays.items(), axs):
        pcm = plot.projected_map(
            da, ax=ax, show_stats=False, add_colorbar=False, **cmap_params
        )

        if split_model_title:
            ax.set_title(format_model_title(model))
        else:
            ax.set_title(model)

    # --- Hide empty axes ---
    for ax in axs[len(model_dataarrays):]:
        ax.set_visible(False)

    # --- Shared extent and colourbar ---
    set_extent(da_for_kwargs, axs, area)
    fig.colorbar(
        pcm,
        ax=axs.flatten(),
        extend=cmap_params.get("extend", "both"),
        location="right",
        label=f"{da_for_kwargs.attrs.get('long_name', '')} [{da_for_kwargs.attrs.get('units', '')}]",
    )

    return fig

def plot_ensemble(
    da_models,
    da_era5=None,
    subplot_kw={"projection": ccrs.PlateCarree()},
    figsize=None,
    layout="constrained",
    cbar_kwargs=None,
    area=area,
    flip_cmaps=False,
    cmap_params=None,
    cmap_params_median=None,
    cmap_params_bias=None,
    cmap_params_std=None,
    **kwargs,
):
    # Default behaviour
    default_kwargs = {"robust": True, "extend": "both"}
    if da_era5 is None and cbar_kwargs is None:
        cbar_kwargs = {"orientation": "horizontal"}

    # Create figure and axes
    fig, axs = plt.subplots(
        *(1 if da_era5 is None else 2, 2),
        subplot_kw=subplot_kw,
        figsize=figsize,
        layout=layout,
    )
    axs = axs.flatten()
    axs_iter = iter(axs)

    # --- ERA5 plot ---
    if da_era5 is not None:
        ax = next(axs_iter)
        params = cmap_params or xr.plot.utils._determine_cmap_params(
            da_era5.values, **default_kwargs
        )
        plot.projected_map(
            da_era5, ax=ax, show_stats=False, cbar_kwargs=cbar_kwargs, **params
        )
        ax.set_title("(a) ERA5")

    # --- Ensemble Median ---
    ax = next(axs_iter)
    median = da_models.median("model", keep_attrs=True)
    params = cmap_params_median or xr.plot.utils._determine_cmap_params(
        da_era5.values, **default_kwargs
    )
    plot.projected_map(
        median, ax=ax, show_stats=False, cbar_kwargs=cbar_kwargs, **params
    )
    ax.set_title("(b) Ensemble Median" if da_era5 is not None else "(a) Ensemble Median")

    # --- Ensemble Bias (Median - ERA5) ---
    if da_era5 is not None:
        ax = next(axs_iter)
        with xr.set_options(keep_attrs=True):
            bias = median - da_era5
        params = cmap_params_bias or xr.plot.utils._determine_cmap_params(
            bias.values, center=0, **default_kwargs
        )
        plot.projected_map(
            bias,
            ax=ax,
            show_stats=False,
            cbar_kwargs=cbar_kwargs,
            **params,
        )
        ax.set_title("(c) Ensemble Median Bias")

    # --- Ensemble Standard Deviation ---
    ax = next(axs_iter)
    std = da_models.std("model", keep_attrs=True)
    params = cmap_params_std or xr.plot.utils._determine_cmap_params(
        std.values, **default_kwargs
    )
    plot.projected_map(
        std,
        ax=ax,
        show_stats=False,
        cbar_kwargs=cbar_kwargs,
        **params,
    )
    ax.set_title("(d) Ensemble Standard Deviation" if da_era5 is not None else "(b) Ensemble Standard Deviation")

    # Set extent
    set_extent(da_models, axs, area)

    return fig

3.2. Plot ensemble maps#

In this section, we invoke the plot_ensemble() function to visualise the SPEI6 change for (a) the ensemble median (defined as the median of the change values for the selected subset of models at each grid cell) and (b) the ensemble spread (calculated as the standard deviation of the change across the selected subset of models).

ANNUAL aggregation

../../_images/3ab2a42d58018ac27b24bafcec384ebcd037bcff3fda7d4a5fb78ae08541977d.png

Fig 1. Spatial patterns of SPEI6 change over the Mediterranean region. Panels show (a) CORDEX ensemble median (median of the selected subset of models at each grid cell) and (b) ensemble spread (standard deviation across the selected models). SPEI6 is standardised over the reference period 1971–2000. SPEI6 change is calculated at each grid point as the difference between the monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000), and then averaged to obtain the ANNUAL mean. Future projections used for the target period follow the RCP8.5 scenario.

3.3. Plot ensemble maps - Seasonal aggregations#

Same as 3.2 section but for: DJF, MAM, JJA and SON

WINTER (DJF)

../../_images/40a594d50d035c502899374cfbee4070edf52415074067773804834f91f92df7.png

Fig 2. Spatial patterns of SPEI6 change over the Mediterranean region. Panels show (a) CORDEX ensemble median (median of the selected subset of models at each grid cell) and (b) ensemble spread (standard deviation across the selected models). SPEI6 is standardised over the reference period 1971–2000. SPEI6 change is calculated at each grid point as the difference between the monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000), and then averaged across the winter to obtain the DJF mean. Future projections used for the target period follow the RCP8.5 scenario.

SPRING (MAM)

../../_images/b0e9c7953932782ae58849bbf8cb897f2e248d09b8ad215d961492005959639e.png

Fig 3. Spatial patterns of SPEI6 change over the Mediterranean region. Panels show (a) CORDEX ensemble median (median of the selected subset of models at each grid cell) and (b) ensemble spread (standard deviation across the selected models). SPEI6 is standardised over the reference period 1971–2000. SPEI6 change is calculated at each grid point as the difference between the monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000), and then averaged across the spring to obtain the MAM mean. Future projections used for the target period follow the RCP8.5 scenario.

SUMMER (JJA)

../../_images/87c0867f8d78ec9d57b00caf967adad7bed569e303d6b8d64d352fb05590d829.png

Fig 4. Spatial patterns of SPEI6 change over the Mediterranean region. Panels show (a) CORDEX ensemble median (median of the selected subset of models at each grid cell) and (b) ensemble spread (standard deviation across the selected models). SPEI6 is standardised over the reference period 1971–2000. SPEI6 change is calculated at each grid point as the difference between the monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000), and then averaged across the summer to obtain the JJA mean. Future projections used for the target period follow the RCP8.5 scenario.

AUTUMN (SON)

../../_images/23772292bbb4c456da3d526e004d2539dba12d2a0e8cb632e64a64c1133b5493.png

Fig 5. Spatial patterns of SPEI6 change over the Mediterranean region. Panels show (a) CORDEX ensemble median (median of the selected subset of models at each grid cell) and (b) ensemble spread (standard deviation across the selected models). SPEI6 is standardised over the reference period 1971–2000. SPEI6 change is calculated at each grid point as the difference between the monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000), and then averaged across the autumn to obtain the SON mean. Future projections used for the target period follow the RCP8.5 scenario.

3.4. Plot model maps#

In this section, we invoke the plot_models() function to visualise the SPEI6 change of every model individually. Note that the model data used in this section maintains its original grid. Only the annual mean is shown here.

../../_images/f0a7c19a791928a89aa245c265ffa0bd4388131d30875236a894c9bbccfcae49.png

Fig 6. Spatial patterns of SPEI6 change for each individual CORDEX model over the Mediterranean region. SPEI6 is standardised over the reference period 1971–2000. For each model, SPEI6 change is calculated at each grid point as the difference between the monthly climatologies of the target period (2026–2055) and the baseline period (1971–2000), and then averaged to obtain the annual mean. Future projections used for the target period follow the RCP8.5 scenario.

3.5. Timeseries#

This section shows the timeseries of SPEI6 anomalies relative to the reference period climatology (same period as the baseline for this notebook). We first use compute_change_4timeseries(), which calculates SPEI6 anomalies for each month of the target period (relative to the baseline monthly climatologies) and aggregates them according to the time_agg parameter. Spatially averaged values are then obtained. The analysis compares the CORDEX ensemble median and individual models and shows the ensemble spread and the slope value of the trend for the CORDEX ensemble median. Only the annual aggregation is shown here.

Show code cell source

Hide code cell source

# Define colors
colors = {"CORDEX": "green"}

# Cutout region
regionalise_kwargs = {
    "lon_slice": slice(area[1], area[3]),
    "lat_slice": slice(area[0], area[2]),
}

# Compute change/anomaly for each model and spatially aggregate
cordex_means = []
for model, ds in model_datasets.items():
    ds_model_change = compute_change_4timeseries(ds, baseline_slice, time_agg)
    ds_model_change["year"] = pd.to_datetime(ds_model_change["year"].values, format='%Y')
    
    ds_model_regional = utils.regionalise(ds_model_change, **regionalise_kwargs)
    cordex_mean = diagnostics.spatial_weighted_mean(ds_model_regional)
    cordex_means.append(cordex_mean)

# Ensemble statistics
cordex_stack = xr.concat(cordex_means, dim="model")
cordex_var = cordex_stack.median("model", keep_attrs=True)
cordex_std = cordex_stack.std("model", keep_attrs=True)

# --- Select ONLY target period ---
cordex_var = cordex_var.sel(year=slice(str(target_slice.start), str(target_slice.stop)))
cordex_std = cordex_std.sel(year=slice(str(target_slice.start), str(target_slice.stop)))

cordex_means = [
    model_data.sel(year=slice(str(target_slice.start), str(target_slice.stop)))
    for model_data in cordex_means
]

# Initialize plot
plt.figure(figsize=(10, 5))
ax = plt.gca()

# Plot individual model time series
for model_data in cordex_means:
    model_values = model_data.values.flatten() if model_data.ndim > 1 else model_data.values
    plt.plot(model_data['year'], model_values, linestyle='-', color="grey", alpha=0.4)

# Add legend entry for individual models
plt.plot([], [], linestyle='-', color='grey', alpha=0.4, label='Individual Models')

# Plot CORDEX median
cordex_var.plot(label="CORDEX Median", color=colors["CORDEX"])

# Fill ±1 std (spread)
plt.fill_between(
    cordex_var['year'],
    cordex_var.values - cordex_std.values,
    cordex_var.values + cordex_std.values,
    color=colors["CORDEX"],
    alpha=0.4,
    label="CORDEX Spread"
)

# Compute trend ONLY for CORDEX (already filtered to target period)
slope_cordex, p_value_cordex = plot_trend(
    cordex_var['year'][~np.isnan(cordex_var.values)],
    cordex_var.values[~np.isnan(cordex_var.values)],
    "CORDEX Median Trend",
    colors["CORDEX"]
)

trend_info = [
    ("CORDEX Median", slope_cordex, p_value_cordex)
]

# Labels and title
plt.xlabel('Year')
plt.ylabel("Fraction of unity")

title = (
    f"SPEI{window} (ref. period {reference_slice.start}-{reference_slice.stop}) | "
    f"Change rel. to {baseline_slice.start}-{baseline_slice.stop}\n"
    f"Target period: ({target_slice.start}-{target_slice.stop}) | {time_agg}"
)
plt.title(title)

plt.legend()
plt.grid(True)

# Format x-axis as years
ax.xaxis.set_major_formatter(mdates.DateFormatter('%Y'))
plt.gcf().autofmt_xdate()

# Annotate trends
plt.gcf().text(
    0.5, -0.015,
    '\n'.join([
        f"{name} Trend: {trend:.4f} (fraction of unity/decade), p-value: {p:.3e}"
        for name, trend, p in trend_info
    ]),
    fontsize=12,
    ha='center',
    va='center',
    bbox=dict(boxstyle="square,pad=0.5", edgecolor="black", facecolor="lightgray")
)

plt.ylim((-2, 2))
plt.show()

../../_images/fa76579674b229a2d902abc01d897524e154a7effaa3106003127736d456380b.png

Fig 7. Timeseries of SPEI6 anomalies over the Mediterranean region. Monthly anomalies are calculated relative to the monthly baseline climatologies, aggregated according to the specified timescale, and then spatially averaged over the region. The figure compares the CORDEX ensemble median and individual models, showing the trend for the ensemble median, as well as the ensemble spread and the corresponding trend slope. Future projections used for the target period follow the RCP8.5 scenario.

3.6. Results summary#

The projected SPEI6 change for the near-term future (2026–2055) indicates an overall decrease in the index (Fig. 1), suggesting a general intensification of drought conditions. However, a large inter-model spread is evident, with substantial disparities not only in the magnitude of change but also in its sign (see Fig. 6). Some models project wetter conditions than the baseline, particularly over regions such as the Balkans, southern France, the Adriatic coast of Italy, and parts of Turkey.
The inter-model spread spatial patterns are sensitive to the seasonal aggregation considered (Figs. 2–5), with the largest inter-model discrepancies occurring in spring.
The timeseries of spatial means over the Mediterranean shows an increase in projected drought conditions for the near-term future period (Fig. 7). For the CORDEX ensemble median, the timeseries shows a slight (but statistically significant) negative trend over this period. Although the trend is weak and the inter-model spread is substantial, most individual models—and the ensemble median—remain consistently below zero. This indicates that drier conditions are projected for the annual aggregation over the region, as anomalies are defined relative to the historical baseline period.
This assessment is broadly consistent with results from the interactive C3S Atlas. The Atlas indicates slightly stronger decreases in SPEI6, which may be partly explained by differences in the standardisation periods used and the larger ensemble of models considered.

3.7. Implications for the users#

The substantial uncertainty in projected drought conditions, particularly at regional and seasonal scales, supports the use of risk-based approaches in water management and drought preparedness, where strategies should be tested against a range of plausible futures rather than relying on a single model or ensemble summary.
The historical bias assessment (CORDEX biases in the SPEI6 drought index over the Mediterranean region) shows that the analysed CORDEX models tend to underestimate drought conditions relative to ERA5. Based on this, the projected intensification of drought in this notebook may be conservative and should be interpreted with this limitation in mind.
Complementary tools such as the interactive C3S Atlas, which include a larger ensemble and additional dimensions such as months, global warming levels, and multiple time horizons, can provide further context for interpreting these results.
Beyond the direct interpretation of results, this notebook emphasises the basis and methodological construction of the SPEI index, supporting its use in operational contexts.

ℹ️ If you want to know more#

Key resources#

Some key resources and further reading were linked throughout this assessment.

The CDS catalogue entries for the data used were:

CORDEX regional climate model data on single levels (daily - Maximum 2m temperature in the last 24 hours, Minimum 2m temperature in the last 24 hours, 2m air temperature and Mean precipitation flux): https://cds.climate.copernicus.eu/datasets/projections-cordex-domains-single-levels?tab=overview
ERA5 monthly averaged data on single levels from 1940 to present (total precipitation): https://cds.climate.copernicus.eu/datasets/reanalysis-era5-single-levels-monthly-means?tab=overview

Code libraries used:

C3S EQC custom functions, c3s_eqc_automatic_quality_control, prepared by B-Open

References#

[1] Vicente-Serrano, S. M., Beguerı́a, S., and López-Moreno, J. I., 2010. A multiscalar drought index sensitive to global warming: The standardized precipitation evapotranspiration index, Journal of Climate, 23, 1696–1718. https://doi.org/10.1175/2009JCLI2909.1

[2] Wilhite, D. A., 2000. Droughts as a natural hazard: concepts and definitions. In: DROUGHT, A Global Assessment, vol. I and II, Routledge Hazards and Disasters Series, Routledge.

[3] Vicente-Serrano, S. M., Beguerı́a, S., Lorenzo-Lacruz, J., Camarero, J. J., López-Moreno, J. I., Azorin-Molina, C., Revuelto, J., Morán-Tejeda, E., and Sanchez-Lorenzo, A., 2012. Performance of drought indices for ecological, agricultural, and hydrological applications, Earth Interactions, 16, https://doi.org/10.1175/2012EI000434.1

[4] Kalisoras, A., Georgoulias, A. K., Akritidis, D., and Zanis, P., 2025. Future Projections in Agricultural Drought Characteristics for Greece Under Different Climate Change Scenarios. Environmental and Earth Sciences Proceedings, 35(1), 29. https://doi.org/10.3390/eesp2025035029

[5] Spinoni, J., Barbosa, P., Bucchignani, E., Cassano, J., Cavazos, T., Cescatti, A., Christensen, J. H., Christensen, O. B., Coppola, E., Evans, J. P., Forzieri, G., Geyer, B., Giorgi, F., Jacob, D., Katzfey, J., Koenigk, T., Laprise, R., Lennard, C. J., Kurnaz, M. L., … Dosio, A., 2021. Global exposure of population and land-use to meteorological droughts under different warming levels and SSPs: A CORDEX-based study. International Journal of Climatology, 41(15), 6825–6853. https://doi.org/10.1002/joc.7302

[6] Iturbide, M., Gutiérrez, J. M., Alves, L. M., Bedia, J., Cerezo-Mota, R., Cimadevilla, E., Cofiño, A. S., Di Luca, A., Faria, S. H., Gorodetskaya, I. V., Hauser, M., Herrera, S., Hennessy, K., Hewitt, H. T., Jones, R. G., Krakovska, S., Manzanas, R., Martínez-Castro, D., Narisma, G. T., Nurhati, I. S., Pinto, I., Seneviratne, S. I., van den Hurk, B., and Vera, C. S., 2020. An update of IPCC climate reference regions for subcontinental analysis of climate model data: definition and aggregated datasets, Earth Syst. Sci. Data, 12, 2959–2970, https://doi.org/10.5194/essd-12-2959-2020

[7] Hargreaves, G. H., and Samani, Z. A., 1985. Reference crop evapotranspiration from temperature, Applied engineering in agriculture, pp. 96–99. https://doi.org/10.13031/2013.26773

[8] Beguerı́a, S., Vicente-Serrano, S. M., Reig, F., and Latorre, B., 2014. Standardized precipitation evapotranspiration index (spei) revisited: Parameter fitting, evapotranspiration models, tools, datasets and drought monitoring, International Journal of Climatology, 34, 3001–3023. https://doi.org/10.1002/joc.3887

Near-future CORDEX projections of the SPEI6 drought index over the Mediterranean region

Contents

5.2.8. Near-future CORDEX projections of the SPEI6 drought index over the Mediterranean region#

🌍 Use case: Calculating drought indices from climate projections data to inform water management#

❓ Quality assessment question#

📢 Quality assessment statement#

📋 Methodology#

📈 Analysis and results#

1. Parameters, requests and functions definition#

1.1. Import packages#

1.2. Define Parameters#

1.3. Define models#

1.4. Define land-sea mask and ERA5 precipitation requests#

1.5. Define model requests#

1.6. Functions to cache#

2. Downloading and processing#

2.1. Download and transform the regridding model#

2.2. Download and transform models#

2.3. Apply land-sea mask and change some attributes#

3. Plot and describe results#

3.1. Define plotting functions#

3.2. Plot ensemble maps#

3.3. Plot ensemble maps - Seasonal aggregations#

3.4. Plot model maps#

3.5. Timeseries#

3.6. Results summary#

3.7. Implications for the users#

ℹ️ If you want to know more#

Key resources#

References#