Introduction

Marine heatwaves (MHWs) are defined as periods when the local sea surface temperature (SST) is significantly higher than typical for the time of year at a specified location. MHWs have occurred throughout the global ocean (Hobday et al. 2016;Holbrook et al. 2019). Typically, MHWs are examined through a local lens. Even when the drivers of marine heatwaves are well known for a particular region (e.g., persistent anticyclonic atmospheric circulation over the North Pacific), the evolution of individual MHWs in these regions has varied considerably (Amaya et al. 2020;Bond et al. 2015;Fewings and Brown 2019).

The motivation to understand the evolution of MHWs is owed to the vulnerability of marine ecosystems to temperature extremes (Smale et al. 2019). MHWs have led to mass mortality in marine invertebrates (Oliver et al. 2017;Garrabou et al. 2009), species range shifts (Mills et al. 2013), habitat destruction including coral bleaching (Hughes et al. 2017), and harmful algal blooms (McCabe et al. 2016). Failure to anticipate the destructive impacts of MHWs leads to fishery management challenges, including changes to the supply chain and loss in value of commercially harvested species (Mills et al. 2013;Pershing et al. 2019;Cheung and Fr ölicher 2020). MHWs can also impact regional atmospheric circulation that perturbs weather patterns over land. Such events have been associated with extreme drought, leading to impacts on agriculture (Williams et al. 2015;Rodriguez 2021) and terrestrial heat extremes (McKinnon and Deser 2018).

By definition, MHWs represent the extremely warm end of the distribution of local sea surface temperature anomalies. Previous studies have used the 90th (Oliver et al. 2018;Hobday et al. 2016) or the 99th (Darmaraki et al. 2019;Frölicher et al. 2018) percentile of the SST distribution to define extremes, where a MHW event is identified when SST exceeds this threshold relative to a long-term seasonal climatology for at least a certain period of time, e.g., 5 days (Hobday et al. 2016).

The distribution of MHWs is influenced by the mean state, natural variability, and long-term anthropogenic change (Frölicher et al. 2018;Oliver et al. 2018). Regions with a large SST variance, for example, in the vicinity of western boundary currents and their extensions, as well as in the equatorial Pacific cold tongue, have the highest MHW intensities globally (Oliver et al. 2018). In addition, extremely long duration MHWs have been linked to modes of interannual to decadal variability in the climate system (Holbrook et al. 2019;Scannell et al. 2016).

Natural variability such as El Niño-Southern Oscillation (ENSO) can impact the presence and persistence of MHWs in the midlatitudes through atmospheric teleconnections from the tropics. For example, anomalies in atmospheric deep convection over the tropics can initiate atmospheric planetaryscale waves that propagate to the midlatitudes where they generate MHWs through changes in local atmospheric conditions, e.g., cloud cover (Hartmann 2015). Large-scale modes of decadal SST variability linked to tropical climate variability, such as the interdecadal Pacific oscillation (Power et al. 1999), can suppress or enhance the likelihood of MHW occurrences depending on the phase and amplitude of the mode (Holbrook et al. 2019;Scannell et al. 2016). They can influence the severity and duration of MHWs by altering the mean strength, direction, and location of ocean currents and heat transport, as well as modulate air-sea heat flux (Perkins- Kirkpatrick et al. 2019;Di Lorenzo and Mantua 2016;Feng et al. 2013).

Interannual and decadal variability within the climate system can be explored using an empirical orthogonal function (EOF) decomposition of climate anomalies. The first few EOF modes have been used to help characterize the spatial structure of some MHWs and their time scales (Di Lorenzo and Mantua 2016). EOFs have also been used to explain the spatial patterns and the long-lived persistence of prominent MHWs (Amaya et al. 2020;Fewings and Brown 2019;Oliver et al. 2018;Di Lorenzo and Mantua 2016). However, EOFs are statistical descriptions of the variability and do not necessarily encapsulate dynamical information. In addition, using a limited number of EOFs to describe the spatiotemporal evolution of MHWs can give an incomplete picture of the evolution of short-lived and smaller-scale MHWs that may not be easily linked to the most dominant modes of climate variability.

Retrospective and contemporaneous studies have generally relied on pointwise metrics (Sen Gupta et al. 2020;Hobday et al. 2018;Oliver et al. 2018), fixed region heat budget analyses (Xu et al. 2018;Oliver et al. 2017;Bond et al. 2015;Chen et al. 2014), or EOFs (Di Lorenzo and Mantua 2016) to characterize the drivers of specific MHW events and to describe their characteristics. These approaches have been successful in explaining the local processes and remote drivers responsible for specific MHWs (Sun et al. 2023a). Here, we expand this view by characterizing the spatiotemporal evolution of MHWs. This description of MHW evolution takes advantage of the 3D evolving field of global SST by detecting and tracking MHWs and characterizing their shape, size, location, duration, and intensity, which may help to identify new patterns in how MHWs evolve and facilitate investigations of MHW dynamics. The 3D here refers to the space and time dimensions (x, y, t) and not spatial dimensions (x, y, z). To do this, we present an object-tracking algorithm called Ocetrac and use it to explore the large-scale spatial connectivity of MHWs as they evolve in time.

Object tracking has been used to describe both atmospheric and oceanic phenomena. For instance, an enhanced watershed method was used to identify hailstorm objects applied to observed gridded radar reflectivity and column-integrated graupel mass estimates from a numerical weather prediction (NWP) model (Gagne et al. 2017). The enhanced watershed method (Lakshmanan et al. 2009) reduces the volume of data that needs to be processed by optimally searching for the local maxima in the storm field and growing the storm object until both area and intensity criteria are met. As with Ocetrac, the watershed object-identification method is parameter sensitive. TempestExtremes, another software package for feature detection, tracking, and analysis of extreme events, uses the algorithmic framework "MapReduce" to first detect candidate extreme events like tropical and extratropical cyclones and atmospheric rivers based on thresholds and other metrics. It also connects these events in time.

