The role of air–sea heat flux for marine heatwaves in the Mediterranean Sea

Recent studies have significantly contributed to understanding physical mechanisms associated with the occurrence of marine heatwaves (MHWs). Building upon prior research, this study investigates the relative role of air–sea heat exchange and oceanic processes during the onset and decline phases of surface MHWs in the Mediterranean Sea based on a joint analysis of remote sensing data and reanalysis outputs over the period 1993–2022. Results show that air–sea heat flux is the major driver in 44 % of the onset and only 17 % of the declining MHW phases. Thus, these findings suggest that oceanic processes play a key role in driving sea surface temperature (SST) anomalies during MHWs, particularly during declines. The role of surface flux becomes more important during warmer months and onset periods. Spatially, the heat flux contribution is greater in the Adriatic and Aegean sub-basins, where it becomes the major driver of most onset phases. Latent heat emerges as the most significant heat flux component in forming the SST evolution across all seasons. Onset and decline phases lasting less than 5 d experience a weaker contribution of heat flux compared to longer phases (lasting 5–10 or more than 10 d). Moreover, an inverse relationship between MHW severity and the contribution of heat flux is observed. At the subsurface, mixed layer shoaling is found over the entire duration of most MHWs, particularly for those of shorter duration. Therefore, the surface cooling right after the peak day is likely not associated with vertical mixing in such cases. These findings suggest that other oceanic processes, potentially horizontal advection, have a key role in modulating SST at the beginning of most MHW declines. In turn, further dissipation of heat is commonly driven by vertical mixing, as indicated by a significant mixed layer deepening after the MHW end day in most cases. This study emphasizes the need to consider subsurface information for future studies of MHWs and highlights the importance of accounting for limitations associated with the definitions employed for MHW phases.

1 Introduction

Marine heatwaves (MHWs) are extreme events, characterized by prolonged periods of anomalously high water temperature, lasting for at least 5 consecutive days (Hobday et al., 2016). These events have gathered increased attention due to their detrimental effects on marine life, especially given the increase observed in their frequency, intensity, and duration over the recent decades, at a global and Mediterranean scale (Oliver et al., 2018; Holbrook et al., 2019; Darmaraki et al., 2019a; Juza et al., 2022; Dayan et al., 2023; Pastor and Khodayar, 2023). Mass mortality events and local extinctions, coral bleaching, and massive shifts of marine species have been extensively reported (Wernberg et al., 2016; Frölicher and Laufkötter, 2018; Smale et al., 2019; Garrabou et al., 2022; Smith et al., 2023), along with socioeconomic impacts on fishery and aquaculture industries (Mills et al., 2013; Cavole, 2016). The intensification of MHW conditions has been attributed mostly to ocean warming (Oliver et al., 2018; Ciappa, 2022), while further intensification is expected in the future (Oliver et al., 2019; Darmaraki et al., 2019b; Plecha and Soares, 2019; Hayashida et al., 2020), driven by anthropogenic forcing and particularly pronounced under high-emission future scenarios (Oliver et al., 2019).

Given these concerns, it is crucial to enhance our understanding of the driving factors behind MHWs at a regional scale, particularly within the framework of exploring predictability options and facilitating marine decision-making (Holbrook et al., 2020; Spillman et al., 2021). Recent research has significantly contributed to identifying physical drivers and MHW-favoring conditions (e.g., Holbrook et al., 2019; Sen Gupta et al., 2020; Oliver et al., 2021; Vogt et al., 2022; Marin et al., 2022). Individual events in the Mediterranean Sea have also been explored, such as the widely known MHW in 2003 in the western Mediterranean Sea (e.g., Sparnocchia et al., 2006; Olita et al., 2007; Bonino et al., 2023), the short-lasting record-breaking MHW in May 2020 in the southeastern Mediterranean Sea (Ibrahim et al., 2021; Denaxa et al., 2022), and the most recent long-lasting MHW in summer 2022 (Marullo et al., 2023; McAdam et al., 2024; Pirro et al., 2024).

