Spatial Variability and Trends of Marine Heat Waves in the Eastern Mediterranean Sea over 39 Years

: Marine heatwaves (MHWs) can cause devastating impacts on marine life. The frequency of MHWs, gauged with respect to historical temperatures, is expected to rise signiﬁcantly as the climate continues to warm. The MHWs intensity and count are pronounced with many parts of the oceans and semi enclosed seas, such as Eastern Mediterranean Sea (EMED). This paper investigates the descriptive spatial variability and trends of MHW events and their main characteristics of the EMED from 1982 to 2020 using Sea Surface Temperature (SST) data obtained from the National Oceanic and Atmospheric Administration Optimum Interpolation ([NOAA] OI SST V2.1). Over the last two decades, we ﬁnd that the mean MHW frequency and duration increased by 40% and 15%, respectively. In the last decade, the shortest signiﬁcant MHW mean duration is 10 days, found in the southern Aegean Sea, while it exceeds 27 days off the Israeli coast. The results demonstrate that the MHW frequency trend increased by 1.2 events per decade between 1982 and 2020, while the MHW cumulative intensity (i cum ) trend increased by 5.4 ◦ C days per decade. During the study period, we discovered that the maximum signiﬁcant MHW SST event was 6.35 ◦ C above the 90th SST climatology threshold, lasted 7 days, and occurred in the year 2020. It was linked to a decrease in wind stress, an increase in air temperature, and an increase in mean sea level pressure.