We use the following definitions in our analysis and description of the spatiotemporal evolution of MHWs (Table 1). Features are individual points where SST is above the locally defined threshold for 1 month. A MHW object is a spatially coherent collection of features. A MHW event is composed of tracked and linked objects. We apply Ocetrac to monthly SST data from 1981 through 2021 to track the evolution of all MHWs globally and examine the distribution of three key MHW metrics (size, intensity, and duration). Four MHWs are used as case studies to describe new insight available with this framework.

Methods

a. Data, preprocessing, and anomaly detection

We analyze the monthly global maps of SST from the 0.258 longitude 3 0.258 latitude gridded Optimum Interpolation SST, version 2.1 (OISSTv2.1), dataset that extends from September 1981 through January 2021. The OISSTv2.1 combines satellite Advanced Very High Resolution Radiometer (AVHRR-only) with observations from ship, buoy, and in situ measurements (including Argo floats and drifters), while accounting for platform differences and using interpolations to fill gaps in the satellite data (Reynolds et al. 2002(Reynolds et al. , 2007)). We create a mask over the Arctic (.658N) and Antarctic (.708S) Oceans to remove data in these regions and to avoid influence from seasonal sea ice and locations where the OISSTv2.1 data are less reliable (Fig. 1).

Using the global maps of SST, we remove the mean, linear trend, and seasonal cycle from September 1981 through January 2021 to compute anomalies (Fig. 2). The total decomposition of the climatological SST is represented as

where the fit SST fit is the linear combination of the mean SST m (Fig. 1a), linear trend SST t , and annual and semiannual harmonics SST s at each grid point. The coefficients of SST fit are found using the least squares regression fit to monthly SST computed over the time period September 1981-August 2021 (473 months). We define SST anomalies SST a as the difference between monthly SST and SST fit such that SST a 5 SST 2 SST fit :

(2)

We remove the trend from SST a so that it has roughly stationary statistics and reflects the common internal ocean and climate dynamics from the beginning to the end of the record. If the long-term trend is not removed, the likelihood of MHWs will increase with time almost everywhere due to the positive SST trend associated with global warming. The trend is the largest in midlatitudes in the subtropical gyres, especially in the northwest Atlantic, western North Pacific, and western South Pacific. We also remove the climatological seasonal cycle from SST a so that we can track MHW events that evolve across multiple seasons or regions with different seasonal cycles. The seasonal cycle of SST is the largest in the subtropics and midlatitudes and weaker in the subpolar and polar regions and the weakest in the deep tropics (Fig. 1c).

Since there are strong spatial variations (Fig. 1b) as well as seasonal cycles in SST a variance, we standardize SST a by dividing by the respective local monthly standard deviation of SST a over the entire period. The resulting standardized anomaly fields SST * a have uniform variance across the globe and in each month. High standard deviations of SST a occur in the eastern equatorial Pacific, western boundary currents, the region connecting the Indian Ocean to the South Atlantic, and in frontal zones with large SST gradients. The subtropics, southern midlatitudes, equatorial Atlantic Ocean, equatorial Indian Ocean, and western tropical Pacific have low standard deviations (Fig. 1b). There is a significant seasonal cycle in SST a variance over most of the global oceans. For example, a stronger SST a variance is observed during springtime near the Gulf Stream, in boreal wintertime in the Pacific cold tongue, in the late summer and early autumn in the subpolar northeast Pacific and northwest Atlantic, and generally during the summertime at polar latitudes (in the online supplemental material 1).

Various methods exist to define SST anomalies and temperature thresholds. Here, we define a given pixel and month to be an MHW candidate if it exceeds the 90th percentile of SST * a in that pixel (Hobday et al. 2016). The results are influenced by the choice of the 90th percentile and the length of the record. With only 40 years of monthly data, our ability to accurately quantify extremes defined by a higher percentile (like the 99th percentile) record is limited by the record length. In this study, we did not assess the sensitivity of our findings to the choice of percentile threshold.

To identify MHWs from the monthly maps of SST * a , we search for candidate MHW points where SST * a exceeds the 90th percentile at each grid point. It may be noted that most of the steps of the detection of MHW candidate points involve choices and that these choices will certainly have a quantitative if not qualitative impact on the results. Not all of the sensitivities associated with these choices are exhaustively documented in this paper, but we discuss some of them at the end of the paper. The purpose here and in the rest of this methods section is to present choices of parameters that yield reasonable results and set the stage for further work.

b. Multiple object tracking

The binary maps highlighting the MHW candidate points produced by the anomaly detection algorithm in section 2a exhibit some small-scale patchiness near the grid scale, which contributes to image noise. To identify the coherent patches a is shown in red and has been divided by its monthly standard deviation. In this example, the monthly standard deviations are similar to 18C (Fig. 1b), so the unitless SSTa * is fairly similar to but not identical to SST a in degrees Celsius. The red circles indicate when SST * a exceeded the 90th percentile of SST * a (shown by the dashed line) computed over the entire period from September 1981 through January 2021.

of anomalies and to simplify processing, these maps are enhanced and cleaned to reduce complexity and remove smaller-scale details and then extracted for features. Smoothing techniques are used to reduce noise and enhance the SST a map quality for object detection. This enhances object features and improves the accuracy of edge detection, making it easier for the object detection algorithm to identify the key patterns and evolution.

The broad overview of how the tracking algorithm works is that features, which are output of the anomaly detection algorithm, are first selected within an image before applying morphological operations. A structuring element, which is a shape used to differentiate objects from others based on their shape or spatial orientation, is then chosen to capture relevant spatial relationships within the image. Morphological operations are applied to the input image using the structuring element, modifying the shape and structure of objects in the image based on the interaction between the image pixels and the structuring element. This is applied iteratively with each iteration refining the image structure.

