Sub-Regional Marine Heat Waves in the Mediterranean Sea From Observations: Long-Term Surface Changes, Sub-Surface and Coastal Responses

Society is facing climate-related challenges and impacts, such as marine heat wave (MHW) events that adversely affect ecosystems, threaten economies and strengthen storms by warming ocean waters. MHWs are substantially increasing in intensity, duration and frequency worldwide, particularly in the Mediterranean Sea, which responds rapidly to climate change. This study proposes a comprehensive analysis of MHWs in the different sub-regions of the Mediterranean, where the strong spatial variability requires focused attention, from surface to sub-surface and from open to coastal oceans. At surface, the MHW indices have dramatically increased over the last four decades from 1982 to 2020, with an unprecedented acceleration rate in recent years in all sub-regions. Besides the sub-regional features of surface MHWs, the propagation of such events into the ocean interior is also examined highlighting sub-regional and seasonal variability in the sub-surface ocean response. The resulting upper-ocean density stratification to these extreme events is enhanced in all sub-regions which would increase the degree of decoupling between surface and deep oceans causing changes in water masses and marine life. Finally, extremely warm events in coastal waters are also addressed through a case study in the Balearic Islands showing their higher intensity and occurrence in near-shore environment as well as the different response from surface to sub-surface that strongly depends on local features. In addition to this study, the Balearic Islands Coastal Observing and Forecasting System (SOCIB) has implemented a smart platform to monitor, visualize and share timely information on sub-regional MHWs, from event detection in real-time to long-term variations in response to global warming, to diverse stakeholders. Society-aligned ocean information at sub-regional scale will support the policy decision-making and the implementation of specific actions at local, national and regional scales, and thus contribute to respond to societal and worldwide environmental challenges.


INTRODUCTION
Marine heat waves (MHW) are extreme warm ocean temperature events during prolonged periods (Hobday et al., 2016). As a consequence of climate change, MHWs are substantially increasing worldwide in frequency, duration and magnitude (Frölicher et al., 2018;Oliver et al., 2018), with harmful impacts on coastal ecosystems and human activities. MHWs are triggered by anomalous atmospheric conditions that are generated or modulated by anthropogenic warming, large-scale climate modes (e.g., El Niño Southern Oscillation and Atlantic Multi-decadal Oscillation), favorable synoptic systems (e.g., highpressure system and wind stress change), remote atmospheric teleconnections (e.g., Rossby waves and Walker circulation) or local processes (air-sea heat flux) (Schaeffer and Roughan, 2017;Holbrook et al., 2019;Gupta et al., 2020). By reducing the cloud cover, enhancing the insolation and/or suppressing the surface wind speed, these drivers increase the solar radiation into the ocean and reduce the latent and sensible heat losses from ocean to atmosphere (Bond et al., 2015;Holbrook et al., 2019;Gupta et al., 2020). Oceanic drivers have been also identified such as remote oceanic teleconnection (e.g., Rossby and Kelvin waves) or local processes affecting the mixed layer temperature (e.g., ocean advection, mesoscale circulation, and vertical processes) (Schaeffer and Roughan, 2017;Elzahaby and Schaeffer, 2019;Holbrook et al., 2019). The impacts of MHWs on the physical ocean may be multiple: mixed layer shoaling, suppression of vertical mixing, reduction of vertical heat transport, reduction of coastal upwelling and Ekman pumping and increase of warm-water horizontal transport (Holbrook et al., 2019). These changes in the physical ocean system affect, in turn, the marine ecosystems.
Several studies on MHW events worldwide have shown their devastating impacts on marine species and ecosystem structure (Wernberg et al., 2013), causing coral bleaching Hughes et al., 2018), declines in seagrass meadows and kelp forests (Wernberg et al., 2013;Arias-Ortiz et al., 2018;Holbrook et al., 2020;Wernberg, 2021), harmful algal blooms , mass mortality of marine organisms (Cerrano et al., 2000;Mills et al., 2013;Short et al., 2015;Jones et al., 2018), irreversible physiological damage, change in species behavior and irreversible redistribution in coastal benthic communities (Johnson et al., 2011;Wernberg et al., 2016;Smale et al., 2019;Verdura et al., 2019). MHWs are emerging as major players in the modification and restructuration of the entire ecosystem, and will severely impact the provision of ecological goods and services in the coming decades (Mills et al., 2013;Smale et al., 2019;Holbrook et al., 2020). MHWs are not restricted to the surface and can propagate deeper in the water column (Elzahaby and Schaeffer, 2019). The depth penetration of such warm events also have negative ecological consequences (Schaeffer and Roughan, 2017) and could also impact the mesopelagic habitats (Santora et al., 2020).
The ocean stratification plays a crucial role in many physical and biogeochemical processes. It directly affects the deepening of the mixed layer that regulates light availability for photosynthesis and nutrient supply from deep to surface oceans, which are major factors in biological productivity. It also influences the ocean ventilation that involves the subduction of surface waters into the ocean interior and the absorption of heat, carbon and oxygen which have fundamental roles in the ocean acidification and the biogeochemical cycle. As a consequence of global warming, globally-averaged upper-ocean thermal stratification has increased due to surface intensification of the warming signal (Capotondi et al., 2012;Rhein et al., 2013;Fu et al., 2016;Moore et al., 2018). Several studies using in situ observations or models have shown the enhanced upper-ocean density stratification, with regional and seasonal patterns (Capotondi et al., 2012;Yamaguchi and Suga, 2019), providing light availability for the phytoplankton community but preventing vertical nutrient supply from sub-surface ocean to euphotic zone (Doney, 2006).
As a semi-enclosed and relatively small basin, including many coastal areas and islands, the Mediterranean Sea is one of the most vulnerable regions to climate change and responds rapidly to global warming (Giorgi, 2006;Lionello and Scarascia, 2018), from the surface to the deep ocean, from coastal to open ocean waters, through oceanographic and environmental changes (Bethoux et al., 1999). The surface Mediterranean warming rate is approximately three or four times higher than the global ocean one (von Schuckmann et al., 2020;Juza and Tintoré, 2021) with basin-averaged trend estimations of sea surface temperature (SST) around 0.041 ± 0.006 • C/year over 1982-2018 (Pisano et al., 2020), 0.037 ± 0.003 • C/year over 1982(von Schuckmann et al., 2020 and 0.038 ± 0.002 • C/year over 1982-2020 (Juza and Tintoré, 2021). The authors also highlighted the strong spatial variability in SST trend estimations. Trend values over 1982-2020 have been estimated at 0.032 and 0.044 ± 0.002 • C/year as averages in the western and eastern Mediterranean, respectively (Juza and Tintoré, 2021), and at 0.033 ± 0.004 • C/year in the Ionian and Levantine Seas (Ibrahim et al., 2021). The rapid sea surface warming trend in this ocean basin has been associated with strong increase in MHW days, particularly during the last two decades since the beginning of the 20th century Ibrahim et al., 2021;Juza and Tintoré, 2021). MHWs in the Mediterranean have been related to large-scale atmospheric anomalies (e.g., persistent high-pressure system) inducing hot and dry weather as the extreme event in summer 2003 (Sparnocchia et al., 2006;Olita et al., 2007;Holbrook et al., 2019). The global warmingenhanced stratification and extreme events may also modify water mass properties in intermediate and deep layers as observed in the Mediterranean by several authors (e.g., Schröder et al., 2006;Schroeder et al., 2016;Houpert et al., 2016;Ozer et al., 2017;Juza et al., 2019;Vargas-Yáñez et al., 2021) and cause mass mortality events . The coastal nearshore ecosystem is highly vulnerable to global warming and extreme temperature events that threaten the biodiversity in the Mediterranean Sea (Cerrano et al., 2000;Garrabou et al., 2009;Bensoussan et al., 2019;Verdura et al., 2019). The climate signal manifests differently from the coastal areas to the open ocean and in the different sub-regions due to the variety of coastal ocean processes such as upwelling, river discharge, coastal dynamics, waves, geomorphology and atmospheric pattern.
The Mediterranean Sea is experiencing strong spatial variations in response to global warming that require specific consideration, in particular to better understand the drivers, mechanisms and consequences of extreme events on the physical, biogeochemical and biological components of the ocean. The main purpose of this study is to comprehensively characterize and analyze MHWs in the Mediterranean Sea at sub-regional scale from surface to sub-surface and from open ocean to coastal near-shore waters, using remote sensing and multi-platform in situ observations. The satellite products deliver the longest observational time series of surface ocean information over the whole basin, the profiling floats give valuable information on sub-surface conditions and the fixed moorings enable a coastal study providing near-shore data. First, the long-term evolution of MHWs characteristics (mean and maximum intensities, mean duration, and frequency) in the Mediterranean Sea is quantitatively estimated and analyzed using the 39-year satellite observation records. The surface ocean response to global warming is also characterized in the different sub-regions to consider the strong spatial variability of the Mediterranean dynamics and to provide suitable estimations and indices for environmental agencies and decision-makings at local, national and regional scales. Then, a 9-year collection of vertical hydrographic profiles from profiling floats is processed to analyze the propagation of sub-regional MHWs into the ocean interior where amount of heat is absorbed affecting directly marine ecosystems, as well as vertical ocean properties and ocean circulation. The global-warming modified stratification in the Mediterranean, which can cause changes in water masses and marine life, is also investigated in the different sub-regions through the depth penetration of MHWs. Finally, the ocean response to extreme temperature events will be also addressed in the coastal environment, using mooring data in the nearshore waters of the Balearic Islands and comparing them with nearby near-shore and off-shore satellite locations. Thus, a comprehensive insight on the behavior of the ocean, enabled by the systematic monitoring over long-term periods, will be provided from local to sub-regional scales, from the surface to the sub-surface and from open ocean to coastal waters. This paper is organized as follows: section "Sub-Regional Approach" introduces the sub-regional approach in the Mediterranean Sea and describes the dynamical sub-regions. Section "Data and Methodology" provides a detailed description of the data used, methodologies applied and diagnostics performed for the computation and analysis of MHWs. Section "Sub-Regional Marine Heat Waves" discusses the sub-regional variability of MHWs at surface over the last four decades from 1982 to 2020. Sections "Depth Penetration of Marine Heat Waves" and "Coastal Marine Heat Waves" address the sub-surface and coastal ocean responses to such extreme events, respectively. Finally, Section "Conclusion and Perspectives" provides conclusions and perspectives.

