Past Variability of Mediterranean Sea Marine Heatwaves

Marine heatwaves (MHWs) are episodes of anomalous warming in the ocean, responsible for widespread impacts on marine ecosystems. For the first time summer MHW variability at surface and three ecosystem‐relevant depths of the upper Mediterranean Sea are assessed here for 1982–2017. We apply a MHW detection algorithm on a hindcast simulation, performed with a high‐resolution, fully coupled regional climate system model. Identified surface events last, on average, 15 days with a mean intensity of 0.6 °C above threshold and a maximum sea surface coverage of around 39%. Subsurface events are seasonally shifted and appear, on average, longer and more intense but less frequent and less extended in space than surface MHWs. We also find significant trends of increase in most MHW properties throughout the period, with severe surface MHWs detected for the first time in 2012, 2015, and 2017. However, MHW spatial ≪ hot spots ≫ are inhomogeneously distributed in surface and deeper layers.


Introduction
Prolonged periods of anomalously warm temperatures in the ocean surface are described as marine heatwaves (MHWs; Frölicher & Laufkötter, 2018;Hobday et al., 2016) and have led to detrimental impacts on marine ecosystems around the world (Smale et al., 2019). For example, they have triggered widespread mortalities and abrupt redistribution of marine species (Garrabou et al., 2001;Wernberg et al., 2013), toxic algal blooms (Cavole et al., 2016), and mass coral bleaching (Hughes et al., 2017) in a matter of weeks or months, with economic losses in fisheries and aquaculture industry (Mills et al., 2013;Oliver et al., 2017). Superimposed on the underlying warming trend of the ocean, MHWs occur regionally from coastal to open ocean, can extend up to thousands of kilometers (Scannell et al., 2016), and may also propagate deeper to the water column (Schaeffer & Roughan, 2017).
To date, our knowledge on past MHWs in the Mediterranean Sea is based mostly on reported mass mortalities of benthic species linked to thermal anomalies, which have been seen to increase since the early 1990s (Coma et al., 2009;Rivetti et al., 2014). One of the first documented, large-scale events occurred in 2003, where surface anomalies of 2-3 • C above normal lasted over a month and resulted in extensive mass mortalities of benthic invertebrates, loss of seagrass meadows, and abrupt changes in community composition (Garrabou et al., 2001). Prior to this, studies focused mostly on local ecological impacts of regional events observed, for instance, in 1994 (Marbà et al., 2015), 1999 (Garrabou et al., 2001), 2006 (Kersting et al., 2013)

10.1029/2019GL082933
Most of these documented episodes affected corraligenous species (e.g., gorgonians, sponges, and seagrass Posidonia oceanica meadows) up to 50 m deep (Rivetti et al., 2014) with rare occasions reported reaching at coast up to 80-160 m (e.g., Arnoux et al., 1992;Rivoire, 1991). Therefore, investigation of subsurface MHWs is equally important to the identification of events at surface. However, most studies so far have analyzed MHWs at the ocean surface, suggesting already an increase in their occurrence, duration, and intensity globally, over the past century (Oliver et al., 2018). Little is known about the evolution of extreme temperatures at depth, since they are harder to observe in large spatial scales and over long time periods. For example, Argo floats during a 2011 MHW in Australia indicated open-ocean temperature anomalies constrained to the mixed layer depth (MLD), while in situ measurements on the shelf revealed their extension down to 100-200 m (Feng et al., 2013). Recently, Jackson et al. (2018) demonstrated that coastal subsurface anomalous warming in the NE Pacific occurred with 1-year lag but persisted for more than 4 years after the surface 2013-2015 MHW was first observed.
Duration of warm anomalies at depth can be longer than at the surface, even reemerging to the MLD in late winter (Deser et al., 2003). Yet it remains unknown how frequently and deep a MHW may penetrate in the open ocean or what are the characteristics of deep events relative to surface MHWs. The only systematic analysis to date is based on coastal in situ data in Australia for 1953-2016 and shows regular warming throughout the water column and a maximum MHW intensity just below the thermocline (Schaeffer & Roughan, 2017). In the case of the Mediterranean Sea, a quantitative assessment using a standarized framework is currently missing for both past surface and subsurface MHWs characteristics. There are, however, indications of an increase in surface MHW frequency and intensity in the basin for the 21st century relative to the past (Darmaraki et al., 2019).
In this study, we apply an existing MHW definition on daily-resolution temperature data between 1982-2017 and describe changes in frequency, duration, intensity, severity, and maximum spatial coverage of past MHWs identified at the surface and selected depths of the upper Mediterranean Sea. To this end, a high-resolution satellite product and a hindcast simulation, performed with a high-resolution, fully coupled, regional climate system model (RCSM), are used.