The standardized anomaly SST * a maps with the MHW candidate points produced by the anomaly detection algorithm in sections 2a and 2b are transformed into a binary image where ones correspond to candidate MHW grid points and zeros correspond to background grid points. Each monthly map is treated as a separate image. Our goal is to identify the groupings of ones that define a MHW object, which meet the defined spatial characteristics in terms of structure and size.

Image processing terminology is defined in Table 1 and illustrated in Fig. 3.

We use mathematical morphology operations from the SciPy multidimensional image processing Python package to remove small, isolated features and to fill small holes within feature clusters. The structuring element S is defined by a quadratic surface with a morphological radius R, which is measured in the number of grid cells, where S 5 x 2 1 y 2 , R 2 :

(3)

Here, x and y define the points within a circle of radius R, where R is defined by the number of grid cells in the map of SST (supplemental material 2). The matrix S is transformed into a two-dimensional binary image with origin at the center of a box with sides of length R. Elements located where S , R 2 have value 1 and elsewhere a value of 0 (Fig. 3). The units of S are in degrees per unit resolution of the grid (e.g., an R of 8 on a 1/48 grid is equal to 28 latitude or longitude).

The term R is chosen to be approximately the same size and shape of the desired objects. The term R is a user-defined parameter and is measured in gridcell units. The other user-defined parameter is P, which is the value of the percentile of object area in square kilometers (km 2 ). The term P gives the minimum size of objects based on a size percentile threshold and is used to filter out objects smaller than the threshold. The structuring element is used to scan over the entire binary image when conducting morphological operations (Gonzalez and Woods 2002). Dilation and erosion are the two most common morphological operations, and both are used here. Erosion eliminates isolated and small features by shrinking features, or in terms of pixels, it removes pixels near the boundaries of the MHWs. Dilation is the opposite of erosion and fills small holes within features or in terms of pixels; dilation adds pixels near the boundaries of features. The precise effect of dilation and erosion, or the number and location of added or removed pixels, depends on the size and shape of the structuring element and thus depends on the radius R.

Erosion and dilation are often performed in succession (Fig. 4). Morphological opening is erosion followed by dilation. Opening is used to eliminate small features while preserving the shape and size of larger features in the image. Alternatively, morphological closing is dilation followed by erosion. Closing fills small holes within features while also preserving the shape and size of other features in the image. Both opening and closing are used to remove small features and smooth the borders of larger features. Generally, these operations can be combined in different sequences to transform an image to isolate specific features. The order and number of morphological operations will lead to different results. Here, we implement a sequence of morphological closing followed by opening, as we found this to effectively clean the feature images and isolate coherent clusters of features, which are now called objects that can be tracked in space and time (Fig. 4; supplemental material 3).

After morphological operations, there are two rounds of object identification and labeling. First, 2D (spatial) objects are detected by labeling connected 2D objects from binary images using Scikit-image's measure module in Python. Zerovalued pixels are background pixels. We define objects when two or more neighboring features with the same value are connected either adjacent or diagonal from each other (e.g., orange pixels in Fig. 3a). The resulting 2D objects are assigned a unique label. This process is repeated for each time step, with each time step as a separate instance of 2D object identification. For each 2D unique object, we use the latitude and longitude coordinates from the Scikit-image's regionprops module to calculate the total object area. Using the distribution of all object areas from September 1981 through January 2021, we calculate the minimum area. Each object's area is then compared to the minimum area value. If an object's area is smaller than the minimum area value, the object is discarded. If an object's area is larger than the minimum area value, the object is retained. The accepted area is the area of the object that has not been discarded. For our purposes, we use the 75th percentile of object area (km 2 ) for the default value of P (Fig. 4). We discuss the sensitivity of MHW characteristics to P in section 3. This first round of object identification eliminates objects smaller than the size percentile threshold. The output of the first round of object identification results in an image that is no longer binary. It also contains the labels for the identified objects.

We next convert an identified object to a binary image. This new binary image does not contain smaller-than-thesize-threshold objects and is less noisy than the original field. Next, using a 3D (x, y, t) centrosymmetric connectivity element, a second set of objects are identified and labeled such that the object label is the same for any two features with similar values that are either adjacent or diagonal to each other and that overlap in space or time. This results in a set of objects with a unique set of IDs that have been tracked sequentially through time. No temporal gaps are allowed. These are 3D objects that are connected in space and time (x, y, t).

No minimum percent overlap is enforced. We allow multiple objects that merge to have the same ID and a single object that splits into multiple objects to retain the ID of the initial object. As a result, any objects that have connectivity at some time in their evolution share an ID. This allows MHWs that are connected through time by the second round of object identification to contain spatially disjoint multiple objects with different labels under the first round FIG. 4. Sequence of morphological operations for closing (dilation I followed by erosion I) and then opening (erosion II followed by dilation II) using a structuring element with (a)-(e) a radius of four grid cells and (f)-(j) a radius of eight grid cells. Orange shading represents the feature area that the morphological operations are performed on. Red stippling in (e) and (j) shows the grid cells identified as potential MHWs before the morphological operations. Green contours outline the final shape of the identified MHW objects. Data shown here are from February 2011. Note again that the first row shows the sequence of morphological operations using R 5 4 grid cells, while the bottom row shows the sequence using R 5 8 grid cells.

of object identification. The connectivity established by Ocetrac is not mechanistic; it is entirely diagnostic and based on the existence of overlap in space and time. Hence, MHWs diagnosed by Ocetrac as a single object can result from a complex combination of mechanisms, which may or may not be coupled (supplemental material 4).

