Frequent marine heatwaves hidden below the surface of the global ocean

Marine heatwaves are extreme warm water events that can cause devastating impacts on ecosystems and have complex socio-economic ramifications. Surface signals and drivers of marine heatwaves have been extensively investigated based on satellite observations, whereas their vertical structure in the global ocean remains unclear. In this study, we identify marine heatwave events in the epipelagic zone (0–200 m) using a four-dimensional spatio-temporal framework based on three ocean reanalysis datasets. We find that only about half of the marine heatwave events have continuous surface signals during their life cycles and around one-third always reside in the subsurface ocean without any imprint on sea surface temperature. The annual number of these subsurface marine heatwave events shows a significant increase in response to subsurface mean-state warming during the past three decades. Our findings reveal the limitation of identifying marine heatwaves solely based on the sea surface temperature and underscore the necessity of subsurface observations for monitoring marine heatwaves. Substantial numbers of marine heatwaves are hidden globally below the sea surface, according to analyses of ocean temperature data.

temperature does not necessarily co-vary with SST.A comprehensive knowledge of MHW statistics below the sea surface is thus pivotal but still lacking at a global scale.Several case studies have indicated various vertical structures of MHWs in different regions caused by their complicated dynamical drivers.Some MHW events are found to be associated with simultaneous warm anomalies across the epipelagic zone, with the amplitude of anomalies either peaking at the sea surface 23 or in the ocean interior [24][25][26][27] .In addition, downward migration of some MHW events causes the fading away of high SST anomalies followed by a subsurface temperature increase [28][29][30] (for example, the well-known 'Blob' in the northeastern Pacific).Last but not least, there are subsurface MHW events with no SST imprint at all during their life cycles 31 .
Learning the vertical structures of MHWs in the epipelagic zone is challenging not only because of the sparsity of subsurface observations but also due to the lack of an appropriate approach to characterizing the spatial structures of MHWs and their temporal evolution.Currently, MHW events are typically identified based on temperature time series Article https://doi.org/10.1038/s41561-023-01325-wfor the ECCO2 whose proportions of the surface, subsurface and mixed MHW events are 52.8%,27.7% and 19.5%, respectively.These quantitative discrepancies among the different reanalysis datasets reflect the differences in their model configurations, temporal coverages (Supplementary Table 2) and data assimilation methods.Nevertheless, all the three reanalysis datasets reveal that a large proportion of MHW events do not have continuous surface signals during their life cycles and the usage of SST is incapable of detecting their presence.

Geographic distribution of MHW events
There is evident spatial heterogeneity for the frequency of different types of MHW event (Fig. 1), which is consistently revealed by the three reanalysis datasets.For the surface MHW events, their frequency shows broad agreement in the geographic pattern with their counterpart derived from the SST in the temporal domain 19,20,34 (Supplementary Fig. 6).They occur more frequently in the western boundary current extensions (WBCEs) and Southern Ocean (SO; Fig. 1a-c), where there is large SST variability due to the advection by the energetic oceanic frontal jets and mesoscale eddies [35][36][37] .In contrast, the subsurface MHW events are mostly distributed in the subtropical gyre interior (SGI; Fig. 1d-f).The geographic distribution of the mixed MHW events is relatively more uniform and resembles a hybrid between those of the surface and subsurface MHW events, exhibiting an increased frequency in the WBCEs, SO and SGI (Fig. 1g-i).The majority of the mixed MHW events have no surface signals at the beginning of their life cycles, and a large fraction at the end (Supplementary Table 3).The former might be caused by the outcropping of warm anomalies via deepening of the surface mixed layer 38 or lifting of isopycnals, while the latter might result from the migration of warm anomalies from the sea surface to the ocean interior via subduction or downward isopycnal deflection 28,29 .