Introduction
In many parts of the world, heatwaves are growing increasingly common and prevalent. It creates significant climatic extremes in oceanic and atmospheric systems, which can have catastrophic ecological impacts and huge socioeconomic effects [1][2][3]. The study of Marine Heat Waves is important for understanding and predicting climate change and its impact on ocean ecosystems. Marine heat wave (MHW), according to [4], is "a prolonged discrete anomalously warm water event that can be described by its duration, intensity, rate of evolution, and spatial coverage". It was suggested that a linear classification scheme be used by [5], where MHWs are defined based on temperature exceedance from local climatology, which will impact the classification. The probability of the occurrence of MHWs is impacted by SST climatological values as described in [2,[6][7][8]. Changes and variability in SST climatological values, as described in [7,8], will have an impact on the definition and probability of occurrence of MHW. These extreme events lead to a widespread mortality of benthic invertebrates [9] and loss of seagrass meadows [10]. Over the last four decades, the international oceanographic community, regulators and policy makers have paid close attention to the EMED, namely because in this region variations in processes of relevance for the global climate change take place over relatively short space and time scale [11][12][13].
In the summer of 2003, one of the first-detected MHWs occurred in the Mediterranean which lasted over a month due to an increase of air temperature, a reduction of wind stress,

Data and Methods of Analysis
The Eastern Mediterranean (EMED) region includes the Ionian and the Levantine basins with two other marginal seas; the Adriatic and the Aegean as shown in Figure 1a. The EMED reaches a depth of over 3000 m off Mersa Matruh at the western part of Egypt and 4000 m to the southeast of Rhodes [20,21] see Figure 1b. The mean depth of the EMED basin is about 2000 m. The shelf is usually extending for no more than 8 km, except off the Nile Delta, where the 200 m contour is traced at 60 km offshore, in the Sicily Strait and the Adriatic Sea. The most important islands in the EMED are Sicily, Malta, Crete, Rhodes, and Cyprus. The study area extending from 30 • to 42 • N and from 13 • to 36.5 • E is shown in Figure 1a,b. This area is covering most of the EMED. In this section the authors will describe the data and methods of analysis used to obtain the results.

Sea Surface Temperature Data
The daily SST data used were obtained from the National Oceanic and Atmospheric Administration Optimum Interpolation Sea Surface Temperature from 1 January 1982 to 31 December 2020 ([NOAA] OI SST V2.1; [22]). The data are freely available on the NOAA website https://www.ncei.noaa.gov/data/sea-surface-temperature-optimum-interpolation/ v2.1/access/avhrr/ (accessed on 10 January 2021). This data set is an interpolation of remotely sensed SSTs from the Advanced Very High Resolution Radiometer (AVHRR) imager into a regular grid of 0.25 • and daily temporal resolution. AVHRR OI data for the EMED Sea were extracted from the global data, providing a 2624-point regularly gridded dataset spanning 14,245 days from 1 January 1982 to 31 December 2020.

Trend and Interannual Variabilty of SST
To demonstrate the spatiotemporal variability in SST over the EMED Sea, an EOF is calculated using the singular value decomposition (SVD) approach [23]. EOF analysis extracts the main mode of SST variability, providing us with a spatial pattern and a timevarying index known as the principal component (PC). The EOF analysis was applied to gridded monthly SST anomalies from 1982 to 2020. As described in [24], the SST anomalies were calculated by subtracting the historical climatological mean value at each grid point from the SST values at the same location (grid). Here, we are interested in interannual and longer time scales variability. So, prior to EOF decomposition, the seasonal cycle and linear trends were removed from the monthly means of SST anomalies at each grid point, and then the SST were normalized by dividing each point time series by its standard deviation. The mean seasonal cycle was estimated by calculating the mean monthly values for each calendar month at each grid point from 1982 to 2020. To obtain de-seasoned monthly values, the mean value of each month was subtracted from all corresponding months in all years. Linear trends of de-seasoned monthly SST anomaly are estimated by using the least squares method [25]. To assess the statistical significance of these trends, the original two-tailed modified Mann-Kendall test at 95% confidence interval was used [26,27]. The authors used a thirteen-month running mean to remove the higher frequency variability in order to highlight the interannual variability, as described in [28]. Refs. [29,30] provide examples of similar data filtration methods. This paper describes the main characteristics of MHWs in the EMED. The authors detect and analyze the most intense MHW events observed in the EMED during the study period, as well as their relationship to atmospheric conditions.

Data and Methods of Analysis
The Eastern Mediterranean (EMED) region includes the Ionian and the Levantine basins with two other marginal seas; the Adriatic and the Aegean as shown in Figure 1a. The EMED reaches a depth of over 3000 m off Mersa Matruh at the western part of Egypt and 4000 m to the southeast of Rhodes [20,21] see Figure 1b. The mean depth of the EMED basin is about 2000 m. The shelf is usually extending for no more than 8 km, except off the Nile Delta, where the 200 m contour is traced at 60 km offshore, in the Sicily Strait and the Adriatic Sea. The most important islands in the EMED are Sicily, Malta, Crete, Rhodes, and Cyprus. The study area extending from 30° to 42° N and from 13° to 36.5° E is shown in Figure 1a,b. This area is covering most of the EMED. In this section the authors will describe the data and methods of analysis used to obtain the results.

Sea Surface Temperature Data
The daily SST data used were obtained from the National Oceanic and Atmospheric Administration Optimum Interpolation Sea Surface Temperature from 1 January 1982 to 31 December 2020 ([NOAA] OI SST V2.1; [22]). The data are freely available on the NOAA It is crucial for our study to determine the warmest period from the whole climatological SST time series, NOAA OI SST V2.1. To reach our goal we divided the SST time series into four periods, (P 1 , 1 January 1982-31 December 1990), (P 2 , 1 January 1991-31 December 2000), (P 3 , 1 January 2001-31 December 2010) and (P 4 , 1 January 2011-31 December 2020). The chosen of the first period P 1 is based on satellite NOAA OI SST V2.1 data availability. Then daily climatology was calculated from the P 1 , P 2 , P 3 and P 4 years employed in the analysis and smoothed with a monthly centred-moving mean.

Marine Heat Waves (MHWs) Calculations
The authors used [5] MHWs, definition, as a discrete prolonged anomalously warmwater events in a location, using the same daily SST database for the 1982-2020 period. The MATLAB toolbox [31] has been used to identify periods of time where temperatures were above the 90th percentile threshold relative to the climatology, and for a duration of at least 5 days. The authors identified the start and end dates and calculated the duration and intensity for each MHW event according to [5]. The climatological mean has been calculated from the daily Sea Surface Temperature (SST) data from 1 January 1982 to 31 December 2020 ([NOAA] OI SST V2.1; [22] following the [5] formula, where, T m is the climatological SST mean in • C calculated over a period from 1 January 1982 to 31 December 2020, to which all values are relative, j is the day of year, y s and y e are the start and end of the climatological base period respectively, and T is the daily SST on day d of year y. Then the 90th percentile threshold based on SST climatology period, T 90 (j) in • C are calculated according to [5] based on climatological SST mean (X) in • C over a period from 1 January 1982 to 31 December 2020 as follow: T 90 (j) = P 90 (X) According to [5] definition, MHWs occur when SSTs exceed the threshold, defined as the 90th percentile of SST variations (T 90 (j)) based on a 39-year climatological period , for at least five consecutive days. At each location and for each MHW, we calculated the event duration (time between start and end dates) and mean cumulative intensity. Associated annual statistics were then calculated, including the frequency of events (i.e., the number of discrete events occurring in each year), where t e, t s : dates on which a MHW begins and ends. And D is the MHW Duration in Days that temperature exceeds the threshold, D = t e − t s . Where P 90 is the 90th percentile and X = {T(y,d)|y s <= y <= y e , j − 5 <= d <= j + 5} Then, the MHWs main characteristics are (i max ) in • C, mean intensity (i mean ) in • C and (i cum ) in ( • C Days): i mean = T(t) − Tm(j)During the MHW eventt s ≤ j ≤ t e and j(t s ) ≤ j ≤ j(t e ) (6) i cum : sum of daily intensity anomalies.
According to [5] during the MHW event: T(t) > T 90 (j), t s is the time, t, while during time (t − 1) presents no MHW event T(t − 1) < T 90 (j), t e is the time, t, where t e > t s T(t) < T 90 (j)

Atmospheric Variables Data
Atmospheric variables were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 [32]. Winds at 10 m elevation (U 10 and V 10 ), air temperature at 2 m elevation, and pressure at mean sea level were used. The dataset has 0.25 • spatial resolution with hourly temporal step. Daily means were computed by averaging hourly data, and a daily climatology was built in the same manner as for SST, using the same period from 1982 to 2020. The wind stress components in east and north directions (τ wx , τ wy ) were calculated according to [33] in (N/m 2 ), where ρ a = 1.2 (kg/m 3 ) is the density of air, C D is the drag coefficient in Equation (8) is computed according to [34] while (W) is the wind speed component (m/s). The authors have computed the wind stress components to track [14][15][16][17][18][19][20] May 2020 MHW events which have the maximum intensity over the whole study period (1982-2020). Figure 2 shows the temporal evolution of daily mean climatology SST ( • C) for P 1 , P 2 , P 3 and P 4 . The figure clearly shows that P 1 has the lowest climatological SST values, with a maximum of about 26.01 • C, a minimum of about 15.07 • C, and a mean value of 19.98 • C (see Table 1). During P 2 , the highest climatological SST was approximately 26.38 • C, and the lowest was approximately 15.25 • C, with a mean of 20.12 • C. The maximum, minimum, and mean statistical values for P 3 were 26.52 • C, 15.44 • C, and 20.48 • C, respectively. During P 4 , the maximum value reached 27.07 • C and the minimum was 16.76 • C with a mean value of 20.92 • C. In summary, the maximum values over the four decades were detected during P 4 in agreement with [19,35,36]. Figure 2 shows that SST values are gradually increasing from one decade to the next. The EMED basin average climatological SST values in Table 1 show that P 4 was the warmest decade by about 0.94 • C within the entire 1982-2020 study period. The authors used anomaly maps to clearly identify cooling and warming regions over the aforementioned time periods.

Description of the SST Data (1982-2020)
As shown in Figure 3a, the P1 SST anomaly map shows a cooling range of −1.00 to −0.60 °C that covers the entire EMED. The P2 east of Crete is cooling to cc. −0.25. (see Figure  3b). Figure 3c shows that the cooling rate for P3 has decreased. The cooling range was between −0.50 to −0.10 °C. The P4 anomaly patterns revealed a warming of about 0.10 to 0.50 °C, over the northern EMED basin except some parts of the north Ionian Sea. This region, which is known as a Bimodal Oscillating System (BiOS), is responsible for decadal reversals of the Ionian basin-wide circulation as described by [37][38][39][40][41]. It is a decadal circulation periodic reversal region from cyclone to anticyclone and vice versa as described by [41]. The BiOS is regarded as a pattern of EMED climate change as mentioned by [40]. In general, the northern EMED basin's SST anomaly for P4 is greater (by 0.25 to 1.00 °C) than the southern one. As shown in Figure 3d, the maximum SST anomaly for P4 was 0.50 °C in the northwest Aegean Sea, while the minimum value was −0.50 °C along the southwestern coast of EMED. This may be attributed to the fact that, the northern part of the EMED basin is occupied by a cyclonic motion such as Rhodes cyclonic gyre as mentioned in [42], while the southern part is dominated by anticyclonic gyres such as Mersa Matruh and Shikmona anticyclonic gyres as described in [37][38][39][40][41][42][43]. Hence, warming in anticyclonic areas is less than in the cyclonic ones due to strong stratification in anticyclonic regions [44,45].
Overall, the patterns of SST anomaly for the four periods show that P4 is the warmest period across the study area. These information highlight some specific hotspot areas with the highest SST anomaly during P4, such as the Aegean, north Levantine, and north Ionian. The authors used anomaly maps to clearly identify cooling and warming regions over the aforementioned time periods.
As shown in Figure 3a, the P 1 SST anomaly map shows a cooling range of −1.00 to −0.60 • C that covers the entire EMED. The P 2 east of Crete is cooling to cc. −0.25. (see Figure 3b). Figure 3c shows that the cooling rate for P 3 has decreased. The cooling range was between −0.50 to −0.10 • C. The P 4 anomaly patterns revealed a warming of about 0.10 to 0.50 • C, over the northern EMED basin except some parts of the north Ionian Sea. This region, which is known as a Bimodal Oscillating System (BiOS), is responsible for decadal reversals of the Ionian basin-wide circulation as described by [37][38][39][40][41]. It is a decadal circulation periodic reversal region from cyclone to anticyclone and vice versa as described by [41]. The BiOS is regarded as a pattern of EMED climate change as mentioned by [40]. In general, the northern EMED basin's SST anomaly for P 4 is greater (by 0.25 to 1.00 • C) than the southern one. As shown in Figure 3d, the maximum SST anomaly for P 4 was 0.50 • C in the northwest Aegean Sea, while the minimum value was −0.50 • C along the southwestern coast of EMED. This may be attributed to the fact that, the northern part of the EMED basin is occupied by a cyclonic motion such as Rhodes cyclonic gyre as mentioned in [42], while the southern part is dominated by anticyclonic gyres such as Mersa Matruh and Shikmona anticyclonic gyres as described in [37][38][39][40][41][42][43]. Hence, warming in anticyclonic areas is less than in the cyclonic ones due to strong stratification in anticyclonic regions [44,45].  Figure 4 depicts a linear trend map after removing SST seasonality from 1982 to 2020. The original two-tailed modified Mann-Kendall test [26,27] was applied to the estimated trend at each grid point, the results indicate a statistically significant trend at the 95% confidence interval was observed over the whole region. High spatial variability of SST linear trend was observed over the EMED, ranging from 0.1 to 0.5 °C per decade. The basin average warming rate was about 0.33 ± 0.04 °C per decade. The higher values were found in the Aegean and Levantine Basin, while the lower values were observed in the Ionian and Tyrrhenian Sea.

Trend and Inter-Annual Variability of SST
The first three EOF modes of SST anomaly explain about 71% of the overall de-seasoned variance of the data. We only discuss the first two modes of EOF since the third mode describes less than 5% of the variance.
The spatial pattern of the first EOF1 mode, which accounts for approximately 50% of the overall deseasoned variance, is a positive anomaly over the entire region of the EMED (Figure 5a), suggesting an in-phase oscillation of the entire basin around the steady-state mean. The Ionian Sea and the Mersa Matruh gyre in the west of the Egyptian coast had the highest SST variability, while the Aegean Sea had the lowest variability, owing to the presence of strong and persistent Etesian winds over the Aegean Sea [35] and the input of Black Sea cold water inflow into the Aegean Sea [28]. The first temporal mode (PC1) reached its peak in the summer of 2003 (Figure 5b), which was the warmest summer over the study period, due to the strongest heat waves that hit the northern Europe. The highest negative peaks were observed in the winters of 1987 and 2019 in agreement with [46]. The lower peak values (i.e., cold periods) in the 1990s could be attributed to Mount Pinatubo's eruption in the Philippines, which caused a cold air-temperature anomaly throughout the Middle East during the winter of 1992 As mentioned by [47], The spatial pattern of the second EOF2 mode (Explains 16% of the total variance) demonstrated a dipole in the study area (Figure 5c), with positive anomalies in the Ionian Sea and negative ones in the Levantine and Aegean basins. The temporal coefficient of the second mode (PC2) alternates between positive and negative patterns over the study period 1982-2020. Overall, the patterns of SST anomaly for the four periods show that P 4 is the warmest period across the study area. These information highlight some specific hotspot areas with the highest SST anomaly during P 4 , such as the Aegean, north Levantine, and north Ionian.
Trend and Inter-Annual Variability of SST Figure 4 depicts a linear trend map after removing SST seasonality from 1982 to 2020. The original two-tailed modified Mann-Kendall test [26,27] was applied to the estimated trend at each grid point, the results indicate a statistically significant trend at the 95% confidence interval was observed over the whole region. High spatial variability of SST linear trend was observed over the EMED, ranging from 0.1 to 0.5 • C per decade. The basin average warming rate was about 0.33 ± 0.04 • C per decade. The higher values were found in the Aegean and Levantine Basin, while the lower values were observed in the Ionian and Tyrrhenian Sea.
The first three EOF modes of SST anomaly explain about 71% of the overall deseasoned variance of the data. We only discuss the first two modes of EOF since the third mode describes less than 5% of the variance.
The spatial pattern of the first EOF1 mode, which accounts for approximately 50% of the overall deseasoned variance, is a positive anomaly over the entire region of the EMED (Figure 5a), suggesting an in-phase oscillation of the entire basin around the steady-state mean. The Ionian Sea and the Mersa Matruh gyre in the west of the Egyptian coast had the highest SST variability, while the Aegean Sea had the lowest variability, owing to the presence of strong and persistent Etesian winds over the Aegean Sea [35] and the input of Black Sea cold water inflow into the Aegean Sea [28]. The first temporal mode (PC1) reached its peak in the summer of 2003 (Figure 5b), which was the warmest summer over the study period, due to the strongest heat waves that hit the northern Europe. The highest negative peaks were observed in the winters of 1987 and 2019 in agreement with [46]. The lower peak values (i.e., cold periods) in the 1990s could be attributed to Mount Pinatubo's eruption in the Philippines, which caused a cold air-temperature anomaly throughout the Middle East during the winter of 1992 As mentioned by [47], The spatial pattern of the second EOF2 mode (Explains 16% of the total variance) demonstrated a dipole in the study area (Figure 5c), with positive anomalies in the Ionian Sea and negative ones in the Levantine and Aegean basins. The temporal coefficient of the second mode (PC2) alternates between positive and negative patterns over the study period 1982-2020.
Sci. Eng. 2021, 9, x FOR PEER REVIEW 8 of 21 Figure 4. Spatial trend map of SST anomaly over the EMED from AVHRR during 1982-2020, after the seasonal cycle was removed at each grid point.

Spatial Description and Trends of the MHWs Main Characteristics from 1982 to 2020
Figure 6a-h depicts the spatial distribution and trends per decade of the EMED MHWs main characteristics, frequency, duration, (icum), and (imax), which are overlaid with no significant values (p > 0.05) for the entire study period.
There was a clear decadal increase in MHW frequency in all EMED locations as demonstrated in Figure 6a,b. The frequency ranged from 1-2.5 counts with a maximum significant (p < 0.05) values > 2 events found in most of the Aegean and Levantine Seas, south of Crete and in front of the eastern coast of Libya. As shown in Figure 6b, the MHW frequency trend values were positively significant (p < 0.05) and ranged between 0 and 2.5 counts per decade. There were no significant (p > 0.05) changes in MHW frequency southeast of Crete, in front of Mersa Matruh city (see Figure 1a), and north of the Ionian Sea region (BiOS), which was previously described as a decadal reversal region in Section 3. The patterns of icum (trend) and mean duration (trend) are nearly identical, as shown in

Spatial Description and Trends of the MHWs Main Characteristics from 1982 to 2020
Figure 6a-h depicts the spatial distribution and trends per decade of the EMED MHWs main characteristics, frequency, duration, (i cum ), and (i max ), which are overlaid with no significant values (p > 0.05) for the entire study period.
There was a clear decadal increase in MHW frequency in all EMED locations as demonstrated in Figure 6a,b. The frequency ranged from 1-2.5 counts with a maximum significant (p < 0.05) values > 2 events found in most of the Aegean and Levantine Seas, south of Crete and in front of the eastern coast of Libya. As shown in Figure 6b, the MHW frequency trend values were positively significant (p < 0.05) and ranged between 0 and 2.5 counts per decade. There were no significant (p > 0.05) changes in MHW frequency southeast of Crete, in front of Mersa Matruh city (see Figure 1a), and north of the Ionian Sea region (BiOS), which was previously described as a decadal reversal region in Section 3. The patterns of i cum (trend) and mean duration (trend) are nearly identical, as shown in Figure 6c The authors found a strong relationship between the MHWs mean frequency and mean duration at each location over the study period (see, Figure 6a,c). The spatial pattern of mean MHW frequency and duration was highly negatively correlated (r = −0.76), which means the frequency was low where duration was long, and vice versa, example of similar analysis can be found in [3].  The authors found a strong relationship between the MHWs mean frequency and mean duration at each location over the study period (see, Figure 6a,c). The spatial pattern of mean MHW frequency and duration was highly negatively correlated (r = −0.76), which means the frequency was low where duration was long, and vice versa, example of similar analysis can be found in [3].

Description of the MHWs Mean Frequency Fields (1982-2020)
Figure 7a-d shows the spatial distribution of the EMED MHWs mean frequency events for the study periods. The P 1 mean MHW events are significantly < 2 counts over the whole EMED domain, Figure 7a. While during P 2 the MHW events have significantly increased up to 3 counts, Figure 7b. We found during P 3 that the Levantine basin has more significant (p < 0.05) MHW events (>3) than the Ionian and Aegean Seas, see Figure 7c. P 4 shows that most significant (p < 0.05) MHW events > 4 when compared to the other periods, Figure 7d.

Description of the MHWs Mean Frequency Fields (1982-2020)
Figure 7a-d shows the spatial distribution of the EMED MHWs mean frequency events for the study periods. The P1 mean MHW events are significantly < 2 counts over the whole EMED domain, Figure 7a. While during P2 the MHW events have significantly increased up to 3 counts, Figure 7b. We found during P3 that the Levantine basin has more significant (p < 0.05) MHW events (>3) than the Ionian and Aegean Seas, see Figure 7c. P4 shows that most significant (p < 0.05) MHW events > 4 when compared to the other periods, Figure 7d.

Description of the MHWs Mean Duration (Days) Fields (1982-2020)
Figure 8a-d describes the spatial distribution of the EMED MHWs mean Duration (days) for the different periods, (P1, P2, P3 and P4). The spatial field for MHW mean duration over P1 is mostly nonsignificant (p > 0.05) and less than 10 days for most of the domain. With an exception in the North Ionian area where the duration is significantly more than 18 days as shown in Figure 8a. For the most of the EMED domain, there was no significant difference in MHW mean duration (p > 0.05) during P2. The MHW mean duration for P2 ranges from 6.71 days to 20.53 days. We noticed during P2 a high significant (p < 0.05) MHW mean duration value > 20 days located off South of Crete, see Figure 8b. P3 and P4 show an increase of significant MHW mean duration values throughout the EMED domain. For P3, the MHW mean duration ranges from 9.07 days to 21.36 days, Figure 8c. The MHW mean duration pattern over the last decade shows greater significant values in the Levantine Sea. The least significant (p < 0.05) MHW mean duration is 10 days found in the South of the Aegean Sea, while it significantly exceeds 27 days off the Israeli coast as shown in Figure 8d. Finally, the findings revealed that the mean duration of MHWs has increased rapidly over the last two decades.

Description of the MHWs Mean Duration (Days) Fields (1982-2020)
Figure 8a-d describes the spatial distribution of the EMED MHWs mean Duration (days) for the different periods, (P 1 , P 2 , P 3 and P 4 ). The spatial field for MHW mean duration over P 1 is mostly nonsignificant (p > 0.05) and less than 10 days for most of the domain. With an exception in the North Ionian area where the duration is significantly more than 18 days as shown in Figure 8a. For the most of the EMED domain, there was no significant difference in MHW mean duration (p > 0.05) during P 2 . The MHW mean duration for P 2 ranges from 6.71 days to 20.53 days. We noticed during P 2 a high significant (p < 0.05) MHW mean duration value > 20 days located off South of Crete, see Figure 8b. P 3 and P 4 show an increase of significant MHW mean duration values throughout the EMED domain. For P 3 , the MHW mean duration ranges from 9.07 days to 21.36 days, Figure 8c. The MHW mean duration pattern over the last decade shows greater significant values in the Levantine Sea. The least significant (p < 0.05) MHW mean duration is 10 days found in the South of the Aegean Sea, while it significantly exceeds 27 days off the Israeli coast as shown in Figure 8d. Finally, the findings revealed that the mean duration of MHWs has increased rapidly over the last two decades.

Description of the MHWs (icum) Fields (1982-2020)
Figure 9a-d presents the spatial distribution of the EMED MHWs mean cumulative intensity (°C days) for the different periods, (P1, P2, P3 and P4). The MHW icum, defined as the integral of intensity over the duration of the event [5], represents the MHW severity as mentioned by [48]. The mean spatial pattern for P1 in Figure 9a shows non-significant values for most of the EMED basin. The highest P1 significant (p < 0.05) value is >30 (°C days) and detected to the south of Italy in the Ionian Sea. Figure 9b demonstrates that the significant (p < 0.05) MHW icum values have increased for P2. It ranges from 9.60 to 32.63 (°C days). The maximum MHW significant icum values for P3 is >30 (°C days) and are located to the South of Cyprus and in the middle of the Ionian Sea, Figure 9c. The map of P4 shows an increase of significant severity, which varies between 14.29 (°C days) and 39.51 (°C days) in the EMED domain. The highest significant MHW icum values > 35 (°C days) can be detected off the Israeli coast near the Shikmona anti-cyclonic gyre location as shown in Figure 9d. The Shikmona anti-cyclonic eddy was located the area between 34° E and 35° E longitudes and 33.5° N and 34.5° N latitudes [43,49,50].

Description of the MHWs (i cum ) Fields (1982-2020)
Figure 9a-d presents the spatial distribution of the EMED MHWs mean cumulative intensity ( • C days) for the different periods, (P 1 , P 2 , P 3 and P 4 ). The MHW i cum , defined as the integral of intensity over the duration of the event [5], represents the MHW severity as mentioned by [48]. The mean spatial pattern for P 1 in Figure 9a shows non-significant values for most of the EMED basin. The highest P 1 significant (p < 0.05) value is >30 ( • C days) and detected to the south of Italy in the Ionian Sea. Figure 9b demonstrates that the significant (p < 0.05) MHW i cum values have increased for P 2 . It ranges from 9.60 to 32.63 ( • C days). The maximum MHW significant i cum values for P 3 is >30 ( • C days) and are located to the South of Cyprus and in the middle of the Ionian Sea, Figure 9c. The map of P 4 shows an increase of significant severity, which varies between 14.29 ( • C days) and 39.51 ( • C days) in the EMED domain. The highest significant MHW i cum values > 35 ( • C days) can be detected off the Israeli coast near the Shikmona anti-cyclonic gyre location as shown in Figure 9d. The Shikmona anti-cyclonic eddy was located the area between 34 • E and 35 • E longitudes and 33.5 • N and 34.5 • N latitudes [43,49,50].

Description of the MHWs (i max ) Fields (1982-2020)
Figure 10a-d shows the spatial distribution of the EMED MHWs i max ( • C) for the study periods. During P 1 , the maximum significant (p < 0.05) mean MHW i max was found throughout most of the Ionian Sea with values > 2 • C, Figure 10a. The significant MHW i max values during P 2 vary between 1.33 • C at southeast Levantine basin and 2.45 • C at northeast Ionian basin near the Otranto strait, Figure 10b. The highest significant value for P 3 is 2.5 • C and is found north of the Aegean Sea, while the lowest significant value (1.31 • C) is found southeast of the Levantine Sea in front of the Egyptian coast. Among the entire study period, the MHW i max spatial pattern for P 4 has the highest significant (p < 0.05). The MHW i max values are significant during P 4 throughout the whole Levantine basin, Figure 10d. Figure 10d shows that the significant MHW i max values range from 1.42 • C off Mersa Matruh to more than 2.5 • C in the northeast Aegean basin.
(°C days). The maximum MHW significant icum values for P3 is >30 (°C days) and are located to the South of Cyprus and in the middle of the Ionian Sea, Figure 9c. The map of P4 shows an increase of significant severity, which varies between 14.29 (°C days) and 39.51 (°C days) in the EMED domain. The highest significant MHW icum values > 35 (°C days) can be detected off the Israeli coast near the Shikmona anti-cyclonic gyre location as shown in Figure 9d. The Shikmona anti-cyclonic eddy was located the area between 34° E and 35° E longitudes and 33.5° N and 34.5° N latitudes [43,49,50].

Description of the MHWs (imax) Fields (1982-2020)
Figure 10a-d shows the spatial distribution of the EMED MHWs imax (°C) for the study periods. During P1, the maximum significant (p < 0.05) mean MHW imax was found throughout most of the Ionian Sea with values > 2 °C, Figure 10a. The significant MHW imax values during P2 vary between 1.33 °C at southeast Levantine basin and 2.45 °C at northeast Ionian basin near the Otranto strait, Figure 10b. The highest significant value for P3 is 2.5 °C and is found north of the Aegean Sea, while the lowest significant value (1.31 °C) is found southeast of the Levantine Sea in front of the Egyptian coast. Among the entire study period, the MHW imax spatial pattern for P4 has the highest significant (p < 0.05). The MHW imax values are significant during P4 throughout the whole Levantine basin, Figure 10d. Figure 10d shows that the significant MHW imax values range from 1.42 °C off Mersa Matruh to more than 2.5 °C in the northeast Aegean basin.

Description of the MHWs (imax) Fields (1982-2020)
Figure 10a-d shows the spatial distribution of the EMED MHWs imax (°C) for the study periods. During P1, the maximum significant (p < 0.05) mean MHW imax was found throughout most of the Ionian Sea with values > 2 °C, Figure 10a. The significant MHW imax values during P2 vary between 1.33 °C at southeast Levantine basin and 2.45 °C at northeast Ionian basin near the Otranto strait, Figure 10b. The highest significant value for P3 is 2.5 °C and is found north of the Aegean Sea, while the lowest significant value (1.31 °C) is found southeast of the Levantine Sea in front of the Egyptian coast. Among the entire study period, the MHW imax spatial pattern for P4 has the highest significant (p < 0.05). The MHW imax values are significant during P4 throughout the whole Levantine basin, Figure 10d. Figure 10d shows that the significant MHW imax values range from 1.42 °C off Mersa Matruh to more than 2.5 °C in the northeast Aegean basin.  Table 2 shows the significant EMED extreme MHW events characteristics (mean duration, imax, imean, and icum) over the study period  with the corresponding MHWonset, MHWend, Longitude and Latitude. The maximum MHW imean was 3.92 °C in 2013, Figure 11a and lasted for 15 days (24 April-08 May) with an imax of 5.06 °C and icum  Table 2 shows the significant EMED extreme MHW events characteristics (mean duration, i max , i mean , and i cum ) over the study period  with the corresponding MHW onset , MHW end , Longitude and Latitude. The maximum MHW i mean was 3.92 • C in 2013, Figure 11a and lasted for 15 days (24 April-08 May) with an i max of 5.06 • C and i cum of 58.77 • C days. The most severe event (i cum of 319.65 • C days) was still ongoing at the time of the analysis with 126 days, i max of 3.41 • C and i mean of 2.54 • C, Figure 11b. The most intense event occurred in 2020, with an i max temperature of 6.35 • C and a duration of 7 days (this event will be discussed in detail in the following subsection). In general, the most extreme MHW significant events were mostly detected in 2020. Table 2. The EMED extreme significant (p < 0.05) MHW events characteristics (mean duration (Days), i max ( • C), i mean ( • C), and i cum ( • C days)) over the study period  with the corresponding MHW onset , MHW end , Longitude and Latitude in degrees where the highest events are highlighted in yellow color.  Figure 11b. The most intense event occurred in 2020, with an imax temperature of 6.35 °C and a duration of 7 days (this event will be discussed in detail in the following subsection). In general, the most extreme MHW significant events were mostly detected in 2020.   The EMED Most Intense [14][15][16][17][18][19][20] May 2020] MHW Event over the Study Period (1982Period ( -2020 In this subsection, the authors will describe in detail the most significant intense event of [14][15][16][17][18][19][20] May 2020. We found the most significant intense MHW event located at longitude~28.500 • E and latitude~34.250 • N south of Rhodes. Two more MHW significant events were discovered to be associated with the most intense main event, according to our analysis. The three significant events occur in the Levantine basin. It is critical in our research to assist in identifying and quantifying the main atmospheric causes of these extreme events. Table 3 demonstrates the EMED most intense [14][15][16][17][18][19][20] May 2020] MHW significant event characteristics, including duration, i max , i mean , i cum , wind stress magnitude, 2m air temperature and mean sea level pressure. The i max for the four events are 6.35, 5.66, 4.74, and 2.95 • C, respectively. Table 3. The EMED most intense (14-20 May 2020) MHW events characteristics (mean duration (Days), i max ( • C), i mean ( • C), i cum ( • C days)), wind stress magnitude (N/m 2 ), 2 m air temperature ( • C) and mean sea level pressure (mb) with the corresponding Longitude and Latitude in degrees.  Figure 12a shows the SST anomaly event map created by subtracting the historical climatological mean value at each grid from the mean temperature [14][15][16][17][18][19][20] May 2020 at the same location (grid). The description of temporal evolution of the SST most intense MHW event based on threshold 90th percentile and SST climatology are shown in Figure 12b. The spatial distribution patterns of the wind stress vectors combined with magnitude, 2 m mean air temperature and mean sea level pressure of ERA5 were chosen at the same time and location as the SST anomaly analysis, as shown in Figure 12c-e. We can see an interesting relationship between the SST anomaly, wind stress, and 2 m mean air temperature fields. We observed a reduction in the magnitude of wind stress fields (less than 0.001 N/m 2 ) at the same locations as the most intense MHW events (see Figure 12c and Table 3). The comparison of wind stress fields and SST anomaly fields revealed that MHW extreme events are related to wind stress reduction, which is consistent with the findings of [14][15][16]18]. Another intriguing parallel is found between the SST anomaly and the 2 m mean air temperature patterns. High air temperature values (>25 • C) were found at the same locations as the MHW extreme events, as shown in Figure 12d and Table 3. Furthermore, we observed a very persistent ridge of high atmospheric pressure >1014 mb over the MHW extreme events locations in the EMED as shown in Figure 12e and Table 3.

MHW
In conclusion, the atmospheric conditions, such as, mean sea level pressure, wind stress, and air temperature, may be linked to MHWs existence, Figure 12c-e.

Discussion
According to our SST climatology analysis (Figures 2 and 3), the maximum values over the four decades were detected during the last decade, which agrees with [19,35,36]. Table 1 have showed that the last decade was the warmest decade by about 0.94 • C during the entire 1982-2020 study period in agreement with [51,52]. During the last decade, the SST anomaly patterns revealed a warming of about 0.10 to 0.50 • C over the northern EMED basin, with the exception of some parts of the north Ionian Sea (Figure 3d). The authors have attributed this to the effect of BiOS which is responsible for decadal reversals of the Ionian basin-wide circulation as mentioned by [37][38][39][40][41]. The maximum SST anomaly for the last decade was 0.50 • C in the northwest Aegean Sea, while the minimum value was −0.50 • C along EMED's southwestern coast (Figure 3d).
The authors have discovered a high spatial variability of SST linear trend over the EMED, ranging from 0.1 to 0.5 • C per decade (Figure 4). The basin average warming rate was about 0.33 ± 0.04 • C per decade. Over the last twenty years, several SST studies for the Mediterranean and the EMED region were published, addressing trends and changes in ocean surface temperature [6,19,28,35,[52][53][54][55][56][57][58]. Ref. [53] have observed a trend of 0.3 • C per decade in the western basin and 0.5 • C per decade in the eastern basin over the period from 1985 to 2006. Similar results were recorded in [28], with an average mean warming rate of around 0.37 • C per decade for the whole Mediterranean basin during 1985-2008, and 0.26 • C per decade and 0.42 • C per decade for the western and eastern sub-basins, respectively. Over a 25-year (1993-2017), the warming rate in the Mediterranean was about 0.36 • C per decade, which was highly correlated with sea level trend, especially in the EMED [52]. According to [35], the annual warming trend for the deseasonalized Mediterranean global averaged SST was 0.36 • C per decade over the period 1982-2016. More recently, studies with longer time series, such as [19], estimated a 0.35 • C per decade warming for the global Mediterranean basin for the 38-year period 1982-2019. Figure 5 have revealed the highest negative peaks were observed in the winters of 1987 and 2019 in accordance with [46], The lower peak values (i.e., cold periods) in the 1990s could be attributed to the eruption of Mount Pinatubo in the Philippines, which caused a cold air-temperature anomaly throughout the Middle East during the winter of 1992, as stated by [47].
We have shown in Figure 6 the spatial distribution and trends per decade of the MHWs main characteristics (mean frequency, duration, (i cum ) and (i max )) overlaid with no significant values (p > 0.05) for the entire study period. According to Figure 6c-f we found the patterns of i cum (trend) and mean duration (trend) were nearly identical. The MHW mean duration has increased significantly (p < 0.05) with values > 17 days with a trend of > 10 days per decade in front of the Egyptian coast in the southeast of the Levantine Sea and off the Israeli coast (Figure 6d). No significant values (p > 0.05) were found in MHW duration at southeast of Crete, off Mersa Matruh city and the north Ionian Sea region. This could be attributed to the accuracy and precision of SST NOAA AVHRR data near coastlines, as described by [59,60]. According to their research, the accuracy and precision of Earth Observation (EO) SST data at the coastline are not well defined. According to Figure 7, the Levantine basin has more significant (p < 0.05) MHW events (>4) in the last decade than the Ionian and Aegean Seas.
In the last decade, the highest significant MHW (p < 0.05) duration > 27 days and i cum values > 35 ( • C days) were detected off the Israeli coast (Figures 8d and 9d), near the Shikmona anti-cyclonic gyre location [43,50]. The MHW i max spatial distribution pattern were significant (p < 0.05) throughout the entire Levantine basin over the last decade as demonstrated in Figure 10d. The significant MHW i max values have ranged from 1.42 • C off Mersa Matruh to more than 2.5 • C in the northeast Aegean basin, as presented in Figure 10d. An extreme MHW event was discovered in 2013, with the maximum MHW i mean temperature of 3.92 • C and lasted 15 days (24 April-08 May) with an i max of 5.06 • C and i cum of 58.77 • C days (Figure 11a and Table 2). While, the most severe event (i cum of 319.65 • C days) was still ongoing at the time of the analysis with 126 days, i max of 3.41 • C and i mean of 2.54 • C, Figure 11b. The most significant intense MHW event was located south of Rhodes at longitude~28.500 • E and latitude~34.250 • N and was associated with another two extreme events located in Levantine basin, that occurred simultaneously in [14][15][16][17][18][19][20] May 2020 with very high intensities (Figure 12a and Table 3). At the same locations as the most intense MHW events, we have observed a reduction in the magnitude of wind stress fields (less than 0.001 N/m 2 ) (see Figure 12c and Table 3). We have concluded that the MHW extreme events are related to the wind stress reduction in agreement with [14][15][16]18].

Conclusions
This study has quantified the SST interannual variability, trends, and spatial distribution of main MHW characteristics in the EMED over 39 years (1982-2020).
Some specific hotspot areas in the Mediterranean Sea as described by [61] such as Aegean, north Levantine and north Ionian have been recorded as the highest SST anomaly patterns during P 4 (2011-2020) which agrees with [18]. During last decade (P 4 (2011-2020)), the northern EMED basin's mean SST anomaly was nearly 1.00 • C higher than the southern one. The minimum SST anomaly value found −0.50 • C and was located on the southwestern coast of EMED in agreement with [18]. This could be due to upwelling associated with the Rhodes cyclonic gyre in the northern EMED basin, as mentioned in [42].
The SST EOF analysis revealed that the Ionian Sea had the highest interannual variability, while the Aegean had the lowest.
Over the entire study period, there was a clear decadal increase in MHW frequency in all EMED locations. The MHW frequency trend values were significantly (p < 0.05) positive. The correlation analysis between the mean MHWs spatial distribution of duration and frequency revealed a strong negative relationship, which was consistent with the findings of [3]. The authors found that the mean MHW frequency and duration increased by 40% and 15%, respectively We have characterized trends and variability of MHWs for the EMED from 1982 to 2020 based on satellite NOAA OI SST V2.1 data and a unified MHW framework. We have shown that between 1982 and 2020, on EMED average, MHW frequency trend increased by 1.2 events per decade and the average MHW i cum trend increased by 5.4 • C days.
The MHW main characteristics patterns over the last decade have shown greater significant values in most of the Levantine Sea. There were no significant (p > 0.05) found in the MHW main characteristics at southeast of Crete, off Mersa Matruh city and the north Ionian Sea region over the entire the study period.