In summary, we describe a new tracking algorithm to detect and follow the evolution of MHWs. A sequence of morphological closing then opening is performed to create smoother and connected images. The results of this tracking algorithm will vary depending on the spatial and temporal resolution of the SST data as well as the choice of morphological radius R and minimum size percentile threshold P. In section 3, we discuss the sensitivities of these choices of R and P, along with useful metrics, including the number of MHW events detected, the average monthly duration of MHW events, the percentage of MHW events composed of multiple objects, and the percentage of the MHW area retained, for characterizing the global spatiotemporal evolution of MHWs. The percentage of the total MHW area retained is the sum of the accepted areas of all retained objects divided by the sum of areas initially identified as regions with MHW conditions. An area value is calculated for each object. The minimum area is calculated using the entire initial set of objects and the minimum size percentile threshold P.

Sensitivity analysis

The representation of MHWs is dependent on the criteria used to define their intensity, size, duration, and shape. This can be influenced by the horizontal and temporal resolution of the SST data, and whether or not the long-term trend in SST is removed. Here, we investigate the sensitivity of MHW metrics to the morphological radius R and minimum size percentile threshold P criteria implemented in Ocetrac. The term P is the percentile used to define the threshold of the smallest area object that is not discarded during object identification. The minimum area is calculated using the entire initial set of identified objects and the minimum size percentile threshold P. The minimum area value corresponds to the percentile of the area based on P. We quantify the effect of these criteria on the number of MHW events detected, average MHW duration, minimum MHW area, and the percent of MHWs with multiple centroids.

As R and P increase, fewer MHWs are detected (Fig. 5a). Large values of R increase the connectedness of features in the binary images, which results in larger objects. A MHW event is composed of one or more objects. Larger R values thus result in fewer but larger MHW events. These well-connected MHWs are also likely to persist for longer than 3 months (Fig. 5e). Each object has a centroid. As such, the number of centroids is also the number of objects making up a MHW. The percentage of MHWs with multiple centroids (in other words, multiple objects making up the MHW) decreases with increasing R (Fig. 5d). Fewer MHWs have multiple centroids when R is larger as a result of increased connectivity among features.

The average monthly duration of MHWs initially increases with R and P for values of P , 70% (Fig. 5b); however, for large R, the average monthly duration peaks for P near 75%. This nonlinear behavior is the result of the decline in the number of MHWs detected as the minimum size percentile increases. A smaller population size decreases the average duration at larger R values (Figs. 5b,e). Duration is most sensitive Recall that R is measured in gridcell units, so R 5 3 grid cells is equivalent to 0.758 longitude 3 0.758 latitude. The value R 5 10 grid cells is equivalent to 2.58 longitude 3 2.58 latitude.

to smoothing radius R, where large radii increase connectivity between neighboring features allowing MHWs to persist for longer periods of time.

Large minimum size percentile thresholds P reduce the percentage of the total MHW area retained (Fig. 5f). Smaller values of P result in a greater percent of the original MHW area retained that results in more MHWs of smaller size (Figs. 5a,c,f). As the size percentile threshold increases, the percent of total MHW area retained quickly declines to less than 50% (Fig. 5f). As P increases, the number of MHWs detected declines with the smallest size events increasing in size (Figs. 5a,c). If P is held constant, the percent of the total MHW area retained (Fig. 5f) decreases and the minimum MHW area (Fig. 5c) increases with increasing smoothing radius R. The larger smoothing radii act to join neighboring features and fill holes within feature clusters. Thus, a large smoothing radii creates larger MHWs, while also decreasing the total number of MHWs detected.

For a demonstration of the sensitivity of an example MHW to the smoothing radius and size percentile threshold, we examine the sensitivity of the 2011 MHW off Western Australia (Fig. 6). The shape and size of the detected objects are noticeably different between radii of 4 and 8 gridcell units, and the results are independent of area threshold P. A smoothing radius of 4 gridcell units produces objects with sharp and jagged edges and interior holes (Figs. 6a,d,g). The object shape difference between an R of 8 and 10 grid cell units is nearly negligible, with the exception of small features disappearing (e.g., Fig. 6b versus Fig. 6c). For example, for an R of 4 gridcell units and an area size percentile threshold P of 65, 5 objects are detected (Fig. 6a). For the same radius but with an area size percentile threshold of 75, 5 objects are still detected (Fig. 6b). For the same radius but with an area size percentile threshold of 90, however, only 3 objects are detected (Fig. 6c). For a radius of 8 grid cell units and area size percentile thresholds of 65, 75, and 90, 4, 3, and 3 objects are detected, FIG. 6. Sensitivity of objects detected from the morphological operations in February 2011. Each panel represents a unique combination of radius and minimum size percentile threshold from 4 to 10 grid cells and 65th-90th percentiles, respectively. Detected objects are outlined in green, brown stippling indicates the grid points where the SST exceeds the 90th percentile, and orange shading represents the filled-in MHW regions to create closed contour objects outlined in green. The OISST grid cell is 0.258 longitude 3 0.258 latitude. Note that R 5 4 gridcell units is equivalent to 1.08 longitude 3 1.08 latitude. The value R 5 8 gridcell units is equivalent to 2.08 longitude 3 2.08 latitude, and R 5 10 gridcell units is equivalent to 2.58 longitude 3 2.58 latitude. respectively (Figs. 6b,e,h). For a radius of 10 gridcell units and area size percentile thresholds of 65, 75, and 90, the number of objects detected is 3, 3, and 2 (Figs. 6c,f,i). As the minimum size percentile threshold P increases, objects disappear when the areas fall below the threshold. The sensitivities to the radius R and size parameters P give insight into how the set of detected MHW events can be impacted by choices made in the application of Ocetrac. Here, we use a radius of 8 as it provides enough detail of the original objects while creating smooth edges. We also choose the 75th percentile for the minimum size percentile threshold as it isolates the well-known MHWs that have occurred in the twenty-first century, including the event of Western Australia in 2011 (Fig. 6e). Different choices of radius and size parameters can be made and may be more suited for tracking specific events.

