Marine heatwaves (MHWs), defined as prolonged periods of anomalously warm sea surface temperatures (SSTs), are increasingly recognized as a major threat to marine ecosystems. These extreme events cause devastating consequences, including mass mortality of marine organisms, coral bleaching, shifts in fish distributions, and alterations to the structure and function of marine ecosystems. The impacts of MHWs are likely to intensify with continued ocean warming. While substantial research has focused on surface MHWs using readily available satellite SST data, a comprehensive understanding of their vertical structure remains elusive. This knowledge gap is critical because the epipelagic zone, the ocean's uppermost layer, which extends to approximately 200 meters and harbors the majority of marine life, does not always exhibit temperature changes that mirror the surface. Many MHWs originate or persist below the surface, making SST-based detection methods incomplete. Several studies have shown this with case studies that demonstrate varied vertical structures of MHWs in different regions, potentially due to complex drivers, such as the downward migration of MHWs and events existing entirely within the subsurface layers. Characterizing the vertical structures of MHWs is challenging due to the limited availability of subsurface oceanographic data and the lack of suitable methodologies for four-dimensional (spatio-temporal) analysis. Current approaches often rely on one-dimensional time series of temperature at individual locations, neglecting the important spatial dynamics of MHW events and leading to an incomplete picture. This study addresses these challenges by applying a four-dimensional framework to investigate the prevalence and characteristics of MHWs across the global ocean.

Previous studies on marine heatwaves have primarily relied on sea surface temperature (SST) data from satellites, offering valuable insights into the spatial extent and frequency of surface events. These studies have documented the increasing frequency and intensity of MHWs globally, directly linking them to anthropogenic climate change and the resulting ocean warming. However, the focus on SST data has limited our understanding of the subsurface characteristics of these events. While some studies have provided regional insights into subsurface MHW behavior, a global analysis using a consistent methodology across different ocean regions is lacking. Past research highlights the varied vertical structures of MHW events, showing that some are confined to the surface layer, while others extend to deeper layers or even reside entirely in the subsurface. The existing methods for identifying MHWs primarily use time series data from individual locations, overlooking the spatial characteristics and four-dimensional dynamics that are crucial for a thorough understanding of their behavior. This study's novel approach directly addresses this limitation by implementing a four-dimensional analysis framework.

This study employs a novel four-dimensional spatio-temporal framework to identify and characterize marine heatwaves (MHWs) in the global ocean's epipelagic zone (0-200m). The methodology utilizes temperature data from three widely used ocean reanalysis datasets: the Global Ocean Physics Reanalysis (GLORYS), the Hybrid Coordinate Ocean Model (HYCOM), and the Estimating the Circulation and Climate of the Ocean Phase II (ECCO2). These datasets were chosen for their broad spatial coverage, relatively long time spans, and varied data assimilation techniques. The data from each dataset were pre-processed using a 1° x 1° horizontal running mean to smooth out high-frequency noise before being interpolated onto a uniform 1° x 1° x 10m grid across the upper 200m of the ocean. The temporal resolutions of the datasets were standardized to daily values. Regions beyond 60°N and 60°S were excluded to minimize the influence of sea ice. The core of the methodology lies in identifying MHWs within a four-dimensional framework. This involved a two-step process: first, identifying potential MHW events based on temperature time series at each grid cell exceeding a seasonally varying 90th percentile threshold; second, smoothing the binary (MHW/non-MHW) result spatially using a nearest neighbor (NN) method based on correlation distance, rather than Euclidean distance, to account for the ocean's inherent spatial anisotropy. A correlation coefficient threshold (rc = 0.5) was used to define the neighborhood for smoothing, with sensitivity analyses performed to assess the impact of this parameter choice. The choice of correlation distance was justified by the fact that the use of the standard spatial distance was inappropriate given the inherent anisotropy of the ocean. Events occupying fewer than 125 grid cells were excluded to ensure the robustness of the methodology. This step aimed to eliminate smaller events which might have been artifacts of the smoothing process. A tracking algorithm was implemented to follow the temporal evolution of each MHW event by linking consecutive snapshots based on their fractional overlap. This algorithm addressed the complexities introduced by MHW splitting and merging. The identified MHWs were then classified into three categories: surface MHWs (continuous surface signals), subsurface MHWs (no surface signals), and mixed MHWs (both surface and subsurface signals). Finally, the horizontal area and intensity of each MHW were calculated at each depth and then temporally averaged over its lifespan. To minimize biases stemming from inconsistencies in long-term temperature trends across different reanalysis datasets compared with observations, the analysis was repeated using synthetic temperature data. In this adjusted dataset, the linear trend from each reanalysis dataset was replaced with the observed trend from a high-resolution global ocean temperature gridded dataset, maintaining the temporal variability from the original reanalysis products.