Model and Simulation
The numerical simulation is performed with the state-of-the-art, fully coupled RCSM model CNRM-RCSM6 (participating in the Med-CORDEX initiative; Ruti et al. (2016)) developed in CNRM (see Nabat et al., 2015;Sevault et al., 2014;Somot et al., 2008, for previous versions). It has a 12-km resolution atmospheric and land component (ALADIN-Climate v6, 91 vertical levels; Daniel et al. (2019)), a 50-km river model (CTRIP; Decharme et al., 2019), an interactive aerosol scheme (TACTIC; Drugé et al. (2019)), and a 6-to 8-km ocean model (NEMOMED12; 75 vertical levels, Beuvier et al., 2012;Hamon et al., 2016). The coupling strategy follows Voldoire et al. (2017) with 1-hr coupling frequency (OASIS3-MCT; Craig et al., 2017). The ocean model domain covers the entire Mediterranean Sea and a small part of the Atlantic as a buffer zone. The Black Sea is parametrized as a freshwater flux combining the total river runoffs of its drainage area and the evaporation minus precipitation budget over the sea. The coupled model is driven by reanalysis at its lateral boundary conditions. We use ERA-Interim (Berrisford et al., 2009) to drive the atmosphere, applying a spectral nudging technique (Colin et al., 2010) and the ORAS4 (Balmaseda et al., 2013) at the oceanic boundary conditions (temperature, salinity, and sea surface height) in the Atlantic part of the domain plus a correction for the ORAS4 sea level following Adloff et al. (2018). The model initial conditions come from the MEDHYMAP data set (Jordà et al., 2017). Before the simulation , 10 years of coupled spin-up has been performed, after a 7-year ocean spin-up and a 78-year spin-up of the river routing model.

Observations
The performance of the model on the yearly extreme sea surface temperature (SST) evolution and the identification of surface MHWs is first evaluated against daily satellite observations between 1982 and 2017. We use the Mediterranean Sea high-resolution (0.04 • grid) L4 data set from Copernicus Marine Service and CNR-ISAC ROME, which provides interpolated remotely sensed surface temperatures from the AVHRR Pathfinder Version 5.2 onto a regular grid (Pisano et al., 2016). Due to difficulties in acquiring long-term and daily in situ time series below surface, subsurface temperature variability of the model was evaluated on interannual timescale, using a 3-D and time-gridded data set based on in situ observations, named here as MEDHYMAP. This is a monthly resolution product which has been already used in the estimation  Imean ( • C) (57.1 ± 6.2) * 10 −2 (65.9 ± 5.3) * 10 −2 (66.9 ± 2.6) * 10 −2 (76.0 ± 1.9) * 10 −2 (62.81 ± 1.67) * 10 −2 Imax trend ( • C/year) 10.8 * 10 −2 73.5 * 10 −3 15.5 * 10 −2 11.9 * 10 −2 12.8 * 10 −2 Note. Shown here also are the MHW domain-averaged mean surface and subsurface MHW frequency, Imean, Duration, Severity (Icum), Surfmax, and Imax and error bar at 95% confidence level for every layer. Linear trends of domain-averaged time series with statistical significance higher than 95% level are indicated in bold. SST = sea surface temperature; MHW = marine heatwave.
of Mediterranean Sea properties and is partly described in Jordà et al. (2017). Prior to any comparisons, observations were first interpolated to the NEMOMED12 grid using the nearest neighbor method.

