Oceanic mesoscale eddies as crucial drivers of global marine heatwaves

Marine heatwaves (MHWs) are prolonged extreme warm water events in the ocean, exerting devastating impacts on marine ecosystems. A comprehensive knowledge of physical processes controlling MHW life cycles is pivotal to improve MHW forecast capacity, yet it is still lacking. Here, we use a historical simulation from a global eddy-resolving climate model with improved representation of MHWs, and show that heat flux convergence by oceanic mesoscale eddies acts as a dominant driver of MHW life cycles over most parts of the global ocean. In particular, the mesoscale eddies make an important contribution to growth and decay of MHWs, whose characteristic spatial scale is comparable or even larger than that of mesoscale eddies. The effect of mesoscale eddies is spatially heterogeneous, becoming more dominant in the western boundary currents and their extensions, the Southern Ocean, as well as the eastern boundary upwelling systems. This study reveals the crucial role of mesoscale eddies in controlling the global MHW life cycles and highlights that using eddy-resolving ocean models is essential, albeit not necessarily fully sufficient, for accurate MHW forecasts.

Marine heatwaves (MHWs) are extreme warm water events in the ocean [1][2][3][4] . They cause severe environmental and socioeconomic impacts, including the loss of biodiversity, reduction of fishery catching rates, damage to aquaculture as well as changes in the behavior of species [5][6][7][8] . Satellite observations have revealed a significant increase in frequency, duration, and intensity of MHWs during the past several decades over most parts of the global ocean 9 , primarily due to the gradual sea surface warming caused by rising greenhouse gas emissions 10 . The growing threat of MHWs on marine ecosystems underscores the imperative for a comprehensive understanding of physical mechanisms responsible for the generation, maintenance, and decay of MHWs in the global ocean 2,[11][12][13] . Such understanding is a prerequisite for establishing a reliable forecast system for MHWs and developing sensible management strategies in a timely manner to alleviate the ecosystem stress and associated socioeconomic ramifications [14][15][16] .
Although there are plenty of case studies focusing on the physical drivers of some major MHW events, such as the well-known Blob 17 in 2014/15 and Blob 2.0 18 in 2019 over the northeastern Pacific, globalscale analysis is still limited. Holbrook et al. 12 established the first global view of the MHW drivers based on the correlation between the observed MHW occurrence and a variety of climate modes. Furthermore, ref. 19 investigated drivers of the most extreme MHWs and found that a large fraction of MHWs in the subtropical ocean coincide with persistent atmospheric high-pressure systems and weakened surface winds. Despite the important role of air-sea interactions in the MHW life cycles implied by these observational studies, a heat budget analysis based on a global ocean-only simulation 13 reveals that the global MHWs are primarily generated by heat flux convergence of oceanic flows, whereas the sea surface heat flux is the main driver of MHW decay. However, the absence of air-sea coupling in the ocean-only simulation causes biases in the simulated sea surface temperature (SST) variability [20][21][22] that may further propagate into the simulated MHWs.
Oceanic mesoscale eddies with a horizontal scale from a few tens to several hundreds of kilometers, manifested in the form of fronts, filaments, and coherent vortices, are the most prominent feature in the upper ocean 23,24 . They account for 70% of oceanic kinetic energy 25,26 and contribute importantly to the SST variability via their induced heat flux convergence [27][28][29] . Yet the effects of mesoscale eddies on the MHW life cycles in the global ocean remain unexplored and largely overlooked. In this study, we use a historical simulation from an eddyresolving global coupled climate model (CGCM) 30 (See "CESM-H" in Methods) to evaluate the role of mesoscale eddies in the global MHW life cycles. As will be demonstrated below, the heat flux convergence by oceanic mesoscale eddies acts as a crucial driver of MHW growth and decay in the global ocean.

