The second study case (peaked on 18 October 2010) centered around the eastern side of the Gulf of Aden (12.875◦ N and 50.875◦ E), occupied 37% of the Red+ (36% in moderate; 1% in strong) area and lasted for 25 days (4 days as a strong MHW and 16 days as a moderate MHW), with a mean intensity of 2.02 ◦C (equivalent to accumulative intensity of 50.5 ◦C-day), as seen in Figure 12. There is a significant correlation between the marine heat wave intensity during this study case and different components of the surface radiation at the 95% significance level, most markedly with SLHF (R = 0.65, *n* = 25), as seen in Figure 13. Similarly, Floss shows a significant correlation with the MHW intensity (R = 0.64, *n* = 25), as seen in Figure 13e.

During this second study case, the chlor\_a concentration is significantly lower than its climatological values over 82% of the time (Figure 14a), partly showing the effect of MHW on the chlor\_a concentration and hence the Red Sea ecology. Moreover, there is a significant correlation between the MHW intensity and chlor\_a concentration during this study case (R = −0.40, *n* = 25) at the 95% significance level, as seen in Figure 14b.

**Figure 10.** Representative of the MHW event (study case 1) showing the sea surface temperature (**a**), intensity [◦C] of marine heat wave event (**b**) and category of that event (**c**) at the peak of the event and the time series during the event's year (**d**), spatially averaged over an area around the maximum intensity grid, which is shown in bubbles.

**Figure 11.** *Cont*.

**Figure 11.** Scatter plots of marine heat wave intensity during the study case 1 (10 July 2001 to 22 August 2001) and surface latent heat flux (**a**), surface sensible heat flux (**b**), surface net thermal radiation (**c**), surface net solar radiation (**d**), and total heat loss to the atmosphere (**e**).

**Figure 12.** Representative of the MHW event (study case 2) showing the sea surface temperature (**a**), intensity [◦C] of marine heat wave event (**b**) and category of that event (**c**) at the peak of the event and the time series during the event's year (**d**), spatially averaged over an area around the maximum intensity grid, which is shown in bubbles.

**Figure 13.** Scatter plots of marine heat wave intensity during the study case 2 (29 September 2010 to 23 October 2010) and surface latent heat flux (**a**), surface sensible heat flux (**b**), surface net thermal radiation (**c**), surface net solar radiation (**d**), and total heat loss to the atmosphere (**e**).

**Figure 14.** Representative of the chlorophyll\_concentration\_in\_sea\_water (chlor\_a) values and its climatological values during the study case 2, which extends from 29 September 2010 to 23 October 2010 (**a**); scatter plots of marine heat wave intensity and chlor\_a concentration during the study case 2 (**b**).

The third study case, which peaked on 18th of July 2019, extended for 24.1% (23% in moderate; 0.5% in strong; 0.4% in severe and 0.2% in extreme) of the study area, most markedly over the strait of Bab al Mandab (Figure 15). Over this strait, this MHW had a mean intensity of 2.99 ◦C over a duration of 26 days (equivalent to accumulative intensity of 77.84 ◦C-day): 3 days as an extreme MHW, 5 days as a severe MHW, 12 days as a strong MHW and 6 days as a moderate MHW. At the 95% level of significance, the marine heat wave intensity during this study case showed a significant correlation only with SNTR (R = −0.45, *n* = 26), as seen in Figure 16c. In the same context, the other studied surface

radiation components (SLHF, SSHF, SNSR and Floss) had an insignificant correlation with the marine heat wave intensity.

**Figure 15.** Representative of the MHW event (study case 3) showing the sea surface temperature (**a**), intensity [◦C] of marine heat wave event (**b**) and category of that event (**c**) at the peak of the event and the time series during the event's year (**d**), spatially averaged over an area around the maximum intensity grid, which is shown in bubbles.

**Figure 16.** *Cont*.

**Figure 16.** Scatter plots of marine heat wave intensity during the study case 3 (27 June 2019 to 22 July 2019) and surface latent heat flux (**a**), surface sensible heat flux (**b**), surface net thermal radiation (**c**), surface net solar radiation (**d**), and total heat loss to the atmosphere (**e**).

During this study case, the link between this MHW event and the chlor\_a concentration is not discussed due to the missing chlor\_a concentration data.

The fourth study case, which peaked on 13 September 2020, extended over 85% (62% in moderate; 23% in strong) of the study area, most markedly around the northern side of the Red Sea (27.625◦ N and 34.875◦ E), as seen in Figure 17. Over the area centered around (27.625◦ N and 34.875◦ E), the MHW duration extended for a 98-day period from 30 August 2020 to 5 December 2020. This MHW had a mean intensity of 2.56 ◦C over a 98-day duration (equivalent to an accumulative intensity of 251.3 ◦C- day): 42 days as a strong MHW and 54 days as a moderate MHW. At the 95% significance level, the correlation between the marine heat wave intensity during this study case and different components of surface radiation, together with the Floss, are significant, most markedly with SNTR (R = −0.65, *n* = 98), as seen in Figure 18.

**Figure 17.** *Cont*.

**Figure 17.** Representative of the MHW event (study case 4) showing the sea surface temperature (**a**), intensity [◦C] of marine heat wave event (**b**) and category of that event (**c**) at the peak of the event and the time series during the event's year (**d**), spatially averaged over an area around the maximum intensity grid, which is shown in bubbles.

**Figure 18.** Scatter plots of marine heat wave intensity during the study case 4 (30 August 2020 to 5 December 2020) and surface latent heat flux (**a**), surface sensible heat flux (**b**), surface net thermal radiation (**c**), surface net solar radiation (**d**), and total heat loss to the atmosphere (**e**).

During this study case, the chlor\_a concentration is significantly lower than its climatological values over 100% of the time (Figure 19a), confirming the effect of MHWs on the Red Sea ecology similarly to study case 2. In the same context and at the 95% significance

level, the MHW intensity and chlor\_a concentration had a significant correlation (R = −0.43, *n* = 98) during this study case, as seen in Figure 19b.

**Figure 19.** Representative of the chlorophyll\_concentration\_in\_sea\_water (chlor\_a) values and its climatological values during the study case 4, which extends from 30 August 2020 to 5 December 2020 (**a**); scatter plots of marine heat wave intensity and chlor\_a concentration during the study case 2 (**b**).

#### *3.4. Main Spatiotemporal Characteristics of Marine MHWs over the Red+*

The annual average of the MHWs intensity over Red+ has significant spatial variability, where the lowest values (<1.2 ◦C) were found along the Saudi Arabia coast at around 18.875◦ N latitude, and the highest values (>2 ◦C) were found to the south of the Gulf of Aqaba at around 29.375◦ N latitude, as described by Figure 20a. On the other hand, the MHWs duration exhibits a markedly spatial pattern and reached its maximum values along the Saudi Arabia coast at around 28.375◦ N latitude and its minimum values in the southern part of the Red Sea at around 15.625◦ N, as seen in Figure 20b. Moreover, the MHW frequency shows a different spatial pattern, where its maximum values were found along the Saudi Arabia coast near 19.375◦ N and its minimum values at the north part of the Gulf of Suez near 29.625◦ N, as seen in Figure 20c. In general, the Red Sea displays a meridional gradient of an increasing annual average of the MHW intensity and MHW duration from north to south. Conversely, the Red Sea displays a meridional gradient of a decreasing annual average of the MHW frequency from north to south.

**Figure 20.** *Cont*.

**Figure 20.** Marine heat wave (MHW) spatial/temporal characteristics over Red+. The latitudinal changes in annual average MHW intensity (**a**), MHW duration (**b**), and MHW frequency (**c**) are also shown.

#### **4. Summary Discussion**

Extreme Red+ sea surface temperature is calculated based on analyses of seasonal mean time series from 1982 to 2020, confirming that the thresholds of warm (cold) ETEs is 30.03 (25.13) ◦C using a 90th (10th) percentile definition. Warm ETEs will be used as an additional thermal restriction to MHWs' threshold values.

The changeable (based on each day of the year) climatological SST mean and 90th percentile threshold was used to identify the MHW events. To avoid the cold season period, warm ETEs were used as an additional thermal restriction to the MHWs' threshold values.

From 1982 to 2020, there are 28 different MHW events that extend for 378 days; 360 of those days are in the moderate category, whereas the other 18 days are in a strong category. On the other hand, the MHW average intensity, annual duration and frequency described a significant increasing trend over Red+ from 1982 to 2020. These results follow the global trend pattern (Oliver et al. [3]). The longest MHWs over Red+ extended for 57 days from 31 August 2020 to 26 October 2020, with a maximum (mean) intensity of 1.76 (1.11) ◦C and a cumulative intensity of 63.78 ◦C-day. For the longest MHW event, the currently identified MHW duration is approximately 3.2 times more than previously detected by Shaltout [6]. This difference in determining the longest duration of MHWs is due to the fact that the current study extends into 2020 and also because different methods were used. In general, the MHW over Red+ is in a moderate category, with a frequency of 1.38 events annually. The average value of an event intensity and duration is 0.8 ◦C and 22.5 days, respectively. Furthermore, the annual average days of MHW over Red+ is 9 days, which agrees with the previous finding of Shaltout [6].