SUB-REGIONAL APPROACH
The present study aims at estimating the evolution of MHW characteristics in the different sub-regions of Mediterranean Sea. The sub-regional approach has been shown to be relevant in the Mediterranean which has strong spatial variability in terms of ocean dynamics, water mass formation, ocean indicators including ocean warming or trophic regimes (D'Ortenzio and Ribera d'Alcalá, 2009;Pinardi et al., 2015;Juza et al., 2016a;Mayot et al., 2016;Juza and Tintoré, 2021). The linear trends with 95% confidence level (according to the methodology from Chelton (1983) using the effective degrees of freedom) of SST over 1982-2020 in the Mediterranean (Figure 1A) also highlight the longitudinal difference and the strong variability inside each sub-basin in surface warming rate with values oscillating between 0.021 and 0.066 • C/year [as also observed in Juza and Tintoré (2021)]. MHWs may also vary in spatial extent, locally confined or broadly distributed in sub-regions (Olita et al., 2007) or ocean basins (Lee et al., 2010), differing in surface and deep layers (Darmaraki et al., 2019), depending on processes that cause and maintain them as well as the geometry of regions (Holbrook et al., 2020). The global average of MHW intensity, duration and frequency increased from 1982 to 2016 with a trend of 0.085 • C, 1.3 days and 0.45 annual event per decade, respectively . In the eastern Mediterranean, the sub-basin average of the MHW frequency trend increased by 1.2 events per decade over 1982-2020 (Ibrahim et al., 2021). The authors also observed local MHW duration and maximum intensity trends reaching values up to 10 days and 0.25 • C per decade, respectively. Maps of annual or period-averaged MHW characteristics in the Mediterranean also highlighted strong spatial variability in MHWs associated with local or sub-regional responses to climate change Darmaraki et al., 2019).
Therefore, in order to properly estimate the long-term variations of the MHW characteristics in the Mediterranean Sea, they are computed in the different sub-regions of the ocean basin, as defined by Manca et al. (2004) (Table 1 and Figure 1A). These sub-regions have been distinguished accordingly to the schematic representation of the upper sub-basin scale dynamics from Millot (1999) for the western Mediterranean (WMed), from Robinson et al. (1991) and Malanotte-Rizzoli et al. (1999) for the eastern Mediterranean (EMed) and from Mosetti (1983) for the Adriatic Sea. Advective processes prevail in some sub-regions. Other subregions are mostly affected by convective mixing, in particular winter convection that plays a major role in vertical water column homogenization such as in the Gulf of Lion, South Adriatic, Cretan Sea and Rhode Gyre (Pinardi et al., 2015;Testor et al., 2018). These sub-regions also allow taking into consideration the spatial variability of the surface ocean warming rate in response to global change ( Figure 1A).

39-Year Satellite Products
Marine heat waves are identified using sea surface temperature (SST) data from reprocessed (REP) and near real-time (NRT) satellite products, which are freely available and distributed by the Copernicus Marine Service 1 , providing the longest continuous time series of SST since 1982. The REP  Smith and Sandwell, 1997) with moorings at Cala Millor (at 39.59 • N -3.4 • E, black circle) and Son Bou (at 39.9 • N -4.06 • E, blue circle), and their closest near-shore (at 39.6 • N -3.4 • E and 39.9 • N -4.06 • E, respectively) and off-shore (at 39.44 • N -3.85 • E and 39.77 • N -3.98 • E, respectively) grid points from satellite product (stars).
(SST_MED_SST_L4_REP_OBSERVATIONS_010_021 2 ) and NRT (SST_MED_SST_L4_NRT_OBSERVATIONS_010_004_a_ V2 3 ) satellite products in the Mediterranean Sea provide optimally-interpolated estimates of SST into a regular grid (1/20 • × 1/20 • and 1/16 • × 1/16 • , respectively) at daily resolution (Buongiorno Nardelli et al., 2013;Pisano et al., 2016) over the 1982-2019 and 2020-present periods, respectively. Given the quality of the NRT product 4 and its good agreement with the REP product in MHW detection over recent years (not shown), it is used to extend the period of study in recent year(s) in order to provide timely ocean temperature information to diverse stakeholders. The daily REP product is also used to build the climatologies that are used to compute anomalies and MHWs. The long-term changes of MHW characteristics are analyzed over the period 1982-2020 using the period 1982-2015 as reference for the daily climatologies (Juza and Tintoré, 2021). Since the proposed study aims at both detecting sub-regional MHWs in the recent and future years and providing information about the long-term variability, the reference period chosen starts as early as possible (in 1982 for the satellite SST) and ends in 2015 (discarding the recent years to be compared between themselves with a common independent reference). The 1982-2015 reference period for the climatologies also allows being aligned with the methodology applied for the sub-regional indicators recently developed in the Mediterranean Sea (Juza and Tintoré, 2021).