Simulated MHWs in the CESM-H
Performance of the CESM-H in simulating the MHWs defined based on SST (see "Definition of MHWs" in Methods) is evaluated against satellite observations (see "Observational products" in Methods) and compared with an ensemble of state-of-the-art coarse-resolution (~1°) CGCMs in the Coupled Model Intercomparison Project Phase 6 (CMIP6) 31 (Supplementary Table 1 and Fig. 1). The CESM-H reproduces the spatial variability of MHW intensity in the observation reasonably well (Fig. 1a, d), with the correlation coefficient between the observed and simulated spatial patterns of MHW intensity reaching 0.81. In contrast, the correlation coefficient decreases to 0.68 between the observation and CMIP6 ensemble mean, mainly due to the absence of strong MHWs in the southern hemisphere western boundary currents and their extensions (WBCEs) and the Southern Ocean in the CMIP6 CGCMs (Fig. 1g).
The CMIP6 ensemble mean underestimates the frequency but overestimates the duration of MHWs (Fig. 1k, l), which has already been noted by the existing literature 9,32 . These biases are evidently alleviated in the CESM-H, especially for the duration of MHWs. Specifically, the globally averaged MHW frequency (duration) in the CMIP6 ensemble mean biases low (high) by 8% (94%) compared to the observational counterpart, whereas this bias is reduced to 3% (59%) in the CESM-H. According to the above comparisons, we conclude that the CESM-H provides a qualitatively reliable simulation of MHWs in the global ocean, outperforming the CMIP6 ensemble mean. This lends support to the fidelity of the CESM-H in simulating the physical drivers of MHW life cycles.