However, a limited number of studies have assessed physical drivers separately for the buildup and decay of MHW events, employing different methodologies and datasets. Schlegel et al. (2021) demonstrated that nearly 50 % of surface MHWs in the northwest Atlantic are heat-flux-driven but less than 20 % decay due to heat flux, suggesting that oceanic processes are mainly responsible for the MHW decline. Marin et al. (2022) investigated upper-ocean MHWs based on global ocean circulation model output and found that heat advection, followed by anomalous air–sea heat flux, explains most of the upper-ocean temperature anomalies during both MHW onset and decline phases. Within this context, the present study utilizes high-resolution observational sea surface temperature (SST) and modeled heat flux data to assess the driving role of air–sea heat exchange during onset and decline phases of surface MHWs in the Mediterranean Sea. Furthermore, it provides insights into the concurrent subsurface conditions by examining the mixed layer dynamics during MHWs.

2 Data and methods

MHWs in this study were identified based on high-resolution gridded satellite SST data in the Mediterranean Sea. Daily SST values from the reprocessed and near-real-time datasets from the Copernicus Marine (product ref. nos. 01 and 02 – Table 1; 0.05° × 0.05° and 0.01° × 0.01° horizontal resolution, respectively) were used to cover the period January 1993–December 2022. To study the net surface heat budget (Q_net), turbulent and radiative surface fluxes were obtained from the ECMWF ERA5 Reanalysis dataset at hourly frequency and 0.25° × 0.25° horizontal resolution (product ref. no. 03; Table 1). Finally, daily values of mixed layer depth (MLD) for the period January 1993–December 2022 were obtained from the Mediterranean Sea Physics Reanalysis and the Mediterranean Sea Physics Analysis and Forecast (product ref. no. 04 and no. 05 – Table 1; 0.042° × 0.042° horizontal resolution). Collocation of SST and MLD with the coarser ERA5 data was performed by using the nearest neighboring value to each ERA5 grid point. The paired values at daily frequency were used in the heat budget analysis.

MHW detection was performed based on the definition and detection methodology of Hobday et al. (2016) using the MATLAB toolbox provided by Zhao and Marin (2019). The reference period used in this study to create the daily climatology required for the event detection is the same as the 30-year study period (1993–2022). To compute the climatology, a time window of 11 d was employed, centered on the day when each daily climatological value was computed. For the computation of the daily threshold time series, the 90th percentile was selected, being widely used in MHW studies, thus allowing for a more direct intercomparison with the literature. Additionally, a 30 d window was applied for smoothing the threshold time series.

Basic properties were computed for each identified event (e.g., start and end day, mean and max intensity, duration). Next, events were split into their onset and decline phase. The onset phase was considered to last from the first day until the day of maximum intensity (I_max) and the decline phase from I_max day until the last day of the event. As in Schlegel et al. (2021), for each phase, an ocean mixed layer heat budget analysis was applied to derive the change in SST attributed to Q_net based on the following equation:

$$\begin{array}{}\text{(1)}& {\mathrm{SST}}_{t_{\mathrm{2}}}^{\prime }-{\mathrm{SST}}_{t_{\mathrm{1}}}^{\prime }={\int }_{t_{\mathrm{1}}}^{t_{\mathrm{2}}}{\displaystyle \frac{Q^{\prime }}{{\mathit{\rho }}_{\mathrm{o}}{c}_{p}h}}\mathrm{d}t+R.\end{array}$$

The left-hand side of Eq. (1) represents the observed change in SST anomaly (${\mathrm{DSST}}_{\mathrm{obs}}^{\prime }$) relative to climatology during a specific phase. Each phase starts at day t₁ and ends at day t₂, which are the start day and I_max day for onset phases or I_max day and end day for decline phases, respectively. On the right-hand side of Eq. (1), Q^′ is the daily anomaly of Q_net. The latter consists of the latent and sensible heat flux as well as net shortwave and net longwave radiation (LH, SH, SWR, and LWR, respectively), as follows:

$$\begin{array}{}\text{(2)}& Q_{\mathrm{net}}=\mathrm{LH}+\mathrm{SH}+\mathrm{SWR}+\mathrm{LWR}.\end{array}$$

It should be noted that a simplified approach was used that does not account for the penetration of solar radiation below the mixed layer. Considering the Jerlov water type IA for relatively clear seawater, 77 % of the solar radiation is expected to be absorbed within the upper 10 m of the ocean based on the solar radiation attenuation equation in Paulson and Simpson (1977). Taking into account the MLD values obtained from the utilized reanalysis dataset, we have considered that the followed approach does not significantly affect our conclusions.