Properties of MHW events
Identifying MHW events in the four-dimensional domain allows characterizing both the magnitude and vertical structure of their time-mean horizontal area A(z) and intensity I(z) ('Properties of MHW events in the four-dimensional domain' in Methods).The A(z) averaged over the global surface MHW events peaks at the sea surface and decreases monotonically with depth to less than one-third of its surface value at 200 m (Fig. 2a).In contrast, the global mean A(z) for the subsurface MHW events is vertically reversed, peaking at 200 m and decreasing monotonically to zero at the sea surface (Fig. 2b).It thus suggests that these subsurface MHW events are likely to extend well below 200 m to the mesopelagic zone.The global mean A(z) for the mixed MHW events is relatively more uniform in the vertical, with its value varying by less than 50% in the upper 200 m (Fig. 2c).
The global mean I(z) for the surface MHW events varies a little in the upper 200 m (Fig. 2d), suggesting that the surface MHW events on average can cause a subsurface warm anomaly as strong as the surface warm anomaly.As to the subsurface MHW events, their global mean I(z) peaks around 100 m, attenuating to both sides (Fig. 2e).Similar is the case for the global mean I(z) of the mixed MHW events (Fig. 2f).We remark that the magnitudes and vertical structures of A(z) and I at a single location 18,32 .Such an identification method can be treated as a projection of MHW events actually occurring in the four-dimensional spatio-temporal domain to the temporal domain, with the information on their spatial structures lost during the projection.In particular, vertical displacement or deformation of MHW events can be properly represented in a four-dimensional domain but may cause misrepresentation of MHW events in the temporal domain.
In this Article, we extend the identification of MHW events from the temporal domain to the four-dimensional domain and apply this identification framework to three ocean reanalysis datasets ('Ocean reanalysis datasets' in Methods) to capture the vertical structures of MHW events in the global ocean.It is understood that there is not always one-to-one correspondence between the MHW events in the observation and reanalysis datasets.Therefore, our study focuses on the statistics of MHW events that are found to show qualitative consistency among the different reanalysis datasets.As will be shown in the next section, around one-third of MHW events in the global ocean have no surface signals at all throughout their life cycles.These subsurface MHW events have occurred increasingly frequently during the past three decades, primarily due to the subsurface mean-state warming.