Physical divers of MHW life cycles in the global ocean
A heat budget analysis in the upper 50 m water column based on the CESM-H's diagnostic output (see "Heat budget analysis" in Methods) is used to quantify contributions to temperature changes during MHW life cycles by different physical processes, including the net sea surface heat flux (NHF), heat flux convergence by mean flows (HFC-M) and mesoscale eddies (HFC-E), and subgrid-scale mixing (MIX). To be consistent with the heat budget analysis, the MHWs are redefined based on the vertical mean temperature in the upper 50 m (denoted as T h i) rather than SST (see "Definition of MHWs" in Methods). As MHWs in this study are defined based on a seasonally varying threshold 1 , it is the anomaly of T h i relative to its climatological mean seasonal cycle (denoted as hT a i henceforth) that is related to MHWs. Accordingly, all the terms in the heat budget are subtracted by their corresponding climatological mean seasonal cycles to quantify their induced changes of hT a i during the MHW life cycles (see "Heat budget analysis" in Methods).
We perform a heat budget analysis over a fixed depth range rather than over the mixed layer, because the latter is difficult to close based on the available model output due to spatio-temporal variations of mixed layer depth (MLD). Moreover, although the heat budget analysis over the mixed layer is dynamically more suitable for analyzing the SST variability than that over a fixed depth range like 0-50 m, it is not necessarily so from biological concerns, because the depth of epipelagic zone, home to a massive number of organisms, may differ from the MLD in many parts of the global ocean 33,34 .
The MHW statistics defined based on SST and T h i generally agree with each other except in the tropical eastern Pacific where the MHWs defined based on T h i are shorter and more frequent than those based on SST ( Supplementary Fig. 1). Such differences may reflect the different temporal variabilities between SST and T h i due to the shallow mixed layer (shallower than 50 m) in this region but should not necessarily be interpreted as deficiencies in defining MHWs based on T h i, as the epipelagic zone here is deeper than the mixed layer 33,34 . In fact, the common use of SST to define MHWs in the existing literature 1,19 is mainly due to its availability in the observation and could be insufficient to measure the thermal stress on marine ecosystems 35,36 .  Fig. 2) averaged at each grid point, respectively. The mesoscale eddies play a crucial role in driving the MHW life cycles over most parts of the global ocean. In particular, the HFC-E accounts for 81% (74%) of hT a i increase (decrease) during the growing (decaying) phase of the MHWs averaged over the global seaice-free ocean (Fig. 3a). The contribution of the HFC-E to the hT a i change during the MHW life cycles varies in space (Fig. 2a, e). It is stronger in the WBCEs and Southern Ocean, consistent with the more energetic mesoscale eddies in these regions 23 .
During the growing phase of MHWs, the HFC-M induced hT a i increase is much smaller than that induced by the HFC-E except in the central-to-eastern equatorial Pacific, where the HFC-M plays a dominant role in driving the MHW growth (Fig. 2b). During the decaying phase of MHWs, the HFC-M still acts to increase hT a i in many parts of the global ocean, especially in the WBCEs (Fig. 2f). The NHF causes hT a i to decrease during the growing phase of MHWs in the WBCEs and Southern Ocean, whereas the opposite is true elsewhere (Fig. 2c). During the decaying phase, the NHF leads to a universal hT a i decrease in the global sea-ice-free ocean (Fig. 2g). The NHF-induced hT a i decrease reflects the sea surface heat flux feedback 37,38 , i.e., the generation of sea surface heat flux anomaly by hT a i that, in turn, damps hT a i itself. As to the MIX that is primarily attributed to the vertical mixing ( Supplementary Fig. 3), it acts to increase hT a i both in the growing and decaying phases of the MHWs in the WBCEs, Southern Ocean and central-to-eastern equatorial Pacific, whereas its contribution to hT a i change is close to zero elsewhere ( Fig. 2d and h).
The spatially heterogeneous effects of the different physical processes on the hT a i changes during the MHW life cycles suggest that the dominant physical drivers of the MHW growth and decay are region-dependent ( Supplementary Fig. 4). In the WBCEs with energetic mesoscale eddies, the HTC-E accounts for average for 97% (89%) of the hT a i increase (decrease) during the growing (decaying) phase of the MHWs (Fig. 3b), acting as the single dominant driver of the MHW life cycles. Similar is the case for the Southern Ocean, with the HFC-E contributing 88% to the increase of hT a i during the growing phase and 87% to the decrease of hT a i during the decaying phase, respectively (Fig. 3c). In contrast, 47% of hT a i increase during the growing phase in the central-to-eastern equatorial Pacific is attributed to the HFC-M (Fig. 3d), implying the association of MHW generation with the sea surface warming during El Niño events via the Bjerknes feedback 12,39 . Nevertheless, the HFC-E plays an important role in driving the MHW growth in this region, accounting for 45% of the hT a i increase during the growing phase. As to the decrease of hT a i during the decaying phase in the central-to-eastern equatorial Pacific, it is mainly ascribed to the HFC-E (70%), consistent with the damping of El Niño events by the HFC-E 40,41 . In the biologically productive eastern boundary upwelling systems 42,43 , the HFC-E contributes 91% to the increase of hT a i during the growing phase, while the decrease of hT a i during the decaying phase is contributed primarily by the HFC-E (74%) and secondarily by the NHF (20%) (Fig. 3e). In the subtropical gyre interior, HFC-E still plays a dominant role (70%) in driving the hT a i increase during the growing phase (Fig. 3f). During the decaying phase, the contribution to the hT a i decrease by the NHF (44%) becomes comparable to that of the HFC-E (60%).
Finally, we re-perform the heat budget analysis by varying the lower bound of the water column from 20 m to 200 m that covers the range of euphotic zone depth over the global ocean 33,34 (Supplementary Figs. 5-7). As expected, contributions of the NHF and MIX to the hT a i change during the MHW life cycles become more important for the shallower water column. Nevertheless, for the range of lower bound considered here, the HFC-E always plays a dominant role, lending strong support to the crucial role of oceanic mesoscale eddies in driving the global MHW life cycles.