During El Niño periods, the MHW intensity tends to reach its maximum value. In the same context, the ISM rarely explains the MHW intensity (the higher values of the MHWs' intensity may occur simultaneously with lower ISMI values and, at another time, may occur with higher ISMI values). Furthermore, the sea surface radiation budget of different components may describe critical conditions of the occurrence of MHW events (more than 7Wm−2 of Floss going to the ocean). As ENSO and ISMI indexes, together with different components of the sea surface radiation budget, cannot completely describe the controls of MHWs, there is another climatic mechanism that, if coupled with them, will describe the MHWs' characteristics. This mechanism merits our consideration and will be discussed in our future work by describing a new climatic index for the study area.

There is a markedly interannual variability for the studied MHW property over Red+, which is followed well by the global properties (Oliver et al. [35]). The pattern of the mean MHW frequency, duration, intensity and SST were positively correlated over Red+, indicating that global warming is the main reason for the positive trends of the MHW frequency, duration and intensity. On a global scale, the average duration and frequency of MHWs were negatively correlated (Oliver et al. [35]), indicating that their relations over Red+ (depending on the current result) do not follow the global patterns. Moreover, higher frequency and intensity values on a global scale appear to relate to El Niño period events (Oliver et al., 2018); however, over Red+, only the intensity is clearly related to El Niño periods. On the other side, the MHWs' intensity has a clearer relation with ISMI than the frequency and duration have with ISMI. This encourages us to find and describe a new valid climatic indicator for the study area in our future study.

For further details about the MHW characteristics, the longest four MHW events were analyzed in depth. These further analyses prove that the chlor\_a concentration has lower values than its climatic values during the MHW events, providing early awareness about the impact of MHWs as a heat stress on various marine sectors while shedding light on decision makers in finding a suitable regional climate policy to cope with global warming issues.

Finally, the spatiotemporal analysis of MHWs confirmed that the spatial distribution of the MHW annual average intensity, duration and frequency had a distinct spatial pattern, which agrees with the previous finding of Genevier et al. [20]. The northern region of Red+ (south of the Gulf of Suez and Aqaba) witnessed the most intense MHW events, which disagrees with the previous finding of Genevier et al. [20]. On the other hand, the eastern coast of the northern region had the greatest number of MHW days, which agrees with the previous finding of Genevier et al. [20]. The northeastern coast of the southern region had the most frequent events, which agrees with Genevier et al. [20]. Generally, the current results agree with the previous finding of Genevier et al. [20] concerning the MHW intensity and duration spatial distribution. The only exception in the spatial pattern of the MHW duration is possibly due to the use of warm ETEs as an additional thermal restriction to MHWs' threshold values.

**Author Contributions:** Conceptualization, M.S. and K.T.; methodology, M.S.; software, M.S. and A.B.; validation, M.S., N.G., N.B. and K.T.; formal analysis, N.G.; investigation, N.B.; resources, A.B.; data curation, M.S. and K.T.; writing original draft preparation, M.S.; writing review and editing, M.S. and K.T.; visualization, M.S. and N.G.; supervision, M.S.; project administration, A.B.; funding acquisition, A.B. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia, grant number 1213 and "The APC was funded by the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia".

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors extend their appreciation to the Deputyship for Research & Innovation, Ministry of Education in Saudi Arabia for funding this research work through the project number 1213.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

## *Article* **Spatiotemporal Variability and Trends of Marine Heat Waves in the Red Sea over 38 Years**

**Bayoumy Mohamed 1,2, Hazem Nagy 1,3,\* and Omneya Ibrahim 1**

**Abstract:** Marine heat waves (MHWs) can have catastrophic consequences for the socio-environmental system. Especially in the Red Sea, which has the world's second longest coral reef system. Here, we investigate the sea surface temperature (SST) variability and trends, as well as the spatiotemporal characteristics of marine heat waves (MHWs) in the Red Sea, using high resolution daily gridded (1/20◦) SST data obtained from the Copernicus Marine Environment Monitoring Service (CMEMS) for the period 1982–2019. Results show that the average warming rate was about 0.342 ± 0.047 ◦C/decade over the entire Red Sea over the whole study period. The Empirical Orthogonal Function (EOF) analysis reveals that the maximum variability is over the central part of the Red Sea, while the minimum variability is in the southernmost part of the Red Sea. Over the last two decades (2000–2019), we have discovered that the average MHW frequency and duration increased by 35% and 67%, respectively. The results illustrate that the MHW frequency and duration trends have increased by 1.17 counts/decade and 1.79 days/decade, respectively, over the study period. The highest annual MHW frequencies were detected in the years 2018, 2019, 2010, and 2017. A strong correlation (R = 0.89) was found between the annual MHW frequency and the annual mean SST.

**Keywords:** Red Sea; marine heat wave; sea surface temperature; duration; frequency

#### **1. Introduction**

Heat waves are characterized as extended periods of unusually hot weather that have become more frequent and longer in recent decades, posing a threat to human health and ecosystems. A similar phenomenon known as a 'Marine Heat Wave' (MHW) has been observed in the oceans, endangering marine ecosystems and productivity [1–4]. The most recent definition of MHWs is discrete periods of abnormally warm sea surface temperatures that can last from five days to months above the 90th percentile threshold of SST climatology [5–7]. Despite extensive knowledge of global SST changes, the study and detection of MHWs and associated climate processes remains limited [2].

The Red Sea has a rich and diverse ecosystem [8]. It has the world's second longest coral reef system [9]. The importance of the Red Sea is highlighted by its high tourism volume and aquaculture, and also its relation to the gas and oil and fishing industries. Furthermore, the Red Sea is a hub for shipping activity, connecting European ports with China and Eastern Asia through Suez Canal.

The Red Sea has a negative water balance, which means that evaporation exceeds precipitation and combined river runoff [10,11]. It is considered to be an arid climate region with little precipitation and runoff, as there is a lack of any permanent rivers that flow into it [12–14].

The impact and consequences of MHW on Red Sea marine life are now critical for decision makers and stakeholders. There are no published MHW studies of the Red Sea. However, some studies in the Rea Sea have discovered a link between extreme heat events

**Citation:** Mohamed, B.; Nagy, H.; Ibrahim, O. Spatiotemporal Variability and Trends of Marine Heat Waves in the Red Sea over 38 Years. *J. Mar. Sci. Eng.* **2021**, *9*, 842. https:// doi.org/10.3390/jmse9080842

Academic Editor: Francisco Pastor Guzman

Received: 15 July 2021 Accepted: 1 August 2021 Published: 4 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

such as ocean heat waves and the spread of the El Nino-Southern Oscillation (ENSO) [15]. In addition, the authors have found several SST studies for the Red Sea region, addressing trends and changes in ocean surface temperature [16–27].

According to [21], the annual mean SST has changed from 27.4 ◦C from 1985 to 1993 to 28.1 ◦C from 1994 to 2007, which was accompanied with significant warming in the mid-1990s. The authors of [15] discovered that the Red Sea maximum SST is warming at a rate of 0.17 ◦C/decade, with the northern Red Sea alone warming from 0.4 to 0.45 ◦C/decade from 1982 to 2015, which exceeds the global trend of 0.11 ◦C/decade [15]. The amplification of recent summer warming and decadal variability of the Red Sea SST was studied in [25]. They found that the significant warming which began in the mid-1990s precisely corresponded with the onset of a positive phase of the Atlantic Multi-Decadal Oscillation (AMO) and a negative phase of the Silk Road Pattern (SRP) in 1997.

In this study, we have investigated the spatiotemporal characteristics and occurrence of MHWs in the Red Sea from 1982 to 2019 using high resolution SST OSTIA (Operational SST and Ice Analysis) data (0.05◦ × 0.05◦). In addition, the SST spatial trend and variability were examined.

Section 2 describes the data and methodology of SST analysis. Section 3 reports the SST climatology, trend, and interannual variability. Section 4 shows the MHW results over the study period. Section 5 presents the discussion and Section 6 provides the conclusions.

#### **2. Data and Methods**

The Red Sea is an almost enclosed marginal tropical sea characterized by an extended shape that separates the African and Asian continents and is connected to the Indian Ocean through the Bab El-Mandab Strait, as shown in Figure 1. It has a moderate size of 4.5 × 106 km3, which is roughly one fifth the size of the Mediterranean Sea. The Bab El-Mandeb strait, which connects the basin to the open ocean, is a small and shallow channel (sill depth of 160 m and the narrowest point is 25 km wide). The study area extending from 12.5◦ N to 30◦ N and from 32◦ E to 44◦ E has an average width of ~220 km and a length of ~2000 km, with a maximum centered depth of 2200 m (Figure 1).

126

#### *2.1. Sea Surface Temperature Climatological Data*

The daily Sea Surface Temperature (SST) data used were freely obtained from the Copernicus Marine Environment Monitoring Service website (CMEMS) (https://marine. copernicus.eu (accessed on 20 April 2021)) from 1 January 1982 to 31 December 2019. The data website is https://resources.marine.copernicus.eu/?option=com\_csw&view= details&product\_id=SST\_GLO\_SST\_L4\_REP\_OBSERVATIONS\_010\_011 (accessed on 23 March 2021). The product name is 'SST\_GLO\_SST\_L4\_REP\_OBSERVATIONS\_010\_011'. The CMEMS OSTIA (Operational SST and Ice Analysis) reprocessed analysis product is a satellite and in situ foundation SST analysis created by the OSTIA [28]. The SST dataset includes daily gap-free maps of sea surface temperature and ice concentration (known as the L4 product) with a horizontal grid resolution of 0.05◦ × 0.05◦ [29]. CMEMS OSTIA data for the Red Sea were extracted from the global data, providing a 16534-point regularly gridded dataset spanning 13,879 days from 1 January 1982 to 31 December 2019.

#### *2.2. Trend and Interannual Variability of SST*

The Empirical Orthogonal Functions (EOFs) analysis [30] was used to evaluate the dominant spatiotemporal characteristics of SST variability from 1982 to 2019, utilizing the monthly de-seasonalized and detrended SST dataset to focus on non-seasonal variability. Thus, before decomposing the EOF, the seasonal cycle and linear trends were removed from the monthly means of SST anomalies at each grid cell, and each point time series was divided by its standard deviation to normalize the SST [31,32]. The mean seasonal cycle was estimated by calculating the mean monthly values for each calendar month at each grid point. SST anomalies were computed by removing the historical climatological mean value at each grid point from the SST values at the same location (grid). The strong seasonal signal was removed from the datasets at each grid cell to produce de-seasonalized maps and time series [33]. This was accomplished by subtracting the climatological mean of each calendar month from the corresponding months in all years. The least squares method [34] was used to estimate linear trends in de-seasonalized monthly sea surface temperature anomalies (SSTA). The Modified Mann–Kendall test [35,36] was used to test the statistical significance of the trends. In this study, the significance level was set to *p* ≤ 0.05, ensuring that all given trends are statistically significant at least at the 95% confidence level. The error of linear trend at the 95% confidence level was calculated using [34] the formula which was mentioned in [37].

#### *2.3. Marine Heat Waves (MHWs) Calculations*

The authors of [5] used the "discrete prolonged anomalously warm-water events in a location" definition of MHWs, using the same daily SST dataset for the 1982–2019 period. The MATLAB toolbox [38] was used to detect periods of time where temperatures were above the 90th percentile threshold compared to the climatology (1982–2019), and with a duration of at least 5 days. The climatological mean has been calculated from the daily Sea Surface Temperature (SST) data from 1 January 1982 to 31 December 2019; the equations of MHW main characteristics are according to [5] and described in [7]. Daily SSTs were spatially averaged over the Red Sea (12.5◦–30◦ N, 32◦–44◦ E) in our study to provide a regional daily SST time series spanning 1982–2019. In order to identify the MHWs in the Red Sea, the MHW definition was then applied to the SST time series.

Every MHW event has a start and end date defined by a mean (imean) and a maximum intensity (imax) in ◦C (mean and spatiotemporal maximum temperature compared to the threshold over the MHW event duration (in days)), as described in [39]. The MHW cumulative intensity (icum) in ◦C days is defined as the spatiotemporal summation of daily temperature compared to the threshold over the event duration as mentioned by [39]. Following the methodology of [4,40], the annual statistics were computed, including the frequency of events (i.e., the number of discrete events occurring each year). The MHW total number of days in each year can be used to examine how the frequency of occurrences in counts has changed over time, using the metrics for each location [5]. According to [5,7],

the main characteristics of MHWs are maximum intensity (imax) in ◦C, mean intensity (imean) in ◦C, and cumulative intensity (icum) in (◦C days):

$$\mathbf{i}\_{\text{max}} = \text{max}\left(\mathbf{T}(\mathbf{t}) - \mathbf{T}\_{\text{m}}(\mathbf{j})\right) \tag{1}$$

$$\mathbf{i}\_{\text{mean}} = \overline{\mathbf{T}(\mathbf{t}) - \mathbf{T}\_{\text{m}}(\mathbf{j})} \tag{2}$$

$$\mathbf{i}\_{\text{cum}} = \int\_{\mathbf{t}\_{\text{s}}}^{\mathbf{t}\_{\text{s}}} \left( \mathbf{T}(\mathbf{t}) - \mathbf{T}\_{\text{m}} \left( \mathbf{j} \right) \right) dt \tag{3}$$

where Tm is the climatological OSTIA SST mean calculated over the period from 1 January 1982 to 31 December 2019, to which all values are relative; j is the day of year; and ts and te are the start and end dates of the MHW event, respectively.

#### **3. SST Climatology, Trend, and Interannual Variability**

Figure 2a depicts the temporal evolution of the daily mean climatology of Red Sea SST (◦C) from 1982 to 2019. The highest temperatures (>29 ◦C) were in summer and early autumn (i.e., from July to October), while the lowest temperatures (<25 ◦C) were in winter and early spring (i.e., from January to March). The seasonal temperature difference between summer and winter is about 6 ◦C, which is consistent with [22]. Figure 2b shows the monthly de-seasonalized SST temporal trend over the study period. The higher deseasonalized SST values were noticed in the last 4 years of the study period (i.e., 2016–2019). The Red Sea temporal trend was about 0.32 ± 0.03 ◦C/decade.

**Figure 2.** *Cont*.

**Figure 2.** (**a**) The temporal evolution of the daily basin averaged SST (**b**) the de-seasonalized monthly SST time series (denoted by the continuous red line) and the SST trend (denoted by the continuous black line) from 1982 to 2019.

The spatial trend map of the de-seasonalized SSTA from 1982 to 2019 is shown in Figure 3. At the 95% confidence interval, a statistically significant trend was observed over the entire region. The linear trend over the Red Sea is not uniform, varying from 0.1 to 0.47 ◦C/decade. The basin averaged warming trend rate was about 0.342 ± 0.047 ◦C/decade. The highest trend was found in the deep area of the Red Sea (i.e., depth >~1200 m), while lower values were found in the Gulfs of Suez and Aqaba and the southern part of the Red Sea.

**Figure 3.** Spatial trend map of the Red Sea SST anomaly from 1982 to 2019, after the seasonal cycle was removed at each grid point.

The first two EOF modes of SST account for 73% of the non-seasonal variance. The spatial pattern of the first EOF (EOF1) (accounting for approximately 52% of the nonseasonal variance) shows a positive anomaly over the entire Red Sea (Figure 4a), showing an in-phase oscillation. The highest variability was found over central and north part of Red Sea, while the lowest variability was detected over the Gulfs of Suez and the Aqaba and the southernmost region of the Red Sea. Strong and continuous northeasterly monsoon winds from the Arabian Sea redirected by mountains into southeasterly winds along the axis of the Red Sea [41] may have a thermal effect that is responsible for the lowest variability over the southern Red Sea. The lower variability in the Gulfs of Suez and Aqaba is due to strong northwesterly winds that penetrate the western edge of the Sinai Peninsula, blow southwards, and sometimes reach the southern end of the Red Sea. Its temporal coefficient PC1 shows a strong interannual variability, reaching its peak in the summers of 1995, 1997, and 2019. The highest negative peaks (i.e., cold periods) were observed in the winters of 1984, 1992, and 2015. The second mode (accounting for approximately 21% of the non-seasonal variance) is out-of-phase variability of the SST in the Red Sea, with positive anomaly values in the south and negative ones in the north. Over the study period of 1982–2019, the temporal coefficient of the second mode (PC2) fluctuates between positive and negative patterns.

**Figure 4.** *Cont*.

**Figure 4.** (**a**) Spatial distribution for the first EOF (EOF1); (**b**) Spatial distribution for the second EOF (EOF2); (**c**) Temporal evolution for EOF (PC1) in normalized variance units; (**d**) Temporal evolution for EOF (PC2) in normalized variance units; for the Red Sea SSTA over the period from 1982 to 2019.

#### **4. Main Characteristics and Trends of Marine Heat Waves (MHWs)**

The mean and spatial trend of the main characteristics of MHWs (frequency, days, duration, and intensities) over the entire study period are shown in Figure 5. The trend maps are superimposed with non-significant values (*p* > 0.05). Figure 5a demonstrates that the MHW frequency varied from 1 to 2.5 counts, with maximum values of >2 counts found in the coastal areas of Sudan, Eritrea, Yemen and Saudi Arabia. The MHW frequency trend values were statistically significant (*p* < 0.05) over the whole domain. The MHW frequency trend ranged between 0 and 2 counts/decade. The highest MHW frequency trend (>1.8 counts/decade) was detected in the deep area (i.e., depth > 1000 m) between Eritrea and south coast of Saudi Arabia, as shown in Figures 1 and 5b. The MHW total days spatial patterns in the study period (1982–2019) varied between 15 and 30 days, with maximum values of >25 days found inside the Gulfs of Suez and Aqaba (Figure 5c). This could be due to the Gulfs of Suez and Aqaba being well resolved with the high-resolution OSTIA data. In general, the northern Red Sea has higher values of MHW total days than the southern region. The MHW total days trend values were statistically significant (*p* < 0.05) across the whole domain (Figure 5d). We discovered a similarity among the *p*-values between the MHW frequency and total days trend fields (Figure 5b,d). The MHW mean duration ranged from 5 to 15 days, with maximum values of >12 days north of the Red Sea for the Gulfs of Suez and Aqaba (Figure 5e). The maximum MHW duration spatial trend significant values (*p* < 0.05; >8 days/decade) were found in the middle of the Red Sea (i.e., between ~lat 18◦ E and 21◦ E) (Figure 5f). Both the MHW duration and the icum spatial patterns revealed similarities between the two fields (Figure 5e,g). The MHW icum values ranged from 10 to 25 (◦C days), with maximum values of >22 (◦C days) north of the Red Sea, excluding the Gulfs of Suez and Aqaba (i.e., between ~lat 22◦ E and 28◦ E) (Figure 5g). The maximum MHW icum significant trend values (*p* < 0.05; >10 ◦C days/decade) were discovered in the same areas as the MHW duration trend (Figure 5f,h). The MHW imax showed values of >1.8 ◦C throughout most of the Red Sea except the southernmost part beside Bab El-Mandab (Figure 5i). The MHW maximum imax significant trend values (*p* < 0.05; >2.5 ◦C) were found in the middle region of the Red Sea (i.e., between ~lat 19◦ E

and 21◦ E) and in the north of the Red Sea close to the Gulf of Suez (i.e., between ~lat 26◦E and 28◦ E) off the Egyptian coast (Figure 5j).

**Figure 5.** *Cont*.

**Figure 5.** *Cont*.

**Figure 5.** The main characteristics and trends of marine heat waves in the Red Sea over the study period (1982–2019) (**a**) mean MHW frequency (counts); (**b**) MHW frequency trend (counts/decade); (**c**) mean MHW total days (days); (**d**) MHW days trend (days/decade); (**e**) MHW duration (days); (**f**) MHW duration trend (days/decade); (**g**) MHW (icum) (◦C days); (**h**) MHW (icum) trend (◦C days/decade); (**i**) mean MHW (imax) (◦C); (**j**) MHW (imax) trend (◦C/decade); the black dots indicate that the change is not significantly different (*p* > 0.05).

#### *4.1. Marine Heat Waves (MHWs) Temporal Variation*

In this section, the authors present the temporal basin averaged annual MHW main characteristics (frequency, duration, and icum) from 1982 to 2019. The goal was to investigate temporal changes in Red Sea MHW characteristics such as frequency, duration, and icum, as well as their relationship to SST (Figure 6). The MHW matrices for the Red Sea region show a clear ascending trend over the last decade, as shown in Figure 6.

From 1982 to 1994, the annual mean MHW frequency was <1, while from 1994 to 2019, the frequency of MHWs increased significantly, reaching a peak in 2018 (7 counts) (Figure 6a). These findings are consistent with the monthly SST time series shown in Figure 2b, from 1982 to 2019. We noticed that no MHWs were detected in 1984. We discovered a strong correlation (R = 0.89) between the annual mean SST and MHW frequency, (Figure 6b).

The results show a general rise in annual MHW duration from approximately 6 days to 21 days in 2018 (Figure 6c), while icum ranged from 8 ◦C days to 30 ◦C days (Figure 6e). The highest annual MHW duration and icum were found during 2018 (~21 days and 30 ◦C days). We discovered that the average frequency and duration of MHWs increased by 35% and 67%, respectively, over the last two decades (2000–2019).

The correlation coefficient between MHW duration and annual mean SST was 0.64, while the correlation coefficient between MHW icum and SST was 0.72.

**Figure 6.** Bars of regional averaged annual means of the MHW (**a**) frequency, (**c**) duration, and (**e**) cumulative intensity (icum), and scatter plots of annual mean SST and MHW (**b**) frequency, (**d**) duration, and (**f**) cumulative intensity (icum), over the Red Sea over the period from 1982 to 2019, with the black lines representing the best-fit linear curve.

#### *4.2. Examples of MHW Events*

In this section, we show two example MHW events discovered in our MHW analysis between 1982 and 2019. In terms of duration and imax, we chose two MHW events in 2019. Table 1 shows examples of the Red Sea MHW event characteristics (mean duration, imax, imean, and icum) over the study period (1982–2019) with the corresponding MHWonset (start date), MHWend (end date), and longitude and latitude.

**Table 1.** Examples of the Red Sea MHW event characteristics (mean duration (days), imax ( ◦C), imean ( ◦C), and icum ( ◦C days)) over the study period (1982–2019) with the corresponding MHWonset, MHWend, and longitude and latitude (in degrees).

The long duration event (236 days) with icum 333.94 ◦C days was detected in the Gulf of Suez during the period from 10 May 2019 to 31 December 2019, with an imax of 2.94 ◦C and imean of 1.42 ◦C (Figure 7a,b). A high MHW imax event of 4.21 ◦C was noticed in 2019 off the Sudanese coast (i.e., between ~lat 18.5◦ E and 20.2◦ E) (Figure 7c,d) and lasted for 66 days (17 May–21 July) with an imean of 2.15 ◦C and icum of 142.06 ◦C days.

**Figure 7.** *Cont*.

**Figure 7.** Examples of the Red Sea MHW events over the study period (1982–2019) chosen in the year 2019, The SST climatology from OSTIA SST for the detection of MHWs (blue color), 90th percentile MHW threshold (green color), and SST time series (black color) for each MHW. The area filled with pink color shows the period accompanied with the identified MHW. (**a**) MHW duration event (236 days) (10 May 2019 to 31 December 2019); (**b**) Spatial mean SST anomaly distribution ◦C (10 May 2019 to 31 December 2019); (**c**) MHW maximum intensity event (imax 4.21 ◦C) (17 May 2019–21 July 2019); (**d**) Spatial mean SST anomaly distribution ◦C (17 May 2019–21 July 2019).

#### **5. Discussion**

Regarding to our SST climatology and trends analysis (Figures 2 and 3), the highest temperatures (>29 ◦C, Figure 2a) were found in summer and early autumn. The lowest temperatures (<25 ◦C, Figure 2a) were recorded in winter and early spring. We found the seasonal temperature difference between summer and winter to be approximately 6 ◦C (Figure 2a) in agreement with [22]. According to Figure 2b, the Red Sea temporal SST trend was approximately 0.32 ± 0.03 ◦C/decade from 1982 to 2019, which was higher than [26], possibly due to the use of longer time series and high-resolution OSTIA data.

Many efforts have been made to investigate the Red Sea SST trends using NOAA OI SST coarse resolution data. The authors of [24] showed that August and February are the warmest and coolest months in the Red Sea, with a significant trend of 0.5 ◦C and 0.3 ◦C/decade, respectively. According to [15], the overall warming trend of the Red Sea maximum SST was approximately 0.17 ◦C/decade, while the northern Red Sea alone warmed by between 0.4 and 0.45 ◦C/decade from 1982 to 2015. The authors of [42] examined the warming of the Red Sea SST from 1982 to 2016 and found an intense warming trend of 0.29 ◦C/decade. Recently, [26] demonstrated that the SST warming trend ranged from 0.1 ◦C/decade in the south to 0.4 ◦C/decade in the north Red Sea between 1982 and 2017. Our analysis revealed a high spatial variability of SST trends over the Red Sea (Figure 3) in the period from 1982 to 2019, varying from 0.10 to 0.47 ◦C/decade. The basin averaged warming rate was approximately 0.342 ± 0.047 ◦C/decade. The highest trend was found in the deep area of the Red Sea (i.e., depth >~1200 m), while lower values were found in the Gulf of Suez, the Gulf of Aqaba, and the southern part of the Red Sea.

The EOF analysis (Figure 4) indicated that the highest SST variability was found over the central and northern part of the Red Sea, while the lowest variability was found over the Gulfs of Suez and Aqaba, as well as the southernmost part of the Red Sea. This could be attributed to the persistence of monsoon winds from the Arabian Sea in the southern Red Sea and northwesterly winds from the Sinai Peninsula in the northern Red Sea [41].

In Figure 5, the authors investigated the mean and spatial trend of the main characteristics of MHWs (frequency, days, duration, and intensities) for the Red Sea over the study period (1982–2019). The maximum MHW frequency values of >2 counts were found in the coastal areas of Sudan, Eritrea, Yemen, and Saudi Arabia. The MHW frequency trend values were statistically significant (*p* < 0.05) over the whole Red Sea. According to Figures 1 and 5b, the deep area (i.e., depth > 1000 m) between Eritrea and the south coast of Saudi Arabia has the highest MHW frequency trend (>1.8 counts/decade). The MHW total days' spatial patterns (Figure 5c) show maximum values of more than 25 days in the Gulfs of Suez and Aqaba. This may be due to the fact that the Gulfs of Suez and Aqaba are well resolved using high-resolution OSTIA data. Figure 5b,d show a similarity among the *p*-values between the MHW frequency and total days trend fields. We found the maximum MHW duration spatial trend significant values (*p* < 0.05; >8 days/decade) to be in the middle of the Red Sea (i.e., between ~lat 18◦ E and 21◦ E), as demonstrated by Figure 5f. Maximum MHW icum values of >22 (◦C days) were found in the northern Red Sea, excluding the Gulfs of Suez and Aqaba (i.e., between lat 22◦ E and 28◦ E) (Figure 5g). The MHW imax shows high values (>1.8 ◦C) throughout most of the Red Sea (Figure 5i), with the exception of the southernmost part near Bab El-Mandab. The MHW matrices for the Red Sea region show a clear ascending trend over the last decade (Figure 6). The highest annual MHW duration and icum were found during the year 2018 (~21 days and 30 ◦C days). The correlation coefficients between the annual mean SST, MHW duration, and MHW icum were 0.64 and 0.72, respectively.

The MHW duration event of 236 days was detected in the Gulf of Suez in 2019 from 10 May 2019 to 31 December 2019 (Figure 7a,b). The high MHW imax event of 4.21 ◦C was discovered in May 2019 off the Sudanese coast (i.e., between ~lat 18.5◦ E and 20.2◦ E), (Figure 7c,d).

#### **6. Conclusions**

This study provides a comprehensive analysis of SST variability, trends, and spatiotemporal patterns of marine heat waves in the Red Sea from 1982 to 2019. To the authors' knowledge, this is the first study to examine marine heat waves and SST trends in the Red Sea using high-resolution data. A statistically significant SST trend was observed over the Red Sea. The EOF analysis showed that the highest interannual variability was detected over the central and northern part of the Red Sea, while the lowest variability was found over the Gulfs of Suez and Aqaba and the southernmost part of the Red Sea.

We found that all MHW characteristic (frequency, duration, and cumulative intensity (icum)) trends increased over the study period. The highest annual MHW frequencies were detected in the years 2018, 2019, 2010, and 2017. There were no MHWs found during the year 1984. A strong correlation was found between the annual MHW frequency and the annual mean SST. Over the period from 2000 to 2019, we discovered that the average MHW frequency and duration increased by 35% and 67%, respectively. Further research is needed to investigate the causes of MHWs and their effects on marine life in the Red Sea. We believe that our study is just the beginning of a series of studies on this hot topic.

**Author Contributions:** Conceptualization, O.I., B.M. and H.N.; methodology, H.N., B.M. and O.I.; formal analysis, O.I., H.N. and B.M.; investigation, O.I., B.M. and H.N.; resources, B.M., O.I. and H.N.; data curation, B.M., H.N. and O.I.; writing—original draft preparation, O.I., B.M. and H.N.; writing—review and editing, H.N., O.I. and B.M.; visualization, B.M., O.I. and H.N.; supervision, H.N., B.M. and O.I. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to express their appreciation to the two anonymous reviewers for their constructive comments.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**

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

**Omneya Ibrahim 1, Bayoumy Mohamed 1,2and Hazem Nagy 1,3,\***

**Abstract:** Marine heatwaves (MHWs) can cause devastating impacts on marine life. The frequency of MHWs, gauged with respect to historical temperatures, is expected to rise significantly 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 find that the mean MHW frequency and duration increased by 40% and 15%, respectively. In the last decade, the shortest significant 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 (icum) trend increased by 5.4 ◦C days per decade. During the study period, we discovered that the maximum significant 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.

**Keywords:** marine heatwaves; sea surface temperature; Eastern Mediterranean Sea; maximum intensity; wind stress; mean sea level pressure

#### **1. 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–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–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–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,

**Citation:** Ibrahim, O.; Mohamed, B.; Nagy, H. Spatial Variability and Trends of Marine Heat Waves in the Eastern Mediterranean Sea over 39 Years. *J. Mar. Sci. Eng.* **2021**, *9*, 643. https://doi.org/10.3390/jmse9060643

Academic Editor: Francisco Pastor Guzman

Received: 21 May 2021 Accepted: 9 June 2021 Published: 10 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

and air-sea exchanges [14–16]. That phenomenon led to an anomalous SST warming in the EMED area during the heatwave of 2003, on the order of +5 ◦C above climatology [17].

Reference [18] investigated the future evolution (1976–2100) of SST and marine heatwaves in the Mediterranean Sea using a local, climatological 99th percentile threshold based on historical-climate daily SST (1976–2005). The dataset 1982–2012 revealed an annual frequency of 0.8 events per year, with MHWs lasting for a maximum of 1.5 months between July and September, while covering a maximum of 90% of the Mediterranean Sea surface. The longest and most severe event of that period corresponded to the MHW of 2003, with the highest mean intensity and maximum event coverage. The authors have concluded that the MHWs will become stronger and more intense in response to increasing greenhouse gas forcing, particular near the end of the century.

Unfortunately, very few studies on the MHW have been published in the EMED, and the topic is still being debated, as mentioned by [19].

In this study 39 years (1982–2020) of SST data have been analyzed to demonstrate the SST trend and interannual variability based on Empirical Orthogonal Function (EOF). 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.

#### **2. 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.

#### *2.1. 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.

#### *2.2. 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.

**Figure 1.** (**a**) The Eastern Mediterranean Sea with the main geographic features (e.g., Basins) where, the solid red square is the position of Mersa Matruh City, the solid blue circle is the position of Crete, the solid green circle is the position of Rhodes, the solid yellow circle is the position of Cyprus and the solid continuous red lines delineate the boundaries between EMED basins. (**b**) Bathymetry according to GEBCO bathymetric dataset (https://www.gebco.net/ (accessed on 20 March 2021)) with 10 m as the minimum depth.

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, (P1, 1 January 1982–31 December 1990), (P2, 1 January 1991–31 December 2000), (P3, 1 January 2001–31 December 2010) and (P4, 1 January 2011–31 December 2020). The chosen of the first period P1 is based on satellite NOAA OI SST V2.1 data availability. Then daily climatology was calculated from the P1, P2, P3 and P4 years employed in the analysis and smoothed with a monthly centred-moving mean.

#### *2.3. 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,

$$\text{T}\_{\text{m}}(\text{j}) = \sum\_{\text{y}=\text{y}\_{\text{s}}}^{\text{y}\_{\text{e}}} \sum\_{d=j-5}^{j+5} \frac{T\left(y, d\right)}{\text{11 (y}\_{\text{e}} - \text{y}\_{\text{s}} + 1\right)}\tag{1}$$

where, Tm 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, ys and ye 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, T90 (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:

$$\mathcal{T}\_{\mathcal{90}}\text{ (j)} = \mathcal{P}\_{\mathcal{90}}\text{ (\%)}\tag{2}$$

According to [5] definition, MHWs occur when SSTs exceed the threshold, defined as the 90th percentile of SST variations (T90 (j)) based on a 39-year climatological period (1982–2020), 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),

$$\mathbf{t}\_{\mathbf{c}} - \mathbf{t}\_{\mathbf{s}} \ge \mathbf{5} \tag{3}$$

where te, ts: dates on which a MHW begins and ends. And D is the MHW Duration in Days that temperature exceeds the threshold,D=te − ts.

Where P90 is the 90th percentile and

$$\mathbf{X} = \langle \mathbf{T}(\mathbf{y}, \mathbf{d}) \mid \mathbf{y}\_{\circ} \text{ <=} \\ \mathbf{y} \prec= \mathbf{y}\_{\circ}, \text{ j} \\ - \ \mathbf{5} \prec= \mathbf{d} \ \lnot= \mathbf{j} + \mathbf{5} \rangle \tag{4}$$

Then, the MHWs main characteristics are (imax) in ◦C, mean intensity (imean) in ◦C and (icum) in (◦C Days):

$$\mathbf{i}\_{\text{max}} = \max \left( \mathbf{T}(\mathbf{t}) - \mathbf{T}\_{\text{m}} \text{ (j)} \right) \tag{5}$$

$$\mathbf{i\_{mean}} = T(t) - \mathrm{Tm}(j)\mathbf{During the MHW event} \\ \mathbf{t\_{e}} \le j \le t\_{\mathbf{e}} \text{ and } j(\mathbf{t\_{e}}) \le j \le j(\mathbf{t\_{e}}) \quad (6)$$

$$\mathbf{i\_{cum}} = \int\_{\mathbf{t\_s}}^{\mathbf{t\_{s-1}}} \left( T(\mathbf{t}) - T\_\mathbf{m}(\mathbf{j}) \right) \, dt \tag{7}$$

icum: sum of daily intensity anomalies.

According to [5] during the MHW event:

$$\mathbf{T(t) > T\_{90} \text{ (j), t.s. is the time, t.s.}}$$

while during time (t − 1) presents no MHW event

$$\mathbf{T(t-1) \mathbf{t\_{s}}$$

$$\text{T(t)} < \text{T}\_{90} \text{ (i)}$$

#### *2.4. Atmospheric Variables Data*

Atmospheric variables were obtained from the European Centre for Medium-Range Weather Forecasts (ECMWF) ERA5 [32]. Winds at 10 m elevation (U10 and V10), 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/m2),

$$
\pi = \mathbb{C}\_{\text{D}} \text{ } \rho\_{\text{a}} \mid \mathbb{W} \mid \text{ } \* \text{ W} \tag{8}
$$

where ρa = 1.2 (kg/m3) is the density of air, CD 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–20 May 2020 MHW events which have the maximum intensity over the whole study period (1982–2020).

#### **3. Description of the SST Data (1982–2020)**

Figure 2 shows the temporal evolution of daily mean climatology SST (◦C) for P1, P2, P3 and P4. The figure clearly shows that P1 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 P2, 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 P3 were 26.52 ◦C, 15.44 ◦C, and 20.48 ◦C, respectively. During P4, the maximum value reached 27.07 ◦C and the minimum was 16.76 ◦C with a mean value of 20.92 ◦C.

**Table 1.** The basin average, minimum and maximum climatological SST ◦C during P1, P2, P3 and P4.

In summary, the maximum values over the four decades were detected during P4 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 P4 was the warmest decade by about 0.94 ◦C within the entire 1982–2020 study period.

**Figure 2.** Temporal evolution of the basin averaged plotted SST in (◦C) where P1 (1982–1990) in continuous blue line, P2 (1991–2000) in continuous green line, P3 (2001–2010) in continuous cyan line and P4 (2011–2020) in continuous red line.

The authors used anomaly maps to clearly identify cooling and warming regions over the aforementioned time periods.

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–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–43]. Hence, warming in anticyclonic areas is less than in the cyclonic ones due to strong stratification in anticyclonic regions [44,45].

**Figure 3.** The EMED SST climatological anomaly maps where (**a**) P1 (1982–1990), (**b**) P2 (1991–2000), (**c**) P3 (2001–2010), and (**d**) P4 (2011–2020). Anomaly maps are calculated by subtracting the climatological mean value over (1982–2020) from each period value (P1, P2, P3 and P4) at the same location (grid).

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.

#### *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.

**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.

**Figure 5.** *Cont.*

**Figure 5.** (**a**) Spatial distribution for the principal (EOF) component (PC1). (**b**) Temporal evolution for EOF (PC1) in normalised variance units. (**c**) Spatial distribution for the second (EOF) component (PC2). (**d**) Temporal evolution for EOF (PC2) in normalised variance units; for the EMED SST anomaly covering the period from 1982 to 2020. In (**b**,**d**), the deseasoned monthly PCs time series is shown in continuous black lines, and the low passed using a 13-month running mean is shown in continuous red lines.

#### **4. 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 Figure 6c–f. The MHW mean duration 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) are found in MHW duration at southeast of Crete, off Mersa Matruh city and the north Ionian Sea region. The maximum significant (*p* < 0.05) MHW icum values found in the Levantine Sea > 30 ◦C Days with a trend >15 ◦C Days per decade (Figure 6e,f). The MHW imax patterns reveal maximum significant values > 1.5◦ with a trend > 0.25 ◦C throughout most of the Levantine Sea with an exception of south east corner and off Mersa Matruh city (see Figure 6g,h).

**Figure 6.** The Marine heatwave properties over the 1982–2020 period, averaged time series of (**a**) mean MHW frequency (Count), (**b**) MHW frequency trend (Count/Decade), (**c**) mean MHW duration (Days), (**d**) MHW duration trend (Days/Decade), (**e**) mean MHW (icum) (◦C Days), (**f**) MHW (icum) trend (◦C Days/Decade), (**g**) mean MHW (imax) (◦C), (**h**) MHW (imax) trend (◦C/Decade), and the black dots indicate the change is not significantly different (*p* > 0.05). Where, the solid red square is the position of Mersa Matruh City, the solid blue circle is the position of Crete, the solid green circle is the position of Rhodes, and the solid yellow circle is the position of Cyprus.

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].

#### *4.1. 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.

**Figure 7.** The EMED MHWs mean frequency (Count) for (**a**) P1 (1982–1990), (**b**) P2 (1991–2000), (**c**) P3 (2001–2010) and (**d**) P4 (2011–2020), and the black dots indicate the change is not significantly different (*p* > 0.05).

#### *4.2. 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.

**Figure 8.** The EMED MHWs mean duration (days) for (**a**) P1 (1982–1990), (**b**) P2 (1991–2000), (**c**) P3 (2001–2010) and (**d**) P4 (2011–2020), and the black dots indicate the change is not significantly different (*p* > 0.05).

#### *4.3. 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].

#### *4.4. 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.

**Figure 9.** The EMED MHWs mean cumulative intensity (◦C days) for (**a**) P1 (1982–1990), (**b**) P2 (1991–2000), (**c**) P3 (2001–2010) and (**d**) P4 (2011–2020), and the black dots indicate the change is not significantly different (*p* > 0.05).

**Figure 10.** The EMED MHWs imax ( ◦C) for (**a**) P1 (1982–1990), (**b**) P2 (1991–2000), (**c**) P3 (2001–2010) and (**d**) P4 (2011–2020), and the black dots indicate the change is not significantly different (*p* > 0.05).

#### **5. Extreme MHWs Events in the EMED over 1982–2020**

Table 2 shows the significant EMED extreme MHW events characteristics (mean duration, imax, imean, and icum) over the study period (1982–2020) 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 of 58.77 ◦C days. The most severe event (icum of 319.65 ◦C days) was still ongoing at the time of the analysis with 126 days, imax of 3.41 ◦C and imean of 2.54 ◦C, 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.

**Table 2.** The EMED extreme significant (*p* < 0.05) MHW events characteristics (mean duration (Days), imax ( ◦C), imean ( ◦C), and icum ( ◦C days)) over the study period (1982–2020) with the corresponding MHWonset, MHWend, Longitude and Latitude in degrees where the highest events are highlighted in yellow color.

**Figure 11.** The EMED extreme MHW significant events over the study period (1982–2020), The SST climatology from NOAA OI SST for the detection of MHWs (blue), 90th percentile MHW threshold (green), and SST time series (black) for each MHW. The pink filled area indicates the period associated with the identified MHW, where (**a**) MHW extreme highest mean intensity event (imean 3.92 ◦C) (24-04-2013 to 08-05-2013), (**b**) MHW extreme highest cumulative intensity event (icum 319.65 ◦C days) (28-08-2020—ongoing).

#### *The EMED Most Intense [14–20 May 2020] MHW Event over the Study Period (1982–2020)*

In this subsection, the authors will describe in detail the most significant intense event of 14–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–20 May 2020] MHW significant event characteristics, including duration, imax, imean, icum, wind stress magnitude, 2m air temperature and mean sea level pressure. The imax 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), imax ( ◦C), imean ( ◦C), icum ( ◦C days)), wind stress magnitude (N/m2), 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–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/m2) 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–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. 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.