MHW events in the four-dimensional domain
A MHW event in the four-dimensional domain is defined as a discrete prolonged anomalously warm water event whose spatial structure continuously evolves with time (Supplementary Fig. 1).This can be treated as an extension of the definition of MHW events in the temporal domain 18 .Identification of four-dimensional MHW events involves first constructing snapshots of MHWs with compact spatial structures and then tracking their temporal evolution ('Identifying MHW events in the four-dimensional domain' in Methods).
The identification approach is applied to the three ocean reanalysis datasets including the Global Ocean Physics Reanalysis (GLORYS), Hybrid Coordinate Ocean Model (HYCOM) and Estimating the Circulation and Climate of the Ocean, Phase II (ECCO2).Temperature data from different ocean reanalysis datasets are first horizontally smoothed using a 1° × 1° running mean and then interpolated onto a 1° × 1° × 10 m uniform grid over the upper 200 m covering the epipelagic zone 22 .It should be noted that the MHW events occupying less than 125 grid cells (corresponding to a volume from 1.6 × 10 4 km 3 at the equator to 0.8 × 10 4 km 3 at 60°) are excluded in the analysis due to the limitation of our identification method ('Identifying MHW events in the four-dimensional domain' in Methods).These small MHW events in space are likely to have less adverse biological effects than the larger ones in general.The fixed grid cell number criterion, compared to a fixed volume criterion, could avoid excessive discard of MHW events in higher-latitude regions where the oceanic processes including MHW events tend to be horizontally smaller in size 33 , partially due to the smaller Rossby deformation radius there.Nevertheless, replacing this fixed grid cell number criterion with a fixed volume criterion of 1.3 × 10 4 km 2 (the volume of 125 grid cells at 30°) does not have a major impact on results presented below (Supplementary Table 1 and Supplementary Figs.2-5).Finally, the identified MHW events are classified into different kinds according to their vertical structures.Specifically, a MHW event is classified as a surface MHW event if it has continuous surface signals during its life cycle, as a subsurface MHW event if it has no surface signals at all during its life cycle or as a mixed MHW event otherwise.
The long-term mean number of MHW events per year globally (60° S-60° N) differs among the three reanalysis datasets, being 309 for the GLORYS, 374 for the HYCOM and 219 for the ECCO2 (Table 1).For the GLORYS, the surface, subsurface and mixed MHW events account for 46.9%, 31.2% and 21.9% of the total MHW events, respectively.The HYCOM has a relatively lower proportion of the surface MHW events (32.6%) but higher proportions of the subsurface (42.4%) and mixed MHW events (25.0%), compared with the GLORYS.The opposite is true Article https://doi.org/10.1038/s41561-023-01325-w(z) derived from the three reanalysis datasets are qualitatively consistent, lending support to their reliability.
We next examine the differences in A(z) and I(z) among different regions including the WBCEs, SO, SGI and equatorial Pacific (EP).The regional mean A(z) for the three kinds of MHW events in the individual regions qualitatively resembles their global mean counterparts in terms of the vertical structure (Figs.2a-c and 3a-c), whereas the magnitude of the regional mean A(z) exhibits evident inter-region difference.b,c,e,f,h,i, Same as a,d,g, but for the HYCOM (b,e,h) and ECCO2 (c,f,i).A MHW event is allocated to some bin if its time-mean geometric centre during the life cycle is located within that bin.The horizontal area of a MHW event is first computed at each depth at each time and then temporally averaged over its life cycle to yield A(z).The intensity of an event, measured as the horizontal mean temperature anomaly relative to the climatological seasonal cycle, is first computed at each depth at each time and then temporally averaged over its life cycle to yield I(z).The pink, blue and purple lines correspond to the results derived from the GLORYS, HYCOM and ECCO2.

Article
https://doi.org/10.1038/s41561-023-01325-w All the three kinds of MHW events are considerably larger in the EP than in the other regions, which may reflect the central role of large-scale air-sea coupling in generating the temperature variation and MHW events in the EP 19,39 .In addition, the surface MHW events in the SGI have a larger size than those in the WBCEs and SO, especially near the sea surface.This is consistent with the recent finding 21 that a large fraction of the surface MHW events in the SGI are associated with the large-scale persistent atmospheric high-pressure systems.Similar to A(z), the vertical structures of the regional mean I(z) for the three kinds of MHW events in the individual regions are also qualitatively consistent with their global mean counterparts, except for the surface MHW events in the EP (Figs. 2d-f and 3d-f).Far from being vertically uniform in the upper 200 m, the regional mean I(z) for the surface MHW events in the EP peaks around 50 m and attenuates rapidly downwards.This might be partially due to the much shallower thermocline in the EP especially during the El Niño events than in the other regions 40 .As to the magnitude of the regional mean I(z), its inter-region difference is less evident compared to that of A(z) yet still noticeable.On one hand, the surface MHW events are most intense in the WBCEs, consistent with the prominent SST variability in these regions 35,36 .On the other hand, the EP has most intense subsurface and mixed MHW events.This might partially result from the strong vertical temperature gradients across the thermocline there so that the vertical displacement of isotherms can generate the large temperature anomaly 41 .