To compute heat flux anomalies, first a daily climatology was computed following the same methodology as with SST for the period 1993–2022. Daily heat flux anomalies were then constructed relative to the mean climatological value over the phase duration (e.g., as in Fewings and Brown, 2019). The time integral of Q^′ divided by the product of the constant values ρ (seawater density), c_p (specific heat capacity), and h (mixed layer thickness) represents the part of the ${\mathrm{DSST}}_{\mathrm{obs}}^{\prime }$ explained by Q^′ during this phase (DSST_Q). The second term of the right-hand side in Eq. (1) stands for any contribution to ${\mathrm{DSST}}_{\mathrm{obs}}^{\prime }$ from other mechanisms affecting SST: horizontal advection, vertical mixing processes, horizontal diffusion of heat flux, and radiative heat loss below the mixed layer. Therefore, the role of air–sea heat flux in this study is assessed in relation to a single residual term (i.e., the non-heat-flux terms merged into a single one), representing the cumulative effect of all other (oceanic) factors influencing the SST tendency during an MHW phase.

The contribution of Q^′ in driving an MHW onset or decline phase N is then quantified through the following proportion of change:

$$\begin{array}{}\text{(3)}& P\left(N\right)={\displaystyle \frac{{\mathrm{DSST}}_Q\left(N\right)}{{\mathrm{DSST}}_{\mathrm{obs}}^{\prime }\left(N\right)}}.\end{array}$$

Therefore, a positive heat flux contribution value during an MHW phase indicates a favoring role of heat flux in the corresponding SST evolution, i.e., a warming (cooling) effect of heat flux during onset (decline). Analogously, a negative contribution value during either an onset or a decline phase indicates that heat flux opposes the corresponding SST tendency. Finally, to examine the evolution of MLD during MHWs, time series of cumulative MLD anomalies (MLDA) were constructed for onset and decline phases. These time series were computed by adding, at each MHW day, the daily anomaly of MLD of the previous day in order to account for longer timescales associated with the mixed layer evolution (as in Schlegel et al., 2021). To explore the correlation between MLDA and SST anomalies (SSTA), Pearson correlation coefficients were computed for onset and decline phases separately.

3 Results

3.1 MHW detection

Properties of MHWs exhibit high variability throughout the Mediterranean Sea (Fig. 1a–c). The northwestern part of the basin, along with the northern Adriatic and northern Aegean seas, experienced on average the highest event intensity over the period 1993–2022, exceeding 2.5 and 2 °C, respectively (Fig. 1b). Events tend to last longer in the eastern part of the basin (Aegean and Levantine seas) and in the central-western Mediterranean Sea, while the shortest durations are mostly found in the Ionian and Alboran seas (Fig. 1c). The mean event frequency over the study period shows some similarities to the mean intensity spatial distribution, suggesting that the most (least) intense and most (least) frequent MHWs are encountered in the northernmost (southernmost) flanks of the Mediterranean Sea (Fig. 1a, b). Results generally agree with the literature on MHWs in the Mediterranean Sea during the recent decades (Darmaraki et al., 2019a; Ibrahim et al., 2021; Juza et al., 2022; Dayan et al., 2023), despite some differences in event detection methods (e.g., choice for percentile-based threshold, accounting for MHW spatial extent) or study periods (choice of climatological period, period for event detection).

Figure 1 (a–c) Mean annual frequency, intensity, and duration of MHWs for the period 1993–2022 (a, b, c, respectively). (d–f) Linear trends of annual values of frequency, intensity, and duration (d, e, f, respectively). Black dots superimposed on trend fields correspond to statistically nonsignificant trends (Mann–Kendall test, 95 % confidence level).

MHW frequency has been increasing over the past 30 years, with a trend of 1 ± 0.6 events per decade for the entire basin (Fig. 1d). Also, MHW duration has increased, particularly in the eastern basin, with trend values locally exceeding 10 d per decade (1.7 ± 1.5 d per decade for the entire basin) (Fig. 1f). Notably, intensity has not increased over the entire basin during the study period (0.08 ± 0.2 °C per decade for the basin, not significant at the 95 % confidence level) (Fig. 1e). The northernmost regions, which are characterized by higher MHW intensity, present small decreasing MHW intensity trends (though not statistically significant in most cases), which is in agreement with Dayan et al. (2023) and Ibrahim et al. (2021) for the eastern basin.