The sensitivity analysis reveals that the choice of parameters R and P influences the basic characteristics of MHWs such as the number, duration, and size. Our choices of radius and minimum area size percentile threshold used in this study were determined by our aim to have approximately 20 MHWs per year (approx. 800 from 1982 to 2020), a minimum area roughly the size of Alaska (approximately 2 3 10 6 km 2 ), and lasting on average 3 months, allowing us to focus on the more extensive and long-lasting MHWs (Holbrook et al. 2019).

Metrics

Ocetrac allows for the characterization of discrete MHWs in time and space. We define a set of measures that are computed over the lifetime of each event at monthly increments (Table 2). To describe the intensity within the MHW, we use the entire SST a field within the object contour (green outlines in Fig. 6) to calculate the mean, maximum, and cumulative intensity. The MHW anomalies are summed over the area and duration of the event to calculate the cumulative intensity. Degree heating weeks (8C weeks) are commonly used to study the impacts of coral bleaching in tropical reef ecosystems (Kayanne 2017;Eakin et al. 2010). The cumulative intensity (8C km 2 months) provides a measure of accumulated heating over the lifetime of the MHW and can be informative when assessing the time, space, and temperature dependence of ecological impacts related to MHWs.

Area is an important descriptor of MHWs. The total area is defined as the sum of gridbox areas contained within each object and takes into consideration grid resolution and latitude. Since the MHW with multiple objects can contain several centroids, we also compute the area for each object within the MHW. Given that MHWs evolve in space over their lifetime, it is informative to find the total MHW area as the sum of unique grid points contained within the MHW over its duration. The mean and maximum areas over time are also computed for each MHW.

The distributions of MHW duration and area are heavytailed, meaning that short-lived or small-area events occur more frequently than long-lasting or large area events (Fig. 7). The largest MHW encompassed the 2013-17 northeast (NE) Pacific "The Blob," impacting a total area of 1.94 3 10 8 km 2 and persisting for 60 months. The MHW off Western Australia has a total area and duration covering 1.53 3 10 8 km 2 for 47 months (Table 3). The Gulf of Maine and Mediterranean Sea MHWs were closer to the global average duration (2.99 months) and average total area (5.74 3 10 6 km 2 ) of all 813 MHWs detected from September 1981 through January 2021.

The maximum MHW intensity has a positively skewed distribution with a mean of 2.68C, maximum of 9.18C, and minimum of 0.28C (Fig. 7). The 2013-17 northeast Pacific The Blob had a maximum SST a of 7.138C, which is larger than the 2009-11 Western Australia (6.08C), 2012 Gulf of Maine (5.88C), and 2003 Mediterranean Sea (3.68C) MHWs, although the maximum intensities of all four MHWs were above average (Fig. 7a, Table 3).

The measures presented in Table 2 are useful to describe MHWs and characterize their evolutions in both time and space. In the following section, we use Ocetrac to detect and follow four well-known MHWs occurring during the twentyfirst century, including the 2013-17 northeast Pacific (Bond et al. 2015;Di Lorenzo and Mantua 2016), 2009-11 Western Australia (Pearce and Feng 2013), 2012 Gulf of Maine (Mills et al. 2013), and 2003Mediterranean Sea MHWs (Black et al. 2004;Sparnocchia et al. 2006).

Case studies

Ocetrac provides a global dataset of MHW spatiotemporal metrics that can be used to explore how past events evolved (Table 3). We provide a summary of how Ocetrac defines a MHW event. Ocetrac first identifies grid cells that are in MHW conditions at each time step and groups nearby cells into spatially distinct objects, each assigned a unique ID. However, at this stage, the algorithm does not recognize whether these objects represent a single, continuous MHW event across multiple time steps. In the next step, a 3D connectivity element is used to determine if these objects are part of the same event.

This method can result in discontiguous objects that are linked together as being part of a single MHW event. This highlights a limitation of using Ocetrac for predictive studies, as the full MHW event is only identified retrospectively.

Here, we explore these recent events and determine 1) if their representation using Ocetrac is consistent with the past literature and 2) if there is anything new that can be learned about MHWs by taking into consideration their spatial and temporal connectivity. We focus on four events that had major impacts on both socioeconomic and ecological systems and that sample from unique geographic regions in both the tropics and midlatitudes. The way these case study MHW events were identified was first by using Ocetrac to track MHW events using the entire data period. We looked for previously identified MHWs in the Ocetrac set of MHWs. For instance, Ocetrac identifies the Gulf of Maine MHW and the Mediterranean Sea events, which are regionally confined. We also examined the Ocetrac results for the 2010/11 event off the west coast of Australia as defined in the literature and found that event in Ocetrac. However, we found the Ocetrac-defined MHW event to be longer and larger. Similarly, using the set of Ocetrac-defined events, we found The Blob to be part of a larger and longer MHW event. 7. Distribution of (a) maximum intensity (mean 5 2.68C, minimum 5 0.28C, maximum 5 9.18C), (b) duration (mean 5 2.99 months, minimum 5 1 month, maximum 5 60 months), and (c) total area (mean 5 5.74 3 10 6 km 2 , minimum 5 6.87 3 10 5 km 2 , maximum 5 1.94 3 10 8 km 2 ) for 813 MHWs detected between September 1981 and January 2021. Named MHWs are indicated by the colored dots using definitions in Table 3.

a. Northeast Pacific