Response of MHW events to anthropogenic climate change
Ocean mean-state warming due to rising greenhouse gas emissions has been demonstrated to cause an increasing frequency of MHW events identified in the temporal domain based on the SST 15,17 .Here we examine the effects of the mean-state warming on the MHW events identified in the four-dimensional domain.It should be noted that although all three reanalysis datasets exhibit the mean-state warming throughout the upper 200 m, their warming rates differ quantitatively from the observed values (Fig. 4a and Supplementary Fig. 7), which may bias the long-term changes of the MHW events 42 .To eliminate such biases, MHW events are re-identified based on the synthetic data by replacing the long-term linear trend in the reanalysis datasets with that derived from the observations ('Synthetic temperature data' in Methods).
The annual numbers of the surface, subsurface and mixed MHW events in the global ocean all exhibit significant (P value < 0.05) positive trends during the past three decades (Fig. 4b-d).The slope of the linear trend is largest for the surface MHW events, varying from 46 to 60 annual events per decade among the three reanalysis datasets.The subsurface MHW events and mixed MHW events have weaker but still significant (P value < 0.05) trends in their annual numbers with the slope of the linear trend varying from 20 to 49 annual events per decade and 21 to 42 annual events per decade, respectively.All these trends of the annual MHW event numbers result primarily from the mean-state warming, with the changing temperature variability playing a minor or negligible role (Supplementary Fig. 8).As the mean-state warming is fastest at the sea surface and becomes slower with depth (Fig. 4a), this explains why the long-term change of the annual number is largest for the surface MHW events but smaller for the subsurface and mixed MHW events.Finally, it should be noted that the trends of the annual numbers of global MHW events in the individual synthetic reanalysis datasets differ quantitatively from each other, even though these datasets share the same mean-state warming derived from the observation (Fig. 4b-d).Such differences imply the differed sensitivity of the MHW events to the mean-state warming among the different reanalysis datasets, which might be partially related to their different high-order temperature statistics (for example, variance and skewness) 43 .

Further applications and future directions
In this study, by constructing a framework for identifying MHW events in the four-dimensional spatio-temporal domain, we discover great numbers of MHW events hidden below the sea surface.These events either reside in the subsurface layer sometime during their life cycles or never have a SST imprint from their emergence to disappearance.Even for the MHW events with continuous surface signals, they are not confined to the thin surface layer but can extend to the entire epipelagic zone and further downward with comparable magnitudes of their induced surface and subsurface warm anomalies.
Although the MHW identification in this study is confined to the upper 200 m, our identification framework can be applied to any depth range.Extending the lower bound to 400 m leads to an even larger fraction of subsurface MHW events in the global ocean (Supplementary Table 4).In addition, by discarding MHW events occupying less than 125 grid cells, we may not be able to capture the MHW events caused by mesoscale eddies.It is possible to identify smaller MHW events in space, by using a smaller grid cell number threshold associated with carefully tuning the parameters in the identification method.Our identification framework provides a flexible tool and can easily adapt to various applications.
The discovery of frequent subsurface MHW events suggest the deficiency of identifying MHW events solely based on SST and the desirability of subsurface observations.Currently, the capacity to identify the subsurface MHW events is strongly limited by the sparsity of long-term high-frequency observations below the sea surface.Most of observational three-dimensional gridded temperature data in the Fig. 3 | Spatial variability of MHW properties of different kinds of MHW events.a-c, Regional mean vertical profile of the time-mean horizontal area A(z) of the surface (a), subsurface (b) and mixed MHW events (c) in the EP (orange), SGI (light blue), WBCEs (navy blue) and SO (green).The vertical profiles in the different regions are normalized by their respective maximum values (coloured numbers) to facilitate illustration.d-f, Same as a-c but for the time-mean intensity I(z).A MHW event is allocated to some region if its time-mean geometric centre during the life cycle is located within that region.The regional mean is first computed for MHW events in the individual reanalysis datasets and then averaged over the three reanalysis datasets.g, Domain of different regions: EP (orange), SGI (light blue), WBCEs (navy blue) and SO (green).

Article
https://doi.org/10.1038/s41561-023-01325-wglobal ocean are available at a monthly timescale [44][45][46] , but insufficient to directly resolve MHW events.There have been some attempts to use these monthly temperature data as a proxy to estimate the statistics of the surface MHW events 20,47 .Although the results seem promising, more efforts are needed to improve the estimation methods.Another potential is to use satellite-measured sea surface height to infer the subsurface MHW events, as these events are typically associated with descending isopycnals 29,30 that may lead to rising sea surface height through the hydrostatic balance.Enhanced subsurface observation systems along with more advanced methods are imperative for monitoring and forecasting the subsurface MHW events, which will be a prerequisite of making sensible management strategies to alleviate MHW-caused ecosystem stress and associated socio-economic ramifications.