MHW Detection Method
MHW characteristics are investigated at the sea surface using satellite observations (approximately micrometer thick) and the model's first layer depth (∼1 m thick) and also at 23-m (∼4 m thick), 41-m (∼6 m thick), and at 55-m (∼8.5 m thick) depth, using daily simulated seawater potential temperatures. We focus our analysis on these layers, as they cover the depth range where thermal stress-related mass mortalities of Mediterranean marine species have been regularly observed in the past (see Table S1 in the supporting information). For this study, no a priori connection was assumed between surface and deeper events.
Individual events are detected using a MHW identification framework, defined by Darmaraki et al. (2019). First, a depth-dependent threshold is calculated for each grid point over the 1982-2012 period. Locally, this threshold is the 30-year average of the yearly 99th percentile computed from daily temperatures, hereby called SST 99Q at surface or T 99Q at depth. Abnormally warm days are then identified when the grid point temperature is higher than the local threshold for five consecutive days or more. The final spatiotemporal extent of the MHW, however, is only considered when a minimum of 20% of the basin is affected. The area covered by MHWs though is considered non-contiguous. Therefore, our detection algorithm targets the detection of large-scale, long-lasting, summer MHW events, a choice which is different from Hobday et al. (2016) definition (see Darmaraki et al., 2019, for a more in-depth discussion of the differences).
For each surface or deeper event a set of metrics, adapted from Hobday et al. (2016), is then computed, such as the frequency (number of events per year), duration (number of days above threshold), mean (Imean), and maximum (Imax) intensity (spatiotemporal mean and maximum temperature anomaly with respect to threshold over the event duration). Further, we consider the maximum event coverage (surfmax) and the cumulative intensity Icum (spatiotemporal sum of daily temperature anomalies over the event duration), which represents the severity of each MHW.

Temperature Evolution
The model's behavior in the extreme Mediterranean Sea SST between 1982 and 2017 was first evaluated. Simulated surface temperature threshold patterns (Figures 1a and 1b, top) are well correlated (Corr. Coeff = 0.9) relative to the satellite observations, albeit a slightly warm bias. In addition, the model appeared to capture well the basin-mean SST 99Q interannual variability (Corr. Coeff = 0.9; Figure 1, bottom), despite an underestimation of the positive observational trend (Table 1). On average, the southeast Mediterranean, Ionian, and Tyrrhenian Seas show the highest SST 99Q (∼28-30 • C), in contrast to the northwestern Mediterranean Basin and Aegean and Alboran Seas that manifest the lowest SST 99Q (∼ 24-27 • C) of the period.
To increase our confidence on the model's performance at depth, we then compare its yearly mean surface and deep temperature variability with MEDHYMAP data set ( Figure S1 and Table S2). Simulated surface temperatures this time show a smaller positive bias, a better trend, and, in general, a better agreement with MEDHYMAP data set than with the satellite-based product. The results of the simulation are also in agreement with the mean SST trends proposed by Pastor et al. (2019), Nabat et al. (2014), and Kennedy et al. (2011) for the Mediterranean Sea between 1982 and 2012. Simulated temperatures at depth are well represented, particularly in terms of basin-mean values, interannual variability, and trend (see Table S2). We note here that evaluation of the model's interannual variability at depth with MEDHYMAP data set (monthly resolution) is performed at the absence of long-term, daily, in situ data that would be necessary for a proper validation of MHWs at depth. Nevertheless, it implies a generally consistent behavior of the model in deeper layers.
In particular, the equivalent basin-mean T 99Q time series (Figure 1, bottom) reveal a warming of about 0.02 • C/year at every layer, which decreases slightly with depth (Table 1). In terms of deep temperature threshold patterns (Figures 1c-1e, top), a northwest-southeast, cold (15-24 • C)-warm (25-32 • C) pattern is more evident in the T 99Q map of 23 m, whereas extreme temperatures between 15 and 20 • C span the whole basin at 41-and 55-m depths, except for the southeast Mediterranean Sea, where T 99Q is more pronounced.

Surface MHW Characteristics
A total of 29 events was identified in the observations over 1982-2017 (Table 1 and Figure 2a). On average, they lasted around 20 days, showed an Imean of about 0.6 • C, and covered a maximum of 44% of the Mediterranean Basin. Although fewer in number, the characteristics of the simulated surface events are within the observed range, except for the mean intensity that is slightly higher (Figure 2b). Their trends, however, are underestimated with respect to the observed ones. In fact, all the observed trends at the surface are statistically significant at 95%, except for the severity, on the contrary to the simulations, where only Imax trend is found significant. Despite that, the main observed MHW characteristics are well captured by the model.