A MHW, colloquially referred to as The Blob, in the northeast Pacific was notorious for its unusually large scale, its persistence, and the magnitude of its temperature anomaly (Bond et al. 2015). In the literature, the MHW anomalies that developed in late 2013 were connected to the warm SSTs in the western tropical Pacific months prior through the excitement of atmospheric Rossby waves that weakened the mean state of atmospheric circulation over the North Pacific (Capotondi et al. 2019;Hartmann 2015;Shi et al. 2019). This resulted in an exceptionally high ridge of atmospheric pressure through the winter of 2014 that weakened surface wind speeds, lowered the rates of turbulent heat loss from the ocean to the atmosphere, and reduced the normal Ekman transport of cold water from the north (Bond et al. 2015). Offshore SST anomalies formed during the boreal winter of 2013/14. The NE Pacific phase of the event was initiated during atmospheric conditions resembling the North Pacific Oscillation (NPO), which were linked to the development of warm conditions in the western-central equatorial Pacific, in turn resulting in atmospheric teleconnections that produced a deepened Aleutian low. A deepened Aleutian low is associated with more southerly winds along the U.S. West Coast, reducing the strength of the northerly climatological winds. This resulted in weakened offshore Ekman transport and warmer temperatures along the coast, as suggested by Di Lorenzo and Mantua (2016). In addition, atmospheric anomalies in the northeast Pacific in 2013-15 were likely dynamically linked through atmospheric variability and thermodynamic coupling in addition to North Pacific decadal SST variability (Tseng et al. 2017;Di Lorenzo and Mantua 2016;Lee et al. 2015). The link between the initial NPO-like sea level pressure anomalies in the North Pacific and ENSO development in the equatorial Pacific has also been recently examined by Capotondi et al. (2022) who discuss a decadal dynamical mode of variability, termed the North Pacific-central Pacific ("NP-CP") mode, which is a component of the Pacific decadal oscillation (PDO) according to Newman et al. (2016). This mode of variability, which is associated with SST precursors of Blob-type events in the northeast Pacific, includes the development of ENSO precursors like the North Pacific meridional mode, favoring the initiation of central equatorial Pacific warming, which in turn increases the persistence of northeast Pacific MHWs.

We use Ocetrac to explore the spatiotemporal connectivity of anomalies during 2013-18 that contains The Blob (Fig. 8a, supplemental material 5). The first signature of the 2013-18 event occurs just south of the Gulf of Alaska as described by Bond et al. (2015) in late 2013. The SST anomaly was confined to the western and northeast Pacific through late 2014. At the same time, SST anomalies in the Indian Ocean were above average for most of 2014, which was a factor in the failed development of a major El Niño event in 2014/15 (Dong and McPhaden 2018;McPhaden 2015). The warm background SSTs likely enabled the MHW to grow in the Indian Ocean and persist through 2015. As defined by Ocetrac, in this MHW event, SST a in the Indian Ocean is linked with SST a in the North Pacific. The North Pacific portion of this mega MHW resembled the spatial pattern of the positive PDO in winter 2015 that extended from the Gulf of Alaska to the central/eastern tropical Pacific. Di Lorenzo and Mantua (2016) showed that the weak El Niño of 2014/15 provided the Aleutian low with enough variability to drive this PDO-like expression of SST anomalies. This variability, along with increased heat content in the tropical Pacific, was an important precursor to the development of one of the most powerful El Niño events on record in 2015/16. Individual snapshots of the monthly evolution of the objects contained within this event demonstrate its global reach (supplemental material 6).

b. Gulf of Maine

The Gulf of Maine MHW in 2012 covered the ocean from Cape Hatteras, North Carolina, to Iceland and in the Labrador Sea (Fig. 8b; Mills et al. 2013). A northward meridional shift in the atmospheric jet stream over North America during the late autumn and early winters of 2011/12 stabilized atmospheric high pressure over the western North Atlantic (Chen et al. 2014). This led to an overall reduction in surface wind speeds and higher than normal air humidity and temperature, which acted to inhibit turbulent heat loss from the ocean to the atmosphere and increase water column stratification (Chen et al. 2014). As a result, SSTs systematically warmed over the continental shelf from November 2011 through at least June 2012 (Chen et al. 2014). Anomalous warming in the spring of 2012 was attributed to large-scale atmospheric variability during the winter of 2011/12, whereas local advective heat flux played a secondary role to cool SSTs (Chen et al. 2014(Chen et al. , 2015)).

The results from Ocetrac show that the Gulf of Maine MHW was a regional event that was confined to the northwest Atlantic (Fig. 9; supplemental material 6). The center of action was centered offshore of Newfoundland with maximum cumulative intensities occurring in the Gulf of Maine, Gulf of St. Lawrence, and part of the Labrador Sea (Fig. 7b). The MHW, which began in April 2012, persisted for 9 months and covered a total ocean area of 7.91 3 10 6 km 2 with a maximum intensity of 5.828C (Table 3). Scannell et al. (2016) also tracked the 2012 Gulf of Maine MHW using 28 latitude 3 28 longitude resolution monthly detrended SST for 3 months, between June and August 2012, and found its area to be 7.60 3 10 6 km 2 with a maximum intensity exceeding 38C. Scannell et al. (2016) used the same detrending method, and warming events were detected by mapping contours of a particular anomaly level, with the default threshold being 1 standard deviation above the regional mean. They also defined the threshold in terms of standard deviation units. Scannell et al. (2016) defined MHWs based on their probability of occurrence. Here, we define a minimum size criterion, which they did not enforce. Scannell et al. (2016) also showed that the likelihood of a MHW of this size is enhanced during the negative phase of the North Atlantic Oscillation (NAO) and positive phase of the Atlantic multidecadal oscillation (AMO), with the AMO being more dominant. The AMO had been positive since the early 1990s, and the NAO took a negative excursion in 2012. The resulting relationship between natural modes of SST variability and MHW size may have favored the large-scale nature of the 2012 warm anomalies.

c. West coast of Australia