3.2 The role of heat flux during MHW onset and decline

Surface heat flux contributes to the observed surface warming during the onset phase of the majority of the events (92 %) detected in the Mediterranean Sea within 1993–2022 (Fig. 2a). In particular, in 44 % of the events, heat flux plays a major role during the development of MHWs in terms of explaining more than half of the observed change in SST anomaly (Fig. 2b). During the rest of the events (8 %), heat flux opposes the surface warming over the onset phases; therefore, other mechanisms compensate for their cooling effect and are responsible for the observed warming in such cases.

Figure 2 (a) Contribution of heat flux to the observed change in SSTA during MHW onset and decline phases for the period 1993–2022. (b) Percentage of events primarily driven by heat flux (i.e., with more than half of the observed warming or cooling attributed to heat flux during onset or decline) for the entire Mediterranean Sea and subregions mapped in (c). (d–e) Box plots for the contribution of heat flux during onset (red) and decline (blue) for the entire Mediterranean Sea and the subregions. (f–g) Box plots for heat flux contribution per season for onset (red) and decline (blue). (h–i) Box plots per season for the contribution of heat flux components during onset (red) and decline (blue). From left to right: latent heat flux, sensible heat flux, net shortwave radiation, net longwave radiation. Note: a positive heat flux contribution during onset (decline) means that heat flux will warm (cool) the sea surface; boxes in box plots define the interquartile range from the 25th up to the 75th percentile, and whisker bars correspond to values falling within 1.5 times the interquartile range.

Air–sea heat exchange contributes to surface cooling during decline phases in a much smaller percentage of MHWs (58 %), while it is the primary driver in only 17 % of the decline phases (Fig. 2a, b). In other words, almost half (42 %) of the observed MHWs in the basin decay under non-favorable heat flux conditions (i.e., heat flux opposing the SST decrease), while oceanic processes are the dominant driver of most (83 %) MHW declines.

Heat flux exchange explains the MHW onset in all examined subregions to a great extent (Fig. 2d). The Adriatic and Aegean seas stand out with most of the observed heating during onset attributed to heat flux in more than half of the identified events (Fig. 2b). In six out of the seven subregions, the percentage of events primarily driven by heat flux during their onset ranges 39 %–53 %. The Alboran Sea presents a much lower percentage of events attributed to heat flux (22 %), potentially suggesting an enhanced role of advective processes in regulating SST in the area, likely associated with its proximity to the Strait of Gibraltar. In contrast, the corresponding percentages for the decline phases show less spatial variability across all subregions, with only 15 %–20 % of the examined declines found to be mainly driven by heat flux exchange (Fig. 2b).

Seasonal analysis was performed considering winter to be the period from December to February, while events spanning different seasons were assigned to the season when their intensity maximizes. Most of the events were found to occur during summer, followed by spring, autumn, and then winter. As expected, the magnitude of SST anomalies during MHWs varies among the seasons (e.g., Thoral et al., 2022), but the contribution of heat flux to their formation is consistently weaker during decline compared to onset phases throughout the year (Fig. 2f, g).

Results suggest a greater contribution of heat flux exchange to the MHW evolution within warmer seasons (from an ocean perspective, i.e., summer and autumn) (Fig. 2f, g). Autumn shows the highest percentage of events driven by Q_net in both onset and decline. LH flux mainly regulates the contribution of Q_net during both phases and throughout the year, followed by SWR and then SH, while LWR exhibits less clear behavior in all seasons (Fig. 2h, i). With the exception of SH, the relative contribution of each heat flux component presents higher variability during decline compared to onset phases across all seasons. Specifically, approximately an equal number of decline cases are associated with positive and negative contributions of the heat flux components to the observed cooling, indicating a less predictable role of Q_net during decline compared to onset periods throughout the year (Fig. 2i). The lowest percentage of heat-flux-driven events is observed during the winter onset and spring decline phases. Particularly for spring, heat flux during most declines acts against the observed surface cooling (mainly through suppressed LH losses), suggesting that MHW dissipation in spring is commonly driven by oceanic factors.

