# Future evolution of Marine Heatwaves in the Mediterranean Sea

Sofia Darmaraki, Samuel Somot, Florence Sevault, Pierre Nabat, William

David Cabos Narvaez, Leone Cavicchia, Vladimir Djurdjevic, Laurent Li,

Gianmaria Sannino, Dmitry V. Sein

► To cite this version:

Sofia Darmaraki, Samuel Somot, Florence Sevault, Pierre Nabat, William David Cabos Narvaez, et al.. Future evolution of Marine Heatwaves in the Mediterranean Sea. Climate Dynamics, 2019, 53, pp.1371-1392. 10.1007/s00382-019-04661-z. insu-03668353

HAL Id: insu-03668353

https://insu.hal.science/insu-03668353

Submitted on 15 May 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Future evolution of Marine Heatwaves in the Mediterranean Sea

Sofia Darmaraki¹ · Samuel Somot¹ · Florence Sevault¹ · Pierre Nabat¹ · William David Cabos Narvaez² · Leone Cavicchia³ · Vladimir Djurdjevic⁴ · Laurent Li⁵ · Gianmaria Sannino⁶ · Dmitry V. Sein⁷,⁸

Received: 28 June 2018 / Accepted: 1 February 2019 / Published online: 21 February 2019
© The Author(s) 2019

## Abstract

Extreme ocean warming events, known as marine heatwaves (MHWs), have been observed to perturb significantly marine ecosystems and fisheries around the world. Here, we propose a detection method for long-lasting and large-scale summer MHWs, using a local, climatological 99th percentile threshold, based on present-climate (1976–2005) daily SST. To assess their future evolution in the Mediterranean Sea we use, for the first time, a dedicated ensemble of fully-coupled Regional Climate System Models from the Med-CORDEX initiative and a multi-scenario approach. The models appear to simulate well MHW properties during historical period, despite biases in mean and extreme SST. In response to increasing greenhouse gas forcing, the events become stronger and more intense under RCP4.5 and RCP8.5 than RCP2.6. By 2100 and under RCP8.5, simulations project at least one long-lasting MHW every year, up to three months longer, about 4 times more intense and 42 times more severe than present-day events. They are expected to occur from June-October and to affect at peak the entire basin. Their evolution is found to occur mainly due to an increase in the mean SST, but increased daily SST variability also plays a noticeable role. Until the mid-21st century, MHW characteristics rise independently of the choice of the emission scenario, the influence of which becomes more evident by the end of the period. Further analysis reveals different climate change responses in certain configurations, more likely linked to their driving global climate model rather than to the individual model biases.

**Keywords** Marine Heatwaves · Mediterranean Sea · Coupled regional climate models · Future scenario · Extreme ocean temperatures · Med-CORDEX · Climate change · Climate simulations

## 1 Introduction

Episodes of large-scale warm temperature anomalies in the ocean may prompt substantial disruptions to marine ecosystems (Frölicher and Laufkötter 2018; Hobday et al. 2016) and major implications for fisheries as well (Mills et al. 2013). Known as marine heatwaves (MHW), these extreme events describe abrupt but prolonged periods of high sea surface temperatures (SST) (Scannell et al. 2016) that can occur anywhere, at any time, with the potential to propagate deeper to the water column (Schaeffer and Roughan 2017). They have received little attention until improved

observational systems revealed adverse consequences emanating from them. Their occurrence is likely to intensify under continued anthropogenic warming (Frölicher et al. 2018; Oliver et al. 2018a), engendering the need for a more comprehensive examination of their spatiotemporal distribution and underlying physical causes.

In the Mediterranean area, a well-known “Hot Spot” region for climate change (Giorgi 2006), the annual mean basin SST by the end of the 21st century is expected to increase from + 1.5 °C to + 3 °C relative to present-day levels, depending on the greenhouse gas (GHG) emission scenario (Somot et al. 2006; Mariotti et al. 2015; Adloff et al. 2015). This significant rise in SST is expected to accelerate future MHW occurrence, in congruence with projections for GHG-induced heat stress intensification of 200–500% throughout the region (Diffenbaugh et al. 2007). The Mediterranean area’s sensitivity to increased GHG forcing is mainly attributed to a significant mean warming and increased interannual warm-season variability, along