Role of mesoscale eddies in driving the life cycles of MHWs with large spatial scales
The HFC-E induced hT a i change occurs primarily at the oceanic mesoscales ( Supplementary Figs. 8, 9). Correspondingly, the HFC-E should be most efficient in driving the growth and decay of the MHWs whose characteristic spatial scale is comparable to that of mesoscale eddies. However, the correlation between the velocity and temperature anomalies induced by mesoscale eddies can generate HFC-E that is coherent at spatial scales larger than that of mesoscale eddies ( Supplementary Figs. 8b, 9b; See "Heat budget analysis" in Methods). Such correlation can arise from the baroclinic instability 44,45 , frontogenesis/frontolysis 46,47 , and turbulent thermal wind 48,49 . Therefore, the HFC-E may also play a role in driving the life cycles of MHWs with a characteristic spatial scale larger than that of mesoscale eddies. To demonstrate this point, we redefine MHWs based on the large-scale T h i (denoted as T ) that filters out the mesoscale perturbations and refer to them as the large-scale MHWs henceforth. Then we quantify the effects of different physical processes on the changes of T a during the life cycles of large-scale MHWs (See "Heat budget analysis" in Methods).
Contributions of the HFC-E to the T a changes during the growing and decaying phases of large-scale MHWs are systematically weaker than its counterparts for the all-scale MHWs ( Supplementary  Fig. 10), whereas the contributions of the HFC-M, NHF, and MIX are affected to a less extent by the spatially low-pass filtering (Figs. 2, 4). Correspondingly, the dominant physical drivers of the large-scale MHW life cycles in the global sea-ice-free ocean are taken over by the HFC-M and NHF for the growing phase and by the NHF alone for the decaying phase (Fig. 4). Nevertheless, the HFC-E still plays an important role in the regional large-scale MHW life cycles, contributing primarily to the T a changes during the growing and decaying phases in the WBCEs, the T a increase during the growing phase in the eastern boundary upwelling systems, and the T a decrease during the decaying phase in the central-to-eastern equatorial Pacific ( Fig. 4 and Supplementary Fig. 10).

Discussion
Our results reveal the crucial role of mesoscale eddies in driving the growth and decay of global MHWs, which is largely overlooked in the existing literature. As these mesoscale eddies are not resolved by coarse-resolution CGCMs, it may partially account for the less frequent MHWs in the CMIP6 CGCMs than the observation (Fig. 1b, h, k). In particular, the MHW frequency in the observations and CESM-H is locally increased in the WBCEs and Southern Ocean with active mesoscale eddies, whereas such increase is largely absent in the CMIP6 CGCMs (Fig. 1b, e, h). Moreover, as mesoscale eddies have a relatively shorter time scale than that of mean flows, the MHWs generated by the HFC-E are expected to last for a shorter period than those generated by the HFC-M. This may explain the, on average longer duration of MHWs in the CMIP6 CGCMs than the CESM-H and observations (Fig. 1c, f, i, l). Nevertheless, we remark that the overly small MHW number and overly long MHW duration still persist in the CESM-H albeit alleviated. Such biases in the CESM-H may partially result from the unresolved heat flux convergence by submesoscale eddies 50 and deficiencies of vertical mixing parameterization 51 that do not account for the effects of surface waves 52,53 and Langmuir turbulence 54 .
The mesoscale eddies, also known as the geostrophic turbulence 55 , are essentially chaotic with limited predictability. This imposes a strong restriction on the forecast capacity of MHWs driven by mesoscale eddies and is consistent with the recent finding 16 that the forecast skills of MHWs by numerical models are evidently lower in the eddy-rich regions than elsewhere 12,16,56 . For these mesoscale eddydriven MHWs, it is more feasible to predict their statistics instead of individual characteristics. It has been well recognized that mesoscale eddies exhibit evident variabilities on multiple time scales regulated by changes in large-scale oceanic and atmospheric circulations [57][58][59][60] . How the variabilities of mesoscale eddies may affect the statistics of MHWs, remain unexplored but will be pivotal for proactive decisionmaking.

Observational products
The observational SST comes from the National Oceanic and Atmospheric Administration (NOAA) Optimum Interaction Sea Surface Temperature V2.0 high-resolution (OISST), which is derived from the advanced very high-resolution radiometer (AVHRR). The daily SST is provided on a 0.25°× 0.25°spatial grid. The data from Jan 1982-Dec 2021 are used for analysis.