Results show that the air–sea interaction, with a dominant role of LH flux, plays a major role in the development of nearly half (44 %) of the MHWs in the Mediterranean Sea. This finding suggests that oceanic processes play a key role during 56 % of the onset cases. A further weakened role of heat flux is found during decline periods (being the major contributor in only 17 % of declines), indicating that MHW decay is also primarily driven by oceanic processes.

3.3 Links with mixed layer depth and MHW characteristics

To gain insight into subsurface conditions during MHWs, we examined how MLD evolves in relation to SST. Negative (positive) correlation between MLDA and SSTA that was found for a large number of events during onset (decline) implies mixed layer shoaling over the entire MHW duration in these cases (Fig. 3a – top). While an MHW event develops, a reduction in MLD is commonly expected, as surface warming may strengthen the stratification of the water column (D'Ortenzio and Prieur, 2012). Given that heat flux is found to contribute to the warming phases in most cases, the concurrent mixed layer shoaling found during most onsets may be interpreted as an effect of the warming driven by the atmosphere. Nevertheless, a thinner-than-usual mixed layer may also exist before the event occurrence and act as a pre-conditioning factor (e.g., Lee et al., 2023). The greatest MLDA–SSTA correlation is found during spring and summer events, while no significant correlation is observed in winter (not shown), as the deeper mixed layer during colder months is expected to be less responsive to surface SST variations (D'Ortenzio and Prieur, 2012). High positive MLDA–SSTA correlation observed during declines implies that the mixed layer in such cases continues to shoal while SST decreases (Fig. 3a), similar to the mixed layer analysis by Schlegel et al. (2021) for the northwest Atlantic. This finding further supports the fact that oceanic processes play a significant role during MHW decline periods, as also indicated by the weaker contribution of heat flux during declines. It also suggests that the surface cooling occurring right after the peak intensity day is likely not due to mixing in the vertical (further discussed below).

Figure 3 (a) Heatmap relating the net heat flux contribution to the Pearson correlation coefficient (CC) between time series of SST anomalies and time series of mixed layer depth cumulative anomalies computed for each MHW onset and decline phase (upper and lower heatmap, respectively) for the period 1993–2022. Colors correspond to the number of events falling in each bin. (b) Box plots for the distribution of CC values for onset and decline phases lasting less than 5 d, between 5–10 d, and more than 10 d (short, medium, and long duration phases, respectively). (c) Percentage change (%) in MLD over a fixed 7 d period after the MHW end day compared to MLD during the onset period. Note: a positive heat flux contribution during onset (decline) means that heat flux will warm (cool) the sea surface. Negative (positive) CC values close to −1 (+1) during onset (decline) correspond to a reduction of the mixed layer depth while SST increases (decreases). Boxes in box plots define the interquartile range from the 25th up to the 75th percentile, and whisker bars correspond to values falling within 1.5 times the interquartile range.

To gain a better understanding of the relationship between subsurface conditions and the contribution of heat flux exchange at the air–sea interface during the evolution of MHWs, this information is integrated into Fig. 3a. Whereas during onset, MHWs are largely driven by heat flux exchange and most of them are accompanied by mixed layer shoaling, there are onset cases where MLDA is strongly positively correlated with SSTA, indicating that surface warming evolves while the mixed layer deepens (Fig. 3a – top). Regarding decline periods, almost equal numbers of cases show positive and negative heat flux contributions (Fig. 3a – bottom). While a decrease in MLD is evident in most decline phases (as indicated by correlation coefficient close to 1), a significant MLDA–SSTA correlation is absent in a considerable number of decline phases. Additionally, there are cases during declines where a high negative correlation between MLDA and SSTA is observed, indicating mixed layer deepening while SST decreases. Such cases are encountered when heat flux contributes to surface cooling (Fig. 3a – bottom), suggesting that vertical mixing (most probably wind-induced) works in the same direction.