Marine Heat Waves
Marine heat waves occur when ocean temperatures are much warmer than normal during a prolonged period. Hobday et al. (2016) characterize as marine heat spike (MHS) SST that exceeds Frontiers in Marine Science | www.frontiersin.org  Hobday et al. (2016), two successive MHW events with 2-day or less time break will be considered as a continuous event. This also allows discarding the unrealistic jumps in regional SST time series due to sparse erroneous daily interpolated data in the NRT product. Finally, as recommended by the authors to extract useful climatologies from inherently variable data, the daily climatological mean and threshold time series are smoothed using a 30-day moving window. In addition, MHW indices have been defined to estimate long-term changes in MHW characteristics over the last four decades from 1982 to 2020. In particular, for each year, the following MHW metrics are computed: mean and maximum intensities above the threshold, mean duration and number of discrete events. The seasonal component of MHW characteristics will be also addressed applying the same approach to each season of every year independently over the whole period. The seasons are defined as follows: winter (January, February, and March), spring (April, May, and June), summer (July, August, and September) and autumn (October, November, and December). The spatio-temporal variations of the above-mentioned MHW indices will be comprehensively analyzed. Moreover, some key numbers of MHW total days will be also provided as interesting indicators for ecosystem stressor, although this metric is an aggregation of MHW duration and frequency.

Web-Based Visualization Tool
The needs for reliable, available and timely information about climate change-related events are increasing for the coastal communities and those dependent on coastal regions. The development of smart platforms, which can rapidly communicate information and warning to users, is largely enabled by the Copernicus Marine Service which provides full, free and open access to historical and NRT data. In this context, the Balearic Islands Coastal Observation and Forecasting System (SOCIB) has implemented a visualization tool to provide continuous information on the ocean temperature state and variability at sub-regional scale, detecting MHWs through the daily monitoring of SST (see Section "Daily Monitoring"), following their evolution through annual indices (see Section "Annual Indices") and their temporal variations since 1982 (see Section "From Seasonal to Decadal Variability"). This user-friendly interface is available at https://apps.socib.es/subregmed-marineheatwaves/ and is automatically updated every day for the daily monitoring and every year for the annual MHW indices.

9-Year Profiling Floats
To analyze the depth extent of MHWs, quality-controlled hydrographic profiles from profiling floats are extracted from the EN4.2.1 dataset 5 (Gouretski and Reseghetti, 2010;Good et al., 2013), which is freely available and distributed by Met-Office, providing temperature and salinity profiles in the Mediterranean Sea from 5 to 2,000 m depth. The Mediterranean basin is globally covered by profiling floats that sample all sub-regions since 2012 (not shown), with an increasing number of profiles over the period 2012-2020 varying between 2,397 and 10,683 profiles per year (Figure 2A). Hence, the hydrographic profiles over 2012-2020 are used to investigate the MHW vertical propagation into the ocean interior.

SeaDataNet Climatology and Anomaly Profiles
Temperature (T), salinity (S), and density (ρ) anomaly profiles (Ta, Sa, and ρa, respectively) are computed using the SeaDataNet 6 SDC_MED_CLIM_TS_V2 product 7 , in particular the monthly climatology of T and S for the Mediterranean Sea over the period 1985-2018. This product is interpolated into a regular grid of 0.125 • × 0.125 • spatial resolution with 92 vertical levels from surface to 5,000 m depth. By linear interpolation, the climatological T, S, and ρ profiles (Tclim, Sclim, and ρclim, respectively) are computed for every in situ profile. Then, Ta, Sa, and ρa are computed such as Ta(z) = T(z) − Tclim(z), Sa(z) = S(z) − Sclim(z), and ρa(z) = ρ(z)-ρclim(z), where z is the depth from the climatology.

Marine Heat Wave Anomaly Profiles
Since the climatology profiles from SeaDataNet are built using the 1985-2018 period, local MHWs are also identified from satellite products using the 1985-2018 reference period, following the methodology from Hobday et al. (2016) described in Section "Marine Heat Waves, " over the period 2012-2020 for their joint analysis with profiling floats. Then, T, S, and ρ profiles are identified as MHW profiles when coinciding in time and space with MHW as detected at surface from the satellite products, more precisely, when the nearest time (day) and location (longitude and latitude) grid point of satellite SST has been locally identified as MHW.
To remove inconsistent profiles, additional tests have been applied to reject MHW profiles with negative surface temperature anomaly that would contradict the positive anomaly which is indicated by the satellite SST field (Elzahaby and Schaeffer, 2019). 1,705 out of 6,973 MHW profiles have been rejected. Between 193 and 1,672 MHW profiles per year were extracted from 2012 to 2020, representing between 3.6 and 16.4% of the annual total number of profiles in the Mediterranean. The total of MHW anomaly profiles which are available over 2012-2020 in the whole Mediterranean ( Figure 2B) represent approximately 9% of all profiles. Moreover, the satellite SST used to detect MHW events and the most surfacing temperature measurements in profiles are highly correlated (0.99) with a mean absolute error (0.38 • C) that is likely due to differences in instrumentation, product resolution and more importantly depth of measurement (68% of the shallowest measurements along the selected MHW profiles being at 5 m, 28% at 10 m and 3% at 25 m).
The sub-regional distributions of Ta, Sa, and ρa profiles are characterized by their median to estimate the vertical propagation of MHW. This statistical approach is appropriate since the distributions of ocean hydrographic properties are not Gaussian (Vargas-Yáñez et al., 2020;Juza and Tintoré, 2021). Note that the median estimations are computed when a minimum of three profiles are available in the sub-regions.

Stratification Changes
Several methodologies have been established to define the upperocean stratification such as the Stratification Index (Herrmann et al., 2010;Somot et al., 2018), Convective Resistance (Frajka-Williams et al., 2014) or density difference between the surface and fixed sub-surface depth Capotondi et al., 2012;Fu et al., 2016). In this study, in order to have a more widely applicable indicator of change in the upper-ocean temperature, salinity and density stratifications in various dynamical sub-regions, the vertical gradients of Ta, Sa, and ρa between the surface and sub-surface are used to characterize the upper-ocean stratification changes induced by MHW events in the Mediterranean at subregional scale. More precisely, the upper-ocean stratification is defined here as the difference between the surface and 150 m depth, which is indicative of the degree of decoupling between the surface and deeper oceans, with important consequences for many biogeochemical processes (e.g., nutrient fluxes between the euphotic zone and ocean interior) (Capotondi et al., 2012;Fu et al., 2016). In the Mediterranean Sea, the 150 m depth is representative of the separation between the surface to intermediate water masses (Robinson et al., 2001;Vargas-Yáñez et al., 2021). The upper 150 m of the ocean is also the layer which strongly and rapidly interacts with the atmosphere and where primary productivity is observed (Lavigne et al., 2015). Although the 150 m depth is below the mixed layer in strongly stratified sub-regions and within the mixed layer in weakly stratified sub-regions, the vertical gradients between the surface and 150 m are useful indicators of changes in the upper ocean in response to global warming and extremely warm events that are modifying ocean properties, affecting marine ecosystems and intensifying storms. The upper-ocean stratification changes in T, S, and ρ are computed here for every MHW anomaly profile as being the difference between the shallowest and 150 m depth measurements along the profile such as T stratification =Ta 150m −ta shallowest , S stratification = Sa 150m −Sa shallowest , and ρ stratification = ρa 150m −ρa shallowest .