A major, unprecedented MHW occurred in late February 2011 off the coast of Western Australia (Pearce and Feng 2013). An important driver of this MHW was the fast phase transition from central Pacific El Niño in 2009/10 to La Niña in 2010/11 that was in part driven by strong easterly wind stress over Indonesia and the western end of the Pacific caused by warm SSTs in the Indian Ocean (Kim et al. 2011). Easterly wind anomalies in the western Tropical Pacific and over Indonesia excited an eastward upwelling oceanic equatorial Kelvin wave that quickly terminated warming associated with El Niño in 2009/10 (Kim et al. 2011;Kug and Kang 2006;Yoo et al. 2010). An extraordinary La Niña quickly ensued, which increased SSTs and sea level in the western tropical Pacific and off the northwest coast of Australia. High steric height anomalies forced a stronger than normal poleward flowing Leeuwin Current (Feng et al. 2013). In addition, northerly wind anomalies associated with low sea level pressure anomalies off of the coast of Western Australia helped to intensify the Leeuwin Current and reduce turbulent heat loss from the ocean (Feng et al. 2013). The poleward advection of warm water contributed to two-thirds of the warming, while positive air-sea heat flux anomalies into the ocean accounted for approximately the other one-third of the warming (Benthuysen et al. 2020). The anomalous air-sea heat flux in February 2011 acted to reinforce the MHW rather than damp the warming effects from La Niña (Feng et al. 2013). The exceptional MHW that resulted along Australia's western coast was dubbed "Ningaloo Niño" for its resemblance to other coupled ocean-atmosphere phenomena in the Pacific (El Niño) and Atlantic (Benguela Niño) (Feng et al. 2013). After the peak warming in March 2011 along the coast, positive sea level and SST anomalies propagated offshore following the propagation of mesoscale eddies (Benthuysen et al. 2014).

Indian Ocean SSTs during the following summers of 2012 and 2013 remained anomalously warm off Western Australia (Caputi et al. 2014) The persistence of anomalies was part of an increasing trend of Ningaloo Niño conditions since the late 1990s, where the shift to a negative PDO phase occurred with the strong La Niña that followed the 1997/98 El Niño (Feng et al. 2013;Sen Gupta and McNeil 2012). The trend was driven in part by a change to the negative phase of the interdecadal Pacific oscillation (IPO) and enhanced ENSO variance, and the former sustains positive heat content anomalies off Western Australia and favors cyclonic wind anomalies that reduce the prevailing alongshore southerly winds and enhance poleward heat transport by the Leeuwin Current (Feng et al. 2013). Further coupling between the alongshore winds and coastal SST has been shown to amplify Ningaloo Niño events (Kataoka et al. 2014).

Ocetrac identifies the initiation of a MHW off the west coast of Australia in 2008. It was detected later near the east coast of Australia in April 2010 and then off the west coast of Australia in October 2010 (supplemental materials 8 and 9). This earlier and large spatial footprint results from the connectivity established by Ocetrac via a physical overlap of warm SSTa * in space and time. Again, the connectivity between different components of this MHW, as detected by Ocetrac, is not necessarily mechanistic, although it may be worth exploring, in future studies, possible common causes.

d. Mediterranean Sea

During the summer of 2003, Western Europe experienced its worst heatwave in over 500 years, which caused excessive morbidity throughout the region, especially in hard hit France (Luterbacher et al. 2004;Valleron and Boumendil 2004). The extremely hot conditions over land from May through August stemmed from a persistent anticyclonic circulation centered over northern France that reduced cloud cover and precipitation (Black et al. 2004;Grazzini and Viterbo 2003). Although short-lived, the anomalous atmospheric anomalies quickly warmed SSTs in the central Mediterranean Sea in May before affecting the entire basin by July, with the exception of the Aegean Sea (Grazzini and Viterbo 2003). The Mediterranean Sea MHW warmed passively as a result of increased surface air temperatures, reduced surface wind speeds, and a reduction of all components of the upward heat flux (evaporation, longwave radiation, and sensible heat) (Olita et al. 2006). Upward heat fluxes increased in the winter and spring of 2003 and decreased in summer 2003. The MHW dissipated abruptly in late August to early September when strong westerly winds cooled surface air temperatures and induced wind-driven turbulent mixing that cooled SSTs (Sparnocchia et al. 2006).

Remote forcing from the northward shift and intensification of the intertropical convergence zone over West Africa, as well as Rossby waves emanating from tropical America that intensified the Azores anticyclone, contributed to the unusual atmospheric conditions driving the 2003 Mediterranean Sea MHW (Black et al. 2004). Decadal fluctuations in North Atlantic SSTs and the thermohaline circulation are known to influence European weather over long time scales. During 2003, the AMO index was positive and associated with elevated air temperatures and reduced wind stress over western Europe (Sutton and Hodson 2005).

The Mediterranean Sea MHW in Ocetrac during the summer of 2003 started in June and persisted through August (supplemental materials 10 and 11). Owing to the nature of the semienclosed region, MHW anomalies in the Mediterranean Sea did not connect with those in the Atlantic and had only one centroid per month. This meant that the MHW was highly localized with maximum anomalies over 48C and a total surface area of 1.78 3 10 6 km 2 , where the maximum cumulative anomalies occurred in the central and western regions of the basin (Table 3, Fig. 8d). The 2003 Mediterranean Sea MHW was the smallest size event of the four case studies examined here and however was intense enough to decimate rocky benthic macroinvertebrate species (Table 3; Garrabou et al. 2009).

Conclusions

Here, we discuss a MHW tracking algorithm called Ocetrac that can be used to characterize the spatiotemporal evolution of MHWs globally. This software tool Ocetrac highlights the spatial connectivity and temporal behavior of MHWs. Using Ocetrac, we characterize the spatial patterns and evolution of some of the most dangerous MHWs of the twenty-first century together with other MHWs identified from global monthly SST observations. A summary of our approach in this study is as follows:

1) Preprocess global SSTs to exclude the long-term warming trend and define anomalies SST a with respect to the local climatology. Anomalies are then standardized by the local monthly standard deviation of SST a over the entire climatological period. The climatological period is defined by a 40-yr period. The long-term warming trend is removed to examine MHW as consistently rare extreme events driven by internal climate variability while excluding the effects of global warming. 2) Extract the binary image of candidate MHW points, called features, where the standardized anomaly SST * a exceeds the 90th percentile. We then define the minimum radius R and minimum area size threshold P. Using morphological operations, opening and closing, the spatial structure is simplified giving SST * a features. We then identify MHW objects as features that are larger than the minimum area size threshold P.

3) Track MHWs in space and time using a 3D (x, y, t) centrosymmetric connectivity element, keeping track of multiple objects of a given MHW as the MHW splits and merges.

Other MHW tracking schemes exist, including the tracking scheme presented in Sun et al. (2023a), which was later extended to track MHW volumes (Sun et al. 2023b). The tracking scheme presented in Sun et al. (2023a) constructs MHW snapshots and tracks MHWs in the time domain, where candidate MHW points are detected using the Hobday et al. (2016) method. Sun et al. (2023a) used daily SST and found the raw MHWs to be spatially incoherent and addresses this by applying a K-nearest neighbor (KNN) algorithm, using a great circle distance to identify the K-nearest grid cells and filtering out grid cells where less than half of its K-nearest grid cells are identified as MHWs. This smoothing procedure has a similar impact as the Ocetrac tracking algorithm procedure of using morphological operations and a structuring element. There are several ways to treat the splitting and merging of MHW events in the time domain. The time domain linkage is where the method presented in Sun et al. (2023a) and Ocetrac diverge. In Sun et al. (2023a), the time domain linkage is determined by the fraction of overlapped domain between MHW units at two consecutive time steps, which is a user-defined parameter a. The results of this tracking algorithm depend on their K-parameter and to a smaller extent a (Sun et al. 2023a). In Ocetrac, no minimum overlap in the time domain is enforced. The labels of the objects in each MHW during the first stage of object identification, which is in the spatial domain, can be extracted.

We demonstrate the usefulness of Ocetrac in following the evolutions of four well-known MHWs in the Pacific, Indian, and Atlantic Oceans and Mediterranean Sea. The advantage of using Ocetrac globally, rather than within a specified region, is that it captures the large-scale and dynamically linked connections between remote SST anomalies that connect seemingly disconnected MHWs. In combination with dynamical studies, Ocetrac can provide a tool to better understand the origin of MHWs and their evolution. However, it is important to note again that the results of the tracking algorithm are sensitive to the choice of the radius R and minimum size percentile threshold P, which are currently user-defined parameters. In the case of the Ocetrac-tracked NE Pacific heatwave, the globally distributed objects outside of the North Pacific are identified as linked to the NE Pacific heatwave. In this example, the user could consider changing P and R and the structure of the connectivity element, which is currently user-defined. This flexibility could be addressed in a future version of the tracking software package. A future version of the tracking software package will also incorporate a submodule to calculate the measures of the characteristics of MHWs, such as the 90th percentile intensity and 90th percentile area. Other measures that will be published in the next stage of the package include the center of mass and centroid displacement, ratio of convex hull area and area of heatwaves, and deformation among others could be developed as measures to understand the compactness of the shape and dispersion and the change of these characteristics over time.

To a large extent, our interpretation of extreme events is dependent on how thresholds are defined. In many circumstances, extreme events are determined based on the space and time scales of their impacts and associated risks. For example, extreme flooding events are often classified by their extent and frequency in terms of their potential for damage (ten Veldhuis 2011). It is therefore useful to consider MHWs as temperature variance outside the normal range of thermal tolerance to native species. However, here, we remove the long-term 40-yr warming trend in order to better isolate the behavior of SST variance above the trend to be able to describe the spatiotemporal connectedness of MHWs and because we define MHWs as consistently rare extreme events that are driven by atmospheric variability. The 40-yr trend also removes some multidecadal variability. There are other applications where the inclusion of the trend would be appropriate as well as different ways to define the warming trend, such as a linear trend. However, global warming is nonlinear, so removing a linear trend may not be the best way to remove the forced response to global warming. When we retain the long-term warming trend, a greater proportion of ocean surface area experiences a MHW, which would increase the MHW intensity, duration, and size. If we were to use this tracking algorithm without detrending, we would identify and track only a few long-lasting and spatially expansive MHWs toward the end of the record. This is because this tracking algorithm relies on the binary maps that indicate candidate MHW points defined by the user-defined baseline and userdefined threshold. The use of this tracking algorithm can then be generalized to any baseline and any threshold that produce binary maps that indicate the candidate points of MHWs. The use case of this tracking algorithm is also not limited to just tracking MHWs but can also be extended to track other phenomena that could have coherent spatial structure, for instance, extremes in local sea level.

We also explore the sensitivity of Ocetrac to the resolution of gridded observational data, ranging from eddy-permitting (0.258) to very coarse (28). The overall large-scale spatial patterns agree well among the different resolutions. We expect this by manipulating R and P as demonstrated in Figs. 4a-d that the tracking algorithm could be implemented for daily (and even finer temporal scale) data. Morphological smoothing removes the artifacts of image noise and is completed by two grayscale morphological operations: grayscale erosion and dilation. In the case of daily MHWs, we anticipate that the morphological operations would need to be applied multiple times to be able to track the objects. However, visualizing and quantifying the spatiotemporal connectivity of MHWs in sea surface temperature forecasts using Ocetrac could enhance the usability of sea surface temperature forecasts.