**Figure 12.** The EMED most significant intense (14–20 May 2020) MHW events over the study period (1982–2020), where, (**a**) MHW mean SST anomaly ◦C. (**b**) SST climatology from NOAA OI SST for the detection of MHWs (blue), 90th percentile MHW threshold (green), and SST time series (black) for each MHW. The pink filled area indicates the period associated with the identified MHW event. (**c**) Wind stress vectors combined with magnitude for the MHW in N/m2 based on ERA5 hourly data. (**d**) 2 m mean air temperature ◦C based on ERA5 hourly data. (**e**) Mean sea level pressure (mb) based on ERA5 hourly data.

#### **6. 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–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–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, (icum) and (imax)) overlaid with no significant values (*p* > 0.05) for the entire study period. According to Figure 6c–f we found the patterns of icum (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 icum 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 imax 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 imax 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 imean temperature of 3.92 ◦C and lasted 15 days (24 April–08 May) with an imax of 5.06 ◦C and icum of 58.77 ◦C days (Figure 11a and Table 2). While, the most severe event (icum of 319.65 ◦C days) was still ongoing at the time of the analysis with 126 days, imax of 3.41 ◦C and imean 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–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/m2) (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–16,18].

#### **7. 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 P4 (2011–2020) which agrees with [18]. During last decade (P4 (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 icum 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.