Online content
Any methods, additional references, Nature Portfolio reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41561-023-01325-w.

Ocean reanalysis datasets
Ocean temperature data in three common reanalysis datasets are used in this study to identify MHW events.The first is the daily 1/12° Global Ocean Physics Reanalysis products (GLORYS) 48 covering the period 1993-2020.The second is the Hybrid Coordinate Ocean Model (HYCOM) 49 providing three-hour snapshots on a 0.08° × 0.08° grid from 1993 to 2012.The third is the Jet Propulsion Laboratory of National Aeronautics and Space Administration (NASA/JPL) Estimating the Circulation and Climate of the Ocean (ECCO), Phase II (ECCO2) project 50 archived as a three-day mean field on a 0.25° × 0.25° grid from 1992 to 2020.
Temperature data from different ocean reanalysis datasets are first horizontally smoothed using a 1° × 1° running mean and then interpolated onto a 1° × 1° × 10 m uniform grid over the upper 200 m covering the epipelagic zone 22 .Temporal resolutions of the different datasets are uniformized to be daily either via interpolation or bin average.The data outside 60° N-60° S are discarded to remove the effects of sea ice.

Identifying MHW events in the four-dimensional domain
In this study, we propose to identify MHW events in the four-dimensional domain defined as a discrete prolonged anomalously warm water event whose spatial structure continuously evolves with time.'Discrete' means that a MHW event should have well-defined start and end times and be compact in the spatial domain.'Prolonged' means a sufficiently long duration, that is, at least five consecutive days 18 .'Anomalously warm' means that the temperature of a MHW event should be well above the baseline climatology.We require the spatial structure of a MHW event to 'continuously evolve with time' so that the snapshots at consecutive times can be tracked.The identification procedure is a generalization of the method proposed by Sun et al. 33 , consisting of constructing MHW snapshots in three-dimensional space and then tracking their temporal evolution.
Constructing snapshots of MHWs with a compact spatial structure at each time involves the following two steps 33 .The first step is to identify MHWs in the temporal domain based on the time series of daily temperature at each grid cell, that is, a discrete period of at least five consecutive days with the temperature above a predefined seasonally varying 90th percentile θ 18 .The value of θ is calculated over a baseline period chosen as the entire period of the ocean reanalysis datasets, that is, 1993-2020 for the GLORYS, 1993-2012 for the HYCOM and 1992-2020 for the ECCO2.However, the identification of MHWs in the temporal domain does not account for the spatial dependence of temperature and leads to spatially incoherent structures of MHWs.This deficiency is remedied in the second step by smoothing the binary variable (MHW or not) over the three-dimensional spatial domain.Here we use the Nearest-Neighbour (NN) method 51,52 for smoothing.Specifically, a grid cell is determined as a MHW if the majority of its nearest grid cells belong to MHWs according to the raw MHW snapshots obtained in the first step and not a MHW otherwise.Supplementary Fig. 9 shows an example of MHW snapshots before and after the NN smoothing.
To determine the neighbourhood, some sort of distance measure between two grid cells has to be defined.Although the spatial distance is the commonly used choice 53 and adopted by Sun et al. 33 , it is not appropriate for the problem at hand due to the strong anisotropy of the ocean 54 .Instead, we use the correlation-based distance 53 , that is, two grid cells are considered to have smaller distance if the correlation coefficient r between their time series of temperature anomalies relative to the climatological seasonal cycle is higher.It should be noted that the rationality of using the spatial distance is largely rooted from the First Law of Geography 55 that variables at spatially proximate locations tend to be more correlated to each other.The correlation-based distance is thus a more effective distance measure than the spatial distance.
The threshold of r (denoted as r c ) to define the neighbourhood is a tuning parameter affecting the performance of the NN smoothing.
An overly large r c may lead to insufficient neighbours for smoothing and is incapable of acquiring spatially coherent MHW events, while a too-small r c may cause excessive smoothing and distort the spatial structure of major MHW events evidently.As there are no pre-trained statistical methods to find the optimal threshold, we tentatively chose r c = 0.5 in this study based on visual inspection of hundreds of MHW snapshots.Sensitivity analysis suggests that the proportions of different kinds of MHW events vary only slightly for r c ranging from 0.4 to 0.6, and these variations are much smaller than the inter-dataset variations (Table 1 and Supplementary Table 5).Moreover, to remove the sporadic high correlation between remote locations that does not encourage the spatial coherent structure of MHW events, we confine the neighbourhood to a 10° × 10° × 200 m cube centred at the grid cell in question when performing the NN smoothing.This also lowers the computational burden for searching the neighbourhood.Finally, as the NN smoothing may misrepresent small-volume MHWs, we discard MHW events occupying less than 125 grid cells, corresponding to a latitudinally varying volume criterion from 1.6 × 10 4 km 3 (at the equator) to 0.8 × 10 4 km 3 (at 60°).As the value of this grid cell number threshold increases, fewer MHW events are retained as expected yet the proportions of different kinds of MHW events keep relatively stable (Supplementary Table 6).
To evaluate effects of the smoothing strategy on the identified MHW events, we replace the correlation-based NN smoothing with another kind of NN smoothing that defines the neighbourhood as a 5° × 5° × 50 m cube centred at the grid cell in question.Using the 5° × 5° × 50 m cubic neighbourhood leads to an even larger proportion of subsurface and mixed MHW events but a smaller proportion of surface MHW events (Supplementary Table 7).Therefore, the correlation-based NN smoothing does not overestimate the proportions of subsurface and mixed MHW events.Furthermore, as the neighbourhood in the correlation-based NN smoothing is not affected by the strong anisotropy of the ocean and statistically more meaningful than the 5° × 5° × 50 m cubic neighbourhood, the results derived from the correlation-based NN smoothing should be quantitatively more reliable.
We track the evolution of a MHW event by linking their consecutive snapshots.This can be done by generalizing the MHW tracking method in a two-dimensional spatial domain (that is, at the sea surface) proposed by Sun et al. 33 to a three-dimensional spatial domain.Under the premise that the spatial structure of a MHW event should continuously evolve with time, MHWs at two consecutive snapshots are linked if the fraction of their overlapped domain (FOD) is larger than 0.5, following where Ω t i (Ω t+1 j ) represents the set of grid cells occupied by the ith (jth) MHW on the tth (t + 1st) snapshot denoted as MHW t,i (MHW t+1, j ), the absolute value symbol represents the cardinality of a set (that is, the number of grid cells) and ∩ represents the intersection of sets.
For a MHW event without splitting or merging during its life cycle, it can be tracked from emergence to disappearance without ambiguity.Here the start (end) of a MHW event is defined as the time when its volume first (last) exceeds 125 grid cells.On one hand, MHW t,i is assigned as the start of a MHW event if it has no associated predecessor with FOD > 0.5.On the other hand, MHW t,i is assigned as the end of a MHW event if there is no successor with FOD > 0.5.However, splitting and merging processes of MHWs complicate the above tracking procedure by breaking the one-to-one linkage between MHWs at consecutive snapshots.Here we adopt the strategy proposed by Sun et al. 33 to handle the splitting and merging processes.Simply put, all the parts of a split MHW event are treated to belong to a single event and tracked as a whole, whereas the merged MHW events are treated as different MHW events and tracked separately.One advantage of this strategy is that https://doi.org/10.1038/s41561-023-01325-w it conserves the number of MHW events during the MHW splitting or merging.Supplementary Fig. 10 illustrates the tracking of a MHW event under the different scenarios, that is, splitting, merging and otherwise.

