Frequent marine heatwaves hidden below the surface of the global ocean

The study addresses how marine heatwaves (MHWs)—typically identified and analyzed using sea surface temperature (SST)—manifest below the ocean surface and how often such events lack surface signatures. While MHWs drive abrupt ecological and socio-economic impacts, most prior global assessments have relied on satellite SST due to observational convenience, leaving the vertical structure in the epipelagic zone (0–200 m) poorly characterized. Case studies indicate diverse vertical structures, including subsurface intensification, downward migration, and purely subsurface events with no SST imprint. Traditional time-series-based methods at single locations lose crucial spatial and vertical information and can misrepresent events undergoing vertical displacement or deformation. This study extends MHW identification to a four-dimensional spatio-temporal framework to quantify the frequency, distribution, properties, and long-term changes of surface, subsurface, and mixed MHWs globally.

Prior work has documented numerous high-impact MHWs and linked increasing frequency and intensity to anthropogenic warming. Most global studies identify MHWs from SST time series, reporting trends and drivers (for example, atmospheric forcing, ocean advection, mesoscale variability) and ecological consequences (mass mortalities, coral bleaching, shifts in fish distributions). Regional and process studies have revealed various vertical structures of MHWs, including events with peak anomalies at the surface or subsurface, downward migration from surface to interior (such as the Northeast Pacific “Blob”), and events confined to subsurface layers without SST anomalies. However, a coherent global description of subsurface MHW occurrence and properties has been lacking due to observational sparsity and methodological limitations of 1D time-series approaches. The present work situates itself in this gap by providing a global, four-dimensional characterization across datasets.

Data: Three ocean reanalysis datasets were used: GLORYS (daily 1/12°, 1993–2020), HYCOM (3-hourly, 0.08°, 1993–2012), and ECCO2 (3-day mean, 0.25°, 1992–2020). Temperatures were horizontally smoothed with a 1°×1° running mean, interpolated to a uniform 1°×1°×10 m grid over 0–200 m, and re-timed to daily resolution by interpolation or bin averaging. Data poleward of 60°N/S were excluded. Identification framework: Four-dimensional MHW events are defined as discrete, prolonged (≥5 consecutive days), anomalously warm entities whose compact three-dimensional spatial structure evolves continuously in time. Step 1 (temporal detection): At each grid cell and depth, MHWs are identified where daily temperature exceeds a seasonally varying 90th percentile threshold computed over each dataset’s full period. Step 2 (spatial coherence via NN smoothing): The raw binary MHW field is smoothed in 3D using a correlation-based nearest-neighbour approach. Grid cells are deemed MHW if the majority of neighbours (defined within a 10°×10°×200 m cube) have high correlation (r≥0.5) in their temperature anomaly time series; r is a tuned parameter (tested 0.4–0.6). This correlation-based neighbourhood mitigates ocean anisotropy versus simple spatial distance. Sensitivity tests using a fixed 5°×5°×50 m cube produced even more subsurface/mixed events, indicating the correlation-based approach does not overestimate subsurface incidence. Size filter: Events occupying <125 grid cells (volume ~1.6×10^4 km^3 at equator to 0.8×10^4 km^3 at 60°) were excluded to avoid artifacts and very small features. Tracking: Consecutive 3D snapshots are linked if the fraction of overlapped domain (FOD) exceeds 0.5. Splitting and merging are handled by assigning all split parts to a single event and tracking merged events separately, conserving event counts. Classification by vertical signature: Events with continuous SST (surface) signals throughout their life cycle are classed as surface MHWs; those with no surface signal at any time as subsurface MHWs; others as mixed. Event properties: For each event, time-mean horizontal area A(z) and intensity I(z) (horizontal mean temperature anomaly relative to the climatological seasonal cycle) are computed at each depth; intensity at depths not occupied by the event is estimated by horizontal extrapolation from the nearest occupied depth to ensure vertical continuity (as zero-padding biases I(z) low). Anthropogenic trend adjustment: To assess long-term changes, reanalysis linear temperature trends at each grid cell and depth were replaced with observed trends from an independent ocean temperature analysis to create synthetic temperature data, preserving variability while aligning mean-state warming with observations. MHWs were re-identified in these synthetic datasets to evaluate trends in annual event counts. Regional analyses: Geographic distributions and vertical property profiles were summarized for key regions (western boundary current extensions, Southern Ocean, subtropical gyre interiors, equatorial Pacific) using 20°×10° bins and event geometric centers.