Subsurface MHW Characteristics
Relative to surface simulated events, frequency of subsurface MHWs reduces as we go deeper in the water column (see Table 1, Figure 2, and also Tables S5-S7). Conversely, their ensemble-mean duration and severity monotonically increase with depth, while their Imean and Imax increase from surface to 41 m and decrease at 55 m and vice versa for the maximum event coverage. From the equivalent trends only the average subsurface Imean showed a monotonic decrease from 23 to 55 m. It is worth noting that the ensemble-mean Icum is higher at depth, even though there are no deep events as severe as the surface MHWs of 2003, 2012, and 2015 (Table 1). This mean behavior may seem counterintuitive, but it could be explained by the longer, on average, durations of events at depth (30-40 days) than at the surface (15-20 days). The highest (lowest) mean duration of the period is found for the group of events at 55-m depth (surface) and the highest (lowest) mean maximum event coverage at surface (41-m depth).
Individually, the most notable MHWs at depth occur in different years from those at surface (Figures 2c-2e). For example, the subsurface MHW 2014 was the longest (52 days), most severe (32 * 10 6 • C·days·km −2 ) and spatially largest (∼43%) event of the period at 23 m. Similarly, the subsurface MHW 2001 was the longest  (2015), covering, however, at peak half of its area.

MHW Seasonality
Examination of MHW timing revealed a seasonality in their occurrence, following the MLD deepening (see Figures 2f and S2). Surface MHWs were detected between July and September. Subsurface events at 23 m followed between September and October, at 41 m between October and November, and at 55-m depth from mid-October to mid-December. Whether subsurface MHWs are always concomitant to the identified surface events in the same year is a hypothesis that was not investigated here. However, it is worth noting that most of the time a surface MHW was followed by a higher-or lower-intensity subsurface event and only 6 (4) years exhibited MHWs only at surface (depth) (see Tables S3-S7).

MHW Spatial Distribution
The spatial distribution of the average event Imean during 1982-2017 seems to vary slightly between the different depths (Figure 3, top). In general, most parts of the Mediterranean Sea show Imean between ∼0.3 and ∼0.9 • C at every layer. North Mediterranean regions display the most pronounced surface Imean (∼1 • C) both in the model and the observations. In deeper layers they only sustain a local MHW intensification (∼1-1.7 • C) that is displaced southward with depth. In fact, the highest Imean at depth (∼1-2 • C) emerges in the Levantine, Ionian, and Southwestern Basins.
By contrast, the corresponding average duration patterns per MHW of the period reveal marked differences with depth (Figure 3, bottom). At surface, mean event duration ranges between 8 and 13 days in most areas, except for the Levantine Basin (∼20-70 days) and some parts of the Aegean (∼20 days), whose duration is underestimated by the model. In deeper layers, average MHW duration appears to increase progressively from 20 days at 23 m (in the Levantine Basin and some regions around the northwestern Mediterranean, Ionian, and Adriatic Seas) to 20-50 days almost everywhere at 41 and 55 m but far from boundary currents.
It is interesting to note that as far as the Western Mediterranean Sea is concerned, from the surface up to 23-m depth, it is characterized by the most intense but infrequent events. This could imply stronger impacts for the marine ecosystems of the region, which are not adapted to such rare and intense events. There are also specific spots in the Alboran, Ionian, and Tyrrhenian Seas, which seem to sustain MHWs between 8 and 13 days and others that are under MHWs for 60-167 days.