Properties of MHW events in the four-dimensional domain
Properties to characterize a MHW event in this study include the time-mean horizontal area A(z) and the time-mean intensity I(z).The horizontal area occupied by a MHW event is first computed at each depth at each time and then temporally averaged during its life cycle to yield A(z).The intensity of an event, measured as the horizontal mean temperature anomaly relative to the climatological seasonal cycle (T clim ), is first computed at each depth at each time and then temporally averaged over its life cycle to yield I(z).
It should be noted that a MHW event does not necessarily occupy the entire water column (0-200 m) during its life cycle.Meaningful values have to be prescribed to the intensity and horizontal area of some MHW event at the depths not occupied by that MHW event to compute A(z) and I(z) and their mean values over multiple MHW events.It is natural to set the horizontal area of some MHW event as zero at the depths not occupied by that MHW event, as its value cannot be positive there by definition and must be non-negative to be meaningful, making the zero as the sole choice.However, application of such zero-padding to the intensity of some MHW event will cause an undesirable jump discontinuity of the intensity at the interface between the water column occupied by that MHW event and otherwise.This is because MHW events are identified based on the exceedance of the threshold θ yet their intensity is measured as the temperature anomaly relative to T clim 18 , in which case the lower bound of the intensity is not zero but a positive quantity ∆ = θ − T clim .To remedy this problem, the intensity at the depths where some MHW event does not extend is computed by averaging the temperature anomaly relative to T clim over a horizontal domain extrapolated from the horizontal domain at the nearest depth occupied by that MHW event.In this case, the intensity varies continuously with depth but does not necessarily decrease to zero outside the vertical extent of a MHW event.Supplementary Fig. 11 compares the global mean I(z) computed based on the zero-padding and the extrapolation methods.Although the global mean I(z) from the zero-padding biases low compared to that from the extrapolation especially near the sea surface for subsurface MHW events as expected, their vertical structures are qualitatively consistent with each other.Therefore, major conclusions of the global mean I(z) are not substantially affected by the ways to compute the intensity at the depths not occupied by some MHW event.