MHW onset and decline phases were also examined in relation to their duration. They were grouped into short, medium, and long duration, considering phases lasting less than 5 d, from 5 to 10 d, and more than 10 d, respectively (Fig. 4a). Shorter onset or decline periods, being the most prevalent category, slightly overshadow the contribution of heat flux during longer-lasting MHW phases (Fig. 4a–c vs. Fig. 2a). In particular, compared to medium followed by long durations, short durations (both for onset and decline) tend to exhibit a smaller contribution of heat flux exchanges triggering the SST evolution (Fig. 4a–c). For long-duration onset phases, the contribution of heat flux surpasses the contribution of oceanic factors, as the former explains more than half of the surface warming in more than half of the events (Fig. 4c). Similarly, it is mostly during shorter declines that heat flux does not present a systematically positive or negative contribution to the SST decrease, while a shift towards higher contributions of heat flux is found for longer phases (Fig. 4a–c). Importantly, the relationship between MLDA and SSTA becomes less clear during onset and decline phases of medium and, in turn, long duration, as indicated by weaker correlation (Fig. 3b). In the case of decline periods, mixed layer shoaling is more frequently observed during short-lasting declines, suggesting a continuation of mixed layer shoaling already present during the onset in such cases.

Figure 4 (a–c) Contribution of heat flux to the observed change in SSTA during MHW onset and decline phases for the period 1993–2022 for phases lasting less than 5 d (a), between 5–10 d (b), and more than 10 d (c) (short, medium, and long duration phases, respectively). (d–f) Same as left but for different MHW severity categories: moderate (d), strong (e), and severe/extreme (f) based on the categorization scheme by Hobday et al. (2018).

These findings suggest that the definition followed for an MHW evolution phase needs to be cautiously taken into account when interpreting results, especially when surface and subsurface conditions are examined for the same phases. This is due to the time needed for a surface warming signal to penetrate below as well as due to the longer timescales associated with processes at deeper layers. Notably, on the end day of a surface MHW, which is the end day of the decline phase (as defined in this study), SST is still above the event detection threshold and will fall below this value by the following day. Therefore, events further dissipate after the decline phase ends. Given such considerations and depending on the specific purpose of a study, different approaches for defining a decline period may be followed. For instance, Darmaraki (2019) considered the entire MHW duration to be the event development phase and treated the following period as decline, and Schlegel et al. (2021) used the same definitions as the present study, while Marin et al. (2022) considered a dynamic criterion for defining the decay period, including a larger portion of the temperature change during the event dissipation after the MHW end day.

On these grounds, we additionally examine how MLD evolves after the MHW end day to shed light on what follows the progressively increasing stratification we observed during MHWs. To this aim, we compute for each event the percentage change in the mean MLD between the onset period and a fixed 7 d period following the MHW end day (Fig. 3c). Results reveal the existence of a deeper mixed layer after the decline of most events (83.5 %) compared to the onset period. In these cases, the magnitude of change in MLD is also significantly greater, occasionally exceeding 100 %, compared to instances when the mixed layer is thinner during this period than during the onset (changing the predetermined length of the examined post-decline period does not significantly alter these conclusions). This analysis shows that a significant mixed layer deepening most likely occurs after the end of MHW decline phases in the Mediterranean Sea, suggesting that vertical mixing eventually contributes to the heat dissipation.

Considering the continuation of mixed layer shoaling found mostly during shorter declines, results suggest that the oceanic factor regulating the SST decrease after the peak intensity day of most events is probably heat advection, while vertical mixing in these cases becomes important over the following days. In line with our findings, Marin et al. (2022) found that the principal driver of the upper-ocean temperature changes during MHW onset and decline phases in midlatitudes is horizontal heat advection. However, their results cannot be directly intercompared to ours, as they use modeled depth-integrated upper-ocean temperatures (thus not detecting surface events), they focus only on extreme cases, and they employ a different definition of the onset and decay period.

Finally, results were examined in relation to MHW intensity and severity to account for the perspective of MHW extremity in terms of absolute and normalized SST anomalies, respectively. Events were split based on their severity category following the categorization scheme proposed by Hobday et al. (2018) (definitions included therein). Results suggest a smaller contribution of heat flux in the evolution of more severe events (Fig. 4d–f). Half of the moderate events were found to be mainly driven by heat flux, while this percentage decreases for strong, severe, and extreme events, as illustrated in Fig. 4d–f. This may be associated with the enhanced role of heat flux during longer events, as was previously discussed, and the fact that events of higher severity categories tend to present a shorter duration. Similar results were found when examining the contribution of air–sea heat flux in relation to the mean intensity of MHWs (not included). Nevertheless, further investigations are needed to unravel how our methodological choices affect these findings.