The authors have showed the significant EMED extreme MHW events characteristics (mean duration, imax, imean, and icum) over the study period (1982–2020). During these extreme events, there was a similarity in 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/m2) at the same locations as the most intense MHW events. High air temperature values (>25 ◦C) have been linked to the locations of MHW extreme events. We concluded that atmospheric conditions such as mean sea level pressure, wind stress, and air temperature could be related to the existence of MHW.

In the future, we should focus more on the anthropogenic factors that contribute to the increasing severity and recurrence of MHW events. Furthermore, possible MHW links to different climatic and atmospheric indices [13,62–64] such as North Atlantic Oscillations (NAO), East Atlantic (EA), East Atlantic–West Russia (EAW), Mediterranean Oscillation (MO) and North-Sea Caspian Pattern (NCP). There is a lot of speculation about what caused the weather pattern we saw during this most severe MHW event. Was it just a coincidental occurrence? More information about the detection, characterization, impact assessment, and prediction of MHWs could be obtained.

**Author Contributions:** Conceptualization, O.I., B.M. and H.N.; methodology, H.N. and O.I.; formal analysis, O.I., H.N. and B.M., investigation, O.I., B.M. and H.N.; resources, O.I., B.M. and H.N.; data curation, H.N. and O.I.; writing—original draft preparation, O.I., B.M. and H.N.; writing—review and editing H.N., O.I. and B.M.; visualization, O.I., B.M. and H.N.; supervision, H.N. and O.I.; project administration, funding acquisition. This research received no external funding. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Not applicable.