Definition of MHWs
An MHW is defined as an event with at least five contiguous days of a given temperature index above its seasonally varying 90th percentile calculated over a baseline period 1   The geographic distributions of MHW statistics computed by using different temperature indices and baseline periods are qualitatively consistent with each other (Supplementary Fig. 1). As greenhouse warming is insignificant before the 1950s, such consistency implies that by now, the anthropogenic climate changes have not altered the physical processes underpinning the MHWs substantially, lending supports that the results derived for the period 1920-1934 are also representative of the present-day situation.
The start and end times of an MHW are typically defined as the time when the temperature index rises above and declines below its threshold (denoted as t s0 and t e0 hereinafter), respectively 1 (Supplementary Fig. 2). However, we remark that such defined start and end times may not be suitable for analyzing the physical drivers of MHW life cycles. For instance, the temperature index may have already increased by an evident amount before it exceeds the threshold and it is the physical process driving this increase that is mainly responsible for the MHW generation. Therefore, we define the start time t s of an MHW as the local minimum point just before t s0 and its end time t e as the local minimum point just after t e0 (Supplementary Fig. 2). The peaking time t p of an MHW is defined as the maximum point of the temperature index anomaly relative to its climatological mean seasonal cycle. The period between the start and peak times of an MHW is defined as its growing phase, while the period between the peak and end times is defined as its decaying phase.

Heat budget analysis
To evaluate the physical drivers of MHW life cycles, a heat budget in the upper 50 m is performed: where T is the temperature, u = ðu,v,wÞ is the three-dimensional velocity, ∇ is the three-dimensional gradient operator, Q NHF is the surface heat flux into the ocean minus the solar radiation penetrated out of the layer base at 50m,ρ 0 = 1026 kg Á m À3 is the reference seawater density, C p = 4000Jðkg Á KÞ À1 is the heat capacity of seawater, H = 50 m is the layer thickness, HMIX is the subgrid-scale horizontal mixing, VMIX is the subgrid-scale vertical mixing parameterized by a Kprofile parameterization 51 the overbar denotes the mean flow signals obtained by a 3 × 3 horizontal running mean, the prime denotes the mesoscale eddy field computed as the perturbation from the 3 × 3 horizontal running mean, and the angle brackets denote the vertical average in the upper 50-m layer. The term on the left-hand side of Eq.
(1) is the temperature tendency. The first À∇ Á u T À Á and second terms À∇ Á u 0 T' ð ÞÀ∇ Á uT' ð ÞÀ∇ Á u 0 T À Á on the right-hand side are the heat flux convergence by mean flows and mesoscale eddies (HFC-M and HFC-E), respectively. The third term is the temperature change caused by the net sea surface heat flux into the upper 50-m water column (NHF). All the terms in Eq. (1) can be explicitly computed based on the CESM-H's diagnostic output. As it is the anomaly of T h i(denoted as hT a i) relative to its climatological mean seasonal cycle that is related to MHWs by definition 1 , the climatological mean seasonal cycles are subtracted from all the terms in Eq. (1). Then individual terms in Eq. (1) are integrated over the growing or decaying phase of MHWs to quantify the contributions of different physical processes to the changes of hT a i. To analyze the physical drivers of MHW life cycles at a given grid point (Fig. 2), change of hT a i and its decomposition into components contributed by different processes are averaged over all the MHWs at that grid point. For the physical drivers of MHW life cycles in a given region (Fig. 3), change of hT a i and its contributions by different processes are weighted average over all the MHWs at the grid points within that region, with the grid area taken as the weight.
It should be noted that mesoscale eddies can cause temperature change with a spatial scale larger than that of mesoscale eddies via the correlation between u 0 and T 0 . This can be shown by horizontally averaging Eq. (1) using the 3 × 3 horizontal running mean: ∂T ∂t * +

Code availability
The iHESP version of CESM-H code is available at ZENODO via https:// doi.org/10.5281/zenodo.3637771. The code used to analyze these data and generate the results presented in the study can be obtained from https://github.com/cecbian/MHW. The Matlab2022b is used for plotting. The MATLAB code of MHW distinguish is obtained from https:// github.com/ZijieZhaoMMHW/m_mhw1.0 62 .