Fixed Moorings
As part of the systematic and continuous beach monitoring program at SOCIB, Nortek AWACs (Acoustics Wave And Current) have been deployed along the coasts of the Balearic Islands since 2011 (Tintoré et al., 2013), in particular at the Cala Millor and Son Bou beaches ( Figure 1B). Cala Millor is located along the north-eastern coast of Mallorca Island. It is a 1.7 km long micro-tidal semi-embayed anthropized beach. In terms of coastal morphodynamics, it is an intermediate beach with highly variable longitudinal sinuous-parallel bars and troughs, related to wave-driven currents dynamics Gómez-Pujol et al., 2011) and a system of three submerged canyons that merge at 20 m depths. The beach face is open to the east and well exposed to north-northwest to southeast incoming waves (Tintoré et al., 2009;Enríquez et al., 2017). Son Bou beach is located along the southern coast of Menorca Island. It is a 4 km micro-tidal long-shore natural beach that comprises a complex system of small morphodynamic cells separated by natural tombolos, shallow rocky outcrops and a sand-barred lagoon. It is open to the south and well exposed to southwest to southeast incoming waves (Enríquez et al., 2019).
At both beaches, two AWACs have been installed nearshore at 17 m depth, both being equipped with a temperature sensor, measuring hourly near-bottom water temperatures 8 . For the two monitoring site time series, daily means from in situ data are computed over the period 2012-2020. Then, the daily 90th percentile climatologies are built using the 2012-2020 reference period and used for the MHW identification. Although the in situ timeseries are shorter than the recommended 30year minimum for the characterization of MHWs (Hobday et al., 2016), there is a benefit of using in situ data as a complement of satellite products since they provide a more accurate representation of the thermal characteristics in the near-shore environment (Schlegel et al., 2017a). The calculation of MHWs from in situ data using their own climatologies allows quantifying the amount they differ from their localities (Schlegel et al., 2017b).

Marine Heat Wave in Near-Shore Coastal Waters
Local MHWs are identified from the in situ measurement time series at Cala Millor and Son Bou beaches over the period 2012-2020 with respect to the 2012-2020 reference period, following the same methodology as described in Section "Marine Heat Waves." For the joint analysis with satellite products, local satellite data are extracted over 2012-2020 from the closest grid point to the monitoring site and from the off-shore grid point ( Figure 1B). Satellite MHWs are then detected from the near-shore and off-shore time series using their respective 90th percentile climatologies over the 2021-2020 reference period as performed for the mooring.

Daily Monitoring
Through the developed monitoring system, daily SST, SST anomaly (SSTA) and MHS are displayed every day for the last available day providing the last ocean temperature state in the Mediterranean (Figure 3). SSTA and MHS are computed using the daily climatologies of local mean and threshold, respectively, established over the period 1982-2015. The MHS value is equal to 0 when the SST value is smaller than the threshold. The positive value, which indicates the MHS event, is the difference between the SST and threshold values. Time series of daily SST averaged in the different sub-regions of the Mediterranean are displayed for the current year, as well as the daily mean and threshold climatologies. Every day, the time series is updated extending to the last available day. When the SST value becomes higher than the threshold, a yellow shadow is drawn to automatically warn about the MHS event. The year 2020 shows important episodes of MHWs in both the WMed, Adriatic Sea and EMed, with particularly longduration (several months in some sub-regions) and intense magnitude (up to 4 • C) (Figure 4). In the WMed, the subregions experienced very strong and long-duration MHWs in particular in winter, spring and summer 2020. The north and middle WMed suffered MHWs from January/February to May 2020 with one or 2-3 successive events, with at least one lasting more than 60 days (up to 2 • C above the threshold), as well as several events reaching 1-2 • C in end-June/early-July, August and September 2020 ( Figure 4A for the Gulf of Lion [DF2]). The south WMed experienced 40-day duration events in February-March 2020 (up to 0.5 • C), as well as strong events in May, July, and August 2020 (with magnitude higher than 2 • C in the Alboran Sea). In the Tyrrhenian Sea, very long-duration MHWs lasting 68-83 days occurred from January to March 2020, as well as several-week duration events in April/May, August, and September 2020. In the Adriatic Sea and EMed, all sub-regions experienced intense, long and frequent MHW events, with an unprecedented intense and long-duration event in the eastern EMed. The Adriatic Sea, north-eastern and middle-east Ionian Sea also suffered several MHWs especially from January to March 2020 with continuous 70-90-day events and intensity up to 0.5 • C, and several events in summer and autumn 2020 up to 1 • C ( Figure 4B for the north Adriatic [DJ1]). The south Ionian Sea also experienced important MHW events in May (48-day duration), August-September (48-day duration), and October (20-day duration) 2020, while the north-western and middlewestern Ionian Sea suffered MHWs in August-September 2020 during 45 and 34 days, respectively. In the Aegean Sea, very long-duration and intense MHWs occurred from January to March and from September/October to November/December 2020 (lasting 84-99 days with intensity exceeding 1 • C) as well as additional events in spring/summer. Further south, in the Cretan Passage, a very intense event reaching 3 • C above the threshold occurred in May 2020 and a 73-day duration event in September-October 2020. Finally, the Levantine basin suffered a very intense MHW in May (up to 4 • C) as well as a very intense and very long-duration event of 125-130 days from end-August to end-December 2020 ( Figure 4C for the south-east Levantine basin [DL4]).

Annual Indices
Annual MHW characteristics (mean and maximum intensities, mean duration and frequency) are computed from 1982 to 2020. The 2D maps of the MHW indices in the Mediterranean are shown for the year 2020 (Figures 5A-D) and highlight the strong sub-regional variability in the ocean basin, with major differences between the western and eastern sub-basins. Concerning the temperature characteristics, the annual mean intensity in the Mediterranean Sea in 2020 varying between 0.23 and 1.39 • C above the threshold with maxima reached in the Levantine basin, Aegean Sea, north Adriatic Sea, and north WMed (Figure 5A). The annual maximum intensity in 2020 oscillates between 0.57 and 5.16 • C above the threshold with maxima found in the EMed, particularly in south of Crete Island, Rhode gyre, and north Cyprus ( Figure 5B). Regarding the time component, the mean duration averaged in 2020 ranges between 6.1 and 96.5 days with maxima experienced in the Gulf of Lion, north Aegean Sea and Levantine basin ( Figure 5C). The number of events in 2020 varies between 2 and 14 with high numbers found in the WMed (Figure 5D). The total days of MHWs oscillates between 20 and 304, which correspond to 0.5 and 10 months, respectively, over the year 2020 (not shown). Note that the extremely high maximum of total days has been found locally along the Greek coast of the Peloponnese.
The 2D maps of the annual MHW indices averaged over the whole period 1982-2020 (Figures 5E-H) also show a strong spatial variability in the Mediterranean Sea highlighting longitudinal and latitudinal differences in MHW characteristics. The MHW mean and maximum intensities over 1982-2020 (Figures 5E,F) shows values ranging over 0.32-0.89 and 0.88-2.26 • C, respectively, and strong differences between the western and eastern Mediterranean Sea and between the northern and southern sub-basins. For both temperature indices, the maxima are found in the northern WMed (Gulf of Lion and Ligurian Sea) and Adriatic Sea (particularly in the northern sub-region); high values are also detected in the Aegean Sea. The minima are found in the southern EMed (Ionian Sea and Levantine subbasin) for the mean intensity and in the southern Mediterranean (southern WMed, southern-middle Ionian Sea, and southern Levantine sub-basin) for the maximum intensity. The MHW mean duration over 1982-2020 ( Figure 5G) oscillates between 7.7 and 22.6 days with maximum values found in the eastern Mediterranean (Levantine basin) as well as in several coastal areas of the northern WMed, Adriatic Sea and EMed (particularly in the Aegean Sea). The minimum values of MHW duration are identified in the Alboran Sea and southern Ionian Sea. The MHW frequency over 1982-2020 ( Figure 5H) varies between 1.4 and 5 annual events with maximum values in the northern WMed as well as in many places of the WMed, Adriatic Sea, and EMed. The minimum values of MHW frequency are found in the southern EMed. The annual MHW total days averaged over 1982-2020 ( Figure 6A) oscillates between 18 and 51 days with maxima detected in the northern WMed, Adriatic Sea, Aegean Sea, Levantine basin as well as in coastal areas (e.g., along the western, southern and north-eastern Greek coasts) and minima in the Alboran Sea and southern Ionian Sea.
For every MHW parameter, the sub-regionally integrated time series are also computed from 1982 to 2020 to follow their longterm evolution in response to global change. The results are presented and discussed in the next section.