**Acknowledgments:** The authors would like to thank Nadia Pinardi for her constructive remarks. The authors also would like to thank Khaled Elhusseiny, MSc Translation, York University and Doha Shahin, Language Interpreters Instructor at Conestoga College for their diligent proofreading of this article. Finally, we would like to express our sincere gratitude to Luke Marsden from UNIS in Norway for his wonderful collaboration in editing our manuscript. Special thanks to the three anonymous reviewers for their constructive comments, which helped to improve the manuscript.

**Conflicts of Interest:** No conflict of interest.

#### **References**

## *Article* **CMIP5-Based Projection of Decadal and Seasonal Sea Surface Temperature Variations in East China Shelf Seas**

**Huiqiang Lu 1, Chuan Xie 2, Cuicui Zhang 1,\* and Jingsheng Zhai 1,\***

**Abstract:** The East China Shelf Seas, comprising the Bohai Sea, the Yellow Sea, and the shelf region of East China Sea, play significant roles among the shelf seas of the Western North Pacific Ocean. The projection of sea surface temperature (SST) changes in these regions is a hot research topic in marine science. However, this is a very difficult task due to the lack of available long-term projection data. Recently, with the high development of simulation technology based on numerical models, the model intercomparison projects, e.g., Phase 5 of the Climate Model Intercomparison Project (CMIP5), have become important ways of understanding climate changes. CMIP5 provides multiple models that can be used to estimate SST changes by 2100 under different representative concentration pathways (RCPs). This paper developed a CMIP5-based SST investigation framework for the projection of decadal and seasonal variation of SST in East China Shelf Seas by 2100. Since the simulation results of CMIP5 models may have degrees of errors, this paper uses hydrological observation data from World Ocean Atlas 2018 (WOA18) for model validation and correction. This paper selects seven representative ones including ACCESS1.3, CCSM4, FIO-ESM, CESM1-CAM5, CMCC-CMS, NorESM1-ME, and Max Planck Institute Earth System Model of medium resolution (MPI-ESM-MR). The decadal and seasonal SST changes in the next 100 years (2030, 2060, 2090) are investigated by comparing with the present analysis in 2010. The experimental results demonstrate that SST will increase significantly by 2100: the decadal SST will increase by about 1.55 ◦C, while the seasonal SST will increase by 1.03–1.95 ◦C.

