The added value of Med-CORDEX coupled high-resolution regional climate models in representing sea surface temperature and marine heatwaves in the Mediterranean Sea

Marine heatwaves (MHWs) pose significant threats to Mediterranean marine ecosystems and coastal economies, with frequency and severity projected to increase under future climate change. While Med-CORDEX fully-coupled Regional Climate System Models (RCSMs) offer enhanced resolution and improved representation of local processes relative to their parent Global Climate Models (GCMs), a systematic assessment of their added value for sea surface temperature (SST) and MHW properties has been lacking. This study quantifies the added value of Med-CORDEX RCSMs over the Mediterranean basin, evaluating their capacity to correct GCM biases and improve the representation of SST and MHW probability distributions. Results show that added value is scale-dependent and metric-specific. RCSMs generally improve the SST spatial pattern and the shape and upper tail of its temporal distribution, but mostly fail to correct Mediterranean basin-averaged errors in the mean, standard deviation, 90th percentile and linear trend. SST improvements are most consistent along coastlines and in semi-enclosed areas, where fine-scale ocean-atmosphere interactions and topographic constraints are best resolved by RCSMs. For MHW duration, downscaling provides consistent and spatially widespread improvements across nearly all models, driven by a better representation of short-lived events. For MHW intensity, added value is model-dependent and not systematic: while the majority of RCSMs improve this metric, some models exhibit deterioration linked to model-specific features. These results demonstrate that higher horizontal resolution is a necessary but not sufficient condition for improved MHW representation, and that simultaneous advances in other model components are required to fully exploit the potential of regional downscaling. In particular, within the Med-CORDEX ensemble, inter-model differences in MHW added value appear to be correlated with improvements in the thickness of the first layer.

1 Introduction

The Mediterranean region has long been recognized as one of the most prominent and vulnerable climate change hotspots of the 21st century (Giorgi, 2006; Lazoglou et al., 2024). This semi-enclosed basin is experiencing temperature increases at rates exceeding the global average, with surface temperatures approximately 1.5 °C above late 19th-century levels (IPCC, 2023). In the 21st century, Mediterranean temperatures are projected to increase 20 % faster than the global average, particularly in summer (Lionello and Scarascia, 2018), intensifying pressures across key sectors as water resources, ecosystems, food production, human health, and regional security (Cramer et al., 2018; Smith et al., 2021, 2025). Marine heatwaves (MHWs, Hobday et al., 2016), i.e., prolonged periods of unusually high sea surface temperatures (SSTs), have emerged as critical manifestations of Mediterranean climate change as their frequency strongly increased from the 1980s (Darmaraki et al., 2019a; Dayan et al., 2023; Hamdeno and Alvera-Azcaráte, 2023). Such events are projected to increase their occurrence, severity and areal extension compared to today's events (Darmaraki et al., 2019b), with geographically varying impacts on marine ecosystems and coastal communities (Garrabou et al., 2022).