• Large hidden subsurface incidence: Across datasets, only about half of MHW events exhibit continuous surface signatures; roughly one-third are entirely subsurface with no SST imprint throughout their life cycles. - Event frequencies by dataset (global, 60°S–60°N, mean annual counts): GLORYS ~308.5; HYCOM ~374.4; ECCO2 ~218.8. Proportions (surface/subsurface/mixed): GLORYS 46.9%/31.2%/21.9%; HYCOM 32.6%/42.4%/25.0%; ECCO2 52.8%/27.7%/19.5%. - Geographic patterns: Surface MHWs are most frequent in western boundary current extensions and the Southern Ocean, consistent with high SST variability regions. Subsurface MHWs occur predominantly in subtropical gyre interiors. Mixed events show a hybrid, relatively uniform distribution with elevated frequencies in WBCEs, SO, and SGI. Many mixed events lack surface signals at onset and/or at termination, consistent with processes like mixed-layer deepening, isopycnal lifting, subduction, or downward isopycnal deflection. - Vertical structure of properties: For surface events, mean time-averaged horizontal area A(z) peaks at the surface and decreases with depth; intensity I(z) is comparatively uniform in the upper 200 m, indicating subsurface anomalies of similar magnitude to surface anomalies. For subsurface events, A(z) increases with depth (zero at surface, peak at 200 m), suggesting extension into the mesopelagic; I(z) peaks around ~100 m and attenuates above and below. Mixed events exhibit relatively uniform A(z) and a subsurface-peaked I(z). These patterns are qualitatively consistent across GLORYS, HYCOM, and ECCO2. - Long-term changes under anthropogenic warming (synthetic datasets with observed trends): Annual counts of surface, subsurface, and mixed MHW events all show significant positive trends over the past three decades (P<0.05). Trend slopes (events per decade) span: surface ≈46–60; subsurface ≈20–49; mixed ≈21–42, varying by dataset. Mean-state warming is the primary driver of these increases, with changes in variability playing minor roles. Faster surface warming relative to depth explains larger trends for surface events. Despite imposing the same observed trends, sensitivity of counts to warming differs among reanalyses, likely reflecting differing higher-order statistics.

By extending MHW identification to a four-dimensional framework, the study reveals that a substantial fraction of MHWs are undetectable from SST alone, directly addressing the limitation of traditional 1D time-series approaches. The prevalence of subsurface and mixed events, their distinct regional distributions, and their vertical property profiles underscore the importance of subsurface processes (advection, isopycnal motions, mixed-layer dynamics, subduction) in shaping MHWs. The significant increases in annual counts across all event types are primarily attributable to mean-state warming, consistent with theoretical expectations and prior SST-based findings, but newly quantified here for subsurface-resident events. These results emphasize that ecosystem and fisheries management strategies based solely on SST will miss a large share of impactful warm anomalies, particularly in subtropical gyres, and that subsurface monitoring and modeling are essential for comprehensive detection, attribution, and projection. The cross-dataset qualitative consistency supports robustness, while quantitative differences highlight sensitivity to background variability and reanalysis specifics.

The study provides the first global, four-dimensional census of MHWs in the epipelagic ocean, demonstrating that many events either spend part of their life cycles below the surface or remain entirely subsurface without SST signatures. It quantifies regional patterns, vertical structures of area and intensity, and significant anthropogenically driven increases in event counts for surface, subsurface, and mixed categories, predominantly driven by mean-state warming. Methodologically, it introduces a robust 4D identification and tracking framework adaptable to different depth ranges and resolutions. Future directions include extending analyses deeper (for example, to 400 m and beyond, which further increases the subsurface fraction), refining detection of smaller-scale events (such as those associated with mesoscale eddies) by tuning size thresholds and smoothing parameters, leveraging monthly 3D datasets and statistical methods for proxy estimation, and integrating satellite sea surface height to infer subsurface events. Expansion and enhancement of subsurface observing systems and advanced inference/forecasting tools are crucial to monitor and manage MHW impacts effectively.

• Observational sparsity below the surface limits direct validation; reliance on reanalyses introduces model- and assimilation-dependent biases. - Long-term temperature trends in reanalyses differ from observations; a synthetic trend-correction was applied, yet datasets still display differing sensitivities due to varying higher-order statistics. - Events smaller than 125 grid cells were excluded, potentially omitting mesoscale-eddy-associated MHWs and biasing counts toward larger events. - Analysis confined to 0–200 m and 60°S–60°N; extending depth or latitude could change proportions and distributions. - Classification depends on smoothing parameters (for example, correlation threshold r and neighbourhood extent) and tracking criteria (FOD>0.5), which, while sensitivity-tested, may influence quantitative results. - No one-to-one correspondence between observed and reanalysis events; focus is on consistent statistics rather than exact event matching. - Monthly gridded observational 3D datasets cannot directly resolve MHWs, limiting independent subsurface validation.