From Seasonal to Decadal Variability
The sub-regional variability from seasonal to interannual and decadal temporal scales is analyzed through Hövmoller representations of the annual MHW indices from 1982 to 2020 as well as through integrated diagrams over the period 1982-2020 (Figures 7, 8). The sub-regional decadal variability of the annual MHW indices is also highlighted comparing the four decades 1982-1991, 1992-2001, 2002-2011, and 2012-2020 and computing the linear trends with 95% confidence level over the entire period 1982-2020 in the different sub-regions. The level of significance and confidence intervals are calculated using the effective number of degrees of freedom, according to the methodology from Chelton (1983). The proposed subperiod division follows the temporal evolution of SST values that are gradually increasing from one decade to the next (Ibrahim et al., 2021) and the increase of the annual MHW total days over the sub-periods (Figures 6B,C). Averaged over the Mediterranean Sea, the mean SSTs from the CMEMS satellite products (described in Section "39-Year Satellite Products") over the decades 1982-1991, 1992-2001, 2002-2011, and 2012-2020 increase with values of 19.06, 19.4, 19.85, and 20.22 • C, respectively. The annual MHW total days averaged over the whole basin also increase from 1982 to 2020 with mean values of 5.7, 12.8, 35.1, and 75.9 days over the four sub-periods, respectively. This decadal increase of SST and MHW total days is also observed in the WMed, Adriatic Sea and EMed sub-basins (not shown). The four sub-periods analysis also allows highlighting the accelerating intensification rate in MHW intensity, duration and frequency. Associated estimations are provided in Tables 2, 3.
Regarding the temperature MHW characteristics during the period 1982-2020, the sub-regionally integrated mean and maximum intensities reach annual mean values up to 1.54 • C (in north WMed in DF3 in 2003) and 3.98 • C (in south WMed in DS1 in 2015), respectively (Figures 7A,E). Both MHW indices have strong temporal variability over 1982-2020 experiencing remarkable events, interannual variability and substantial increase over the period, particularly in recent years. Among the most extreme events, the MHWs which were related to the heat waves in 2003 and 2006 in most of sub-regions (e.g., Sparnocchia et al., 2006;Olita et al., 2007;Garrabou et al., 2009;Hobday et al., 2018;Smale et al., 2019;Holbrook et al., 2020) are well marked. The mean (maximum) intensity is higher than 1 • C (2.5 • C) in most of the WMed, Adriatic Sea and north-middle Ionian Sea in 2003 and in the northern WMed in 2006. In 2013, as also observed by Ibrahim et al. (2021), the EMed experienced particularly intense MHW, in particular the north-eastern Levantine with mean (maximum) intensity larger than 1 • C (2.5 • C). The year 2015 also shows very high mean and maximum intensity values in several sub-regions of the WMed including the Alboran Sea [as also observed in Darmaraki et al. (2019) at basin scale], as the years 2018 and 2019 in the northern WMed, Adriatic Sea and eastern EMed. More recently, in 2020, as addressed above and also observed by Ibrahim et al. (2021) in the EMed, MHWs with maximum intensity exceeding 2.5 • C have been detected in the sub-regions of the WMed (such as the Alboran Sea and Gulf of Lion) and EMed (such as the northern and north-eastern Levantine basin, southern Aegean Sea and Cretan Passage). It is also worth noting that the north Adriatic suffered very intense events with MHW maximum intensity exceeding 2.5 • C many years over the last two decades. Although the temperature MHW characteristics and their changes have strong sub-regional variability with maxima in trend values reached in the Adriatic Sea and northern WMed, all sub-regions are experiencing increase in MHW mean and maximum intensities. The significant trend values with 95% confidence level over the period 1982-2020, which are observed in the Adriatic Sea and most of the WMed and EMed sub-regions, range over 0.24-0.52 and 1.03-2.13 • C, respectively (Table 2 and Figures 7D,H). No significant values are found in MHW mean and maximum intensities in some sub-regions of the EMed (north-eastern, southern and south-eastern Levantine basin, southern Aegean Sea, middle-eastern Ionian Sea) and WMed (Sardinia Channel and Algero-Provençal sub-region) (Figures 7D,H). The increase of temperature characteristics of MHWs has accelerated in recent years. Comparing the four decades, the mean values through the sub-periods increase in the 28 sub-regions of the Mediterranean Sea (Table 2 and Figures 7C,G). More precisely, the sub-regional values of annual mean and maximum intensities are much higher during the two last decades (2002-2011 and 2012-2020) than the two first ones (1982-1991 and 1992-2001) in the WMed, Adriatic Sea and EMed. The subregional maxima of the annual mean intensity are found in  -periods 1982-1991, 1992-2001, 2002-2011, and 2012-2020, as well as over the whole period 1982-2020 and linear trends with 95% confidence level (only) over 1982-2020 of sub-regional MHW indices in the WMed, Adriatic Sea, and EMed.
the WMed and Adriatic Sea during the decade 2002-2011 due to the extremely warm and long-duration event in 2003 in these sub-basins. Otherwise, the last decade is the sub-period when the annual mean intensities are the largest in the EMed and the maximum intensities are the highest in the WMed, Adriatic Sea and EMed. Seasonal variability is also highlighted through estimations of seasonal means over the period 1982-2020 and the four sub-periods in both MHW temperature indices (Table 3 and Figures 7B,F). The seasonal mean and maximum intensities are the largest in spring (with sub-regional values ranging over 0.42-0.84 and 0.84-1.70 • C, respectively) and the smallest in winter (with values over 0.20-0.38 and 0.44-0.76 • C, respectively) in the WMed, Adriatic Sea and EMed. The difference of sub-regional seasonal means between the most recent (2012-2020) and the oldest (1982)(1983)(1984)(1985)(1986)(1987)(1988)(1989)(1990)(1991) decades of the study are higher in summer for the WMed and Adriatic Sea and in spring for the EMed. The maximum differences between the two sub-periods in mean and maximum intensities are found for the EMed in spring with values up to 0.47 and 1.14 • C, respectively.
Concerning the MHW occurrence, the sub-regionally integrated mean duration and frequency reach annual values up to 44.7 days (in EMed in DL3 in 2015) and11.8 events (in EMed in DH2 in 2010) (Figures 8A,E). Both related-time variables show extreme events, interannual variability and significant increase over the period 1982-2020. Among the years with the longest annual mean duration of MHWs exceeding 1 month are the year 2003 in the WMed sub-regions (such as the Balearic Sea and Algerian sub-regions) affected by the extreme event in summer, the year 2015 in the southern Levantine sub-basin, and the recent years from 2018 to 2020 in the eastern EMed, particularly in the Aegean Sea and Levantine basin. The spatial variability is also highlighted through the 28 sub-regional estimates. All sub-regions experienced positive increases over the 39-year period in mean duration and frequency. The MHW mean duration trend values are statistically significant at the 95% level in the Adriatic Sea and most of the sub-regions of the WMed and EMed ranging from 4.8 to 14.9 days over 1982-2020 (Table 2 and Figures 8D,H). There are no significant changes in MHW mean duration in some sub-regions of the EMed and WMed, as for the sub-regional MHW mean and maximum intensity trend values (Figures 8D,H). The sub-regional MHW frequency trend estimations are positively significant at 95% confidence level in all sub-regions of the Mediterranean Sea with values varying from 4.2 to 7.1 annual events over 1982-2020 and maxima reached in the EMed (Table 2 and Figures 8D,H). These estimates lead to a positive trend of sub-regional annual total days of 55.0-125.5 ± 0.2 days over the 39-year period (not shown). The four sub-periods comparison also highlights the increase over the last four decades from 1982 to 2020 and the strong acceleration rate in the last decade when the mean duration and number of events are much higher in all sub-regions than during the three previous decades (Table 2 and Figures 8C,G). The analysis of seasonal variations over the whole period (Table 3 and Figures 8B,F) shows that the MHW duration reaches its maximal sub-regional mean over 1982-2020 in winter in the WMed and Adriatic Sea (11.0 and 10.2 days, respectively) and in autumn in the EMed (12.2 days). The differences of sub-regional mean duration between the 2012-2020 and 1982-1991 sub-periods are the highest in winter in the WMed and Adriatic Sea (10.6 and 12.4 days, respectively) and in autumn in the EMed (17.1 days). The seasonal means over the whole period and differences between the most recent and the oldest sub-periods of discrete numbers of events tend to be the smallest in autumn-winter and the highest in spring-summer with values oscillating over 0.4-0.8 and 0.2-1.8 event per year, respectively.
All annual MHW indices (mean and maximum intensities, mean duration and number of events) are substantially increasing over the period 1982-2020 with positive linear trends which are significant at 95% confidence level in many sub-regions of the Mediterranean Sea (Table 2 and Figures 7D,H, 8D,H). Although these significant increases affect most of the sub-regions, the Mediterranean experiences a strong sub-regional variability in MHW intensity, duration and frequency with significant subregional trends that are the highest in the Adriatic Sea and WMed for the intensity and the largest in the EMed for the duration and frequency. Additionally, these results show an acceleration in the increase of duration and number of events in the last decade (2012-2020) leading to an accelerated increasing of MHW total days in all Mediterranean sub-regions with significant local trend values ranging from 24 to 177 days over the period 1982-2020 and with higher acceleration in the Adriatic Sea and EMed than in the WMed (Figures 6B,C). In the WMed, the highest trends in MHW mean and maximum intensities, mean duration and events are found in the north WMed including the Gulf of Lion [DF2] and Ligurian Sea [DF3]. In the Adriatic Sea, the highest trends in mean and maximum intensities are found in the North Adriatic [DJ1]. The EMed is the place of the highest trends in duration and number of events, particularly in coastal areas (not shown).
These annual indices at sub-regional scale allow following the long-term evolution in intensity and occurrence of MHWs on an annual basis. It is worth mentioning that these estimations are annual and sub-regional averages, thus they include but do not reflect local event (in time and/or in space). Local estimations would be much higher with even more alarming results than the estimates obtained in this study.