**Keywords:** decadal and seasonal SST variation; East China Shelf Seas; CMIP5; WOA18

#### **1. Introduction**

The climate is changing. Our Earth is warming up. Many agree that climate change may be one of the greatest threats facing our planet. Ocean warming, which contributes much to global warming, has become a more serious problem. Recent observation-based estimates and model simulation results show that ocean warming is accelerating and will continue in the next one hundred years [1]. Since global warming may lead to many ecological problems, more efforts need to be made in the assessment and projection of the warming rate of the oceans in the future. The sea surface temperature (SST), which is an important physical parameter of oceans, can reflect the effect of climate change. The estimation of SST variation has become a hot research topic in marine science.

The East China Shelf Seas, which consist of the Bohai Sea, the Yellow Sea, and the shelf region of East China Sea, have been recognized as the most significant marginal seas in the Western North Pacific Ocean (WNPO) [2]. The Bohai Sea is a shallow semi-enclosed sea with an average depth of only 18 m, and the changes of its water temperature have a greater influence on its ecosystem than the physical forcing from the external oceans [3]. The Yellow Sea is a wide and shallow sea, where the depth in most regions is less than 50 m. Its water column can be seriously affected by the atmospheric conditions such as heating,

**Citation:** Lu, H.; Xie, C.; Zhang, C.; Zhai, J. CMIP5-Based Projection of Decadal and Seasonal Sea Surface Temperature Variations in East China Shelf Seas. *J. Mar. Sci. Eng.* **2021**, *9*, 367. https://doi.org/10.3390/ jmse9040367

Academic Editor: Kyung-Ae Park

Received: 31 January 2021 Accepted: 23 March 2021 Published: 30 March 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

cooling, and wind stress within it than from the open seas [4]. The East China Sea, which is located between the largest continent and the largest ocean, has been largely effected by the climatic forcing from both the high-latitude Northern Hemisphere (East Asian Monsoon System) and the tropic ocean (Kuroshio Current (KC)) and generates characteristic climate patterns with strong horizontal and vertical temperature gradients [5–7]. These East China Shelf Seas are very sensitive to global warming because of their shallow waters [8]. The SST increase over these regions is about 0.8 to 2 ◦C per century, which is nearly twice the global average increase of SST [9]. The SST increases can not only affect the metabolic rates of marine organisms, but also influence other oceanic states, such as local currents [3]. Therefore, we need to analyze the SST variation, especially to project the long-term spatial and temporal variations in the future. However, this is a very difficult task because of the lack of effective long-term projection data [8].

Recently, with the rapid development of ocean simulations, model intercomparison projects have become important ways to investigate climate changes. Among them, the Coupled Climate Model Intercomparison Project plays an important role [10]. Phase 5 of the Climate Model Intercomparison Project (CMIP5) is an international collaboration framework, which provides a multimodel context to help understand the responses of climate models to a common forcing with the aim of promoting the climate model projection and assessment for the Fifth Assessment Report (AR5) of the Intergovernmental Panel on Climate Change (IPCC) [11]. CMIP5 contains multiple models, which can be used to project climate changes and sea level rise in the future under different climate change scenarios or representative concentration pathways (RCPs) including RCP2.7, RCP4.5, RCP6.0, and RCP8.5. These scenarios correspond to the peak of the atmosphere radiative imbalance of 2.6, 4.5, 6.0, and 8.5 W/m2 by 2100, respectively [12]. Following the previous CMIPs, CMIP5 is a new experimental framework, which can be widely applied for the analysis of decadal and seasonal climate changes in the future.

In the literature, several CMIP5-based SST projection works have been developed. For example, Zhou and Ying [13] analyzed the interannual variability of the SST over the Pacific in the historical simulation and future analysis under RCP4.5 and RCP8.5. Qu and Huang [14] investigated the decadal variability of the tropical Indian Ocean SST–South Asian High (SAH) relation, as well as its response to global warming. Qin and Xie [15] studied the connections between the precipitation extremes during 1953– 2002 in the dry and wet regions of China and the SST in the eastern tropical Pacific Ocean based on two sets of observation data, 17 CMIP5 models, and nine regional climate models. Tachibana et al. [16] examined the western Indian Ocean SST biases among the CMIP5 models and found that the multimodel ensemble mean SST biases over the western equatorial Indian Ocean are warmer than the observations during the summer monsoon season. Song et al. [17] evaluated 18 CMIP5 models according to their capability of simulating the SST annual cycle in the eastern equatorial Pacific. Xu et al. [18] tested a number of previously proposed mechanisms responsible for the southeastern tropical Atlantic SST bias based on CMIP5 models. Zhao and Zhang [19] analyzed the impacts of SST warming in the tropical Indian Ocean on the projected change in summer rainfall over Central Asia based on historical and RCP8.5 experiments. Kucharski and Joshi [20] evaluated the teleconnection from the tropical South Atlantic SST anomalies to the Indian monsoon based on observations and CMIP5 model data. Levine et al. [21] examined the extent and impact of cold SST biases developing in the northern Arabian Sea in the CMIP5 multi-model ensemble. Langehaug et al. [22] investigated the projection of SST in the Nordic Seas and Barents Sea using initialized hindcast simulations performed with three climate models (MPI-ESM-LR, CNRM-CM5, IPSL-CM5-LR). For the Chinese coastal seas, a few works have also been presented to study the SST changes using CMIP5 models. For example, Huang et al. [23] evaluated the capacities of 17 selected CMIP5 models on the historical SST simulation in the South China Sea and projected the SST changes in the 21st Century under RCP2.6, RCP 4.5, and RCP8.5, respectively. Tan et al. [24] evaluated the variation trend of SST over offshore China in the 21st Century based on the selected CMIP5