Synthetic temperature data
Mean-state warming is suggested to play a central role in the increasing frequency of MHW events during the past several decades 20,42 .However, the long-term changes of temperature in the reanalysis datasets do not agree with the observed value derived from the Institute of Atmospheric Physics global ocean temperature gridded product 44 , which may bias the long-term changes of MHW events in the reanalysis datasets.To suppress the bias, we replace the linear trend of temperature at each depth at each grid cell in the reanalysis datasets with their observed counterparts.The synthetic temperature data can be treated as a hybrid of long-term temperature changes in the observation and temperature variability in the reanalysis datasets.

ig. 1 |
Geographic distribution of different kinds of MHW events.a,d,g, Climatological mean annual number of the surface (a), subsurface (d) and mixed MHW events (g) occurring in the individual 20° × 10° (longitude × latitude) half-overlapping bins derived from the GLORYS.

Fig. 2 |
Fig. 2 | Vertical distribution of MHW properties of different kinds of MHW events.a-c, Global mean vertical profile of the time-mean horizontal area A(z) of the surface (a), subsurface (b) and mixed MHW events (c).d-f, Same as a-c but for the time-mean intensity I(z).The horizontal area of a MHW event is first computed at each depth at each time and then temporally averaged over

Fig. 4 |
Fig. 4 | Long-term changes of different kinds of MHW events under anthropogenic climate changes.a, Time series of monthly mean temperature averaged over the global ocean (60° N-60° S) derived from the observation.Here the time-mean temperature at each depth is subtracted from the corresponding time series to facilitate the illustration.b, Time series of annual number of the global surface, subsurface and mixed MHW events in the GLORYS.c,d, Same as b but for the HYCOM (c) and ECCO2 (d).Slopes of the fitted linear trends and their corresponding 95% confidence intervals are shown in the top left corner in b, c and d.The pink, blue and purple lines correspond to the results for the surface, subsurface and mixed MHW events, respectively.