Discussion
Despite the model's mean warm bias and underestimated temperature trend with respect to satellite observations, MHW identification at the sea surface appeared consistent. In particular, almost every event identified up to 2013 in the two data sets correspond to years with documented regional, thermal stress-related mass mortalities of marine species in the basin, whereas events identified after 2014 have not yet been discussed in the literature. Inconsistencies with respect to observed events may emerge due to the 20% minimum spatial threshold, which can inhibit the identification of smaller-scale MHWs. On the other hand, there may also be events identified here that were not observed in real time.
For the first time to our knowledge, Mediterranean MHWs were also investigated at different depths. Individually, they can be as long as and as intense as surface MHWs but not as frequent and large in terms of spatial coverage. The excess heat at depth seems to accumulate in smaller regions, which could also explain the higher, on average, Imean and Imax values displayed at depth relative to the surface. The deep layer disconnection from the air-sea exchanges allows for the long-term preservation of the heat content, due to the low-frequency variability of the ocean interior, especially if no strong mixing occurs. Therefore, the longest MHW durations are found at 55-m depth. As far as Imean is concerned, its increase with depth could be locally linked to spatial variations of the thermocline depth that ranges between 10 and 30 m in the model, as proposed by Schaeffer and Roughan (2017).
In principle, for this study, surface MHWs were not associated with subsurface events. However, on a given year, the latter most of the times was concomitant to the former, especially toward the end of the period (Tables S3-S7), following the seasonal MLD deepening. If a connection is to be assumed, then surface temperature anomalies developed during the summer could propagate into the water column when it is weakly stratified. In particular, we hypothesize that mechanical mixing due to wind could partly explain their progressive propagation at 23 and 41 m during autumn and 55 m in December (see Figure S2). However, complex 3-D heat budget computations are required to verify this theory, which are out of the scope of this study. We also speculate that wind conditions may explain how less intense deep events follow more intense surface MHWs (Imean > 0.7 • C), as well as the opposite situation when surface MHW Imean is below 0.7 • C (see Figure S3), explaining that way the patterns of Figures 2a-2e. In the first case, for a given year, moderate wind speeds could decrease heat at the sea surface, reducing surface MHW intensities, transporting, however, the energy below surface, thus intensifying MHWs at depth. On the other hand, years with less intense or no deep MHWs at all following a surface event could be related to two factors: (1) strong stratification in the upper ocean, due to low wind, which intensifies surface MHWs to the expense of deeper events and (2) strong wind speed during the autumn and winter season that does not allow summer MHWs to propagate later deeper in the water column. The cooling effect of surface latent heat flux in that case would damp warm anomalies. This is in line with Sparnocchia et al. (2006) and Schaeffer and Roughan (2017), which have already underlined the role of wind in modulating MHW properties.
Although our study highlights that intense (>T 99Q ) temperature anomalies occur within the MLD at large scales and away from the coast, they may also fall below MLD at local scales (e.g., shelf) or at large scales but with more attenuated anomalies (<T 99Q ). A first approach in the evolution of MHW characteristics with depth (up to 320 m) is provided in Figure S4. For a better understanding of MHWs at depth, however, further analysis of the heat content of specific events is required, which is out of the scope of this study.

Conclusions
We have identified MHWs at surface and at three ecosystem-relevant depths of the Mediterranean Sea between 1982 and 2017, by applying a MHW detection algorithm on a hindcast simulation performed with the latest version of the CNRM-RCSM6. Comparisons with a high-resolution satellite product show that the model is able to reproduce well-extreme Mediterranean Sea surface features of that period. Despite an underestimation of their mean characteristics with respect to observed events, individual event identification was consistent between the model, the observations, and the existing literature. Following a significant increase in extreme SST with time, simulated surface MHWs, on average, last around 15 days, with a mean intensity of 0.6 • C and a (maximum) minimum spatial coverage of around (39) 20% of the basin. MHWs become longer, more severe, and spatially extended with time, and MHWs of 2003, 2012, and 2015 are identified as the most severe surface events of the period.
We have also investigated subsurface MHWs and their evolution at 23, 41, and 55 m. Despite a reduction of extreme temperatures with depth, subsurface events are characterized, on average, by higher mean and maximum intensity and severity relative to surface MHWs, probably due to their longer durations. However, they display, on average, lower frequency and spatial coverage. Their occurrence appears to follow the seasonal MLD deepening, while the most severe deep events were identified in 2001 and 2014. North Mediterranean areas seem more vulnerable to higher mean MHW intensities at surface, whereas southern regions seem more prone to them as we go deeper in the water column. Long-lasting events on the other hand appear to affect mostly the Levantine region at surface and spread all over the basin but far from boundary currents in deeper layers. For this study no relation was assumed between the surface and subsurface MHWs. Dedicated analysis on specific events is further required to better understand the mechanisms behind the vertical extension and heat content evolution of a MHW from the surface to depth. Horizon 2020 research and innovation programme under the Marie Sklodowska-Curie Grant Agreement 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. We would like to thank Dr. Gabriel Jordá for producing and kindly sharing with us the MEDHYMAP data set. This work is also a part of the Med-CORDEX initiative (www.medcordex.eu) and HyMeX programme (www.hymex.org). Model data used in this study are available online (www.medcordex.eu) and observational data through the Copernicus Marine Environment Monitoring Service (http://marine. copernicus.eu/).