Reliable climate projections are critical to inform society about the implications of climate extremes. However, Global Climate Models (GCMs) from the Coupled Model Inter-comparison Project (CMIP) initiative (Taylor et al., 2012; Eyring et al., 2016) typically operate at coarse resolutions (1° ocean) that cannot adequately resolve the complex topography, coastlines, strait dynamics, regional wind systems, and fine-scale ocean-atmosphere feedbacks characteristic of the Mediterranean basin. Capturing these local features and processes is essential for accurately representing present and future SST and MHW characteristics in the region. The Mediterranean Coordinated Regional Downscaling Experiment (Med-CORDEX, https://med-cordex.github.io/, last access: 31 August 2026) represents the largest coordinated multi-model effort to address these limitations through fully coupled, high-resolution Regional Climate System Models (RCSMs) operating at typically 0.1° in the ocean and 0.1–0.5° in the atmosphere, with structural enhancements specifically designed to improve process representation in the Mediterranean (Ruti et al., 2016; Somot et al., 2018). While these improvements offer the potential for a better performance, it is crucial to quantitatively assess whether the increased computational complexity translates into actual skill improvement over the driving models. Such improvement is quantified as added value. Added value is not uniform across all variables, regions, or scales. Analysis of EURO-CORDEX (Jacob et al., 2014, 2020) and Med-CORDEX downscaling experiments revealed that large-scale and time-averaged fields over smooth terrain may show minimal improvement (Rummukainen, 2016), while regions with complex orography and coastal contrasts demonstrate clear added value in process representation (Byun et al., 2023), precipitation (Ciarlo` et al., 2021; Soares and Cardoso, 2018; Careto et al., 2022), wind speed (Herrmann et al., 2011; Molina et al., 2023, 2024), air temperature (Torma et al., 2015; Cardoso and Soares, 2022) and climate extremes (Fantini et al., 2018; Iles et al., 2020). Despite extensive research on Mediterranean climate change (Darmaraki et al., 2019b; Soto-Navarro et al., 2020; Parras-Berrocal et al., 2024, 2025) and the increasing use of Med-CORDEX simulations, a systematic assessment of their added value for SST and MHW properties remains lacking. Addressing this gap is particularly important given the Mediterranean's status as a MHW hotspot (Denaxa et al., 2025) and the crucial role of SST for such extremes. Beyond the Med-CORDEX activity, Mishra et al. (2023) already showed the benefits of enhancing horizontal resolution for oceanic variables in HighResMIP simulations in the Mediterranean region. At the global scale, Pilo et al. (2019) demonstrated the added value of higher resolution simulations, which exhibited weaker biases in MHW properties compared to their parent models. Furthermore, the gain of higher horizontal resolution is mainly focused on the duration of MHWs (Guo et al., 2019; Capotondi et al., 2024), but biases in intensity and frequency still persist (Pilo et al., 2019). By resolving fine-scale processes that GCMs cannot capture, Med-CORDEX simulations have a high potential to improve the representation of SSTs and MHWs in the Mediterranean Sea.

The objective of this study is to quantify the added value of Med-CORDEX models for Mediterranean SST and derived MHW properties (duration and intensity). To our knowledge, this represents the first systematic multi-model added value assessment based on coordinated GCM-RCSMs pairs for oceanic variables, not only in the Mediterranean but across all ocean regions worldwide. We examined the capacity of RCSMs to (i) correct GCM biases and (ii) improve the spatiotemporal representation of the probability distribution function (PDF) of these variables over a 39-year period. This assessment evaluates the effectiveness of downscaling for marine extremes and informs future modeling priorities. The article is structured as follows: Sect. 2 details the datasets and methodology of the analysis; Sect. 3 presents and discusses the main results; Sect. 4 summarizes the main findings and limitations, concluding with key recommendations.

2 Data and Methods

2.1 Data

This study utilizes daily SST fields from both RCSM simulations and the corresponding GCM simulations. These SST outputs, along with their corresponding MHW properties, are compared to daily SST observations and derived MHW characteristics, sourced from the High Resolution L4 Sea Surface Temperature Reprocessed for the Mediterranean Sea, hereafter defined as OBS (see Table S1 in the Supplement for details, https://doi.org/10.48670/moi-00173, Pisano et al., 2016; Embury et al., 2024). OBS SST are available from 1982 and represent the foundation temperature at 1/16°, i.e., the surface ocean temperature free of the diurnal cycle. In contrast, models provide a daily SST average which incorporates the diurnal cycle. However, as shown in Sect. 2.3, this methodological difference is unlikely to be the primary driver of the model-observation discrepancies.

This study analyzes an ensemble of seven fully coupled Med-CORDEX RCSMs driven by CMIP5 and CMIP6 GCMs. The specific configurations (abbreviated as CMCC, AWI, UNIBEL, LMD-CNRM, LMD-IPSL, LMD-MPI, CNRM-RCSM4, CNRM-RCSM6, and GUF) feature high oceanic (6–30 km) and atmospheric (12–50 km) resolutions, enhanced vertical resolution, high-frequency atmosphere-ocean coupling, and dedicated parameterizations of vertical ocean dynamics and air-sea fluxes optimized for the Mediterranean Sea (Ruti et al., 2016). All RCSMs cover the entire Mediterranean basin and a small portion of the adjacent North Atlantic. Full model details are provided in Table 1.

Cavicchia et al. (2015)Conte et al. (2020)Sevault et al. (2014)Sein et al. (2015)Parras-Berrocal et al. (2020)Djurdjevic and Rajkovic (2008)L’Hévéder et al. (2013)Hamouda et al. (2023)Sevault (2024)Scoccimarro et al. (2011)Voldoire et al. (2013)Giorgetta et al. (2013)Giorgetta et al. (2013)Voldoire et al. (2013)Dufresne et al. (2013)Giorgetta et al. (2013)Döscher et al. (2022)Séférian et al. (2019)

Table 1 Characteristics of the Med-CORDEX RCSMs utilized in this study.

Simulations are evaluated over the 1982–2020 period by concatenating historical outputs with early scenario realizations: 2006–2020 from the Representative Concentration Pathway 8.5 (RCP8.5) for CMIP5 GCMs and associated RCSMs; 2015–2020 from the Shared Socio-economic Pathway 5-8.5 (SSP5-8.5) for CMIP6 GCMs and associated RCSMs. We select this period to maximize overlap with high-resolution OBS, improving the statistical robustness of the results while minimizing the use of data from future scenario simulations. This approach is standard practice (Ciarlo` et al., 2021; Sevault et al., 2014; Scoccimarro et al., 2024; Simon et al., 2025), as radiative forcing divergence in early scenario years is minimal while internal variability and model characteristics dominate the signal (Hawkins and Sutton, 2009; Lehner et al., 2020).

2.2 Methods

2.2.1 Assessment of SST and MHW statistics

The capacity of RCSMs to correct biases present in their parent GCMs was assessed by comparing a set of SST statistics and MHW properties across OBS, RCSMs, and GCMs. The SST statistics considered are the mean, standard deviation (STD), 90th percentile, linear trend, and STD of SST residuals. SST residuals were defined as anomalies obtained by removing both the long-term linear trend and the climatological seasonal cycle from each grid point. The considered MHW properties are duration, mean intensity, and number of events. All metrics were computed at each grid point over the full 1982–2020 reference period and then spatially averaged over the Mediterranean domain. GCM biases were defined as the difference between GCM and OBS values, while RCSM increments were calculated as the difference between each RCSM and its driving GCM. Spatial maps of error improvement were constructed by identifying grid points where the absolute RCSM error is smaller than the absolute GCM error for a given SST statistic, providing a spatial assessment of where the RCSMs reduce or amplify GCM errors.

2.2.2 Added Value index

We defined added value as the ability of an RCSM to reproduce observed PDFs more accurately than its driving GCM. Since regional models aim to capture local dynamics smoothed out in global simulations, improving the representation of the observed PDF structure at each simulated location is critical. To quantify this improvement, we employ the Added Value index (AVi) developed by Ciarlo` et al. (2021), computed as:

$$\begin{array}{}\text{(1)}& \mathrm{AVi}=D_{\mathrm{GCM}}-D_{\mathrm{RCSM}}\end{array}$$

where D_GCM and D_RCSM are the so-called “Relative Probability Differences” for GCMs and RCSMs, calculated as

$$\begin{array}{}\text{(2)}& D={\displaystyle \frac{{\sum }_{i=\mathrm{1}}^n\left|(N_{\mathrm{m}}^i-N_{\mathrm{o}}^i)\mathrm{\Delta }v\right|}{{\sum }_{i=\mathrm{1}}^n{N}_{\mathrm{o}}^i\mathrm{\Delta }v}}\end{array}$$

where $N_{\mathrm{m}}^i$ and $N_{\mathrm{o}}^i$ denote the number of data points in the model and observations, respectively, falling into bin i, and Δv is the bin width of the variable. Positive AVi values indicate that the RCSM reproduces the observed PDF more accurately than the driving GCM, while negative values indicate that the GCM provides a closer match to observations. As the definition of Δv influences the calculation of the PDF, and thus the AVi, this parameter is defined through a standard rule. The optimal Δv for each distribution is estimated through the Freedman-Diaconis rule (Freedman and Diaconis, 1981), applied to OBS, as:

$$\begin{array}{}\text{(3)}& \mathrm{\Delta }v={\displaystyle \frac{\mathrm{2}\times \mathrm{IQR}}{n^{\mathrm{1}/\mathrm{3}}}}\end{array}$$

where n is the number of data points considered and IQR is the inter-quartile range. Following Ciarlo` et al. (2021), we apply a correction for bins that are empty in the GCM but populated in OBS and RCSM. To prevent artificial negative AVi scores in these cases, those bins contribute zero to D_RCSM, ensuring the RCSM adds value for capturing events missed by the GCM. If the RCSM fails to populate a bin where both the GCM and OBS show events, the same rule is applied symmetrically. Sensitivity tests confirmed that AVi values are not substantially affected by this correction; in particular, area-averaged AVi rankings and signs across models are preserved for all metrics. At each grid point, the statistical significance of AVi was assessed via bootstrap randomization. A reference distribution of 1000 synthetic AVi values was generated by randomly drawing two new GCM and RCSM samples from their union, simulating the null hypothesis of no systematic difference between global and regional models. The original AVi was defined significant when it fell outside the 5th–95th percentile range of this distribution. The area-averaged AVi was then computed over the Mediterranean domain, with non-significant grid points set to zero.

For each model family, OBS and GCM were bilinearly regridded on the correspondent RCSM grid to facilitate a fair comparison between regional and global models. Interpolating GCM values to a finer resolution may produce unrealistic outputs, while interpolating OBS and RCSM values to a coarser resolution can compromise the representation of the PDF (Prein et al., 2016). However, to effectively illustrate the added value of downscaling, it is crucial to demonstrate that resolving local-scale processes provides more information than simply interpolating a coarse, non-resolving grid onto a finer grid. AVi is computed independently at each point of the RCSM native grid and displayed as spatial maps. Multi-model AVi mean maps are obtained by interpolating individual AVi fields by nearest neighbor onto the NEMOMED8 grid (native to the CNRM-RCSM4 and LMD models) prior to averaging. Grid points where fewer than 6 out of 9 simulations (66 %) agree in sign with the multi-model mean are masked. Global PDFs of the selected parameters are shown, aggregating all data from AVi-significant grid points across OBS, GCMs and RCSMs, with each PDF normalized by the total number of data points. A corresponding AVi value is derived from these aggregated PDFs, providing a measure of added value that integrates the performance of all RCSMs and GCMs across the significant domain.

We analyzed the added value of four parameters: median-removed SST anomalies, upper tail SST anomalies, MHW duration and intensity. For each grid point, median-removed SST anomalies are obtained by subtracting the median from the SST distribution, focusing the analysis on its shape and excluding the influence of background model biases. Upper-tail SST anomalies are computed by subtracting the 90th percentile from all values above this percentile, isolating very high SST values relative to the local climatology, without any persistence requirement. MHWs are defined based on the criteria established by Hobday et al. (2016). An MHW event occurs when SST exceeds the 90th percentile climatology calculated over a 30-year daily climatology (1982–2011) for at least 5 consecutive days, with a maximum gap of 2 d. We focus on two MHW characteristics: duration, defined as the number of consecutive days of each event, and mean intensity, defined as the average SST anomaly above the 90th-percentile threshold computed over the reference climatology. It should be noted that defining MHWs relative to each model's own climatology effectively removes mean-state biases from all data sources, acting as an implicit bias correction. Combined with the median-removed and extreme SST, our analysis therefore captures only the variability component of SST, excluding potential added value associated with improvements in the mean state. This has implications for the interpretation of results and is discussed further.

2.3 Sensitivity to the observational reference

The OBS SST product provides the longest high-resolution SST record available for the Mediterranean Sea, but it is based on nighttime retrievals, whereas daily model outputs represent averages over the full diurnal cycle. To assess whether this difference in temporal sampling affects our results, we compared the satellite-derived nighttime L4 SST used in this study against two other independent products: the High Resolution Diurnal Subskin Sea Surface Temperature Analysis (https://doi.org/10.48670/moi-00170, Marullo et al., 2014; Buongiorno Nardelli et al., 2013) and the CMEMS Mediterranean Sea Physics Reanalysis (https://doi.org/10.48670/mds-00375, Escudier et al. (2021)), which provides hourly reanalyzed SST fields. Details of these products can be found in Table S1. The inter-product differences in SST (Fig. S1a, c in the Supplement) and extreme SST (Fig. S1b,d) PDFs over 2019–2022 are substantially smaller than the RCSM-GCM differences across the ensemble, demonstrating that the diagnosed added value is robust with respect to OBS.

Figure 1 RCSM increment versus GCM bias for area-averaged SST (a) mean, (b) STD, (c) 90th percentile, (d) trend, (e) STD of SST residuals over 1982–2020. The grey line indicates perfect error correction. Blue, red, and orange shading denote skillful error reduction, deterioration, and overcompensation (RCSM error opposite in sign but larger in magnitude than the GCM error), respectively. For the definition of the Mediterranean domain see Fig. 2.

3 Results

3.1 SST

3.1.1 Evaluation of SST statistics

Figure 1 assesses whether RCSMs tend to correct area-averaged SST errors present in GCMs (blue areas) or amplify existing ones (red and orange areas). For the definition of the Mediterranean domain, see Fig. 2. RCSMs struggle to correct the area-averaged statistics. Regarding the mean state, only a subset of 4 models effectively reduces the biases present in their driving GCMs (Fig. 1a). Conversely, the majority of models amplify the existing cold bias in the driving GCMs. Spatial maps of GCM and RCSM biases (Fig. S2) confirm that this tendency is broadly distributed across the basin. However, 7 out of 9 RCSMs improve the spatial structure of the mean SST field relative to their driving GCMs, as shown by higher pattern correlations with observations (Fig. S2). Area-averaged bias improvement is connected to the capacity of RCSMs to reduce the bias across large portions of the basin, with the best-performing models improving more than 69 % of grid points (Fig. S3). Yet, all RCSMs show skillful bias reduction in key regions, particularly along the coastlines and in semi-enclosed areas such as the Adriatic Sea, the Alboran Sea, and the Aegean Sea, where the added value of downscaling is most evident (Fig. S3). Seasonally, winter is the most favorable period for bias reduction (Fig. S4d), whereas other seasons show limited improvement or further deterioration (Fig. S4a, b, c). This suggests that while RCSMs better resolve winter mixing and convection processes, they may not adequately capture the basin-wide heat fluxes or advective inputs that dominate other seasons.

Similarly, the standard deviation (STD), here used as a measure of total variability, is rarely improved (Fig. 1b). Seven regional simulations feature a lower STD than their parent GCMs, worsening the signal in most cases. As with the mean bias, RCSMs tend to degrade the STD over large portions of the basin, while showing localized improvement along the coastlines and in semi-enclosed areas (Fig. S5). For most models, autumn and winter are the seasons in which RCSMs most effectively reduce GCM STD errors (Fig. S4g, h), while in spring and summer they more commonly amplify the existing errors (Fig. S4e, f).

The representation of the 90th percentile is generally degraded (Fig. 1c), with most RCSMs improving fewer than 50 % of grid points (Fig. S6). Most models underestimate the observed 90th percentile more severely than their driving GCMs, consistent with the cold mean bias and reduced variability described above. In summer and autumn, only two models improved the representation of the 90th percentile (Fig. S4j, k), while the majority of RCSMs show correcting skills in winter and spring (Fig. S4i, l).

Most RCSMs underestimate warming relative to their driving GCMs, except LMD-MPI (Fig. 1d). Unlike the mean bias and STD, trend errors show no systematic spatial structure, with RCSMs failing to outperform their driving GCMs even in coastal and semi-enclosed areas (Fig. S7), suggesting a large-scale origin for this deterioration. Issues in representing temperature trends are not unique to Med-CORDEX models. EURO-CORDEX studies have documented similar failures in reproducing long-term surface air temperature trends, largely attributed to misrepresented aerosol evolution and resulting biases in surface solar radiation (Nabat et al., 2014; Bartók et al., 2017; Boé et al., 2020; Schumacher et al., 2024).

The variability of SST residuals is underestimated in CMIP5 GCMs, and slightly overestimated in CMIP6 GCMs (Fig. 1e). However, about half of the RCSMs manage to correct this bias, except AWI, LMD-CNRM, LMD-MPI and GUF, which amplify the variability bias of their GCMs. As noted for the mean bias and STD, spatial improvement is again concentrated along the coastlines and in semi-enclosed areas (Fig. S8).

The inability to improve area-averaged statistics indicates that large-scale SST biases in the Mediterranean are inherited from the driving GCMs and are not corrected by dynamical downscaling. This finding aligns with previous studies showing that downscaling adds value at small scales but is constrained by the large-scale errors of the driving models (Veljovic et al., 2010; Di Luca et al., 2013; Moreno-Chamarro et al., 2022).

Figure 2 AVi for median-removed SST in 1982–2020. (a) Multi-model AVi mean; grid points where fewer than 6 out of 9 simulations agree in sign with the mean are masked. (b) Aggregated PDFs of median-removed SST at AVi-significant grid points across all GCMs, RCSMs and OBS, with the corresponding AVi value (see Sect. 2.2.2). Individual model AVi maps for (c) CMCC, (d) AWI, (e) UNIBEL, (f) LMD-CNRM, (g) LMD-IPSL, (h) LMD-MPI, (i) CNRM-RCSM4, (j) CNRM-RCSM6, and (k) GUF. Positive (negative) values indicate that the RCSM provides a more (less) accurate representation of the observed distribution relative to the parent GCM; non-significant grid points are shown in white. (l) Area-averaged AVi computed over the Mediterranean domain (non-significant grid points set to zero), with red (blue) indicating positive (negative) area-averaged AVi. The percentage of positive and negative AVi grid points in the domain is reported in red (blue) when the positive (negative) fraction dominates.

3.1.2 AVi for SST

AVi for the median-removed SST (Fig. 2c–k) assesses the skill of RCSMs in representing SST PDFs at each grid point, independently of mean biases. Overall, there is a general improvement in the representation of the median-removed SST distribution across RCSMs. The multi-model mean (Fig. 2a) reveals a spatially distinct pattern: improvements are most pronounced in the eastern and central Mediterranean, including the Tyrrhenian Sea, while degradation prevails in the remaining part of the western basin. Six simulations exhibit positive basin-average AVi, ranging from +0.03 to +0.22, with improvements over more than 54 % of grid points (Fig. 2e, f, g, i, j, k). In contrast, AWI, LMD-MPI, and CMCC show predominantly negative AVi, with degradation affecting more than half of the basin grid points (Fig. 2c,d,h), particularly in the western Mediterranean, resulting in basin-average AVi values of −0.02, −0.04, and −0.08, respectively. Given this spatial heterogeneity, the mostly positive added value is not clearly reflected in the aggregated PDFs (Fig. 2b).

The positive added value in the central and eastern Mediterranean stems from RCSMs successfully simulating a broader SST PDF than their parent GCMs, bringing their distributions closer to observations (Fig. S9). Conversely, the Gulf of Lion remains a more nuanced case. SST variability in this region is strongly controlled by intense deep convection, which is highly sensitive to winds, buoyancy fluxes, stratification, and vertical mixing (Demirov and Pinardi, 2002; Ahumada and Cruzado, 2007; Macias et al., 2018). Although GCMs simulate larger SST variability that overlaps more closely with the observed distribution leading to negative AVi values (Fig. S10), this should be interpreted with caution: GCMs lack the fine-scale two-dimensional circulation, deep water mass formation, and wind-driven air-sea exchanges required to realistically represent deep convection in this region. The narrower SST distributions produced by RCSMs may therefore reflect a physically more consistent representation of deep convection, in which winter SST is set by the temperature of the deep layers, thereby limiting the range of surface variability. This interpretation is supported by the PDF of CNRM-RCSM6 (Fig. S10h), which shows a clear reduction on the cold tail of the distribution, consistent with a more realistic suppression of extreme cold events during convective episodes. This interpretation underscores a limitation of PDF-based added value metrics for SST: improved overlap with observed distributions does not always imply better physical representation, especially in regions where GCM performance may reflect error compensation rather than process fidelity.

Figure 3 As in Fig. 2 but for the SST upper tail distribution (SSTs above the 90th percentile) in 1982–2020.

Regarding warm SSTs (above the 90th percentile), the multi-model mean (Fig. 3a) highlights small common areas of improvements in the central Mediterranean, Tyrrhenian Sea, and Ligurian Sea, with negative values confined to the far eastern and western margins of the basin. Individual models exhibit substantially different spatial patterns of added and lost value (Fig. 3c–k), indicating a large inter-model spread in the representation of extreme SSTs. The aggregated PDFs (Fig. 3b) reveal that RCSMs generally stretch the upper tail of the SST distribution, bringing it closer to observations. However, at the individual model level, improvements arise either from an enhanced upper tail or from a constrained one (Fig. S11), both of which bring the simulated distribution closer to the observed SSTs. Six models exhibit zero (LMD-CNRM and UNIBEL) to positive (+0.14, CMCC) basin-average AVi, improving between 30 % to 81 % of grid points while limiting degradation to less than 31 %. The three models with negative AVi, ranging from −0.04 (LMD-MPI) to −0.07 (GUF and AWI), show degradation across 38 % to 57 % of grid points. Negative AVi values persist in the Alboran Sea and eastern Levantine basin, regions heavily influenced by advective dynamics and stratification (Pinardi and Masetti, 2000; Bergamasco and Malanotte-Rizzoli, 2010; Sánchez-Garrido and Nadal, 2022). This is unexpected given the known difficulty of GCMs to represent such dynamics, which would typically favor RCSM added value. The reason why the downscaling process does not lead to improved extreme SST representation remains an open question.

The added value analysis highlights a critical distinction: downscaling adds value in the representation of the shape of the SST distribution and upper tail SSTs, but mostly fails to correct biases in area-averaged statistics. However, improvements in SST statistics are present at the local level, especially along the coastlines and in semi-enclosed basins.

Figure 4 As in Fig. 1 for area-averaged MHW properties (1982–2020): (a) duration, (b) mean intensity, and (c) number of events. MHWs are calculated using a fixed baseline in 1982–2011 (see Sect. 2.2.2).

3.2 MHW properties

3.2.1 Mean errors

Here we evaluate the capacity of the RCSMs to correct the mean biases in MHW properties present in their driving GCMs. All GCMs considerably overestimate the duration of MHWs (Fig. 4a). However, downscaling markedly improves this metric: almost all RCSMs substantially reduce the GCM bias. This improvement is particularly pronounced in summer, autumn, and winter (Fig. S12b, c, d), when the magnitude of the GCM bias is strongest. Spring presents a more nuanced picture: the GCM bias is smaller in this season, and not all RCSMs succeed in further refining the signal (Fig. S12a). RCSMs broadly reduce the strong positive duration biases present in their driving GCMs across the Mediterranean basin (Fig. S13), with the spatial pattern of residual biases in the RCSMs largely mirroring that of the driving GCMs, albeit at a systematically lower magnitude. The only exception is AWI, which amplifies the GCM duration bias. Since these improvements are spatially ubiquitous and do not show systematic dependence on sub-basin location or local dynamical features, they can be attributed to a systematic enhancement brought by the downscaling process.

For MHW mean intensity, GCMs generally underestimate the observed values, but the ability of RCSMs to correct this bias varies considerably among models (Fig. 4b). Interestingly, the two CMIP6 GCMs accurately capture the mean MHW intensity, deviating from the general CMIP5 behaviour, consistent with the improved MHW skill of CMIP6 models reported in the literature (Qiu et al., 2021). Winter stands out as the only season in which every RCSM improves the parent GCM signal (Fig. S12h). In contrast, during the other seasons, several RCSMs tend to amplify the GCM bias (Fig. S12e, f, g). The spatial distribution of RCSM intensity biases (Fig. S14) highlights recurring hotspots of deterioration across both RCSMs and GCMs, particularly in the northern Adriatic Sea, the Gulf of Lion, and the Alboran Sea regions. Unlike the ubiquitous improvement seen for MHW duration, the general spatial pattern of intensity biases is preserved from GCMs to RCSMs, without a systematic reduction in magnitude. This suggests that what drives the consistent improvement in MHW duration does not extend to intensity, pointing to a more complex and model-dependent set of processes governing the latter.

Regarding the number of events, RCSMs generally perform well, improving the area-averaged bias relative to their driving GCMs (Fig. 4c). GCMs tend to underestimate the total number of MHWs, likely as a consequence of their tendency to overestimate event duration.

Figure 5 As in Fig. 2 but for MHW durations in 1982–2020. MHWs are calculated using a fixed baseline in 1982–2011 (see Sect. 2.2.2).

3.2.2 AVi for MHW Duration

MHW duration shows substantial and widespread positive AVi. The multi-model mean reveals that improvements are particularly robust over the Adriatic Sea, parts of the central Mediterranean, and along the coasts of Gulf of Lion and Ligurian Sea (Fig. 5a), while inter-model agreement is weak in the western Mediterranean and the easternmost Levantine basin. The positive AVi largely arises from the improved representation of short-lived MHW events, which are underestimated in the GCMs but more accurately captured by the RCSMs (Fig. 5b). Furthermore, GCMs tend to overestimate long-lived MHW events, while RCSMs better capture the upper tail of the duration distribution. This behaviour is consistent across all AVi-positive models (Fig. S15a, c, d, e, f, g, i). Seven out of nine models exhibit positive basin-average AVi, ranging from +0.03 to +0.09, although the spatial distribution of improvements varies between models (Fig. 5c–k). Positive AVi models improve between 17 % and 36 % of grid points, with degradation remaining below 8 %. Notably, the ability of RCSMs to improve MHW duration appears partly conditioned by the bias of their driving GCM: models driven by GCMs with a large bias (CMCC-CM, MPI-ESM-LR and IPSL-CM5A-MR, Fig. S13d, h, l) tend to exhibit larger positive AVi, while those driven by GCMs with a small bias (CNRM-CM5 and CNRM-ESM2-1, Fig. S13f, j, r) show more limited improvement.

This result reflects the broad improvement in SST variability discussed above and agrees with prior findings (Pilo et al., 2019; Guo et al., 2019; Capotondi et al., 2024), showing that improved performance in MHW duration can be achieved by enhancing horizontal resolution, which plays a key role in better resolving the temporal build-up and dissipation of MHW events through enhanced mesoscale dynamics. RCSMs increase the realism of local and regional dynamics, producing SST fields with sharper gradients and more accurate spatial variability that better capture the conditions leading to the onset and decay of short-lived events. In the Mediterranean, MHW termination is frequently driven by regional wind systems which trigger rapid surface cooling through enhanced heat loss (Bonino et al., 2025). These mesoscale wind patterns are inherently better captured at higher resolution (Herrmann et al., 2011; Obermann et al., 2018), enabling RCSMs to more realistically simulate the atmospheric forcing episodes responsible for MHW decay. Beyond horizontal resolution, air-sea flux parameterizations and the thickness of the uppermost ocean layer further modulate the ocean's response time to such forcing, and are therefore also expected to influence MHW duration.

A notable exception is AWI, which exhibits negative area-averaged AVi for duration, degrading about 12 % of grid points while improving only 6 %, suggesting that increased resolution alone does not universally guarantee better performance and that other aspects of the RCSM configuration must also be considered. AWI is the only model in the ensemble with a coarser uppermost ocean layer, and, unlike the other RCSMs, it uses a globally coupled ocean component (Sein et al., 2015; Parras-Berrocal et al., 2020). Both factors may contribute to more damped SST temporal variability and a reduced ability to resolve the rapid onset and decay of short-lived MHW events, though disentangling their respective roles is beyond the scope of this study.

Figure 6 As in Fig. 2 but for MHW mean intensities in 1982–2020. MHWs are calculated using a fixed baseline in 1982–2011 (see Sect. 2.2.2).

Figure 7 Area-Averaged AVi for MHW intensity across the ensemble for different analysis configurations. Red (blue) shading denotes positive (negative) values, with color intensity proportional to magnitude.

3.2.3 AVi for MHW Intensity

MHW mean intensity exhibits a more variable picture across the ensemble compared to duration. The multi-model mean reveals that the Gulf of Lion is the only region with robust improvement, while inter-model agreement remains weak elsewhere (Fig. 6a). A notable exception is CMCC, which displays consistently strong performance across the whole basin, particularly in the central and the western Mediterranean (Fig. 6c), yielding the largest area-averaged AVi (+0.57), improving the signal over 94 % of the grid points. Five additional models (LMD-IPSL, UNIBEL, LMD-MPI, CNRM-RCSM4, CNRM-RCSM6) achieve positive area-averaged AVi (from +0.04 to +0.17) with areas of improvement between 22 % and 50 %, although the latter are spatially fragmented and accompanied by areas of deterioration (Fig. 6e, g, h, i, j). The remaining models (AWI, LMD-CNRM, GUF) show negative area-averaged AVi, ranging from −0.06 to −0.18, with extensive areas of deterioration (Fig. 6d, f, k) ranging from 35 % to 44 % of the Mediterranean domain. The aggregated PDFs show that RCSMs improve the representation of intensities by slightly shifting the distribution towards warmer values (Fig. 6b). Such behaviour characterizes AVi-positive models (Fig. S16a, c, e, f, g, h). In contrast, AWI and LMD-CNRM show a systematic shift toward cooler intensities compared to their parent GCMs (Fig. S16b, d), while GUF exhibits a broader distribution than both its parent model and OBS (Fig. S16i), leading to an overestimation of strong intensities rather than a simple mean shift.

Basin-average AVi values for MHW maximum intensity, defined as the peak SST anomaly within each event, closely mirror those obtained for mean intensity, suggesting that added value is not metric-specific but instead reflects broader discrepancies in the magnitude of SST departures from climatological baselines during extreme events (Fig. 7). Analyses over two additional sub-periods (1982–2005 and 2006–2020) yield results consistent with those obtained for the full 1982–2020 interval, with a general decrease in basin-average AVi magnitude but minimal change in sign or relative model ranking, demonstrating the robustness of these findings to the choice of evaluation and baseline periods (Fig. 7). Most RCSMs underestimate warming relative to their driving GCMs (Fig. 1d), reducing warm SST anomalies during late-period MHW events. Accurate representation of mean-state evolution is indeed essential for correctly characterizing present and future MHW properties (Frölicher et al., 2018; Schlegel et al., 2019; Oliver et al., 2021; Martínez et al., 2023). Removing the long-term trend generally reduces AVi values but does not alter their signs, indicating that trend biases modulate the magnitude of added value without changing the overall interpretation (Fig. 7). Similarly, applying a 15-year shifting climatology deteriorates basin-average AVi across nearly all models yet does not qualitatively alter model behaviour: AWI, LMD-CNRM, and GUF remain negative-AVi models, confirming that their degradation is not solely attributable to trend misrepresentation.

The representation of SST residuals' variability controls the magnitude of transient thermal anomalies. Negative‑AVi RCSMs degrade this variability: AWI and LMD-CNRM have larger negative biases, whereas GUF shows an amplified positive bias (Fig. 1e). We believe that the reasons for the degradation of SST residuals' variability and MHW intensity are model-dependent, as they strongly affect only a few models of the Med-CORDEX ensemble. For AWI, a key possible limitation is its coarse vertical discretization: its surface layer (16 m) is thicker than that of its parent GCM (12 m) (Table 1). Such a thick surface layer likely damps short-term temporal fluctuations, smooths stratification gradients, and weakens warm anomalies, whereas resolving the fine upper-ocean structure is crucial for simulating MHW processes (Capotondi et al., 2024; Darmaraki et al., 2024). This is supported by UNIBEL, forced by the same GCM but using sigma coordinates with a much thinner first layer (1.8 m); despite coarser horizontal resolution (30 vs. 7–25 km), it attains positive added value in both duration and intensity, showing that thickness of the first layer can be as important as horizontal resolution for MHW characteristics. GUF mainly reflects an inherited problem: its parent EC-Earth3-Veg already overestimates the variability of SST residuals (Fig. 1e) and has the broadest MHW intensity distribution among all GCMs (Fig. S16i). Such variability could have been further amplified by the RCSM, consistent with Kim et al. (2023). Hence, RCSM performance for MHW intensity is highly sensitive to the quality of the parent model and the variability transmitted across scales. For LMD-CNRM, the negative added value for MHW intensity appears rooted in the atmospheric component. LMD-CNRM and CNRM-RCSM4 share the same ocean configuration but differ in: (i) atmospheric model (stretched-grid LMDz4 vs. ALADIN5.2) and (ii) spectral nudging of atmospheric fields (absent vs. implemented) (L’Hévéder et al., 2013; Sevault et al., 2014). Although LMDz4-regional has finer horizontal resolution, it employs fewer vertical levels and exhibits weaker SST variability. The larger variability in SST residuals in CNRM-RCSM4 compared to LMD-CNRM (0.39 vs. 0.34 over 2006–2020) may be partly attributable to spectral nudging, as CNRM-CM5 itself exhibits the largest SST residual variability among the driving GCMs over the same period (0.46), suggesting that this signal is inherited from the parent model rather than generated by the downscaling process. However, the exact role of spectral nudging in shaping extreme events in downscaling simulations remains an open question requiring targeted future work.

These results highlight that the added value for MHW intensity is generally positive although not uniform across the ensemble nor systematically linked to increased horizontal resolution. Improvements in intensity are thus governed by model-specific processes that are not straightforwardly improved by resolution alone, consistently with Pilo et al. (2019).

4 Discussion and Conclusions

Med-CORDEX downscaling simulations offer a powerful framework for assessing climate-change impacts in the Mediterranean Sea, a region where MHWs pose significant risks to ecosystems and human activities and reliable high-resolution simulations are crucial for defining the characteristics of future marine thermal extremes. Notably, Med-CORDEX represents one of the few coordinated multi-model ensembles worldwide to include a fully coupled ocean component, making it a unique tool for studying air-sea interactions and marine extremes at regional scale. Exploiting this ensemble approach rather than relying on a single model or a single pair of runs allows us to derive more robust conclusions about the added value of regional coupling, by separating systematic signals from model-specific behaviour. Here we ask whether Med-CORDEX RCSMs systematically improve the representation of SST and MHW properties relative to their driving GCMs, and under which conditions and at which scales this added value emerges. Identifying where and why regional downscaling succeeds or fails represents a step towards improving future generations of regional climate models.

For SST, the added value of dynamical downscaling is scale-dependent and metric-specific. Most RCSMs fail to correct area-averaged biases in the mean, standard deviation, 90th percentile, and trend, reflecting the inability of dynamical downscaling to reduce large-scale errors inherited from the driving GCMs. Nevertheless, 7 out of 9 RCSMs improve the spatial structure of the mean SST field relative to their driving GCMs and all RCSMs show error improvement localized along the coastlines and in semi-enclosed areas such as the Adriatic, Alboran, and Aegean Seas, regardless of their basin-averaged performance. Together, these results confirm that the added value of dynamical downscaling for SST is primarily a regional-scale phenomenon, concentrated where fine-scale ocean-atmosphere interactions and topographic constraints are most influential. RCSMs generally improve the representation of the shape of the SST distribution, particularly for median-removed SST in the central and eastern Mediterranean, where they simulate a broader PDF closer to observations. However, the Gulf of Lion stands out as a region of consistent negative AVi, where regional downscaling suppresses SST variability, though this may reflect a more physically consistent representation of deep convection rather than a genuine deterioration. For extreme SSTs, RCSMs generally stretch the upper tail of the distribution towards observations, with 6 out of 9 models showing positive or zero area-averaged AVi, the latter cases still featuring a larger fraction of improving than degrading grid points. Improvement remains elusive in the Alboran Sea and eastern Levantine basin for reasons that remain unclear.

For MHW duration, dynamical downscaling provides consistent and spatially widespread added value. Almost all RCSMs substantially reduce the strong positive duration biases of their parent GCMs, with improvements particularly robust over the Adriatic Sea. This added value stems primarily from the improved representation of short-lived events. The ubiquitous nature of these improvements across the ensemble points to a general enhancement brought by the downscaling process, consistent with previous studies showing that increased horizontal resolution and better representation of mesoscale dynamics improve MHW duration (Pilo et al., 2019; Guo et al., 2019; Capotondi et al., 2024). Improvements in MHW duration reinforce the usefulness of Med-CORDEX RCSMs for impact and ecosystem studies in which the persistence of thermal stress is critical. However, exceptions confirm that enhanced horizontal resolution alone is insufficient if other aspects of the model configuration are degraded.

The picture is more variable for MHW intensity. The majority of simulations improve the representation of this metric by reducing the underestimation found in the parent GCMs, yet three RCSMs exhibit clear deterioration. Accurately capturing MHW intensity requires a faithful representation of intra-seasonal SST variability, which governs the magnitude of transient thermal anomalies. In the models showing negative added value, this variability is plausibly linked to model-specific factors, though attributing specific failures to specific model components would require dedicated sensitivity experiments. Consequently, higher horizontal resolution is a necessary but not sufficient condition for improved intensity representation, and concurrent advances in other model components are needed.

We note that RCSM added value for MHW duration tends to be larger when the driving GCM has a larger bias, which is expected by construction. This relationship is not observed for MHW intensity, where added value appears instead to be governed by model-specific configuration factors. This bias-dependence for MHW duration, in any case, concerns basin-averaged metrics and does not reflect an improved representation of the local and regional-scale processes (dynamics of coastal and semi-enclosed areas) that RCSMs are designed to resolve and where their added value remains concentrated.

A key methodological consideration concerns the implicit bias correction introduced by the definition of MHWs. By computing thresholds relative to each model's own climatology, mean-state biases in the spatio-temporal structure of SST are removed for all data sources. This limits the detectability of added value to the variability component of SST. This may partly explain why added value is more readily identified for SST anomalies than for MHW characteristics. Future studies could complement the classic MHW definition with absolute threshold methods to recover the mean-state dimension of added value.

To further explore the relationship between model configuration and added value, we correlated the area-averaged AVi for each of the four metrics with RCSM improvement factors in key configurations: ocean horizontal resolution, uppermost layer thickness, atmospheric resolution, and coupling frequency, each defined as the ratio of the GCM to the RCSM value for these specific aspects (see Table 1). Of the sixteen correlations tested, only two are statistically significant: AVi for MHW duration (r=0.85, p<0.01) and intensity (r=0.77, p=0.01) both correlate with uppermost layer thickness improvement, suggesting that a thinner surface layer is a key factor in improving MHW representation. The absence of significant correlations with horizontal resolution does not contradict the literature but suggests that, within the Med-CORDEX ensemble, inter-model variance in MHW added value is better explained by the uppermost layer thickness. More broadly, the largely non-significant correlations point to either a predominant role of 1D processes and physical parameterizations, or to errors in GCM boundary conditions overwhelming the RCSMs' structural improvements.

The present study focuses exclusively on SST and derived MHW properties, without examining the atmospheric drivers that govern SST variability and MHW onset. Recent work has shown that Mediterranean MHWs are primarily driven by anomalies in surface heat fluxes, in particular reduced latent heat loss associated with weak winds, enhanced near-surface humidity, and reduced cloud cover associated with enhanced surface downwelling radiation (Bonino et al., 2025; Garfinkel et al., 2026; Paredes-Fortuny et al., 2025). RCSMs improve the representation of atmospheric fronts and mesoscale wind systems (Herrmann et al., 2011; Obermann et al., 2018), which is expected to translate into better representation of surface latent heat fluxes and evaporation (Hamon et al., 2016; Herrmann and Somot, 2008), key drivers of MHW onset and decay. However, a process-based assessment of how such atmospheric processes influence MHWs across the Med-CORDEX ensemble was beyond the scope of this study but represents an important direction for future work.

Substantial differences among Med-CORDEX RCSM configurations (boundary forcing, physics, coupling, vertical and horizontal resolution) introduce structural uncertainty that complicates the attribution of added value to individual model features. Furthermore, most RCSMs are driven by a single GCM, which further limits inter-model comparability, as distinct large-scale biases and trends inherited from the driving GCM are imprinted onto the regional simulations. Future work on added value quantification would benefit from coordinated multi-model ensembles driven by a common GCM and a more strict protocol, enabling clearer identification of the mechanisms through which regional downscaling enhances or limits model skill.