DEPTH PENETRATION OF MARINE HEAT WAVES Depth Propagation
The penetration of MHW into the ocean sub-surface is addressed through the analysis of vertical hydrographic profiles from profiling floats in the Mediterranean over the period 2012-2020 (see Section "Marine Heat Wave Profiles"), when the oceanic basin is relatively well-sampled at seasonal and annual scales in the different sub-regions ( Figure 2B). However, the Sicily Strait (DI3) and north Adriatic (DJ1) sub-regions are not addressed in this section due to lack of profiles and insufficient coverage ( Figure 2B).
The medians of the sub-regional distributions of Ta, Sa, and ρa MHW profiles over the period 2012-2020 as a function of depth (Figures 9A,C,E) highlight the vertical penetration of MHWs. Ta from MHW profiles are positive at surface and sub-surface in all sub-regions with higher values at surface (with sub-regional values oscillating between 0.64 and 2.82 • C) decreasing with depth ( Figure 9A). The thermal propagation in depth shows a strong spatial variability over the period. Ta are more intense and reach deeper layer in the EMed than in the WMed. While the medians of the Ta   ( Figure 9E), with negative values higher at surface (ranging from −0.01 to −0.95 kg/m 3 ) than at sub-surface, and reaching 50-300 m depth in the WMed and 100-800 m depth in the EMed.
In order to identify the ocean processes that could be involved in the vertical propagation of MHWs at sub-regional scale, the seasonal distributions of MHW profiles are also examined. The uneven spatio-temporal coverage of the profiling float network is challenging when analyzing such data ( Figure 2B). Some caution has to be taken, in particular for the study at seasonal scale with a reduced number of profiles in the sub-regions (Figures 10D,H,L,P).
The seasonal distributions of MHW profile anomalies show seasonal variations in the vertical extent of MHWs (Figure 10). The Ta of MHW profiles have a strong seasonal variability in its vertical structure (Figures 10A,E,I,M) Figure 2B) and the northeastern corner of the Crete Island are areas of downwelling during the summer and autumn seasons (Bakun and Agostini, 2001) that may contribute to the vertical propagation. Downwelling processes in late-autumn/early-winter could also explain the depth penetration until 200 m in the African side of the Alboran Sea (Bakun and Agostini, 2001). The Sa of MHW profiles have not well-marked seasonality although variations are found (Figures 10B,F,J,N). The ρa vertical distributions also follow the Ta patterns at seasonal scales (Figures 10C,G,K,O), with more intense values in summer and autumn (with maximum median values of −1.16 and −1.43 kg/m 3 , respectively) in the upper 30 and 100 m layers, respectively, and with deeper vertical extent in intermediate and deep layers in winter, when vertical mixing drives the weak stratification enabling MHWs to extend deeper than the surface mixed layer (Schaeffer and Roughan, 2017). Downwelling could also contribute to the vertical propagation of MHWs (as in the southern coast of the Mediterranean) while upwelling would prevent them (such as in Gulf of Lion, eastern Aegean Sea, and off the west coast of Greece). The sub-regions where various upward and downward ocean processes may occur (e.g., convection, downwelling, and upwelling) would require further analyses at local scale.

Ocean Stratification Changes
The changes in T, S, and ρ upper-ocean stratifications are characterized here by the difference between the shallowest and 150 m depth changes (Figures 9B,D,F) over the whole period 2012-2020 and per season over the period. As observed above (Figures 9, 10), most of the T, S, and ρ changes in MHW profiles occur in the upper 150 m layer. Thus, vertical gradients between the surface and 150 m depth appear appropriate to assess the impact of MHWs on the upper-ocean stratification changes in the Mediterranean Sea. The temperature stratification is negative (with sub-regional values over the period varying between −0.27 and −1.4 • C) and associated with positive Ta for all seasons and sub-regions ( Figure 9B). The T increase is higher at surface than in sub-surface leading to the temperature stratification increase (with higher sub-regional values in summer and autumn reaching −3.50 and −3.70 • C, respectively). In the EMed (WMed), the salinity stratification over 2012-2020 is generally negative (positive) reaching sub-regional maximum of −0.26 (0.17) with associated positive (negative) Sa (Figure 9D). The S increase (decrease) in EMed (WMed) is higher at surface than in subsurface reducing (enhancing) the salinity stratification. The density stratification, which is induced by both temperature and salinity stratifications, is positive (with sub-regional values over 2012-2020 ranging over 0.02-0.39 kg/m 3 ) and combined with negative ρa values for all seasons and sub-regions (Figure 9F), indicating that the density decrease is higher at surface than at 150 m depth enhancing the upper-ocean stratification. As for the temperature stratification, the seasonal upper-ocean density stratifications are higher in summer and autumn with values reaching 0.96 and 0.95 kg/m 3 , respectively.
In the WMed, the surface T intensification is associated with S decrease of which the weak wind stress may be the major cause of evaporation decrease through latent heat-flux . The stronger T increase and S decrease at surface have a cumulated effect on the ρ decrease. In the EMed, the T increase is associated with S increase that may be induced by reduction in fresh water inputs (e.g., precipitation or wind conditions). At surface, the S increase does not compensate the T increase effect on ρ decrease. In both WMed and EMed, the ρ decrease induced by MHWs is more important at surface than below enhancing the ocean stratification, in particular in summer and autumn, which correspond to periods of preconditioning phase for deep convection . The decrease in density of surface and intermediate layers as well as the related increase of ocean stratification hinder any further outflow sinking (Schröder et al., 2006). Such changes modify the hydrographic properties of the Mediterranean Outflow Water flowing into the Atlantic Ocean (Somot et al., 2006) of which heat and salt contents are in turn modified (Curry et al., 2003) impacting its thermohaline circulation.