The analysis of three ocean reanalysis datasets (GLORYS, HYCOM, and ECCO2) revealed a previously underestimated prevalence of subsurface marine heatwaves (MHWs). Across the datasets, the proportion of MHWs that lacked continuous surface signals ranged from 27.7% to 42.4%. This is a substantial fraction, implying that significant MHW activity occurs below the surface and is undetectable using surface temperature data alone. The geographic distribution of MHW types showed considerable heterogeneity, with surface MHWs predominantly located in western boundary current extensions and the Southern Ocean, while subsurface MHWs were concentrated in subtropical gyre interiors. Mixed MHWs showed a more uniform distribution, exhibiting features of both surface and subsurface events. Further analysis of MHW properties revealed distinct vertical structures. Surface MHWs exhibited a peak in horizontal area at the surface, decreasing with depth. Conversely, subsurface MHWs showed the opposite, with the largest area at 200m and decreasing towards the surface. Mixed MHWs displayed a more uniform vertical area distribution. The intensity (temperature anomaly) of surface MHWs showed relative consistency across depths, indicating that subsurface anomalies were as strong as those at the surface. Subsurface and mixed MHWs demonstrated intensity peaks at intermediate depths (around 100m), reflecting their subsurface nature. Analysis of the long-term trend in MHW frequency indicated a significant increase in all MHW types over the past three decades. The increase was most pronounced for surface MHWs, with trends ranging from 46 to 60 events per decade across the datasets, reflecting the faster warming rate at the surface. Subsurface and mixed events also showed significant increases, albeit smaller in magnitude. These increases are largely attributed to mean-state warming, as analyses using synthetic data controlled for this factor corroborated the results. However, the quantitative differences in the increase across the three datasets suggest differing sensitivity of MHW events to warming based on higher-order temperature statistics.

The findings of this study challenge the traditional reliance on SST data for MHW detection, revealing a significant number of subsurface events previously unseen. The large proportion of subsurface and mixed MHWs underscores the limitations of current monitoring methods based solely on surface temperature observations. The significant increase in subsurface MHWs over the last three decades, directly linked to subsurface warming, emphasizes the consequences of anthropogenic climate change extending beyond the sea surface. This has profound implications for marine ecosystems, as these subsurface events can have substantial biological effects, despite not being directly detectable at the surface. The observed geographic distribution of MHW types reflects the complex oceanographic processes driving these events. The greater frequency in western boundary current extensions and the Southern Ocean points to the importance of ocean currents and eddies in MHW formation. The concentration of subsurface events in subtropical gyre interiors suggests distinct mechanisms responsible for the generation and maintenance of these hidden events. The study’s finding of similar intensities between surface and subsurface warm anomalies emphasizes the potential impact of subsurface MHWs on marine organisms, even in the absence of surface manifestations. Future research should prioritize the development of improved MHW detection methods capable of utilizing subsurface data, incorporating four-dimensional analysis, and acknowledging the complexities of oceanographic processes that govern the generation and evolution of subsurface heatwaves.

This research demonstrates the widespread occurrence of subsurface marine heatwaves, significantly altering our understanding of these extreme events and their impacts on marine ecosystems. The study highlights the limitations of solely using sea surface temperature for detection and underscores the necessity of integrating subsurface observations and four-dimensional analysis into future MHW monitoring and prediction strategies. Future studies should focus on improving subsurface observation capabilities and developing advanced detection methods to fully capture the impacts of these previously hidden events. The findings emphasize the urgent need for a more comprehensive and accurate monitoring system for marine heatwaves that accounts for their complex vertical structures.

The study's reliance on ocean reanalysis data introduces potential limitations due to model biases and uncertainties inherent in the data assimilation techniques. The choice of a specific correlation coefficient threshold for the nearest neighbor smoothing algorithm and the minimum grid cell size for MHW identification might influence the precise quantification of the different MHW types, although sensitivity tests were conducted to evaluate the robustness of the approach. While the study aimed to address some biases in long-term temperature trends by using synthetic data, it is important to acknowledge that residual model-related biases in temperature variability might still exist. Future research with direct measurements of subsurface temperatures could improve the accuracy and precision of the findings. The study is primarily limited to the epipelagic zone (0-200m), which could be extended to examine deeper water layers for a more comprehensive view of subsurface MHWs.