models under RCP4.5. Song et al. [25] assessed the monthly, seasonal, and interannual SSTs in the China seas over 1960–2002 using five representative CMIP5 models.

Although existing works have made some achievements in the study of the SST changes in the East China Shelf Seas, there are still some problems unsolved. Firstly, existing works mainly focused on the study of annual and decadal variations of the SST in the East China Shelf Seas without investigating the seasonal variations, which may affect ocean organisms and the ecological environment more seriously (e.g., the changes of species distribution and the move up of the phenophase). Secondly, existing works lack precise ocean observation data for model validation. Existing works mainly used the low-resolution observation data HadISSTfrom the Hadley Center in the U.K. (1◦ × 1◦). It is not accurate enough for the evaluation of high-resolution CMIP5 models (e.g., the Max Planck Institute Earth System Model of medium resolution (MPI-ESM-MR) of 0.1◦ × 0.1◦). Moreover, existing works only used the observation data for model validation without establishing model correction or simulation result modification, which is a fundamental procedure to make projections more accurate. Thirdly, the models used in existing works are not very appropriate. Several high-resolution models (e.g., the MPI-ESM-MR) were not fully evolved in the existing works. Fourthly, the study region of existing works was very rough. Existing works mainly referred to the Chinese coastal seas as the rectangle region of 0◦ N–45◦ N, 100◦ E–140◦ E, which is a very coarse region to study on the SST changes in the East China Shelf Seas.

To solve these problems, this paper establishes an SST analysis framework, which can be used for the projection of both decadal and seasonal SST variation in the East China Shelf Seas using high-resolution CMIP5 data. This paper investigates both the decadal and seasonal SST variation. For the second problem, we introduce the high-resolution hydrological observation data WOA18 (0.25◦ × 0.25◦) from World Ocean Atlas 2018 [26] for model validation and calculate the error map of each of the CMIP5 model data for the simulation result modification. This strategy ensures the high accuracy of CMIP5 model data. To solve the third problem, we select the best model MPI-ESM-MR with the highest resolution (0.1◦ × 0.1◦) from seven representative models (ACCESS1.3, CCSM4, FIO-ESM, CESM1-CAM5, CMCC-CMS, NorESM1-ME, and MPI-ESM-MR). MPI-ESM-MR not only has low errors, but its high resolution also guarantees the specific analysis of local climate changes. To solve the fourth problem, we utilize the bathymetric depth map with a finer grid to define the study area (as shown in Figure 1). In this paper, the decadal and seasonal SST variation analysis is performed by comparing the SST simulation result on 2030, 2060, and 2090 with the present analysis in 2010 under RCP4.5. We use RCP4.5 as a representative scenario because it is a medium-mitigation emission scenario that stabilizes direct radiative forcing at 4.5 W/m2 (650 ppm CO2 equivalent) by 2100 [27,28]. Compared to RCP4.5, RCP2.6 is a mitigation scenario leading to a very low forcing level (at 2.6 W/m2 by 2100), while RCP8.5 is a scenario with very high greenhouse gas emissions (at 8.5 W/m2 by 2100). RCP4.5 is neither too high, nor too low for projection. To the best of our knowledge, this is the first work to investigate both decadal and seasonal SST variations in the East China Shelf Seas by 2100 using high-resolution CMIP5 data.

This paper is organized as follows: Section 2 introduces the data and the study area. The model validation and simulation result modification are performed by comparing the CMIP5 model data with WOA18 data in Section 3. Sections 4 and 5 analyze the projection results of decadal and seasonal SST variations, respectively. Section 6 concludes this paper.

**Figure 1.** Depth map of East China Shelf Seas.

#### **2. Data and Study Region**

In this paper, CMIP5 model data were used for the projection of future SST changes, while WOA18 data were used for model validation, as well as for present analysis. CMIP5 consists of 36 models, which are developed by 16 institutes. In the preprocessing procedure, we selected the most representative models according to the following criteria: (1) for each institute, we selected the latest developed, with the highest resolution model for study; (2) the selected models should provide the simulation results of the SST from 2010 to 2100 continuously. We finally chose seven representative models. Their information is shown in Table 1. From Table 1, we can see that the resolution of Max Planck Institute Earth System Model of medium resolution (MPI-ESM-MR) is much higher than that of other models. MPI-ESM-MR couples general circulation models for the ocean and the atmosphere. It has been widely applied in many climate change experiments for either idealized CO2-only forcing or forcings based on observations and RCP scenarios [29].

**Table 1.** CMIP5 models used in this paper: model name, average horizontal resolution (latitude × longitude), and reference. MPI-ESM-MR, Max Planck Institute Earth System Model of medium resolution.

The annual and seasonal SST data in 2010 were obtained from the World Ocean Atlas 2018 (WOA18) hydrological observation averaged over 2005 to 2017. WOA18 provides both objectively analyzed (1◦ grid) climatological fields of in situ temperature for annual and seasonal compositing periods of the World Ocean. It also includes associated statistical observation data interpolated on 5◦, 1◦, and 0.25◦ grids [26]. We used the highest resolution data of 0.25◦ grid. As shown in Figure 1, the study area is comprised of the Bohai Sea, the Yellow Sea, and the shelf region of East China Sea.

Since each selected CMIP5 model has its inner simulation variations and errors, we calculated the climatological annual and seasonal water temperature in 2010 by the average result over 2006 to 2015. A similar processing was done for 2030, 2060, and 2090. They were investigated by the average result of 2026∼2035, 2056∼2064, and 2086∼2094, respectively.

#### **3. Model Validation and Simulation Result Modification**

As we know, the simulation results of numeric ocean models may contain some degree of errors. Before using the CMIP5 model data, we needed to perform model validation to evaluate whether the simulation results are accurate enough for SST projection. The model validation was performed by comparing the annual and seasonal CMIP5 model data with the hydrological observation data of WOA18 in 2010. The qualitative and quantitative comparison results are shown in Figure 2 and Table 2, respectively. From Figure 2, we can find that most selected CMIP5 models have a similar SST distribution as the WOA18 data. However, different models also vary from each other in the model resolution and simulation accuracy. Among these models, the simulation result of MPI-ESM-MR is better than the others. Its resolution (in 0.1◦) is the highest one (even higher than WOA18), and its simulation result is most similar to WOA18. MPI-ESM-MR can illustrate the local climate impact factors such as the Yellow Sea cold water mass in Summer and Autumn, the Yellow Sea warm current in Winter, the Kuroshio invasion, and the northern flow of the Taiwan warm current around Changjiang shore to the Tsushima Strait. Benefiting from the high resolution, its simulation result is even better than the observation data (WOA18). The numerical results in Table 2 also verify this finding. From Table 2, we can see that both the annual and seasonal simulation errors of MPI-ESM-MR are low enough. That means that MPI-ESM-MR is the best model that can be used for the projection of decadal and seasonal SST variations.

**Table 2.** Comparison of CMIP5 model data with the real observation data from WOA18 in 2010.

After model validation, we performed simulation result modification to make the CMIP5 model data more accurate. Firstly, we calculated the error map of MPI-ESM-MR by comparing its results with WOA18 data for the year and the four seasons of 2010, respectively. Then, we subtracted these errors from the original simulation results to obtain more accurate projection data. Figure 3 illustrates the annual and seasonal error maps of MPI-ESM-MR. From Figure 3a–e, we can see that the error of MPI-ESM-MR is mainly concentrated in the deep trough of the Yellow Sea, where its simulation on the influence of the Yellow Sea warm current is stronger than the observation result. This phenomenon leads to the higher SST simulation in the eastern Yellow Sea in spring and winter and in the northern part of East China Shelf Seas in autumn than the WOA18 data. However, its simulation results of other regions were satisfactory enough for projection.

**Figure 2.** Model validation: the comparison of CMIP5 simulation results with WOA18. The first row to seventh row illustrate the result of seven CMIP5 models: ACCESS1.3, CCSM4, FIO-ESM, CESM1-CAM5, CMCC-CMS, NorESM1-ME, and MPI-ESM-MR, respectively. The last row shows the observation result from WOA18 data. The first to fifth column present the results in the year round, spring, summer, autumn, and winter of 2010, respectively.

**Figure 3.** The error maps of MPI-ESM-MR in 2010: (**a**) annual, (**b**) spring, (**c**) summer, (**d**) autumn, and (**e**) winter.

#### **4. Projection Results of Decadal SST Changes**

We analyzed the decadal SST changes in 2030, 2060, and 2090 compared with the present analysis in 2010. The modified projection results of the seven models and their average results are shown in Table 3 and Figure 4, respectively. From Table 3 and Figure 4, we find that most models demonstrate a significant SST increase in the next 100 years. From Table 3, we also find that the simulation results of MPI-ESM-MR are most similar to the average of all models. Moreover, its resolution is much higher than that of other models. Since the resolution of the other six models is too low to distinguish local climate impact factors such as the Kuroshio system, it is inappropriate to use them to project the SST changes in the future. Therefore, next, we used the modified result of MPI-ESM-MR to analyze the decadal SST changes. From Table 3, we can see that the SST increases from 18.09 ◦C in 2010 to 18.68 ◦C in 2030, 19.78 ◦C in 2060, and 19.65 ◦C in 2090. In particular, a remarkable SST increase was obtained from 2030 to 2060, which was 1.10 ◦C. Till 2090, the SST will increase by 1.55 ◦C. Besides the decadal analysis, we also investigated the annual changing rate of SST from 2010 to 2090. The changing rate between 2030 and 2010, 2060 and 2030, and 2090 and 2060 is illustrated in Table 3. We can see that this rate is very small from 2010 to 2030, achieves the top value from 2030 to 2060, and becomes low from 2060 to 2090. This suggests that the SST increase becomes stable from 2060 to 2090.