**Electronic supplementary material** The online version of this article (https://doi.org/10.1007/s00382-019-04661-z) contains supplementary material, which is available to authorized users.

<figure><img src="image_2.png" /></figure>

Sofia Darmaraki
sofia.darmaraki@meteo.fr

Extended author information available on the last page of the article
with a reduction in precipitation (Giorgi 2006). A recent study has already identified significant increases in MHWs globally over the last century, including the Mediterranean Sea (Oliver et al. 2018a).

In fact, one of the first-detected MHWs worldwide occurred in the Mediterranean in the summer of 2003: Surface anomalies of 2–3 °C above climatological mean lasted for over a month due to significant increases in air-temperature and a reduction of wind stress and air-sea exchanges (Grazzini and Viterbo 2003; Sparnocchia et al. 2006; Olita et al. 2007). These factors seem to have also triggered an anomalous SST warming in the eastern Mediterranean area during the heatwave of 2007, at the order of + 5 °C above climatology (Mavrakis and Tsiros 2018). Since then, numerous studies have explored the modulating factors behind individual events around the world. For instance, a combination of local oceanic and large-scale atmospheric forcing was suggested for the Australian MHW of 2011 (Feng et al. 2013; Benthuysen et al. 2014) and the persistent, multi-year (2014–2016) “Pacific Blob” (Bond et al. 2015; Di Lorenzo and Mantua 2016). Other events have been attributed to mainly atmosphere-related drivers, such as the 2012 Atlantic MHW (Chen et al. 2014, 2015) and the extreme marine warming across Tropical Australia (Benthuysen et al. 2018), or to ocean-dominating forcing like the 2015/2016 Tasman Sea MHW (Oliver et al. 2017). The importance of regional influences was further noted in coastal MHWs in South Africa (Schlegel et al. 2017a) and during subsurface MHW intensification around Australia (Schaeffer and Roughan 2017).

As a result of these events, severe impacts on marine ecosystems have been documented worldwide, including biodiversity die-offs and tropicalisation of marine communities (Wernberg et al. 2013, 2016), extensive species migrations (Mills et al. 2013), strandings of marine mammals and seabirds, toxic algal blooms (Cavole et al. 2016) and extensive coral bleaching (Hughes et al. 2017). In the Mediterranean Sea in particular, unprecedented mass mortality events and changes in community composition due to extreme warming were reported in the summers of 1999 (Perez et al. 2000; Cerrano et al. 2000; Garrabou et al. 2001; Linares et al. 2005), 2003 (Garrabou et al. 2009; Schiaparelli et al. 2007; Diaz-Almela et al. 2007; Munari 2011), 2006 (Kersting et al. 2013; Marba and Duarte 2010) and 2008 (Huete-Stauffer et al. 2011; Cebrian et al. 2011), affecting a wide variety of species and taxa (e.g. 80 % of Gorgonian fan colonies and seagrass *Posidonia oceanica*). MHWs can be especially lethal for organisms with reduced mobility that are usually limited to the upper water column; Their severity is determined by both temperature and duration (Galli et al. 2017). Finally, cascading effects have also been observed in fisheries, resulting in huge financial

losses and even economic tensions between nations (Mills et al. 2013; Cavole et al. 2016; Oliver et al. 2017).

However, despite the growing body of MHW-related literature, systematic examination of MHWs as distinct exceptional events with intensity, frequency and duration has only just emerged. Although marine extremes have been investigated before, only a few studies have analysed past trends in extreme ocean temperatures (e.g. Scannell et al. 2016; MacKenzie and Schiedek 2007) and even fewer have dealt with their future evolution. For instance, past trends of extreme SST have been investigated in coastal regions (Lima and Wethey 2012) and through thermal-stress-related coral bleaching records (Lough 2000; Selig et al. 2010; Hughes et al. 2018). Using a more standardised framework, past MHW occurrences have been studied in the Tasman Sea (Oliver et al. 2018b) and the global ocean (Oliver et al. 2018a). For the 21st century, MHW projections have been performed so far on a global scale, with the use of multi-model setups from CMIP5 (Frölicher et al. 2018) and CMIP3 (Hobday and Pecl 2014) and under different GHG emission scenarios. On a regional scale though, ocean extremes have been assessed in Australia (King et al. 2017) and the Tasman Sea (Oliver et al. 2014).

The above-mentioned studies used different definitions for extreme warm temperatures, with some adopting a recent standardised MHW approach proposed by Hobday et al. (2016). The set of dedicated statistical metrics developed in this framework allows for a consistent definition and quantification of the MHW properties. A MHW is now described as a “discrete, prolonged, anomalously warm water event at a particular location”. Using this definition, Schlegel et al. (2017b), for example, identified an increase in MHW frequency around South Africa for the period 1982–2015, while Schaeffer and Roughan (2017) demonstrated subsurface intensification of MHWs in coastal SE Australia between 1953–2016. A linear classification scheme was also proposed by Hobday et al. (2018), where MHWs are defined based on temperature exceedance from local climatology.

In the case of the Mediterranean Sea, however, little is known about past or future MHW trends and their underlying mechanisms. The MHW-related research has mostly been focused on local ecological impacts without systematically assessing MHW occurrence. According to Rivetti et al. (2014) and Coma et al. (2009), most of the mass mortalities documented in the basin were related to positive thermal anomalies in the water column that occurred regionally during the summer. Although they have been reported with increased frequency since the early 1990s, their occurrence has been observed as early as the 1980s. Meanwhile, the evolution of extreme Mediterranean SST in the 21st century has so far been examined in relation to the thermotolerance responses of certain species. For instance, Jordà et al. (2012) used an ensemble of models

<figure><img src="image_3.png" /></figure>
under the moderately optimistic scenario for GHG emissions A1B and suggested an increased seagrass mortality in the future around the Balearic islands due to a projected rise of the annual maximum SST by 2100. Similarly, Bensoussan et al. (2013) evaluated the thermal-stress related risk of mass mortality in Mediterranean benthic ecosystems for the 21st century, based on the average warming estimated between 2090–2099 and 2000–2010, under the pessimistic future warming scenario A2. Finally, Galli et al. (2017) showed an increase in MHW frequency, severity and depth extension in the basin, assuming exceedances from species-specific thermotolerance thresholds under the high-emission IPCC RCP8.5 scenario. [The A1B and A2 emission scenarios correspond to projections of a likely, mean temperature change of 1.7–4.4 °C and 2.0–5.4 °C respectively by the end of the 21st century (IPCC 2007), whereas RCP2.6, RCP4.5 and RCP8.5 to a likely change of 0.3–1.7 °C, 1.1–2.6 °C and 2.6–4.8 °C respectively, by the end of the period (Kirtman et al. 2013)].

In addition, our understanding of the Mediterranean Sea's response to future climate change to date mostly relies on ensembles of low resolution GCMs (CMIP5) (e.g. Jordà et al. 2012; Mariotti et al. 2015) or on numerical experiments carried out with a single regional ocean model under different emission scenarios (i.e Somot et al. 2006; Bensoussan et al. 2013; Adloff et al. 2015; Galli et al. 2017). Consequently, the various sources of uncertainty related to the choice of the socio-economic scenario, choice of climate model and natural variability have not been properly taken into account by climate change impact studies on Mediterranean Sea ecosystems and maritime activities. Since there is an evident link between distinctive climate anomalies and notable ecosystem effects (e.g in the Mediterranean Sea, Crisci et al. 2011; Bensoussan et al. 2010), it is important to address these uncertainties by considering different possible climate futures through multi-model, multi-scenario set ups when possible.

In this context, the aim of this study is to provide a robust assessment of the future evolution of summer MHWs in the Mediterranean Sea using an ensemble of high-resolution coupled regional climate system models (RCSM), driven by GCMs and a multi-scenario approach (RCP2.6, RCP4.5, RCP8.5). The RCSM's ability to reproduce Mediterranean SST features is first evaluated against satellite data. Then, a MHW spatiotemporal definition, based on SST and on Hobday et al. (2016)'s recommendations, is developed and applied to study the response of extreme thermal events to future climate change. For the first time, changes in summer Mediterranean MHW frequency, duration, intensity and severity are investigated with respect to an envelope of possible futures.

This paper is organised as follows: in Sect. 2 we present the ensemble of RCSMs along with the methodology

proposed for the detection and characterisation of summer MHWs. Model evaluation against observed mean and extreme SST is performed in Sect. 3, using daily SST data. We also describe the future evolution of Mediterranean SST and MHW properties under different greenhouse gas emission scenarios from 1976–2100. A discussion and summary of the results are presented in Sects. 4 and 5.

# 2 Material and methods

## 2.1 Model data and simulations

An ensemble of six coupled RCSMs (CNRM-RCSM4, LMDZ-MED, COSMOMED, ROM, EBU-POM, PROTHEUS) with different Mediterranean configurations is employed in this study. Participant members are provided by six research institutes from the Med-CORDEX initiative [Ruti et al. (2016), https://www.medcordex.eu/] and each simulation will be herein referred to by the name of the corresponding institute, as mentioned in Table 1 (e.g. simulations with the CNRM-RCSM4 model will be referred as CNRM, etc.). Med-CORDEX can be considered as a multi-model follow-up to the CIRCE project (Gualdi et al. 2013), which studied the Mediterranean Sea under a single scenario (A1B) with most of the simulations stopped in 2050.

One novel aspect of the Med-CORDEX ensemble is that all models have a high-resolution oceanic (eddy-resolving) and atmospheric component as well as high coupling frequency (see Table 1). The free air-sea exchanges offered by their high-resolution interface is also an advantage for the MHW representation, which depends on ocean-atmosphere interactions. The domains cover the entire Mediterranean and a small part of the Atlantic, while the Black Sea and Nile river are respectively parametrized or represented with climatologies (except for AWI/GERICS, in which the oceanic component is global and explicitly simulates the Black Sea). Boundary conditions come from 4 different general circulation models of CMIP5. Information about each coupled system is summarised in Table 1. To avoid biasing results towards one or more members of the ensemble, only the realization with the highest resolution is selected for each model.

All the numerical simulations produced daily SST data (3D temperatures were stored at a monthly scale) between 1950–2005 for the historical experiment (HIST) and for 2006–2100 under the Representative Concentration Pathway RCP8.5 (high-emission scenario), RCP4.5 (moderate-emission scenario), RCP2.6 (low-emission scenario) IPCC scenarios. As the models use boundary conditions from CMIP5, which are not in phase with the observed variability, simulation chronology does not represent the actual conditions that correspond to each calendar year. Instead, they are expected
<table><caption>**Table 1** Characteristics of the Med-CORDEX coupled regional climate system models (RCSM) and the simulations used in this study</caption><thead><tr><th>INSTITUTE</th><th>CNRM</th><th>LMD</th><th>CMCC</th><th>AWI/GERICS</th><th>U.BELGRADE</th><th>ENEA</th></tr></thead><tbody><tr><td colspan="7">Model characteristics</td></tr><tr><td>RCSM name</td><td>CNRM-RCSM4</td><td>LMDZ-MED</td><td>COSMOMED</td><td>ROM</td><td>EBU-POM</td><td>PROTHEUS</td></tr><tr><td>Driving GCM</td><td>CNRM-CM5</td><td>IPSL-CM5A-MR</td><td>CMCC-CM</td><td>MPI-ESM-LR</td><td>MPI-ESM-LR</td><td>CNRM-CM5</td></tr><tr><td>Med.Sea Model</td><td>NEMOMED8</td><td>NEMOMED8</td><td>NEMO-MFS</td><td>MPIOM</td><td>POM</td><td>MedMIT8</td></tr><tr><td>Ocean Res.</td><td>9–12 km</td><td>9–12 km</td><td>6–7 km</td><td>10–18 km</td><td>30 km</td><td>13 km</td></tr><tr><td>Num. of z-Levels (Ocean)</td><td>43</td><td>43</td><td>72</td><td>40</td><td>21</td><td>42</td></tr><tr><td>SST (1st layer depth)</td><td>6 m</td><td>6 m</td><td>1.5 m</td><td>16 m</td><td>1.8 m</td><td>10 m</td></tr><tr><td>Timestep (Ocean)</td><td>1200 s</td><td>1200 s</td><td>480 s</td><td>900 s</td><td>360 s</td><td>600 s</td></tr><tr><td>Atmosphere model</td><td>ALADIN-climate</td><td>LMDZ</td><td>CCLM</td><td>REMO</td><td>Eta/NCEP</td><td>RegCM</td></tr><tr><td>Atmosphere Res.</td><td>50 km</td><td>30 km</td><td>50 km</td><td>25 km</td><td>50 km</td><td>30 km</td></tr><tr><td>Coupling frequency</td><td>Daily</td><td>Daily</td><td>80 min</td><td>60 min</td><td>6 min</td><td>6 h</td></tr><tr><td colspan="7">Numerical Simulations</td></tr><tr><td>SPIN UP</td><td>130 years</td><td>40 years</td><td>25 years</td><td>56 years</td><td>5 years</td><td>No Spin Up</td></tr><tr><td>Initial conditions</td><td>MEDATLAS</td><td>MEDATLAS</td><td>MEDATLAS</td><td>MEDATLAS</td><td>MEDATLAS</td><td>MEDATLAS</td></tr><tr><td>HIST</td><td>1950–2005</td><td>1950–2005</td><td>1950–2005</td><td>1950–2005</td><td>1950–2005</td><td>1979–2005</td></tr><tr><td>RCP8.5</td><td>2006–2100</td><td>2006–2100</td><td>2006–2099</td><td>2006–2099</td><td>2006–2100</td><td>–</td></tr><tr><td>RCP4.5</td><td>2006–2100</td><td>2006–2100</td><td>2006–2099</td><td>2006–2099</td><td>–</td><td>2006–2100</td></tr><tr><td>RCP2.6</td><td>2006–2100</td><td>–</td><td>–</td><td>–</td><td>–</td><td>–</td></tr><tr><td>References</td><td>Sevault et al. (2014)</td><td>L'Hévéder et al. (2013)</td><td>Cavicchia et al. (2015)</td><td>Sein et al. (2015)</td><td>Djurdjevic and Rajkovic (2008)</td><td>Artale et al. (2010)</td></tr></tbody></table>

More information on MEDATLAS initial conditions can be found in Rixen et al. (2005)

to represent the climate statistics of each period (e.g. average, standard deviation) well. We use SST instead of deeper layer temperatures, as both the models' behaviour and the MHW identification technique can be evaluated at a larger scale using satellite data. A total of 17 simulations were used from six models with variable resolution (Table 1). For the purposes of our analysis, we define 30-year periods from the HIST run between 1976–2005 (from this moment on referred as HIST), the near future (2021–2050) and the far future (2071–2100).

In the case of ENEA, the HIST run span from 1979–2005 due to different simulation initialization, while CMCC and AWI/GERICS simulations reached 2099. The spin-up strategy of the Med-CORDEX ensemble was not prescribed, therefore it was different for every configuration. The lack of a long spin-up (e.g. U.BELGRAD, ENEA) could be detrimental for temperatures at deeper layers but not so relevant for the SST evolution. For the CNRM model, a constant monthly flux (atmosphere to ocean) correction was applied to minimise identified biases, with no significant influence on the climate change signal. Also, a slightly intense SST

signal in the Alboran Sea was noted in the U.BELGRADE configuration for 2021–2050 under RCP8.5 and is probably linked to the simple representation of the connection between the Mediterranean Sea and the Atlantic Ocean in the model: the open boundary condition, as defined in the POM model, was applied in single model point defined on the strait of Gibraltar, without any buffer zone and with prescribed boundary conditions in the Atlantic Ocean; this is, on the other hand, a common approach in many models. Finally, an error has been recently reported concerning the CNRM-CM5 GCM files that were used as atmospheric lateral boundary conditions for CNRM and ENEA (http://www.umr-cnrm.fr/cmip5/spip.php?article24), but this likely has no significant effect on the long-term climate change signal.

Working from the hypothesis that MHWs are usually confined close to the surface, in this study we consider that the model SST data of the 1st layer depth represent surface temperatures between 1–16 m, depending on the model. We acknowledge, however, that MHWs may penetrate deeper to the water column under certain conditions, but assume for
the time being that SST is a reliable sign of possible harmful conditions for deeper layers.

## 2.2 Reference dataset

In order to evaluate the model's capability to simulate trends in regional extreme thermal events, we first perform comparisons with satellite data (OBS) provided by the Copernicus Marine Service and CNR - ISAC ROME. More specifically, the Mediterranean Sea high-resolution L4 dataset is employed, providing daily, reprocessed SSTs on a 0.04 °grid, an interpolation of remotely sensed SSTs from the Advanced Very High Resolution Radiometer (AVHRR) Pathfinder Version 5.2 (PFV52) onto a regular grid (Pisano et al. 2016). They are obtained over a 30-year period of January 1982 to December 2012 and are used as a reference for the models' performance in the mean and extreme climate in the Mediterranean Sea. With the aim of validating the “present-day” climate, we choose the 30-year period (1976–2005) in the model HIST runs that has the greatest overlap with the observed 30-year period (1982–2012). Prior to performing any calculations and in order to compare the results between the models and observations, we first interpolated every dataset to the NEMOMED8 grid, already in use by 2 RCSMs, by implementing the nearest neighbour method.

## 2.3 Defining marine heatwaves

As for their atmospheric counterparts, there is no universal definition for MHWs. However, certain metrics can be applied to compare different events in space and time. In this research, the qualitative MHW definition proposed by Hobday et al. (2016) is followed. We use it as a baseline for developing a quantitative method that will identify MHWs, namely in the summer months, based on the climatology and the geographical characteristics of the area. Although we recognise that heatwaves in colder months might also be essential for certain species, we choose to focus on extreme events related to the highest annual SSTs, when organisms may be beyond their optima, as seen by previous mass mortality events in the Mediterranean (e.g. MHW of 1999 and 2003).

According to Hobday et al. (2016) a MHW is a “prolonged, anomalously warm water event at a particular location” and it should be defined relative to a 30-year period. In our case, a subset of the HIST experiments (1976–2005) and the 1982–2012 period for the observations are chosen, representing the average climate in the latter half of the 20th century. In order to achieve a homogenised yet area-specific temperature diagnostic, for every year of the reference period (HIST) we first compute the 99th quantile of daily SST ($SST_{99Q}$) for every grid point. Then we average these 30 years of extreme values, constructing a 2D threshold map.

Note that individual threshold maps were created for each dataset separately, accounting for the different model characteristics (e.g SST bias). An “anomalously warm day” at every grid point is then any given day when the local $SST_{99Q}$ threshold is exceeded. However, in order to be classified as a “prolonged” event, we set the minimum duration of a MHW to 5 days, following Hobday et al. (2016). Further, we aim to identify long-lasting events, since most of the previous mass mortalities in the basin occurred during thermal anomalies that lasted for more than 5 days (e.g. (Garrabou et al. 2009; Di Camillo and Cerrano 2015; Cerrano et al. 2000; Cebrian et al. 2011). In addition, the average present-day MHW duration in the basin was found around 10 days (not shown). Therefore, a 3-day or 7-day minimum definition threshold would not change significantly the MHW characteristics in the future (see Sect. 4.)

The discrete nature of MHWs also necessitates a well-defined starting and ending day, but gaps with temperatures close to threshold values can also be found, as a result of day-to-day SST fluctuations. At this point, our definition differs slightly from that of Hobday et al. (2016). More specifically, gaps of up to 4 consecutive days or less are allowed inside a local MHW (considered as warm days). This is true, however, only when both the preceding and following 6-day mean SST of a gap day (including the gap day in each mean) are above the local $SST_{99Q}$. For the cool day “neighbourhood” this would represent a tendency to remain above threshold, even though the SST of that particular cool day might be below limit. This also reflects the fact that minor SST deviations from the threshold cannot impact the overall warm conditions of a MHW. It would most likely not offer either an “essential” relief to organisms, even to the less mobile and perhaps less tolerant species, once a MHW has started. Taking advantage of the default statistical sensitivity of the mean to outliers (in this case cold temperatures), we make the assumption that an event with the potential to interrupt a MHW (e.g wind, current) should cause a considerable drop in daily SST. Therefore, a below-threshold drop in either of the 6-day SST averages would not allow any cool day to merge with a MHW, in the same way that a sequence of five cool days or more would interrupt an event entirely. The 11-day window around the gap day is chosen since the minimum duration of a MHW was set to five days.

The spatial coverage of the MHW is then determined by aggregating grid points that are “activated” in a MHW state every day but are not necessarily contiguous. In common with many atmospheric definitions, a minimum 20% of the Mediterranean surface in $km^2$ was chosen in order to detect large-scale events that may have a broad ecosystem impact but also represent rare occurrence for the average climate conditions of HIST period. We, therefore, opt for prolonged, large-scale and extremely warm ocean temperatures that do not occur on a yearly basis in the 20th century,
with a view of quantifying their evolution in the 21st century under different GHG emission scenarios. The advantage of a percentile-based SST threshold in our case is that spatial patterns are also identified independently from the different extreme temperature levels that characterise sub-basins in the Mediterranean. We acknowledge that the detection method is developed based on subjective choices, and the sensitivity of the climate change results to these changes was also tested (See Sect. 4).

Once a MHW is identified, a subset of MHW metrics defined in Hobday et al. (2016) are used to characterise it. We examine the frequency of MHWs (Annual count of events), and the duration of each event is defined as the time between the first ($t_s$) and last day ($t_e$) for which a minimum of 20% of Mediterranean Sea surface is touched by a MHW. Every event is characterised by a mean and max intensity (mean and spatiotemporal maximum temperature anomaly relative to the threshold over the event duration) and a maximum surface coverage. Finally, its severity is represented by cumulative intensity (spatiotemporal sum of

<figure><img src="image_4.png" /><figcaption>**Fig. 1** Schematic of a MHW based on Hobday et al. (2016). The black line represents daily SST variations of one grid point in a random year and red line is the local threshold ($SST_{99Q}$) based on the 30-year average of yearly 99th quantile of daily SST for that point. The blue line is the daily 30-year climatology for this point. Also shown here also are the starting day ($t_s$) and ending day ($t_e$) above $SST_{99Q}$, gap days and the different measures of daily intensity. MHW metrics refer to the total event duration</figcaption></figure>

<table><caption>**Table 2** Marine heatwave (MHW) set of properties and their description after Hobday et al. (2016)</caption><thead><tr><th>Marine Heatwave Metrics</th><th>Description</th></tr></thead><tbody><tr><td>Frequency</td><td>Number of events occurring per year</td></tr><tr><td>Duration</td><td>(t<sub>e</sub> - t<sub>s</sub>) + 1 (days)</td></tr><tr><td>Mean intensity (I<sub>mean</sub>)</td><td>∫ [∫ (SST<sub>(x,y,t)</sub> - SST<sub>99Q</sub>(x, y))dxdy]dt / ∫ dxdy ∫ dt (°C)</td></tr><tr><td>Max intensity (I<sub>max</sub>)</td><td>max<sub>(x,y,t)</sub>(SST<sub>(x,y,t)</sub> - SST<sub>99Q</sub>(x, y)) (°C)</td></tr><tr><td>Severity (I<sub>cum</sub>)</td><td>∫ [∫ (SST<sub>(x,y,t)</sub> - SST<sub>99Q</sub>(x, y))dxdy]dt (°C days km<sup>2</sup>)</td></tr></tbody></table>

daily temperature anomalies relative to the threshold over the event duration) (Fig. 1, Table 2)

# 3 Results

## 3.1 Model evaluation

The first goal of this paper is to evaluate each models' ability to simulate mean ($SST$) and extreme Mediterranean Sea SST ($SST_{99Q}$) correctly. For this reason, pattern correlations were first performed, using the Pearson product-moment coefficient of linear correlation between two variables. The observed annual mean $SST$ between 1982–2012 (Fig. 2, OBS) shows a NW-SE pattern of cold-warm temperatures ranging from ~ 15 °C to 23 °C, respectively. Similarly, all the models demonstrate a warmer Eastern Mediterranean (EM) between 19 and 23 °C while colder deep water formation areas (e.g. Gulf of Lions, Adriatic) are captured well around 15–17 °C. Despite a multi-model mean (MMM) cold bias of about 0.6 °C, spatial correlations between each model mean 1976-2005 $SST$ and observations are high (MMM ~ 0.94). The lowest bias is found in ENEA and the highest in the CMCC and AWI/GERICS models (see Table 3). Note that satellite provides skin and night-time SST values, whereas the model SST represents averaged daily temperatures of the first few meters of mixed layer depth. Part of the model bias can be therefore explained by this difference in SST.

More complex spatial patterns are revealed when examining the 2D threshold maps used as the basis for defining Mediterranean MHWs (Fig. 3). The highest $SST_{99Q}$ are observed in Central Ionian, Gulf of Gabes, Tyrrhenian Sea and Levantine basin varying from approximately 27–31 °C and the lowest (20–22 °C) in deep water formation areas and the Alboran Sea (Fig. 3 OBS). In general, all the models are able to reproduce these patterns, although this time they share lower spatial correlations with the observations (MMM ~ 0.78). The ENEA model shows a warm bias whereas CMCC, U.BELGRAD and AWI/GERICS show a cold bias larger than 1 °C. The similar behaviour of the latter three could perhaps be related to the common atmospheric component (ECHAM) of their driving GCM. On the whole, the difference between the MMM mean and the extreme
<figure><img src="image_5.png" /><figcaption>**Fig. 2** Yearly SST (°C) for the HIST run of every model (1976–2005) and satellite data during 1982–2012. Note that the HIST run for ENEA is from 1979–2005</figcaption></figure>

basinwide SST is found to be ~ 6.6 °C, in good agreement with the observations (7.1 °C). The corresponding MMM spread is small for the SST but higher for the SST$_{99Q}$ (~ 1 °C).

Further, the domain-averaged timeseries of SST illustrate a warming tendency of both the SST and annual SST$_{99Q}$ (Fig. 4; Table 3). Even though MMM SST and SST$_{99Q}$ obtain similar trends (~ 0.02 °C/year) they seem to underestimate the corresponding observed trends (0.04/0.05 °C/year). Particularly for AWI/GERICS and U.BELGRAD, this is a response likely explained by their common driving GCM (MPI-ESM-LR). On the other hand, the amplitude of inter-annual variability is found similar to the observations for most of the models. On the whole, the observed and most of the model trends are statistically significant at a level of

95% except for certain cases indicated in Table 3. Interestingly, none of the simulations peaked as high as the observations during the exceptional MHW year of 2003 (20.4 °C for SST and 28.4 °C for SST$_{99Q}$). This record basinwide SST$_{99Q}$ value is on average 8.7 °C higher than the average SST of 1982–2012 and 2.8 °C greater than the basin-mean SST$_{99Q}$ of that period.

In terms of MHW properties during 1982–2012 (Table 3), observed MHW frequency is found at 0.8 events per year that last a maximum of 1.5 months and range between July and September. The mean intensity of MHWs varies from 0.3 to 0.9 °C, covering a maximum of 20–90% of the Mediterranean Sea surface, with a maximum intensity of 5.0 °C (2002) and a maximum severity of 8.5 × 10⁷ °C days km². The highest values over this period (except from I$_{max}$) refer
<figure><img src="image_6.png" /><figcaption>**Fig. 3** Individual MHW threshold maps of mean SST<sub>99Q</sub> (°C) computed from the HIST run of every model (1976–2005) and satellite data during 1982–2012. Note that the HIST run for ENEA is from 1979–2005</figcaption></figure>

to the characteristics of the well-known MHW 2003. More specifically, they correspond to a Mediterranean-scale event lasting 48 days (20 July–5 September) by our definition, in line with Grazzini and Viterbo (2003) and Sparnocchia et al. (2006). It seems that mainly the phase of the MHW that was both large-scale and intense was captured here.

On average, the simulated events during HIST are well within the equivalent observed range of every variable. They manifest though a slightly lower annual frequency, a potential for slightly higher maximum durations and starting dates up to early September. They also appear to underestimate the upper level of the I<sub>mean</sub>, I<sub>max</sub> and severity range. In particular, event durations of two months or more are exhibited by LMD, CMCC and AWI/GERICS models, while the ENEA model shows the highest I<sub>max</sub> of 5.3 °C. Maximum severity,

on the other hand, appears closer to the observed values only in the LMD and CMCC models. These configurations also show a MHW maximum spatial coverage above 80%, along with CNRM and AWI/GERICS. In general, the Med-CORDEX ensemble appears to perform well given that this is the first time, to our knowledge, that Mediterranean RCSMs have been evaluated for MHWs properties.

To better understand the ensemble variability of the MHW characteristics in the HIST period, we also combine Intensity-Duration-Frequency (IDF) information for every dataset separately (see Fig. 5). The total number of events of this period are organised in bins of I<sub>mean</sub> (every 0.02 °C) and duration (5-day bins progressively increased to 10-day and 20-day bins). Although some models simulate longer events relative to the observed MHW 2003,

<figure><img src="image_7.png" /></figure>
<figure><img src="image_8.png" /><figcaption>**Fig. 4** Timeseries of area-averaged, yearly $\overline{SST}$ °C (left) and SST$_{99Q}$ °C (right), during HIST for every model and satellite data, represented by a solid line. Trends are indicated in dashed lines. The different simulations are represented by different colors</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Fig. 5** IDF plot; Intensity ($I_{mean}$ in °C), Duration (Days), Frequency (Number of MHW during 1976–2005). $I_{mean}$ is organised in bins of 0.02 °C while duration is in bins of 5, 10 and 20 days. Red box indicates observed characteristics corresponding to the exceptional MHW of 2003</figcaption></figure>

only CMCC exhibited equivalent MHWs in terms of duration and intensity. At the same time, 1–3 events are detected for most classes of $I_{mean}$ and duration, in both the observations and the models. There are only a few cases where 3–7 events appear with $I_{mean}$ below 0.6 °C but follow no specific duration pattern.

## 3.2 Future Mediterranean SST evolution

In this section we analyze projections of $\overline{SST}$ and SST$_{99Q}$ in the 21st century by comparing their evolution against the reference period and under different GHG emission scenarios.
<table><caption>**Table 3** Evaluation of SST and MHW properties during HIST run</caption><thead><tr><th>Characteristics</th><th>CNRM</th><th>LMD</th><th>CMCC</th><th>AWI/GERICS</th><th>U.BELGRAD</th><th>ENEA</th><th>Multi-Model</th><th>OBS</th></tr></thead><tbody><tr><td>SST evaluation</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>SST</td><td>19.11</td><td>19.16</td><td>18.80</td><td>18.69</td><td>19.16</td><td>19.41</td><td>19.06 ± 0.30</td><td>19.70</td></tr><tr><td>SST - OBS Corr.Coeff</td><td>0.97</td><td>0.97</td><td>0.95</td><td>0.97</td><td>0.92</td><td>0.89</td><td>0.94 ± 0.03</td><td></td></tr><tr><td>SST -OBS Bias</td><td>-0.59</td><td>-0.54</td><td>-0.90</td><td>-1.01</td><td>-0.54</td><td>-0.29</td><td>-0.63 ± 0.30</td><td></td></tr><tr><td>SST Timeseries Trend</td><td>0.02</td><td>0.04</td><td>0.02</td><td>0.01*</td><td>0.01</td><td>0.01*</td><td>0.02 ± 0.01</td><td>0.04</td></tr><tr><td>SST<sub>99Q</sub></td><td>26.23</td><td>26.03</td><td>24.52</td><td>24.43</td><td>25.68</td><td>26.97</td><td>25.64 ± 1.00</td><td>26.79</td></tr><tr><td>SST<sub>99Q</sub>-OBS Corr虚拟</td><td>0.88</td><td>0.84</td><td>0.77</td><td>0.84</td><td>0.70</td><td>0.74</td><td>0.79 ± 0.07</td><td></td></tr><tr><td>SST<sub>99Q</sub>-OBS Bias</td><td>-0.56</td><td>-0.76</td><td>-2.27</td><td>-2.36</td><td>-1.11</td><td>0.18</td><td>-1.15 ± 1.00</td><td></td></tr><tr><td>SST<sub>99Q</sub> Timeseries trend</td><td>0.03</td><td>0.04</td><td>0.04</td><td>0.01*</td><td>0.00*</td><td>0.02*</td><td>0.02 ± 0.02</td><td>0.05</td></tr><tr><td>MHW characteristics (HIST)</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>Frequency</td><td>0.7</td><td>0.6</td><td>0.6</td><td>1</td><td>0.7</td><td>0.7</td><td>0.7</td><td>0.8</td></tr><tr><td>Duration</td><td>2-32</td><td>4-67</td><td>1-74</td><td>1-61</td><td>4-38</td><td>4-44</td><td>2.6-52.6</td><td>1-48</td></tr><tr><td>Starting Day</td><td>Jul-Sep</td><td>Jul-Aug</td><td>Jul-Aug</td><td>Jul-Sep</td><td>Jul-Aug</td><td>Jul-Sep</td><td>Jul-Sep</td><td>Jul-Aug</td></tr><tr><td>Ending Day</td><td>Jul-Sep</td><td>Aug-Sep</td><td>Aug-Sep</td><td>Aug-Sep</td><td>Jul-Sep</td><td>Aug-Sep</td><td>Aug-Sep</td><td>Aug-Sep</td></tr><tr><td>I<sub>mean</sub></td><td>0.4-0.8</td><td>0.3-0.7</td><td>0.3-0.9</td><td>0.2-0.6</td><td>0.3-0.6</td><td>0.4-0.6</td><td>0.3-0.7</td><td>0.3-0.9</td></tr><tr><td>I<sub>max</sub></td><td>1.8-5.5</td><td>1.5-4.1</td><td>1.4-4.9</td><td>1.0-3.0</td><td>1.5-3.2</td><td>1.8-5.3</td><td>1.5-4.3</td><td>1.3-5.0</td></tr><tr><td>Severity (Icum)</td><td>0.04-3.6</td><td>0.1-7.8</td><td>0.01-8.8</td><td>0.01-4.3</td><td>0.1-1.6</td><td>0.1-2.6</td><td>0.05-4.8</td><td>0.02-8.5</td></tr><tr><td>Max Surface</td><td>20.8-82.5</td><td>21.5-88.3</td><td>20.1-90.8</td><td>20.3-81.9</td><td>22.5-59.8</td><td>20.7-72.1</td><td>21.0-79.2</td><td>20.1-90.1</td></tr></tbody></table>

Mean annual and threshold SST are indicated with SST (°C) and SST<sub>99Q</sub> (°C) respectively. The Mann-Kendall non-parametric test is used to detect the presence of linear or non-linear monotonic trends (°C/year) in domain-averaged SST timeseries. Trends with statistical significance lower than 95% level are indicated with star. Spatial correlations (Corr.Coeff) and bias with respect to observations are given for each dataset. Also shown here, are the range (min and max) of frequency, duration (days), starting day (calendar month), ending day (calendar month), I<sub>mean</sub> (°C) I<sub>max</sub> (°C), Severity*10<sup>7</sup> (°C days km<sup>2</sup>) and maximum surface coverage(%) of MHWs. The multi-model column indicates the ensemble average values and standard deviation for each variable

During 2021–2050 an increase is found for the domain-averaged ensemble mean SST and SST<sub>99Q</sub> with respect to HIST, around 0.8–1 °C and 1–1.2 °C respectively. While

the mid-21st century anomalies appear almost independent from the greenhouse gas forcing, a more diverse and substantial warming occurs towards 2071–2100 (see Fig. 6).

<figure><img src="image_10.png" /><figcaption>**Fig. 6** Area-average, yearly SST °C (left) and extreme SST<sub>99Q</sub> °C (right) anomalies with respect to HIST. Bold colors represent the</figcaption></figure>

multi-model average and lighter colors are the individual simulations. RCP2.6 scenario has only one simulation (CNRM), HIST run is illustrated in grey and observations in dashed black

<figure><img src="image_11.png" /></figure>
In particular, the multi-model mean $\overline{SST}$ and SST$_{99Q}$ anomalies under RCP8.5 are 3.1 °C and 3.6 °C respectively, exhibiting nearly a doubling of their corresponding RCP4.5 rise. Similarly, the equivalent increase of CNRM $\overline{SST}$ and SST$_{99Q}$ under RCP8.5 is about 3 times as high as that under RCP2.6 for the same period. Individually, however, the highest mean and extreme SST anomalies are demonstrated by the LMD and CMCC models under every scenario and for every period (see Table 4). For both SST indices, the effects of the different emission scenarios become more evident by 2060, with the highest/ intermediate warming occurring for every model under RCP8.5/RCP4.5 and the lowest under the (mono-model) RCP2.6 simulation. In the latter, little or no difference is found between the $\overline{SST}$ and SST$_{99Q}$ rise throughout the century. In contrast, under RCP4.5 and RCP8.5, the multi-model SST$_{99Q}$ increase appears greater than the $\overline{SST}$ rise by 20–25% during 2021–2050 and by 16–18% for 2071–2100 (see Table 4 and discussion section). This implies a higher contribution from $\overline{SST}$ to the warming towards the end of the century.

The spatial distribution of the corresponding anomalies, however, appears inhomogeneous. For 2071–2100, some regions in the Levantine basin, Balearic islands, Tyrrenian Sea, Ionian Sea and North Adriatic Sea exhibit the highest MMM $\overline{SST}$ anomalies in every scenario (Fig. 7). In contrast, the lowest anomalies of that time are located in the Alboran Sea, where cold waters are advected from the

Atlantic, and depending on the scenario they may range from ~ 0.6 °C (RCP2.6) to ~ 2.4 °C (RCP8.5).

Meanwhile, the most pronounced extreme warm anomalies (SST$_{99Q}$) for 2071–2100 under RCP4.5 and RCP8.5 are projected for the NW mediterranean, Tyrrenian Sea, Ionian Sea and some parts of North Levantine basin (Fig. 8). Under RCP2.6 though, the greatest SST$_{99Q}$ anomalies (> 1.2 °C) are more confined towards the Aegean Sea, Adriatic, Tyrrenian Sea and the area around Balearic islands. In addition to the highest SST$_{99Q}$ rise, the Adriatic Sea, Ionian Sea, Tyrrenian Sea, some parts around the Balearic islands and the North Levantine basin display also the greatest $\overline{SST}$ rise, for every scenario during the 2nd half of the 21st century. During 2021–2050, however, they exhibit the highest mean and extreme warming under RCP26 and RCP4.5 but not under RCP8.5. The Alboran Sea and the SE Levantine basin, on the other hand, demonstrate the lowest SST$_{99Q}$ anomalies in every period and every scenario

### 3.3 Future evolution of Mediterranean MHWs

The MHW climate change response is examined here using anomalies. These anomalies are computed for the average MHW characteristics in the future relative to the average MHW characteristics in HIST run, for each sub-period, model and scenario (Table 5).

The multi-model mean reveals an increase in frequency of 0.3–0.4 events/year for every period of RCP4.5/RCP8.5 with the mono-model RCP2.6 simulation showing a

<table><caption>**Table 4** Future Mediterranean-averaged, yearly mean ($\overline{SST}$) and extreme (SST$_{99Q}$) anomalies (with respect to HIST) for the near and far future under different emission scenarios</caption><thead><tr><th></th><th>CNRM</th><th>LMD</th><th>CMCC</th><th>AWI/GERICS</th><th>BELGRAD</th><th>ENEA</th><th>Multi-model</th></tr></thead><tbody><tr><td>RCP85 2021–2050</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><span style="text-decoration: overline;">SST</span></td><td>0.9</td><td>1.3</td><td>1.2</td><td>0.7</td><td>0.7</td><td>—</td><td>1.0 ± 0.3</td></tr><tr><td>SST<sub>99Q</sub></td><td>1.1</td><td>1.7</td><td>1.3</td><td>1.0</td><td>1.0</td><td>—</td><td>1.2 ± 0.3</td></tr><tr><td>RCP85 2071–2100</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><span style="text-decoration: overline;">SST</span></td><td>2.7</td><td>3.8</td><td>3.4</td><td>2.7</td><td>2.7</td><td>—</td><td>3.1 ± 0.5</td></tr><tr><td>SST<sub>99Q</sub></td><td>2.9</td><td>4.5</td><td>4.3</td><td>3.1</td><td>3.1</td><td>—</td><td>3.6 ± 0.7</td></tr><tr><td>RCP45 2021–2050</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><span style="text-decoration: overline;">SST</span></td><td>0.7</td><td>1.2</td><td>1.0</td><td>0.6</td><td>—</td><td>0.6</td><td>0.8 ± 0.4</td></tr><tr><td>SST<sub>99Q</sub></td><td>0.8</td><td>1.4</td><td>1.3</td><td>0.8</td><td>—</td><td>0.7</td><td>1.0 ± 0.5</td></tr><tr><td>RCP45 2071–2100</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><span style="text-decoration: overline;">SST</span></td><td>1.6</td><td>2.1</td><td>2.0</td><td>1.1</td><td>—</td><td>1.2</td><td>1.6 ± 0.8</td></tr><tr><td>SST<sub>99Q</sub></td><td>1.8</td><td>2.6</td><td>2.5</td><td>1.4</td><td>—</td><td>1.3</td><td>1.9 ± 0.9</td></tr><tr><td>RCP26 2021–2050</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><span style="text-decoration: overline;">SST</span></td><td>0.8</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>SST<sub>99Q</sub></td><td>1.0</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>RCP26 2071–2100</td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td><span style="text-decoration: overline;">SST</span></td><td>1.0</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>SST<sub>99Q</sub></td><td>1.0</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

The multi-model column indicates the ensemble average values and standard deviation. Values are in °C
<table><caption>**Table 5** Future response (anomalies with respect to HIST) of MHW mean properties for the 6 RCSMs under RCP8.5, RCP4.5 and RCP2.6, for the near (2021-2050) and far future (2071-2100).</caption><thead><tr><th>Variables</th><th>CNRM</th><th>LMD</th><th>CMCC</th><th>AWI/GERICS</th><th>BELGRAD</th><th>ENEA</th><th>Multi-model</th></tr></thead><tbody><tr><td colspan="8">Marine heatwave characteristics</td></tr><tr><td colspan="8">RCP8.5 (2021–2050)</td></tr><tr><td>Frequency</td><td>0.5</td><td>0.4</td><td>0.5</td><td>0.1</td><td>0.4</td><td>—</td><td>0.4 ± 0.2</td></tr><tr><td>Duration</td><td>32.7</td><td>47.8</td><td>39.3</td><td>39.1</td><td>37.4</td><td>—</td><td>39.2 ± 16.7</td></tr><tr><td>Starting day</td><td>Jul</td><td>Jul</td><td>Jul</td><td>Aug</td><td>Jul</td><td>—</td><td>Jul</td></tr><tr><td>Ending day</td><td>Sep</td><td>Sep</td><td>Sep</td><td>Sep</td><td>Sep</td><td>—</td><td>Sep</td></tr><tr><td>I<sub>mean</sub></td><td>0.2</td><td>0.4</td><td>0.3</td><td>0.2</td><td>0.2</td><td>—</td><td>0.3 ± 0.1</td></tr><tr><td>I<sub>max</sub></td><td>1</td><td>1.4</td><td>1.3</td><td>1.2</td><td>0.9</td><td>—</td><td>1.2 ± 0.5</td></tr><tr><td>Max surface</td><td>23.5</td><td>39.2</td><td>29.4</td><td>36.7</td><td>35.6</td><td>—</td><td>32.9 ± 14.6</td></tr><tr><td>Severity (Icum)</td><td>3.9</td><td>11.6</td><td>7.1</td><td>6.7</td><td>4.7</td><td>—</td><td>6.8 ± 3.8</td></tr><tr><td colspan="8">RCP8.5 (2071–2100)</td></tr><tr><td>Frequency</td><td>0.3</td><td>0.4</td><td>0.5</td><td>0</td><td>0.4</td><td>—</td><td>0.3 ± 0.2</td></tr><tr><td>Duration</td><td>83.8</td><td>105.9</td><td>97.7</td><td>99.6</td><td>83.7</td><td>—</td><td>94.1 ± 9.9</td></tr><tr><td>Starting Day</td><td>Jun</td><td>Jun</td><td>Jun</td><td>Jun</td><td>Jun</td><td>—</td><td>Jun</td></tr><tr><td>Ending Day</td><td>Oct</td><td>Oct</td><td>Oct</td><td>Oct</td><td>Oct</td><td>—</td><td>Oct</td></tr><tr><td>I<sub>mean</sub></td><td>0.9</td><td>1.9</td><td>1.6</td><td>1.3</td><td>1.1</td><td>—</td><td>1.4 ± 0.4</td></tr><tr><td>I<sub>max</sub></td><td>3.1</td><td>4.2</td><td>4.4</td><td>3.8</td><td>3.1</td><td>—</td><td>3.7 ± 0.6</td></tr><tr><td>Max Surface</td><td>51.7</td><td>50.7</td><td>46.1</td><td>58.9</td><td>53.5</td><td>—</td><td>52.2 ± 4.6</td></tr><tr><td>Severity (Icum)</td><td>29.2</td><td>73.1</td><td>63.1</td><td>46</td><td>34.1</td><td>—</td><td>49.1 ± 18.7</td></tr><tr><td colspan="8">RCP4.5 (2021–2050)</td></tr><tr><td>Frequency</td><td>0.6</td><td>0.4</td><td>0.5</td><td>0.1</td><td>—</td><td>0.2</td><td>0.4 ± 0.2</td></tr><tr><td>Duration</td><td>17.7</td><td>43.7</td><td>38.6</td><td>33.9</td><td>—</td><td>20.6</td><td>30.9 ± 16.2</td></tr><tr><td>Starting day</td><td>Aug</td><td>Jul</td><td>Jul</td><td>Jul</td><td>—</td><td>Jul</td><td>Jul</td></tr><tr><td>Ending day</td><td>Sep</td><td>Sep</td><td>Sep</td><td>Sep</td><td>—</td><td>Sep</td><td>Sep</td></tr><tr><td>I<sub>mean</sub></td><td>0.1</td><td>0.3</td><td>0.2</td><td>0.2</td><td>—</td><td>0.2</td><td>0.2 ± 0.1</td></tr><tr><td>I<sub>max</sub></td><td>0.5</td><td>1.2</td><td>1.2</td><td>1.1</td><td>—</td><td>0.9</td><td>1.0 ± 0.5</td></tr><tr><td>Max surface</td><td>13.5</td><td>34.1</td><td>27.9</td><td>33.8</td><td>—</td><td>27.9</td><td>27.4 ± 13.5</td></tr><tr><td>Severity (Icum)</td><td>1.9</td><td>8</td><td>6.7</td><td>4.9</td><td>—</td><td>3.7</td><td>5.0 ± 3.0</td></tr><tr><td colspan="8">RCP4.5 (2071–2100)</td></tr><tr><td>Frequency</td><td>0.4</td><td>0.4</td><td>0.5</td><td>0</td><td>—</td><td>0.2</td><td>0.3 ± 0.2</td></tr><tr><td>Duration</td><td>56.4</td><td>69.8</td><td>67.9</td><td>55.4</td><td>—</td><td>45.4</td><td>59 ± 10.0</td></tr><tr><td>Starting day</td><td>Jul</td><td>Jun</td><td>Jul</td><td>Jul</td><td>—</td><td>Jul</td><td>Jul</td></tr><tr><td>Ending day</td><td>Sep</td><td>Oct</td><td>Oct</td><td>Sep</td><td>—</td><td>Sep</td><td>Sep</td></tr><tr><td>I<sub>mean</sub></td><td>0.4</td><td>0.8</td><td>0.7</td><td>0.4</td><td>—</td><td>0.4</td><td>0.5 ± 0.2</td></tr><tr><td>I<sub>max</sub></td><td>1.6</td><td>2.3</td><td>2.3</td><td>1.8</td><td>—</td><td>1.5</td><td>1.9 ± 0.4</td></tr><tr><td>Max surface</td><td>44.4</td><td>48.5</td><td>42.6</td><td>51.2</td><td>—</td><td>39</td><td>45.1 ± 4.8</td></tr><tr><td>Severity (Icum)</td><td>10.6</td><td>25.3</td><td>23.9</td><td>11.8</td><td>—</td><td>8.5</td><td>16.0 ± 7.7</td></tr><tr><td colspan="8">RCP2.6 (2021–2050)</td></tr><tr><td>Frequency</td><td>0.7</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Duration</td><td>17.2</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Starting day</td><td>Jul</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Ending day</td><td>Sep</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>I<sub>mean</sub></td><td>0.1</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>I<sub>max</sub></td><td>0.5</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Max surface</td><td>15</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Severity (Icum)</td><td>2.4</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td colspan="8">RCP2.6 (2071–2100)</td></tr><tr><td>Frequency</td><td>0.5</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Duration</td><td>30.5</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Starting day</td><td>Jul</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Ending day</td><td>Sep</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

<figure><img src="image_12.png" /></figure>
<table><caption>**Table 5** (continued)</caption><thead><tr><th>Variables</th><th>CNRM</th><th>LMD</th><th>CMCC</th><th>AWI/GERICS</th><th>BELGRAD</th><th>ENEA</th><th>Multi-model</th></tr></thead><tbody><tr><td>I<sub>mean</sub></td><td>0.2</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>I<sub>max</sub></td><td>0.9</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Max surface</td><td>22.7</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr><tr><td>Severity (Icum)</td><td>3.6</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td><td>—</td></tr></tbody></table>

Shown here are the average annual event count (frequency), average MHW duration (days), starting day (calendar month), ending day (calendar month), I<sub>mean</sub> (°C), I<sub>max</sub> (°C), severity (* 10<sup>7</sup> °C days km<sup>2</sup>) and maximum surface coverage (%). The multi-model column indicates the ensemble mean values for each variable and their standard deviation. Only the CNRM simulation is available for the RCP2.6 scenario, 5 simulations for RCP8.5 and 5 simulations for RCP4.5

<figure><img src="image_13.png" /><figcaption>**Fig. 7** Multi-model average anomaly of yearly SST (°C) with respect to the corresponding ensemble mean HIST of each scenario, for the near and far future. The RCP2.6 scenario has only one simulation (CNRM)</figcaption></figure>

slightly greater increase of 0.5–0.7 events/year. In fact, individual simulations of Fig. 9 suggest a shift from a period where years without MHWs were common (1976–2030) to a period with at least one long-lasting MHW every year. More specifically, towards 2071–2100, events can start as early as June and finish as late as October under RCP8.5, whereas for RCP4.5 and RCP2.6, the MHW temporal extent appears between July–September (Fig. 10). It is clear that the higher the radiative forcing, the broader the window of occurrence. For example,

MHWs during 2071–2100 may last on average 3 months longer in RCP8.5 than HIST MHWs (~21.8 days, not shown) but almost 2 months longer in RCP4.5 (see Table 5). This is a MMM increase in the duration, almost double the corresponding increase during 2021–2050 under RCP4.5 (~30.9 days) and more than double that under RCP8.5 (~39.2 days). Even under the optimistic RCP2.6 scenario, MHWs by 2050 may be 17.2 days longer than today and may become 1 month longer at maximum by 2100.
<figure><img src="image_14.png" /><figcaption>**Fig. 8** Multi-model average anomaly of extreme SST<sub>99Q</sub> (°C) with respect to corresponding ensemble mean HIST (1976-2005) of each scenario, for the near and far future. The RCP2.6 scenario has only one simulation (CNRM)</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>**Fig. 9** Annual number of MHWs (Annual Frequency) for RCP8.5 (red) RCP4.5 (blue) RCP2.6 (green) HIST (grey) and observations (dashed black). Bold colors indicate the multi-model mean and shaded zones represent individual MHW events identified by the models. Years without MHWs are also included, with shaded areas reaching 0. RCP2.6 has only 1 simulation (CNRM)</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>**Fig. 10** Annual earliest starting (solid lines) and latest ending (dashed lines) day of MHW events for RCP8.5 (red) RCP4.5 (blue) RCP2.6 (green) HIST (grey) and observations (black). Bold colors indicate multi-model average values while lighter dots represent individual event dates</figcaption></figure>

Long-term projections show analogous changes in the $I_{\text{mean}}$ of future MHWs. They are examined through IDF plots that display the total number of MHWs identified by the ensemble during HIST (1976–2005) run, near and far future (Fig. 11). To avoid imbalances in the present-future comparisons arising from the different sets of models for

<figure><img src="image_17.png" /></figure>
<figure><img src="image_18.png" /><figcaption>**Fig. 11** IDF (I<sub>mean</sub>, Duration, Frequency) plots display the total MHW number of every dataset, for every scenario, over 2021-2050 and 2071-2100. RCP8.5 and RCP4.5 include events from 5 simulations, while RCP2.6 from only 1 (CNRM) simulation. HIST run contains</figcaption></figure>

MHWs from the corresponding set of models each time. The number of MHWs is calculated over each 30 year period. For contrast purposes, the red box depicts the observed characteristics of MHW 2003 in the Mediterranean

RCP4.5 and RCP8.5 (see Table 5), all the simulated future events are pooled for every period and juxtaposed against the corresponding sets of HIST events. Therefore, we show 3 HIST IDF plots, one for each scenario. As previously demonstrated the stronger the emission scenario, the longer the duration and the higher the I<sub>mean</sub> of the events. The MMM I<sub>mean</sub> response appears small during 2021–2050 (+ 0.1°C to + 0.3 °C depending on the scenario) but increases towards

the end of the period with higher radiative forcing. For instance, MHWs show durations of up to 170 days (Fig. 11) in the far future of RCP8.5 and I<sub>mean</sub> of 1.8 °C on average (not shown). For the CNRM model though and under RCP2.6, the corresponding response towards 2071–2100 has doubled compared to the mid-21st century, while it becomes 4.5 times higher under RCP8.5 (see Table 5). Longer-lasting MHWs at the end of the period for RCP4.5 and RCP8.5
explain the lower frequency of occurrences compared to RCP2.6. A similar behaviour to $I_{\text{mean}}$ displays the MMM average response of $I_{\text{max}}$, with the highest anomalies indicated towards 2071–2100 (up to 3.7 °C), whereas during the mid-21st century they range between 0.5 °C (RCP2.6) and 1.2 °C(RCP8.5).

It should be also noted that for RCP4.5 and RCP8.5, events with characteristics similar to the observed exceptional MHW 2003 (Fig. 11, red box) seem to become the new standard over 2021–2050 and even constitute weak occurrences for the distant future of RCP8.5. In the more optimistic RCP2.6, MHWs appear more frequent during 2021–2050 but their number is slightly decreased towards 2071–2100. Their characteristics, however, sustain a lower increase throughout the period compared to RCP4.5 and RCP8.5. For example, the response in duration and $I_{\text{mean}}$ is found close to that projected for CNRM during 2021–2050 and under RCP4.5 and RCP8.5 (see Table 5). Therefore, the possibility for an event like the MHW 2003 to occur regularly still features in a scenario close to the Paris Agreement (RCP2.6).

Yet, the range of the uncertainty in future projections evolves not only in time but also throughout the different models. The severity (Icum) distribution of future MHWs was determined in that sense using Whisker diagrams. In these box plots, a specific Icum index is appointed at each simulated event of every dataset for each period and scenario (see Fig. 12, left). By definition, Icum translates the total spatiotemporal MHW impact into numbers. It features an exponential increase from HIST towards the

end of the century from ~1 × 10⁷ C days km² to about ~50 × 10⁷ C days km² for RCP8.5 (2071–2100-not shown). Moreover, the higher the emission forcing, the higher the rate the ensemble mean Icum response escalates from its mid- to end-of-century values; for example, Icum varies from 5 to 33 × 10⁷ C days km² in RCP4.5 and from 6.8 to 49.1 × 10⁷ C days km² in RCP8.5 (see Table 5). This becomes more evident when comparing the equivalent CNRM severity response under RCP2.6 (2.4–3.6 × 10⁷ C days km²) with the significantly higher response under RCP4.5 and RCP8.5. Although all configurations indicate an abrupt escalation through time, there appears to be a family of models (CMCC and LMD) that share a stronger climate change response. Those models exhibit higher changes in Icum, along with higher $I_{\text{mean}}$, $I_{\text{max}}$, and duration values than the remaining models (see also Table 5 and Sect. 4).

The identified families of MHWs are also associated with a maximum spatial coverage illustrated through box plots in Fig. 12. It is estimated that events may affect a maximum of 40% of the Mediterranean Sea, on average, during HIST but may impact almost 100% of the basin by 2071–2100 under RCP8.5. Notwithstanding the large variability found for the mid-21st century, by 2100 the simulated maximum MHW extent seems to be an unanimous projection from every model and under RCP8.5. Conversely, MHWs under RCP2.6 increase their maximum coverage throughout the period, but towards 2071–2100 events may impact, on average, a maximum of 70% of the Mediterranean Sea.

<figure><img src="image_19.png" /><figcaption>**Fig. 12** Whisker diagram of (left) Severity (Icum) and (right) maximum surface coverage of every observed and simulated MHW during HIST, 2021-2050 and 2071-2100. Box plots illustrate minimum, 25th</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>percentile, median, 75th percentile and maximum values of each variable for a given model, scenario and period</figcaption></figure>
# 4 Discussion

## 4.1 MHW detection method

Several sensitivity tests performed on the MHW detection algorithm using only the CNRM model indicate low levels of uncertainty associated with small perturbations on the initial definition. For example, definitions with a different number of gap days, different minimum duration or minimum MHW spatial extent allowed (e.g. 10%) were tested but did not seem to change significantly the response of future MHW characteristics with respect to HIST run (see Supplementary Table 1S). The use of different quantile thresholds also showed that climate change response of duration, $I_{\text{mean}}$ and $I_{\text{max}}$ with respect to HIST does not differ significantly if a lower/higher threshold than the $SST_{99Q}$ is chosen. However, the severity and maximum spatial coverage appear more sensitive to such changes (see Supplementary Table 1S).

However, certain limitations exist: assuming no spatial connectivity, the detection algorithm provides identification of large-scale (>20%) and long-lasting events but does not consider MHW effects during colder months or spatially smaller events. While it describes surface MHWs in the summer, it can be also applied to deeper layers and/or winter season when availability of data allows it.

## 4.2 Model-observation discrepancies

The discrepancies on mean and extreme Mediterranean Sea temperatures with respect to observations on the models were also evaluated using a shorter but common reference period of 1982-2005 for both datasets. Values of $SST$, $SST_{99Q}$, their trends and pattern correlations did not change considerably. However, the multi-model mean bias was slightly reduced by 28% for $SST$ and by 31% for $SST_{99Q}$ (see Supplementary material Table 2). Moreover, MHW identification appeared consistent despite the different SST layer depth of the observations (~ mm) and the models (~ m).

## 4.3 Model uncertainty

By default, the estimate of the uncertainty is given by the variation of the results across the ensemble members in an opportunistic way (Knutti et al. 2010). Although the models we use have a high-resolution representation of the air-sea interactions, uncertainties are introduced due to their individual biases but also due to the small number of the currently available Med-CORDEX simulations (6 RCSMs). To this purpose, more runs will be added in the future as part of the Med-CORDEX initiative. Despite this limitation,

the RCSM ensemble seems to explore well the spread of SST anomalies predicted by earlier studies based on GCMs. For example, for RCP4.5 we estimate annual area-average $SST$ anomalies from 2006–2100 with respect to HIST from approximately 0.7 to 2.6 °C, depending on the model (Fig. 6, left). This covers a large part of the corresponding anomalies found by Mariotti et al. (2015) for 2006–2100 with respect to 1980–2005 mean, which were between 0.5 °C to 3.5 °C for the CMIP5 ensemble of GCMs under the RCP4.5 scenario. Although our ensemble appears to underestimate the upper limit of this CMIP5 range, this could also reflect a better representation of the Mediterranean Sea dynamics by the regional models. Indeed, at higher resolutions the representation of air-sea interaction also improves (e.g. Akhtar et al. 2018; Roberts et al. 2016; Hewitt et al. 2017). At the same time, our results indicate an intensification of MHWs in the Mediterranean Sea with time, in agreement with the results obtained by Oliver et al. (2018a) and Frölicher et al. (2018), which used different MHW definitions.

Albeit some models have demonstrated lower/higher biases than others, we have chosen not to discard any of the configurations since their weak performance in some indices is not related to any specific behaviour of MHW indices in scenario. This choice also favours the holistic presentation of the uncertainty spectrum, without a considerable impact on the climate change response. More specifically, closer examination of the $SST$ and $SST_{99Q}$ bias effect on the anomalies of the average MHW characteristics in RCP8.5 and RCP4.5 with respect to HIST suggested no particular tendency or outliers affecting the range of the outcome, for any of the periods and scenarios (Supplementary material Fig. 1S, Fig.2S). It is however notable that LMD and CMCC have a tendency to show stronger responses in MHW or SST values. This could be due to the driving GCMs (IPSL-CMA5-MR and CMCC-CM), which demonstrate a higher mean surface temperature change over Europe by 2080 compared to CNRM-CM5 and MPI-ESM-LR, according to McSweeney et al. (2015). In that study, the performance of all the GCMs driving the Med-CORDEX RCSMs was characterised as “satisfactory” for downscaling, except that of LMD, which was found with biases.

## 4.4 MHW evolution and changes in SST

Present-day extreme warming at the order of $SST_{99Q}$ might constitute a rare occurrence for the Mediterranean Sea climate; however, in the future it becomes the new normal. In 2071–2100 in particular, the warming signal is found so high that almost every day from June to October can experience such extreme temperatures. This means that future warming in the Mediterranean Sea is practically able to saturate what is considered today as a severe MHW. The difference between the scenarios lies in the fact that under RCP4.5 and
RCP8.5 anomalous temperatures appear more persistent and widespread and therefore fewer but longer and more intense events occur. By contrast, under RCP2.6, events appear less persistent, and therefore more “breaks” between MHWs may occur (frequency of events is increased), since a significant part of the basin is more likely to fall below the $SST_{99Q}$ threshold (Fig. 9).

Most of the future changes in MHW characteristics were seen to increase following the GHG forcing, yet this raises the question of whether this behavior could be explained by changes in the mean (shift of distribution) or the day-to-day SST variability (distribution flattening/narrowing). As a first indicator, we calculated (for the CNRM model only) the $\overline{SST}$ difference between RCP8.5 (2071-2100) and HIST and added it to the current $SST_{99Q}$ threshold map (see Supplementary Table 1). The resulting climate change response (future-present) of MHW characteristics was much lower than the one found when using the initial $SST_{99Q}$ threshold alone. This signifies that the mean SST change alone can

explain a large part, but not all, of the future changes in MHWs. We estimate that 10–20% of the MHW characteristics are due to changes in day-to-day variability.

To further test our hypothesis, we calculated the multi-model mean ratio $R = \Delta(SST_{99Q_{Scenario} - SST_{99Q_{Hist}}})/\Delta(\overline{SST}_{Scenario} - \overline{SST}_{Hist})$ for every scenario and period (see Fig. 13). In regions where $R > 1$ SST daily variability contributes to the extreme temperature increase and only where $R > 2$, it dominates the mean SST change contribution (distribution flattening). For $R < 1$, a narrowing of SST daily distribution lowers the mean SST signal, which makes the dominant contribution when close to $R=1$. Overall, model results indicate a higher contribution from SST daily variability change in the mid-21st century compared to 2071-2100, when $\overline{SST}$ change becomes more important (Table 4; Fig. 13). During 2021-2050 and for every scenario, basin-mean $R \sim 1.2$ with the Alboran Sea, some coastal parts of the Aegean Sea, Adriatic Sea and SE Levantine basin exhibiting a narrowing ($R < 1$) or a shift

<figure><img src="image_21.png" /><figcaption>**Fig. 13** Multi-model mean ratio $R$ of $\Delta SST_{99Q}$ (°C) over $\Delta \overline{SST}$ (°C) for every period and scenario. Regions where $R > 1/R < 1$ indicate</figcaption></figure>

regions where flattening/narrowing of SST distribution is detected in addition to the mean distribution shift. Where $R \sim 1$ the $\overline{SST}$ increase can be considered as the main factor for MHW changes
(R=1) in the SST distribution. Towards 2071-2100, however, and under RCP4.5 and RCP8.5, basin-average R is between 1 < R < 1.2 and more areas demonstrate a range closer to R = 1.

It is worth noting that narrowing of the SST distribution
is a rare situation that appears in small areas only under
RCP2.6 or around the Alboran Sea under RCP85 for the
mid-21st century. Flattening, on the other hand, appears
more common and could possibly reflect the increase in
day-to-day variability of 2 m-air temperature estimated
by Giorgi (2006) for the Mediterranean area in the future.
Finally, a slightly stronger influence of SST daily variability
is seen in the NW Mediterranean area for every scenario
and period (Fig. 13). The possible explanations for such a
spatial pattern could be a future mixed layer depth shoaling,
as projected by (Adloff et al. 2015). This would mean that
heat fluxes would be able to change the heat content of a
shallower MLD faster, creating that way an increase in daily
SST variability.

5 Conclusions

The main objective of this study is to investigate the future evolution (1976–2100) of SST and marine heatwaves in the Mediterranean Sea, using the best dedicated multi-model ensemble available. Here we examine six Regional Climate System Models from the Med-CORDEX initiative, driven by 4 CMIP5 GCMs under the RCP2.6, RCP4.5 and RCP8.5 scenarios. A quantitative MHW definition and detection method based on SST and on Hobday et al. (2016) approach is developed, targeting large-scale and long-lasting events, mostly in the warmer months. The algorithm uses a climatological 99th percentile threshold based on historical simulations (1976–2005) and takes into account a spatially-varying threshold. It delivers MHW metrics such as frequency, duration, mean and maximum intensity along with severity and maximum spatial coverage.

Spatiotemporal indices under a 1976–2005 (HIST) run
reveal that the Med-CORDEX ensemble simulates the
present MHW characteristics well, although it appears to
underestimate the warming trends of SST and SST<sub>99Q</sub> of that
period with respect to observations from 1982–2012. The
latter dataset yields an annual frequency of 0.8 events/year,
with MHWs lasting a maximum of 1.5 months between July
and September, while covering a maximum of 90% of the
Mediterranean Sea surface. The longest and most severe
event of that period corresponded to the MHW of 2003,
which also demonstrates the highest mean intensity and
maximum event coverage.

Analysis of future evolution shows that differences in
the GHG forcing are reflected mostly towards 2071–2100,
whereas uncertainty for the mid-21st century is dominated

by the model uncertainty. Ensemble means by the end of the
century demonstrate the highest *SST* (3.1 °C) and *SST*<sub>99Q</sub>
(3.6 °C) increase under RCP8.5 and lowest under RCP2.6
(mono-model). The corresponding warming for 2021-2050,
however, is less pronounced under RCP4.5 (~0.8°C/1°C)
and RCP8.5 (~1 °C/1.2 °C). In contrast, basinwide mean
and extreme SST for RCP2.6 (~1 °C) does not differ sig-
nificantly from mid- to end of 21st century.

By 2100, models project at least one long-lasting MHW
occurring every year under RCP8.5 up to 3 months longer,
and about 4 times more intense and 42 times more severe
than today's events. Their occurrence is expected between
June and October, affecting at peak, the entire Mediterranean
basin. In fact, with respect to the HIST run, MMM MHW
frequency increases by a factor of ~ 1.6 for RCP8.5 and
RCP4.5 by 2021–2050 and slightly less than that towards
2071–2100 for both scenario. The equivalent CNRM com-
parison between the scenarios reveals a slightly greater fre-
quency increase during 2071–2100 under RCP2.6 (by fac-
tor of 1.7) than under RCP8.5 and RCP4.5. Multi-model
mean duration, on the other hand, is multiplied by a factor
of 3.7 for RCP4.5 and 5.3 for RCP8.5 during 2071–2100.
MHWs under RCP8.5 may also have an I<sub>mean</sub> 3.9 times as
high as today's event, while the equivalent increase under
RCP4.5 and RCP2.6 is significantly lower (see Table 5).
For 2021–2050, however, there is a higher convergence
in the factor of increase in frequency (~ 1.5x) duration (~
2.4x–2.7x), I<sub>mean</sub> (~ 1.5x) and severity (~ 5x–7x) between
MMM of RCP4.5 and RCP8.5.

In general, MHWs become stronger and more intense in
response to increasing greenhouse gas forcing and especially
towards the end of the century. RCP2.6, however, shows
a slight increase in MHW signatures with time but lower
than RCP4.5 and RCP8.5. Note here that certain models
demonstrate stronger climate change responses than others,
likely due to the choice of the driving GCM rather than to
the individual RCSM biases. Much of the MHW evolution
is found to occur mainly due to an increase in the mean SST,
but an increase in daily SST variability also plays a notice-
able role. Complementary sensitivity tests also prove that a
mean shift in SST distributions alone cannot be responsible
for the futures changes in MHWs.

Overall, the MHW and SST changes predicted for the
21st century will clearly impact the vulnerable Mediterra-
nean Sea ecosystems. What was encountered as widespread
consequences from the MHW 2003 could become the “new
normal”, since our analysis signified that future MHWs
become longer and more intense than this event in the near
future. Especially under RCP8.5 and 2071–2100, MHWs
can become three times longer than the MHW 2003, with
mean intensities three times higher. While RCP8.5 is the
business-as-usual scenario, RCP2.6 is the closest to Paris
agreement limits, which could offer a relative stability in
both the SST increase and MHW evolution in the basin after the mid-21st century. MHWs exert a strong influence not only on marine ecosystems but also on marine-dependent economies and hence society. Therefore, more research is needed towards an improved mechanistic understanding of these events and their underlying physical drivers. In a constantly warming world, this information, along with projections of large-scale future MHW evolution, might help identify regions with a physical predisposition to these extreme occurrences. In combination with biogeochemical studies, more light could be shed on the full extent of the biological system risks related to MHWs.

**Acknowledgements** “We would like to thank the anonymous reviewers for their constructive suggestions. This research was funded by the MARmaED project, which has received funding from the European Union’s Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement No 675997. The result of this publication reflects only the author’s view and the Commission is not responsible for any use that may be made of the information it contains. This work is also a part of the Med-CORDEX initiative (www.medcordex.eu) and HyMeX programme (www.hymex.org). Dmitry V.Sein was supported by the EC Horizon 2020 project PRIMAVERA under the grant agreement 641727 and the state assignment of FASO Russia (theme No.0149-2018-0014)”. V.Djordjevic was partially supported by the Serbian Ministry of Science, Education and Technological Development, under grant No.III43007.

**Open Access** This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

## References

Adloff F, Somot S, Sevault F, Jorda G, Aznar R, Deque M, Herrmann M, Marcos M, Dubois C, Padorno E (2015) Mediterranean Sea response to climate change in an ensemble of twenty first century scenarios. Clim Dyn 45(9-10):2775-2802

Akhtar N, Brauch J, Ahrens B (2018) Climate modeling over the Mediterranean Sea: impact of resolution and ocean coupling. Clim Dyn 51(3):933-948

Artale V, Calmanti S, Carillo A, Dell'Aquila A, Herrmann M, Pisacane G, Ruti PM, Sannino G, Struglia MV, Giorgi F (2010) An atmosphere-ocean regional climate model for the Mediterranean area: assessment of a present climate simulation. Clim Dyn 35(5):721-740

Bensoussan N, Romano JC, Harmelin JG, Garrabou J (2010) High resolution characterization of northwest Mediterranean coastal waters thermal regimes: to better understand responses of benthic communities to climate change. Estuar, Coast Shelf Sci 87(3):431-441

Bensoussan N, Garreau P, Pairaud I, Somot S, Garrabou J (2013) Multidisciplinary approach to assess potential risk of mortality of benthic ecosystems facing climate change in the NW Mediterranean Sea. In: OCEANS 2013 MTS/IEEE, San Diego: An Ocean in Common, pp 1-7

Benthuysen J, Feng M, Zhong L (2014) Spatial patterns of warming off Western Australia during the 2011 Ningaloo Niño: quantifying impacts of remote and local forcing. Cont Shelf Res 91:232-246

Benthuysen JA, Oliver EC, Feng M, Marshall AG (2018) Extreme marine warming across tropical Australia during austral summer 2015-2016. J Geophys Res: Oceans 123(2):1301-1326

Bond NA, Cronin MF, Freeland H, Mantua N (2015) Causes and impacts of the 2014 warm anomaly in the NE Pacific. Geophys Res Lett 42(9):3414-3420

Cavicchia L, Gualdi S, Sanna A, Oddo P, et al. (2015) The regional ocean-atmosphere coupled model COSMO-NEMO\_MFS. CMCC Research Paper (RP0254)

Cavole LM, Demko AM, Diner RE, Giddings A, Koester I, Pagniello CM, Paulsen ML, Ramirez-Valdez A, Schwenck SM, Yen NK (2016) Biological impacts of the 2013-2015 Warm-water Anomaly in the Northeast Pacific: winners, losers, and the future. Oceanography 29(2):273-285

Cebrian E, Uriz MJ, Garrabou J, Ballesteros E (2011) Sponge mass mortalities in a warming mediterranean sea: are cyanobacteria-harboring species worse off? PLoS One 6(6):e20-e211

Cerrano C, Bavestrello G, Bianchi CN, Cattaneo-Vietti R, Bava S, Morganti C, Morri C, Picco P, Sara G, Schiaparelli S (2000) A catastrophic mass-mortality episode of gorgonians and other organisms in the Ligurian Sea (North-western Mediterranean), summer 1999. Ecol Lett 3(4):284-293

Chen K, Gawarkiewicz GG, Lentz SJ, Bane JM (2014) Diagnosing the warming of the Northeastern US Coastal Ocean in 2012: a linkage between the atmospheric jet stream variability and ocean response. J Geophys Res: Oceans 119(1):218-227

Chen K, Gawarkiewicz G, Kwon YO, Zhang WG (2015) The role of atmospheric forcing versus ocean advection during the extreme warming of the Northeast US continental shelf in 2012. J Geophys Res: Oceans 120(6):4324-4339

Coma R, Ribes M, Serrano E, Jiménez E, Salat J, Pascual J (2009) Global warming-enhanced stratification and mass mortality events in the Mediterranean. Proc Natl Acad Sci 106(15):6176-6181

Crisci C, Bensoussan N, Romano JC, Garrabou J (2011) Temperature anomalies and mortality events in marine communities: insights on factors behind differential mortality impacts in the NW Mediterranean. PLoS One 6(9):e23814

Di Camillo CG, Cerrano C (2015) Mass mortality events in the nw Adriatic sea: phase shift from slow-to fast-growing organisms. PloS One 10(5):e0126689

Di Lorenzo E, Mantua N (2016) Multi-year persistence of the 2014/15 North Pacific marine heatwave. Nat Clim Change 6(11):1042-1047

Diaz-Almela E, Marba N, Duarte CM (2007) Consequences of Mediterranean warming events in seagrass (posidonia oceanica) flowering records. Global Change Biol 13(1):224-235

Diffenbaugh NS, Pal JS, Giorgi F, Gao X (2007) Heat stress intensification in the Mediterranean climate change hotspot. Geophysical Research Letters 34(11)

Djurdjevic V, Rajkovic B (2008) Verification of a coupled atmosphere-ocean model using satellite observations over the Adriatic Sea. Ann Geophysicae: Atmos Hydrosphered Space Sci 26:1935-1954

Feng M, McPhaden MJ, Xie SP, Hafner J (2013) La Niña forces unprecedented Leeuwin Current warming in 2011. Sci Rep 3:1277

Frölicher TL, Laufkötter C (2018) Emerging risks from marine heat waves. Nat Commun 9(1):650

Frölicher TL, Fischer EM, Gruber N (2018) Marine heatwaves under global warming. Nature 560(7718):360

Galli G, Solidoro C, Lovato T (2017) Marine heat waves hazard 3D maps, and the risk for low motility organisms in a warming Mediterranean Sea. Front Mar Sci 4:136

Garrabou J, Perez T, Sartoretto S, Harmelin J (2001) Mass mortality event in red coral Corallium rubrum populations in the Provence
region (France, NW Mediterranean). Mar Ecol Progress Ser 217:263-272

Garrabou J, Coma R, Bensoussan N, Bally M, Chevaldonné P, Cigliano M, Diaz D, Harmelin JG, Gambi M, Kersting D (2009) Mass mortality in Northwestern Mediterranean rocky benthic communities: effects of the 2003 heat wave. Global Change Biol 15(5):1090-1103

Giorgi F (2006) Climate change hot-spots. Geophysical research letters 33(8)

Grazzini F, Viterbo P (2003) Record-breaking warm sea surface temperature of the Mediterranean Sea. ECMWF Newslett 98:30-31

Gualdi S, Somot S, Li L, Artale V, Adani M, Bellucci A, Braun A, Calmanti S, Carillo A, Dell'Aquila A (2013) The CIRCE simulations: regional climate change projections with realistic representation of the Mediterranean Sea. Bull Am Meteorol Soc 94(1):65-81

Hewitt HT, Bell MJ, Chassignet EP, Czaja A, Ferreira D, Griffies SM, Hyder P, McClean JL, New AL, Roberts MJ (2017) Will high-resolution global ocean models benefit coupled predictions on short-range to climate timescales? Ocean Model 120:120-136

Hobday A, Oliver E, Gupta AS, Benthuysen J, Burrows M, Donat M, Holbrook N, Moore MT PJ, Wernberg T, Smale D (2018) Categorizing and naming marine heatwaves. Oceanography 31(2):1-13

Hobday AJ, Pecl GT (2014) Identification of global marine hotspots: sentinels for change and vanguards for adaptation action. Rev Fish Biol Fish 24(2):415-425

Hobday AJ, Alexander LV, Perkins SE, Smale DA, Straub SC, Oliver EC, Benthuysen JA, Burrows MT, Donat MG, Feng M (2016) A hierarchical approach to defining marine heatwaves. Progress Oceanogr 141:227-238

Huete-Stauffer C, Vielmini I, Palma M, Navone A, Panzalis P, Vezzulli L, Misic C, Cerrano C (2011) Paramuricea clavata (anthozoa, octocorallia) loss in the marine protected area of tavolarola (sardinia, italy) due to a mass mortality event. Mar Ecol 32:107-116

Hughes TP, Kerry JT, Álvarez-Noriega M, Álvarez-Romero JG, Anderson KD, Baird AH, Babcock RC, Beger M, Bellwood DR, Berkelmans R et al. (2017) Global warming and recurrent mass bleaching of corals. Nature 543(7645):373

Hughes TP, Anderson KD, Connolly SR, Heron SF, Kerry JT, Lough JM, Baird AH, Baum JK, Berumen ML, Bridge TC et al. (2018) Spatial and temporal patterns of mass bleaching of corals in the anthropocene. Science 359(6371):80-83

IPCC (2007) In: Solomon S, Qin D, Manning M, Chen Z, Marquis M, Averyt KB, et al. (eds), Climate Change 2007: The Physical Science Basis. Contribution of Working Group i to the Fourth Assessment Report of the IPCC

Jordà G, Marbà N, Duarte CM (2012) Mediterranean seagrass vulnerable to regional climate warming. Nat Clim Change 2(11):821-824

Kersting DK, Bensoussan N, Linares C (2013) Long-term responses of the endemic reef-builder cladocora caespitosa to Mediterranean warming. PLOS ONE 8:1-12, https://doi.org/10.1371/journal.pone.0070820

King AD, Karoly DJ, Henley BJ (2017) Australian climate extremes at 1.5° C and 2° Cof global warming. Nat Clim Change 7(6):412

Kirtman B, SB P, JA A, GJ B, R B, I C, FJ DR, AM F, M K, GA M, M P, A S, C S, R S, GJ vO, G V, HJ W (2013) 13: Near-term climate change: Projections and predictability. In: Climate Change 2013: The physical science basis. contribution of working group i to the fifth assessment report of the intergovernmental panel on climate change

Knutti R, Abramowitz G, Collins M, Eyring V, Gleckler P, Hewitson B, Mearns L (2010) Good practice guidance paper on assessing and combining multi model climate projections. IPCC expert meeting on assessing and combining multi-model climate projections 1-15

Lima FP, Wethey DS (2012) Three decades of high-resolution coastal sea surface temperatures reveal more than warming. Nat Commun 3:704

Linares C, Coma R, Diaz D, Zabala M, Hereu B, Dantart L (2005) Immediate and delayed effects of a mass mortality event on gorgonian population dynamics and benthic community structure in the NW Mediterranean Sea. Mar Ecol Progress Ser 305:127-137

Lough J (2000) 1997-98: unprecedented thermal stress to coral reefs? Geophys Res Lett 27(23):3901-3904

L'Hévéder B, Li L, Sevault F, Somot S (2013) Interannual variability of deep convection in the Northwestern Mediterranean simulated with a coupled AORCM. Clim Dyn 41(3-4):937-960

MacKenzie BR, Schiedek D (2007) Daily ocean monitoring since the 1860s shows record warming of northern European seas. Global Change Biol 13(7):1335-1347

Marba N, Duarte CM (2010) Mediterranean warming triggers sea-grass (posidonia oceanica) shoot mortality. Global Change Biol 16(8):2366-2375

Mariotti A, Pan Y, Zeng N, Alessandri A (2015) Long-term climate change in the Mediterranean region in the midst of decadal variability. Clim Dyn 44(5-6):1437-1456

Mavrakis AF, Tsiros IX (2018) The abrupt increase in the aegean sea surface temperature during the June 2007 southeast mediterranean heatwave-a marine heatwave event? Weather

McSweeney C, Jones R, Lee RW, Rowell D (2015) Selecting CMIP5 GCMs for downscaling over multiple regions. Clim Dyn 44(11-12):3237-3260

Mills KE, Pershing AJ, Brown CJ, Chen Y, Chiang FS, Holland DS, Lehuta S, Nye JA, Sun JC, Thomas AC et al. (2013) Fisheries management in a changing climate: lessons from the 2012 ocean heat wave in the Northwest Atlantic. Oceanography 26(2):191-195

Munari C (2011) Effects of the 2003 European heatwave on the benthic community of a severe transitional ecosystem (comacchio saltworks, italy). Mar Pollut Bull 62(12):2761-2770

Olita A, Sorgente R, Natale S, Ribotti A, Bonanno A, Patti B (2007) Effects of the 2003 European heatwave on the Central Mediterranean Sea: surface fluxes and the dynamical response. Ocean Sci 3(2):273-289

Oliver EC, Wotherspoon SJ, Chamberlain MA, Holbrook NJ (2014) Projected Tasman Sea extremes in sea surface temperature through the twenty-first century. J Clim 27(5):1980-1998

Oliver EC, Benthuysen JA, Bindoff NL, Hobday AJ, Holbrook NJ, Mundy CN, Perkins-Kirkpatrick SE (2017) The unprecedented 2015/16 Tasman Sea marine heatwave. Nat Commun 8(16):101

Oliver EC, Donat MG, Burrows MT, Moore PJ, Smale DA, Alexander LV, Benthuysen JA, Feng M, Gupta AS, Hobday AJ (2018a) Longer and more frequent marine heatwaves over the past century. Nat Commun 9(1):1324

Oliver EC, Lago V, Hobday AJ, Holbrook NJ, Ling SD, Mundy CN (2018b) Marine heatwaves off eastern Tasmania: trends, interannual variability, and predictability. Progress Oceanogr 161:116-130

Perez T, Garrabou J, Sartoretto S, Harmelin JG, Francour P, Vacelet J (2000) Mortalité massive d'invertébrés marins: un événement sans précédent en Méditerranée nord-occidentale. Comptes Rendus de l'Académie des Sciences-Series III-Sciences de la Vie 323(10):853-865

Pisano A, Nardelli BB, Tronconi C, Santoleri R (2016) The new Mediterranean optimally interpolated pathfinder AVHRR SST Dataset (1982-2012). Rem Sens Environ 176:107-116

Rivetti I, Fraschetti S, Lionello P, Zambianchi E, Boero F (2014) Global warming and mass mortalities of benthic invertebrates in the Mediterranean Sea. PloS One 9(12):e115655

Rixen M, Beckers JM, Levitus S, Antonov J, Boyer T, Maillard C, Fichaut M, Balopoulos E, Iona S, Dooley H et al. (2005) The Western Mediterranean Deep Water: a proxy for climate change. Geophysical Research Letters 32(12)
Roberts MJ, Hewitt HT, Hyder P, Ferreira D, Josey SA, Mizielinski M, Shelly A (2016) Impact of ocean resolution on coupled air-sea fluxes and large-scale climate. Geophysical Research Letters 43(19)

Ruti P, Somot S, Giorgi F, Dubois C, Flaounas E, Obermann A, Dell'Aquila A, Pisacane G, Harzallah A, Lombardi E, et al. (2016) MED-CORDEX initiative for Mediterranean climate studies, bams, https://doi.org/10.1175.Tech. rep., BAMS-D-14-00176.1

Scannell HA, Pershing AJ, Alexander MA, Thomas AC, Mills KE (2016) Frequency of marine heatwaves in the North Atlantic and North Pacific since 1950. Geophys Res Lett 43(5):2069–2076

Schaeffer A, Roughan M (2017) Subsurface intensification of marine heatwaves off southeastern Australia: the role of stratification and local winds. Geophys Res Lett 44(10):5025–5033

Schiaparelli S, Castellano M, Povero P, Sartoni G, Cattaneo-Vietti R (2007) A benthic mucilage event in North-Western Mediterranean Sea and its possible relationships with the summer 2003 European heatwave: short term effects on littoral rocky assemblages. Mar Ecol 28(3):341–353

Schlegel RW, Oliver EC, Perkins-Kirkpatrick S, Kruger A, Smit AJ (2017a) Predominant atmospheric and oceanic patterns during coastal marine heatwaves. Front Mar Sci 4:323

Schlegel RW, Oliver EC, Wernberg T, Smit AJ (2017b) Nearshore and offshore co-occurrence of marine heatwaves and cold-spells. Progress Oceanogr 151:189–205

Sein DV, Mikolajewicz U, Gröger M, Fast I, Cabos W, Pinto JG, Hagemann S, Semmler T, Izquierdo A, Jacob D (2015) Regionally coupled atmosphere-ocean-sea ice-marine biogeochemistry model ROM: 1. description and validation. J Adv Model Earth Syst 7(1):268–304

Selig ER, Casey KS, Bruno JF (2010) New insights into global patterns of ocean temperature anomalies: implications for coral reef health and management. Global Ecol Biogeogr 19(3):397–411

Sevault F, Somot S, Alias A, Dubois C, Lebeaupin-Brossier C, Nabat P, Adloff F, Déqué M, Decharme B (2014) A fully coupled Mediterranean regional climate system model: design and evaluation of the ocean component for the 1980–2012 period. Tellus A: Dyn Meteorol Oceanogr 66(1):23,967

Somot S, Sevault F, Déqué M (2006) Transient climate change scenario simulation of the Mediterranean Sea for the twenty-first century using a high-resolution ocean circulation model. Clim Dyn 27(7–8):851–879

Sparnocchia S, Schiano M, Picco P, Bozzano R, Cappelletti A (2006) The anomalous warming of summer 2003 in the surface layer of the Central Ligurian Sea (Western Mediterranean). Annales Geophysicae 24:443–452

Wernberg T, Smale DA, Tuya F, Thomsen MS, Langlois TJ, De Bettignies T, Bennett S, Rousseaux CS (2013) An extreme climatic event alters marine ecosystem structure in a global biodiversity hotspot. Nat Clim Change 3(1):78–82

Wernberg T, Bennett S, Babcock RC, de Bettignies T, Cure K, Depczynski M, Dufois F, Fromont J, Fulton CJ, Hovey RK et al. (2016) Climate-driven regime shift of a temperate marine ecosystem. Science 353(6295):169–172

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

## Affiliations

Sofia Darmaraki¹ · Samuel Somot¹ · Florence Sevault¹ · Pierre Nabat¹ · William David Cabos Narvaez² ·

Leone Cavicchia³ · Vladimir Djurdjevic⁴ · Laurent Li⁵ · Gianmaria Sannino⁶ · Dmitry V. Sein⁷,⁸

Samuel Somot
samuel.somot@meteo.fr

Florence Sevault
florence.sevault@meteo.fr

Pierre Nabat
pierre.nabat@meteo.fr

William David Cabos Narvaez
william.cabos@uah.es

Leone Cavicchia
leone.cavicchia@unimelb.edu.au

Vladimir Djurdjevic
vdj@ff.bg.ac.rs

Laurent Li
laurent.li@lmd.jussieu.fr

Gianmaria Sannino
gianmaria.sannino@enea.it

Dmitry V. Sein
dmitry.sein@awi.de

¹ CNRM, Centre National de Recherches Météorologiques, UMR 3589, Université de Toulouse, Météo-France, CNRS, 42 Avenue Coriolis, 31057 Toulouse, France

² Department of Physics and Mathematics, UAH, University of Alcalá, Madrid, Spain

³ CMCC, Centro Euro Mediterraneo sui Cambiamenti Climatici, Lecce, Italy

⁴ Faculty of Physics, University of Belgrade, Studentski trg 12, Belgrade, Serbia

⁵ LMD Laboratoire de Meteorology Dynamique, Centre National de la Recherche Scientifique (CNRS), Université Pierre et Marie Curie (Paris 6), Paris, France

⁶ ENEA, via Anguillarese 301, 00123 Rome, Italy

⁷ Alfred-Wegener-Institute for Polar and Marine Research (AWI), 27568 Bremerhaven, Germany

⁸ Shirshov Institute of Oceanology, Russian Academy of Science, 36 Nahimovskiy Prospect, Moscow 117997, Russia

<figure><img src="image_22.png" /></figure>