COASTAL MARINE HEAT WAVES Coastal Near-Shore Response
The ocean response to global warming in the coastal nearshore waters is addressed in the Balearic Islands using the near-shore fixed moorings at Cala Millor (Mallorca Island) and Son Bou (Menorca Island) beaches over the period 2012-2020. MHWs have been identified (as described in Section "Near-Shore Coastal Marine Heat Waves") in the 2012-2020 time series of temperature from the fixed mooring as well as from the nearest satellite grid point to the monitoring site and the satellite off-shore point ( Figure 1B). The near-shore satellite points are located at 1.4 and 0.1 km from the Cala Millor and Son Bou moorings, respectively. The off-shore points, which are in deep waters at a nearly equivalent distance from the continental shelf-break as their corresponding near-shore point, are situated at 43.5 and 15.7 km from the Cala Millor and Son Bou near-shore points, respectively. MHW indices (mean and maximum intensities, mean duration, count and total days) are estimated for the mooring, near-shore and off-shore satellite time series over the common monitored period from 2012 to 2020 ( Table 4).
The in situ temperature time series show the marked seasonal variability with values varying from 12.57 to 27.65 • C at Cala Millor and from 12.85 to 27.78 • C at Son Bou (not shown, Table 4). The remote sensed temperature time series at the nearshore satellite grid point display warmer temperatures than in situ records ranging from 12.81 to 28.44 • C and from 12.76 to 28.3 • C in the two monitoring sites, respectively (not shown, Table 4). At Cala Millor and Son Bou, the daily mean time series of in situ and near-shore satellite SST are highly correlated (0.97 and 0.96, respectively) with a mean absolute error of 0.75 and 0.82 • C, respectively. Differences between datasets are expected since in situ measurements are recorded at 17 m depth, while satellite data provide sea surface estimations; in this sense, deviation of values along the water column may occur, particularly during warmer seasons when the ocean is more stratified (Guillén et al., 2018). Moreover, in situ data are daily averaged from hourly measurements, including the diurnal fluctuations while satellite data are night-time images (Buongiorno Nardelli et al., 2013). Finally, additional differences are intrinsically related to the instrumentation and data processing (e.g., sensor type, calibration, quality-control procedure, correction or postprocessing). The temperature time series at the off-shore satellite point have more amplitude than at the near-shore point with values oscillating between 12.73 and 28.8 • C at Cala Millor and between 12.79 and 28.36 • C at Son Bou over the period 2012-2020. The daily mean time series of near-shore and off-shore Co-occurrence (%) of MHW events between the mooring and near-shore satellite points, and between the near-shore and off-shore satellite points, for both sites of study.
Frontiers in Marine Science | www.frontiersin.org satellite data are fully correlated with a mean absolute error of 0.21 and 0.13 • C at Cala Millor and Son Bou, respectively. At Cala Millor, the total days of detected MHWs over 2012-2020 are higher in the satellite surface data (192 near-shore and 197 off-shore) than in the mooring data at 17 m depth (149) from which fewer events are detected. While the mean intensities of MHWs detected over 2012-2020 at the satellite near-shore and off-shore points have values of around 0.36 and 0.34 • C at surface, respectively, the mean intensity of MHWs identified in the mooring at sub-surface has a smaller value of 0.30 • C. At surface, the maximum intensity of the MHWs detected with the satellite from 2012 to 2020 reaches a warmer value near-shore (1.80 • C) than off-shore (1.30 • C). In the nearshore environment, the maximum intensities reach values of 1.80 and 1.77 • C in satellite surface and mooring sub-surface data, respectively. At Cala Millor, in the near-shore area, the satellite surface point experienced higher MHW total days and maximum intensity over 2012-2020 than the mooring subsurface point. At Son Bou beach, the MHW total days are higher in the mooring in sub-surface (184) than in the satellite surface data (162 near-shore and 174 off-shore). The mean intensity of the MHWs detected at surface from the satellite near-shore and off-shore data (0.39 • C) is warmer than the mean intensity of the MHWs identified at sub-surface in mooring data (0.37 • C). The MHW maximum intensity reached over 2012-2020 is warmer at the mooring (1.82 • C) than at satellite grid points (1.33 • C near-shore and 1.31 • C off-shore). At Son Bou beach, the sub-surface near-shore waters as detected by the mooring show higher MHW total days and maximum intensity than the surface near-shore and off-shore waters as observed by satellites.
At both sites of study, the near-shore satellite point displays MHWs that reach higher maximum intensity than the offshore point, among their respective detected events. The differences between near-shore and off-shore thermal signals suggest that further physical processes that are not involved in large-scale dynamics should be considered to understand the evolution of coastal MHWs, such as local wave and wind action as well as coastal orientation and configuration (e.g., embayment, headland, and bathymetry) that introduce eddies and perturbations to global circulation patterns (Schumann et al., 1995;Schlegel et al., 2017a). Moreover, the MHW characteristics, as detected in the near-shore environment, differ from the mooring location to the closest satellite point. As mentioned above, the difference between in situ data and satellite product may be attributed to the difference in datasets as well as in monitoring depths. Additionally, the thermal behavior of the coastal near-shore waters is different in the two beaches. At Cala Millor (Son Bou) beach, the total days of MHW events and the maximum intensity reached over the period 2012-2020 are higher (smaller) in the near-shore surface waters as detected by satellite than in the sub-surface waters as detected by the mooring ( Table 4). The difference of behavior between the two sites may be attributed to specific local features and physical processes of the coast. In particular, the morphodynamic differences between both study sites and their exposure to different wind and wave at both large and local scales are potentially the drivers of such differences. Thus, the coastal nearshore response to MHW in sub-surface would strongly depend on the location. Figure 11 shows the time-span along the year of MHW at mooring, near-shore and off-shore satellite points at both Cala Millor and Son Bou beaches for single years. At first sight, there is not a clear seasonal variability on MHW events at the different locations in terms of duration and months of detection. Additionally, the occurrence (frequency and length) in MHW events varies between the different locations. To determine the relation of detected MHW events between the coastal mooring, near-shore and corresponding off-shore satellite points, the cooccurrence of events is computed for each study site. The percentage of co-occurrence is calculated for the study timeinterval (Table 4) in such a way that two events at different locations co-occur when they coincide in time at least 20% of the duration of the longest event. This minimum percentage of overlapping allows implicitly including the MHW duration. The values stand for the intercomparison between the off-shore and near-shore points, and between the near-shore and mooring points for each area of study.