Then, we analyzed the spatial SST changes in different regions through the qualitative results shown in Figure 5. Figure 5a,c,e shows the SST projection in 2030, 2060, and 2090 and their changes compared to 2010 in (b), (d), and (f), respectively. From Figure 5, we can find that the highest SST increase rate is obtained in 2060 (in Figure 5d), and the increase becomes stable in 2090 (in Figure 5f). Although the SST increase in 2030 (in Figure 5b) is not as significant as that in 2060, its distribution is inhomogeneous. The significant SST increases mainly concentrate in the east of Bohai (especially east of Qinhuangdao), the Changjiang Estuary, and its adjacent shore of Jiangsu Province. The increment can reach to as high as 1.5 ◦C. From Figure 5d, we can see that compared to 2030, the increase of the SST in 2060 expands to almost all the regions and at an even higher rate. The increment reaches about 3.5 ◦C in the Changjiang Estuary and in the outer shelf of the East China Sea. From Figure 5f, we can see that the SST variation in 2090 is similar to that of 2060 and

becomes stable gradually. One remarkable thing to note is that the SST in the northern and central Yellow Sea is relatively lower than that in the other regions (see Figure 5a,c,e), and their SST increases are not very obvious (see Figure 5b). This is mainly due to the presence of the cold water mass in the Northern and Central Yellow Seas, which has a great influence on the distribution of the SST. The increase of the SST can severely influence marine ecosystems and cause many ecological problems, especially for shelf regions with shallow water. For example, the Brown tide has broken out in the Qinhuangdao coastal area of Bohai recurrently since 2009. The Changjiang Estuary and the shore in Northern Jiangsu Province suffered from the Red Tide and Green Tide, and jellyfish disasters often have also occurred in these areas recently. The increase of the SST will intensify these problems through the following mechanisms. First, it can affect the metabolic rate of marine organisms. Second, it can influence the other oceanic states, such as local currents. Third, it can further affect the water column stratification, substrate structure, photosynthetic light intensity, and nutrient cycling [3]. The species distribution can also be disturbed by the increase of the SST. For example, the warm and cold water fish stocks can be reduced greatly in the regions with significant SST increases, and tropical ocean organisms may move up to the middle and high latitude areas with warm water [37].

**Figure 4.** The decadal SST variations in the next 100 years.

**Table 3.** The decadal SST variations in the next 100 years.

**Table 3.** *Cont*.

**Figure 5.** *Cont*.

**Figure 5.** The projection results of decadal SST changes by MPI-ESM-MR in (**a**) 2030, (**c**) 2060, and (**e**) 2090 and their changes compared to that in 2010 in (**b**) 2030–2010, (**d**) 2060–2010, and (**f**) 2090–2010.

#### **5. Projection Results of Seasonal SST Changes**

We analyzed the seasonal SST variations in the next 100 years. The numerical results of the seven models and their average results are shown in Table 4. Similar to the decadal analysis, we used the modified simulation result of MPI-ESM-MR to investigate the seasonal SST changes because of its high resolution and low projection error. From Table 4, we can see that the SST increases significantly from 2010 to 2090 for all seasons. The increment reaches about 2 ◦C by 2090. Next, we analyzed the spatial seasonal SST variations. The qualitative simulation results are shown in Figure 6. The first to fourth rows show the SST changes from 2010 to 2030, 2060, and 2090 in spring, summer, autumn, and winter, respectively. We found that in summer, the regions with significant SST increases mainly are concentrated in the east of Qinhuangdao in Bohai, the shore of Northern Jiangsu Province in the Yellow Sea, and the Changjiang Estuary in the East China Sea, which is coincident with that in the decadal analysis. This is mainly because the upwelling waters become weak and the Kuroshio invasion becomes stronger in summer. In autumn, the increase of the SST becomes smoother, and the regions with a significant SST increase move from the inner shelf area to the middle and even outer shelf regions in winter. In particular, the outer shelf regions have higher SST increases than other regions. A remarkable demonstration of the SST increase in the Northern Yellow Sea is illustrated in Figure 6f, where the Yellow

Sea warm current plays a more significant role in summer than the Yellow Sea cold water mass. The SST increases in this region will intensify the low oxygen and acidification problems. Similar to the decadal analysis, for all seasons, the SST increase in 2060 (shown in the second column) is more significant than that in 2030 (shown in the first column) and becomes stable in 2090 (shown in the third column). As previously mentioned, the increase of the SST may cause many environmental problems and destroy the marine ecosystem. For example, the increase of the SST in Bohai may lead to the Brown Tide problem, and the SST increases in the shore of Jiangsu Province and the Changjiang Estuary may result in the Red and Green Tide problem and even cause the outbreak of jellyfish disasters in these regions. The increase of the SST in the middle shelf and outer shelf regions of the East China Sea will intensify the marine acidification and ocean hypoxia there. The model validation and modification of this paper is relative simple. We use the difference between the CMIP5 model and the WOA data in 2010 as the model error and substrate these errors from CMIP5 model data to modify their projection results. We should use more sophisticated methods with all the model assumptions and conditions to perform model modification. However, it is not an easy task especially for the future projection with many unknown conditions. We donot know whether these errors will remain constant for the 2010–2100 period. Thus. we include the experimental results without projection modification for comparison in this part. The decadal and seasonal SST changes without modification of model projection results are shown in Table A1 and Table A2, respectively. Comparing Table A1 to Table 3, and Table A2 to Table 4, we find that the bottom part of these tables are the same with each other. That is, the investigation of 35 the SST changes in the future using the differences of 2030, 2060, and 2090 relative to that in 2010 is not influenced by the 36 model validation and modification method in this paper.

**Table 4.** The projection results of seasonal SST variations in the next 100 years using seven CMIP5 models and their average result.

**Figure 6.** Seasonal changes of SST from 2010 to 2030 (the first column (**a**,**d**,**g**,**j**)), 2060 (the second column (**b**,**e**,**h**,**k**)), and 2090 (the third column (**c**,**f**,**i**,**l**)), respectively. The first row to fourth row illustrates the changes in Spring, Summer, Autumn, and Winter, respectively.

#### **6. Conclusions**

This paper investigated both decadal and seasonal SST variations in the East China Shelf Seas by 2100 using CMIP5 models under RCP4.5. Seven representative CMIP5 models were used in this paper, including ACCESS1.3, CCSM4, FIO-ESM, CESM1-CAM5, CMCC-CMS, NorESM1-ME, and MPI-ESM-MR. To make the projection more accurate, we introduced the high-resolution hydrological observation data from WOA18 for model validation and simulation result modification. The best model (MPI-ESM-MR) with highest resolution and low errors was selected for the projection of SST changes. The extensive experimental results demonstrated that both decadal and seasonal SSTs in the East China Shelf Seas will increase significantly in the next one hundred years: the decadal SST will increase by about 1.5 ◦C till 2090, and the seasonal SST will increase by about 1.03 ◦C–1.95 ◦C by 2090. The highest increment was obtained from 2030 to 2060, and it became stable from 2060 to 2090. Although the SST increase in 2030 was not as significant as that in 2060 and 2090, the distribution of SST increases was inhomogeneous. Some local regions had a high SST increase of 1.5 ◦C in the east of Bohai Sea (east of Qinhuangdao), the Changjiang Estuary, and its adjacent shore of Jiangsu Province. Compared to the decadal analysis, the seasonal variation may play a more important role in understanding climate changes. We found that in summer, the significant SST increase mainly concentrated in the eastern area of Bohai, the Changjiang Estuary, and the shore of Jiangsu Province, which is mainly due to the weakness of upwelling, the strengthening of the Kuroshio invasion, and the influence of the Taiwan warm current. In autumn, the increase of SST becomes smoother, and the regions with a significant SST increase move from the inner shelf to the outer shelf regions in winter and spring. The significant SST increase may affect the local marine ecology seriously. It can change the distribution of marine organisms, lead to the moving up of the phenophase, intensify the low oxygen and acidification problems, and cause many environmental problems. It has brought about the Brown Tide in the east of Qinhuangdao in the Bohai Sea, led to the Red and Green Tide in the Changjiang Estuary and the shore of Jiangshu Province, and caused the outbreak of jellyfish disasters in these regions. To the best of our knowledge, this is the first work to investigate both decadal and seasonal SST variations in the East China Shelf Seas using high-resolution CMIP5 model data and observation data.

**Author Contributions:** Conceptualization, C.Z.; Methodology, C.Z., C.X., H.L. and J.Z.; Analysis, C.Z., C.X. and H.L.; Writing and Editing, C.Z., H.L.; Project Administration, C.Z. and J.Z. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the National Natural Science Foundation of China Grant 41806116 and the National Key Research and Development Program of China 2016YFC1401203.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** Data used can be accesses via CMIP5 website.

**Acknowledgments:** H.L. acknowledges the National Natural Science Foundation of China Grant 41806116 and the National Key Research and Development Program of China 2016YFC1401203.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **Appendix A**

**Table A1.** The decadal SST variation in the next 100 years using the raw CMIP5 model data without projection result modification.

The model validation and modification of this paper is relative simple. We use the difference between the CMIP5 model and the WOA data in 2010 as the model error and substrate these errors from CMIP5 model data to modify their projection results. We should use more sophisticated methods with all the model assumptions and conditions to perform model modification. However, it is not an easy task especially for the future projection with many unknown conditions. We donot know whether these errors will remain constant for the 2010–2100 period. Thus. we include the experimental results without projection modification for comparison in this part. The decadal and seasonal SST changes without modification of model projection results are shown in Table A1 and Table A2, respectively. Comparing Table A1 to Table 3, and Table A2 to Table 4, we find that the bottom part of these tables are the same with each other. That is, the investigation of 35 the SST changes in the future using the differences of 2030, 2060, and 2090 relative to that in 2010 is not influenced by the 36 model validation and modification method in this paper.

**Table A2.** The projection results of seasonal SST variations in the next 100 years using seven CMIP5 models and their average result using the raw CMIP5 model data without projection result modification.

#### **References**

MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Journal of Marine Science and Engineering* Editorial Office E-mail<EMAIL_ADDRESS>www.mdpi.com/journal/jmse

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-2601-0