4 Conclusions

This study investigates the role of air–sea heat exchange during MHWs in the Mediterranean Sea using satellite and reanalysis data within 1993–2022. An ocean mixed layer heat budget analysis is performed to derive the change in SST attributed to the net surface heat budget during onset and decline phases.

Air–sea heat flux is the primary driver in 44 % of the onset and only in 17 % of the decline phases in the basin. Heat flux acts in favor of the development of most MHWs across all seasons and especially during warmer months and onset phases. Moreover, its contribution is greater in the Adriatic and Aegean sub-basins, where it becomes the major driver of most onset phases. Among the heat flux components, LH emerges as the most significant contributor to SST anomalies, in line with prior studies (Sen Gupta et al., 2020; Oliver et al., 2021; Schlegel et al., 2021; Marin et al., 2022). Short-lasting onset or decline phases (shorter than 5 d) tend to experience a smaller contribution of heat flux in forming the SST evolution compared to longer phases (lasting 5–10 or more than 10 d). Furthermore, there is an inverse relationship between MHW severity and the contribution of heat flux. As fewer than half of the events are primarily driven by air–sea heat flux, our results imply a key role of oceanic processes in regulating SST during most events and particularly during decline periods.

Examining the mixed layer during MHWs revealed a progressively decreasing MLD over the entire event duration, particularly for shorter-lasting events. In turn, a significant mixed layer deepening was found to occur after the end of the decline period of most events in the basin. In cases of smaller contribution of heat flux (e.g., during shorter compared to longer events or during declines compared to onsets), a stronger correlation between SSTA and MLDA is also found, further supporting the key role of oceanic processes in such cases. These findings suggest that the surface cooling occurring right after the peak intensity day is likely not associated with vertical mixing. Moreover, this potentially suggests that horizontal advection is the oceanic factor playing the most significant role by the time the decline period begins, especially for shorter declines. This hypothesis finds support in the results of Marin et al. (2022), who highlighted the role of horizontal heat advection during the MHW evolution, and Schlegel et al. (2021), who suggested that advection and mixing should drive the MHW decline based on similar indications.

Nevertheless, the authors suggest taking into account potential limitations associated with the definition of MHW phases followed within a study, especially while examining concurrent subsurface conditions. These considerations concern the longer timescales below the sea surface, as well as the complexity associated with long-lasting events which are not expected to be adequately described by two single phases. Specifically, subsequent warming and cooling periods may occur within a long onset or decline phase, complicating their representation by the definition employed in this study for onset and decline. These concerns are part of the broader discussion on challenges related to the lack of a standardized framework for analyzing MHW drivers, as many other methodological approaches (e.g., the choice of integration depths in heat budget analyses) may potentially influence findings on MHW drivers. For this reason, alongside recommending the use of definitions and methods aligned with the specific contexts of individual studies, we highlight the fact that clearly articulating the employed methods within a study is vital for both the precise interpretation of results and meaningful comparisons across different studies on MHW drivers. Despite the aforementioned caveats, this study provides useful insights into the role of surface heat flux and mixed layer dynamics during MHWs in the Mediterranean Sea. Importantly, following the methodology applied by Schlegel et al. (2021) allowed for a fair intercomparison of results among two substantially different regions (Mediterranean Sea vs. northwest Atlantic). The results are surprisingly similar, which is largely attributed to the same methodology employed to assess drivers of SST-defined events and especially the same definition followed for onset and decline periods. The striking similarities in the results for the two regions may also imply inherent characteristics of the interplay of air–sea interaction and oceanic processes during anomalous SST fluctuations over similar latitudinal zones, which is also observed to a certain extent in Marin et al. (2022).

This study highlights the need to consider subsurface information in MHW studies to gain insight into ocean internal dynamics throughout their evolution. In this context, combining observations and ocean reanalysis systems, such as in Dayan et al. (2023), and considering MHW evolution periods aligned with the objectives and specific characteristics of a study appear to constitute a promising direction towards understanding physical drivers as well as improving monitoring and therefore enabling early warning of MHWs.