Co-occurrence
The co-occurrences of registered events decrease when approaching the coast with values at Cala Millor and Son Bou between the offshore and near-shore points of 54 and 70%, respectively, and between the near-shore and mooring points of 33 and 40%, respectively. The decreasing concomitance toward the coast for both sites enhances the relevance of local scale processes on the evolution and propagation of MHWs in coastal areas. As aforementioned, the difference in MHW co-occurrences between the near-shore and off-shore locations as well as inside the near-shore environment may be related to the nature and variability of physical ocean processes, the near-shore water circulation being subject to the influence of morphodynamic features (e.g., bathymetry, shallow waters, coasts, and lands) and local response to external forcing (e.g., wind, wave, and current). At this stage, the atmospheric and oceanic conditions of cooccurrence and no co-occurrence as well as their seasonality have not been identified and would require further analysis supported by long-term observations and numerical simulations.

Discussion
In this section, the MHWs have been are identified using the recent 2012-2020 period reference due to the monitored period of the mooring data. Since SST and MHW events have substantially increased from 1982 to 2020 with a clear acceleration in recent years (as introduced in Section "Sub-Regional Approach" and shown in "Sub-Regional Marine Heat Waves"), the intensity, duration and frequency of MHWs as presently studied in coastal near-shore waters are probably underestimated. Indeed, the threshold climatologies established over the period 2012-2020 differ from the ones over 1982-2015 with warmer values (not shown). Consequently, fewer MHWs are detected using the 9-year climatology and have underestimated temperature and time characteristics. More specifically, at the satellite near-shore and off-shore points, the use of the 2012-2020 FIGURE 11 | Time-span of MHW events from the mooring, near-shore and off-shore satellite points for single years for Cala Millor (blue colors) and Son Bou (gray colors) beaches.
period as reference compared to the 1982-2015 period induces underestimations by 35-42, 33-43, 38-56, and 38-51% for the mean and maximum intensities, mean duration and number of events, respectively. Based on the most intense and longest MHW events given the use of a short-term climatology based on recent years, this section suggests the complex ocean behavior in the near-shore environment, where local effects strongly impact the near-shore ocean response to extreme event.

CONCLUSION AND PERSPECTIVES
In a warming world, society is facing major climate-related challenges, risks and disasters. Among them, marine heat wave (MHW) events have devastating impacts on ecosystems, threaten economies that rely on them, and strengthen severe storms and flooding by warming the ocean. In the whole Mediterranean Sea, MHWs are increasing substantially in intensity, duration and frequency, with strong spatial variability highlighting the sub-regional response to global warming. All sub-regions are experiencing positive trends of all MHW characteristics with unprecedented intensification in recent years. The sub-regional trend estimations with 95% confidence for the MHW mean and maximum intensities, mean duration and event number oscillate over 0.06-013 • C/decade, 0.26-0.55 • C/decade, 1.23-3.82 days/decade and 1.1-1.8 events/decade, respectively.
The most intense and shallow events preferentially occur in summer and autumn, when the water column is the most stratified suggesting predominant air-sea fluxes (anomalous solar radiation), and high values extend at depth in winter when wind-induced vertical mixing weaken the stratification. The analysis of vertical hydrographic profiles during MHW events has showed the depth penetration of such events into the ocean interior as well as its sub-regional and seasonal variations. The upper-ocean thermal stratification is increasing due to the surface intensification of the warming signal leading to strengthen the upper-ocean density stratification. In the WMed, the salinity decrease also contributes to the density decrease in the upper layer while in the EMed, the salinity increase does not counteract the temperature increase. Thus, all sub-regions of the Mediterranean are experiencing strong surface density decrease. The resulting strengthening in upper-ocean density stratification will act as a barrier between the surface and deep water masses -modifying the ocean circulation and preventing the ocean ventilation -, and between the nutrient-poor surface and nutrient-rich sub-surface waters -impacting directly on the oceanic biological activities.
Moreover, the ocean response to climate warming also strongly varies from coastal areas to open ocean. A case study along the coasts of the Balearic Islands has showed the higher occurrence and maximum intensity in near-shore than off-shore areas. It has also highlighted that the local surface and sub-surface responses in near-shore environments may differ and strongly depend on the local features.
Giving access to relevant information on ocean state and variability at sub-regional scale through the development of tailored products to users will support the preparation for and the response to climate change, in particular in coastal zones (Juza and Tintoré, 2021). To support the achievement of international initiatives (e.g., UN SDG14 9 and UN Decade of Ocean Science for Sustainable Development) and to respond to both science and society needs, SOCIB has implemented a visualization tool to detect extreme warm temperature events in real-time and to follow their long-term evolution in response to global warming in the Mediterranean Sea. The event-based monitoring at sub-regional scale is a very valuable tool that can offer targeted information for diverse ocean stakeholders (e.g., general public, recreation, aquaculture, environmental agencies, and policy decision-makings) in real-time using observations or with short-term predictability using model forecasts, such as the 10-day forecasts from the Copernicus Marine Service Mediterranean Sea analysis and forecast model (Clementi et al., 2021).
Providing four-dimensional and continuous datasets, regional and coastal numerical simulations in the Mediterranean Sea (e.g., Copernicus Marine Service MED-MFC in the Mediterranean, Simoncelli et al., 2019;Clementi et al., 2021; SOCIB WMOP in the western Mediterranean, Juza et al., 2016b) are valuable data products that could help to better understand the drivers of MHWs and ocean response to warming signals, and thus to improve their predictability. They could complement observations in poorly or unevenly sampled dynamical regions (e.g., deep convection areas, upwelling and downwelling zones, permanent currents or mesoscale eddies) to better characterize the ocean processes involved in the MHW propagation (both vertical extent and horizontal advection), to estimate their contributions and thus to enhance the predictability depending on the sub-regional dynamics. This model-based analysis could be supported by observational platforms deployed in specific dynamical areas such as the fixed mooring and gliders in the deep convection area in Gulf of Lion (MOOSE network 10 ; Testor et al., 2018), the SOCIB gliders in the Ibiza Channel (Heslop et al., 2012) or the fixed mooring in the Sicily Strait (Schroeder et al., 2017).
The coastal ocean is a crucial area [as internationally recognized through the Coastal Working Group of EuroGOOS (El Serafy, 2020) or the CoastPredict programme 11 as a contribution to the UN Decade of Ocean Science] which is highly sensitive to climate change with major societal, economic and environmental impacts. Given the ocean response to global warming from sub-regional to local scales, long-term time series from additional in situ platforms such as coastal fixed moorings (e.g., the network operated by Puertos del Estado 12 ) and repeated cruise fixed stations (e.g., the RADMED monitoring 9 https://un.org/sustainabledevelopment/sustainable-development-goals/ 10 www.moose-network.fr/ 11 https://www.coastpredict.org/ 12 www.puertos.es program 13 along the Iberian coast) could provide long-term based climatologies and relevant information at local scale in coastal areas. Additionally, in complement to satellite products, high temporal resolution in situ measurements at local scale could provide a more accurate representation of the thermal signal near the coast, capturing diurnal fluctuations and local-scale deviations that affect coastal ecological processes (Leichter et al., 2006;Castillo and Lima, 2010).
Better understanding the ocean behavior in response to extreme events, improving their prediction, fostering their monitoring and facilitating access to ocean information from local to sub-regional scales will help to implement adaptive management, to establish adaptation strategy and to support the marine conservation and sustainable management of the oceans in a warming world.

DATA AVAILABILITY STATEMENT
The datasets used in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
MJ contributed to the conception and design of the study, data processing, interpretation of the data, scientific discussion, and writing of the manuscript. ÀF-M contributed to data processing, interpretation of the data, scientific discussion, and writing section. JT contributed to scientific discussion and revision of the manuscript. All authors contributed to the article and approved the submitted version.

FUNDING
Part of this work was supported by the EuroSea project (www.eurosea.eu), which has received funding from the European Union's Horizon 2020 Research and Innovation Programme under grant agreement No. 862626.