Introduction

The oceans are taking up more than 90% of the excess heat in the climate system and are experiencing a long-term warming trend1,2. Many other changes are related to this ocean warming, including rising sea level3, increased upper ocean stratification4,5, and accelerating ocean circulation6,7. Climate change, climate modes, and the movement of atmospheric circulation cells can affect the development, dynamics, and location of ocean fronts8,9,10.

Ocean fronts separate water masses with different physical and biogeochemical properties, including salinity, temperature, nutrients, and biomass11. Fronts are found everywhere and can be identified in satellite sea surface temperature (SST), sea surface height, and other satellite or in-situ observations12. Fronts are mainly generated by convergence, confluence (when two or more flowing bodies of water join to form a single channel), rotation, and shear (when adjacent currents flow at different or opposite speeds)10. Many types of fronts have been observed and described across regions and processes that include tidal mixing, coastal and equatorial upwelling, boundary currents, subtropical convergence zones, and seasonal ice zones11.

Fronts are among the oceans’ most ecologically important structures and play an important role as principal biogeographical boundaries13,14. The stirring of large-scale currents and mesoscale eddies creates complex filaments of temperature, nutrients, and phytoplankton. These filaments become sharp fronts with strong submesoscale circulation15. By promoting the vertical movement of water masses, fronts facilitate the exchange of nutrients between subsurface and surface waters16,17. Nutrient upwelling enhances phytoplankton growth, forming the foundation of productive marine food webs12,18, while convergence and downwelling of water masses enhances carbon exports19. Studies have shown that predators, such as southern elephant seals in the Southern Ocean20, blue sharks in the Gulf Stream21, and loggerhead sea turtles in the southwestern Atlantic22, can track fronts and eddies as part of their foraging strategies. Therefore, trends in ocean fronts should be considered when investigating changes in ecosystems and commercially important fisheries23.

Several studies have examined trends of fronts in specific regions, such as the long-term increase of frontal frequency in the upwelling systems of California24,25, and central Chile26, and the short-term decreasing frequency of ocean fronts in the California Current System during the 2014/15 North Pacific marine heatwave27. Global studies of fronts suggest an increase in frontal abundance28. However, the changes in frontal activity in the context of warming oceans have not been explored globally. Here we focus on identifying commonalities in the trends in fronts and chlorophyll in ocean warming hotspots, distinguishing this study from other frontal trend studies.

We investigate the 2003–2020 trends of frontal activity in ocean warming hotspots, areas that experience warming greater than 1.48 °C per 100 years, equivalent to the highest 10% of warming values29. The trends are analyzed in the context of oceanographic dynamical regions30 (Fig. 1; see “Methods” section). We use satellite observations of chlorophyll (Chl) as a credible indicator of phytoplankton biomass and a precursor of primary productivity31,32. In this study, we: (1) compute trends in frontal characteristics and Chl in ocean warming hotspots; (2) classify global trends of frontal activity and Chl based on commonality across oceanographic dynamical regions (see Fig. 2 as an example for one hotspot), and (3) explore implications for global fisheries. This analysis is predicated on the assumption that ocean warming hotspots defined by Hobday and Pecl29 continue to exist as areas experiencing the most rapid warming during our study period of 2003 to 2020. We chose hotspots as our study regions because they act as windows into the future, making changes within hotspots useful for discussing potential upcoming changes to other parts of the ocean. However, warming is unlikely to be the sole cause of the observed trends. Our results indicate that there are fewer and weaker fronts associated with reduced phytoplankton in equatorial and subtropical gyre hotspots, while in high-latitude hotspots, stronger and more numerous fronts are associated with increased phytoplankton, and boundary current hotspots have mixed results.

Fig. 1: Locations and names of the twenty ocean warming hotspots analyzed in this study.
figure 1

The background and hotspot colors indicate the classification by oceanographic dynamical regions30: equatorial (red), subtropical gyre (yellow), and boundary currents (blue). Hotspots that are located north of 40°N and not part of any oceanographic dynamical regions were classified as ‘high-latitude’ (magenta). The Indo-China hotspot (orange) was classified as an equatorial and subtropical gyre hotspot as it overlaps with multiple regions, most of it is in equatorial and subtropical gyre regions. The Southern-Brazil Uruguay hotspot is divided into two parts: north (Southern-Brazil Uruguay-N) and south (Southern-Brazil Uruguay-S), which are assigned to different oceanographic dynamical regions. Hotspots shown in gray were not analyzed due to either a lack of data (gray solid) or because they were cooling over our study period (gray dashed; See “Methods” section and Supplementary Table 1).

Full size image
Fig. 2: Full Southern California-Baja hotspot to illustrate frontal detection and frontal characteristics.
figure 2

a snapshot of observed frontal positions (thin black lines) overlapping 8-day averaged sea surface temperature (SST) values for 16 Oct 2017 – 23 Oct 2017; b, d, f, h spatial patterns of trends in SST, frontal frequency (Ffreq), frontal strength (Fstre), and log-transformed chlorophyll (Chl) concentration; c, e, g, i time series averaged over the whole area in SST, Ffreq, Fstre, and Chl and associated linear trends. Linear trend values (per decade) are shown at the bottom-right corner of panels c, e, g, and i, where ‘*’ indicates statistical significance at the 95% confidence level. The gray shading in panels c, e, g, and i highlights the overlapping period (2003–2010), between our study and two previous ones.

Full size image

Results

Equatorial and subtropical gyre hotspots

Five of the six equatorial and subtropical gyre hotspots show concurrent declining trends of frontal activity: frontal frequency (Ffreq), frontal density (Fdens), and frontal strength (Fstre) (not all statistically significant). Trends of regional chlorophyll (RChl, mostly statistically significant) also demonstrate concurrent decline in these hotspots since 2003 (Fig. 3 and Table 1). The equatorial and subtropical gyre hotspots are composed of two equatorial hotspots (Galapagos and Indian Ocean), three subtropical gyre hotspots (Southern California-Baja, Southern-Brazil Uruguay, and Micronesia), and the Indo-China hotspot. The Indo-China hotspot overlaps multiple dynamical regions, but most of it is in the equatorial and subtropical gyre regions. We therefore consider it a hotspot in the equatorial and subtropical gyre region but do not classify it as “equatorial” or “subtropical gyre” separately.

Fig. 3: Equatorial and subtropical gyre hotspots: Linear trends of frontal frequency (Ffreq), frontal density (Fdens), frontal strength (Fstre), and regional log chlorophyll (RChl) anomaly.
figure 3

a Hotspots (colored lines) and the equatorial (red shading) and subtropical gyre (yellow shading) regions. Panels b through e: Anomaly trends of (b) Ffreq, (c) Fdens, (d) Fstre, and (e) RChl. In b–e, data share x-axis labels; standard errors are shown with black bars; statistically significant trends (above the 95% confidence level) are marked with an asterisk.

Full size image
Table 1 Ocean warming hotspot frontal characteristics.
Full size table

The northern part of Southern-Brazil Uruguay hotspot is the only equatorial and subtropical hotspot that experienced an increase in frontal activity (only Fstre trend is significant) matched by significantly increased RChl (Fig. 3 and Supplementary Table 1). As such, all equatorial and subtropical gyre hotspots demonstrate consistent changes of frontal activity and Chl, whether they are increasing or decreasing (Fig. 3).

High-latitude hotspots

Three of the five high-latitude hotspots show increasing trends of frontal activity (not all statistically significant) matched by increasing trends of RChl (all statistically significant; Fig. 4 and Table 1). The North Sea is the only high-latitude hotspot that has both decreasing trends for fronts (Ffreq and Fdens trends are statistically significant) and RChl (Fig. 4 and Table 1). The Sea of Okhotsk hotspot has contradicting trends across its frontal metrics (small and non-significant trends) and RChl (significantly increasing; Fig. 4 and Table 1).

Fig. 4: High-latitude hotspots: Linear trends of frontal frequency (Ffreq), frontal density (Fdens), frontal strength (Fstre), and regional log-transformed chlorophyll (RChl) anomaly.
figure 4

a Hotspots (colored lines). Panels b through e: Anomaly trends of (b) Ffreq, (c) Fdens, (d) Fstre, and (e) RChl. In b–e, data share x-axis labels; standard errors are shown with orange bars; statistically significant trends (above the 95% confidence level) are marked with an asterisk.

Full size image

Boundary current hotspots

Frontal trends in boundary current hotspots vary by region. Half the hotspots (South-Brazil Uruguay-S and Mozambique Channel) show increasing trends of frontal activity across all frontal metrics (not all statistically significant), while the other half (Southeast Australia and Southwest Australia) show decreasing frontal trends (not all statistically significant) over the 2003 to 2020 period of our study (Fig. 5 and Table 1). While all four boundary current hotspots have non-significant trends in Chl, three of the four hotspots show consistent changes in frontal activity and RChl. That is, both frontal activity and RChl increasing or both decreasing (Fig. 5 and Supplementary Table 1).

Fig. 5: Boundary current hotspots: Linear trends of frontal frequency (Ffreq), frontal density (Fdens), frontal strength (Fstre), and regional log chlorophyll (RChl) anomaly.
figure 5

a Hotspots (colored lines), and the “boundary current” (blue shading) regions. Panels b through e: Anomaly trends of (b) Ffreq, (c) Fdens, (d) Fstre, and (e) RChl. In b–e, data share x-axis labels; standard errors are shown with orange bars; statistically significant trends (above the 95% confidence level) are marked with an asterisk.

Full size image

Discussion

Decades of work on fronts has shown that they are critical dynamic features of the marine environment33,34. Many processes, like wind-driven vertical mixing, biogeochemical cycling, and carbon export, are enhanced at fronts13,33,34. Fronts are key features of the climate system and are often favorable environments for phytoplankton13,14. These concepts are consistent with our observations of higher Chl associated with increased frontal activity (frequency, density, and strength) and, conversely, lower Chl linked to decreased frontal activity. Our results in almost all hotspots (13 of 15), which cover many different types of frontal systems show such consistency across all our parameters (Ffreq, Fdens, Fstre, and RChl; Table 1).

Most equatorial and subtropical gyre ocean hotspots except Southern-Brazil Uruguay, have experienced a decline in the frequency, density, and strength of ocean fronts, accompanied by a decrease in RChl (not all significant; Fig. 3; Table 1). Ocean warming enhances stratification at large scales4, impeding the upward movement of colder, nutrient-rich waters from the deep ocean35. Increased stratification has been observed in most of our equatorial and subtropical gyre hotspots35. Furthermore, previous work indicates that there is a decrease in submesoscale activity and a reduction in vertical flux of heat at submesoscales associated with a warming of the upper ocean36. Together, this implies that sharp gradients at the surface may become less pronounced or disappear altogether without those potential drivers of heterogeneity, leading to fewer and weaker fronts and a homogenization of the ocean surface30. Decreased upwelling also leads to decreased nutrient availability and phytoplankton biomass37.

In high-latitude hotspots, our results indicate an increase in frequency, density, and strength of fronts, as well as an increase in RChl (not all significant; Fig. 4 and Table 1). Our current knowledge of frontal trends and links with the climate system at high latitudes is limited and a subject of ongoing research. Changes in nearby sea ice and associated changes in ocean stratification are potential drivers shared by most high-latitude hotspots. In the last 40 years, Arctic sea ice has decreased at a rate of 4.7% per decade38. Ocean fronts are often found near the ice edge39, separating cold polar waters from the warmer open ocean, and can remain like imprints in the ocean behind retreating sea ice. For example, a previous study found an increase in fronts in the Bering Sea across 1985-199640. By combining this frontal increase with our results (2003–2020), the Bering Sea may have already been experiencing increasing trends of frontal activity for forty years. These increasing trends might be one of the reasons for increased RChl and FChl in high-latitude regions. Less or thinner sea ice cover also leads to higher light levels in the underlying ocean41, stimulating the growth of phytoplankton41,42,43.

Our results for boundary current hotspots show a mix in the sign and significance of the frontal and Chl trends depending on the hotspot (Fig. 5). However, Ffreq, Fdens, Fstre, and RChl have matching trend signs in three of the four boundary current hotspots (Fig. 5). Fronts in boundary current regions are often associated with strong gradients between the boundary current and surrounding waters44. Consequently, changes in frontal activity of boundary current hotspots are likely related to changes in the boundary currents themselves. The increasing trends we observed for the Southern Brazil-Uruguay – S and Mozambique Channel hotspots may be correlated with the intensification of western boundary currents, which would strengthen the existing boundary current fronts and produce additional frontal features in surrounding waters. In contrast, the decreasing trends for the two hotspots near Australia may in part reveal weakening ocean currents. Understanding the specific relationship between local ocean currents and their associated fronts would require dynamic studies for each boundary current region.

Each hotspot has a unique dynamic environment, altogether covering a wide range of frontal systems. Some hotspots exhibit contrasting results to other hotpots in the same dynamic environment, such as the northern part of the Southern-Brazil Uruguay hotspot in the equatorial and subtropical gyre group and the North Sea hotspot in the high-latitude group (Table 1). The northern part of Southern-Brazil Uruguay hotspot is a long, narrow area along the coast. Here, increased alongshore winds45 may be a contributor to increased upwelling, frontal activity, and Chl. The North Sea is a shallow shelf sea with a broad connection to the ocean and substantial continental influences, such as tidal mixing, freshwater discharge, and heat flow, from northwestern Europe46. Detailed interpretation requires frontal research from a regional perspective, focusing on these regions with inconsistent trends.

A comparison of chlorophyll trends with multi-year averages shows that they exceed 5% of the mean per decade in six hotspots, all of which exhibit statistically significant RChl trends (Supplementary Table 1). Among them, the maximum rate of change was 51.2% per decade in the South-Brazil Uruguay – N hotspot (Supplementary Table 1). Our chlorophyll trends are based on satellite ocean color observations, which characterize the upper few tens of meters of the water column. In many tropical and subtropical areas, subsurface chlorophyll maxima are common, and may or may not be major contributors to total phytoplankton biomass and productivity. Our analysis is unable to say anything about sub-surface variability, but an expanded biogeochemical Argo network could address this.

To illustrate the potential socio-economic impact of the observed changes of fronts and Chl, we examined annual fisheries landings in ocean warming hotspots and compared them with global population density distribution (Fig. 6). The equatorial and subtropical gyre hotspots support fisheries with a total annual catch of around 5 million tonnes (Fig. 6 and Supplementary Fig. 1), of which the Indo-China hotspot alone accounts for more than 4 million tonnes. Our results show a decrease in frontal activity and RChl in equatorial and subtropical gyre hotspots. Declining phytoplankton biomass may explain part of the projected catch decline in the tropics47. The human population is more highly concentrated in mid- and low-latitude regions, especially between the equator and 45°N (Fig. 6b). A decline in fishery production could be challenging for mid- and low-latitude populations48, especially for the more than 3 billion people dependent on fish for their daily protein intake and micronutrients47.

Fig. 6: Annual global fisheries landings and global population density.
figure 6

a Aggregated annual fisheries landings by latitude; (b) Aggregated population density by latitude; (c) Map of annual global fisheries landings (2010–2015) and global population density (2020) with bounds of the studied hotspots: equatorial and subtropical gyre hotspots (red), high-latitude hotspots (magenta), and boundary current hotspots (blue).

Full size image

Some of the high latitude hotspots with increasing RChl, such as the Sea of Okhotsk and the Bering Sea, are well-known productive areas11. A previous study for the Sea of Okhotsk has shown that phytoplankton, zooplankton, and benthic biomass are concentrated at fronts49. Our five high-latitude hotspots support a total fish catch of around 3 million tonnes annually (Fig. 6; Supplementary Fig. 1). This is less than the equatorial and subtropical gyre hotspots, but the fish catch density in high-latitude hotspots is about double that of the tropics (0.995 tonnes per km2 compared to 0.482 tonnes per km2). Previous work has projected increases in fish catch potential at high latitudes under climate change, partly attributed to sea ice retreat, increased primary productivity, and expanding habitat for lower latitude species48. This may correlate with our observed increasing trends of frontal activity and Chl in most high-latitude hotspots.

The fish caught in our boundary current hotspots is approximately 0.2 million tonnes annually (Supplementary Fig. 1). Despite the fact that the total amount is small in comparison to other dynamical groups of hotspots, it is important locally. For example, the Southeast Australia hotspot, according to Australian Commonwealth fisheries, includes important species groups such as tuna, billfish, scalefish, and shark. In 2020-21, the Australia Eastern Tuna and Billfish Fishery, one of the largest Australian fisheries, accounted for 10% of total Commonwealth fishery production value, while the Southern and Eastern Scalefish and Shark Fishery (Commonwealth Trawl Sector) accounted for 19%50. Longer records of frontal activity and Chl are crucial for anticipating changes in the fisheries of Eastern Australia, given that the current Chl trends observed in the Southeast Australia hotspot are not statistically significant.

In previous studies, both Kahru, et al.24 and Obenour28 addressed trends in frontal frequency for the Ensenada region off the west coast of North America, comparable to our South California-Baja hotspot. For the 2003 to 2010 period for which there is an overlap between our study and Kahru, et al.24 and Obenour28 (Fig. 2), all three studies find that frontal frequency in the Ensenada region is increasing. According to Fig. 2, the increase appears to be caused primarily by a substantial jump in Ffreq in 2010, and the trend appears to have changed since then. However, this comparison is still limited due to the use of different datasets and study areas.

While we have identified cross-hotspot commonality in trends of frontal activity and Chl across dynamical regions, some of these trends are not statistically significant (e.g., RChl trends in boundary current hotspots). A longer time series would help ascertain whether all the observed trends are significant and persist or if some of these signals are linked to natural variability. The trends we observed in hotspots are a potential window into future large-scale changes in other areas of the ocean29. We acknowledge that ocean warming is not spatially uniform, and some regions are cooling on decadal time scales. However, based on recent IPCC data51, the global mean ocean warming rate is currently faster than the 1950-1999 rate on which our hotspots were defined. Warming drives changes in the mechanisms that generate fronts. We, therefore, expect similar but more widespread changes in frontal activity and Chl in decades to come.

This study has documented changes over the past two decades for frontal activity and Chl in ocean warming hotspots. The observed trends show a systematic change in frontal activity and Chl on a global scale. In equatorial and subtropical gyre hotspots, we observed a decline in frontal activity and Chl. In high-latitude hotspots, frontal activity and Chl have generally increased. We thus expect the equatorial and subtropical surface oceans to become more thermally uniform and poor in phytoplankton in the future. Conversely, high-latitude areas have the potential to be more productive through increased frontal dynamics and nutrient fluxes. Hotspots in boundary current regions show mixed results, some increasing and others decreasing. The commonality that we have found between frontal activity and Chl trends in ocean warming hotspots distinguishes this study from other studies reporting frontal trends. Considering the population and fisheries landings in and near ocean hotspots, better projections for coming decades of ocean frontal systems and impacts on phytoplankton are vital for local and global economies.

Methods

Data

To detect sea surface temperature (SST) fronts and to calculate frontal properties in ocean warming hotspots we use the National Aeronautics and Space Administration (NASA) Moderate Resolution Imaging Spectroradiometer (MODIS Aqua) spatially gridded (Level 3) infrared sea surface temperature (SST, °C) data product version R2019.0. 8-day averaged data are available at 4.63 km spatial resolution. The first full year of data was 2003 and as a result, the time frame we selected is the 18 years between 2003-01-01 and 2020-12-31, a total of 828 8-day maps. We applied the frontal detection method only in hotspots to reduce computation needs. To extract chlorophyll (Chl, mg m−3) concentration data in hotspots and along SST fronts, we used the NASA MODIS Aqua Level 3 Mapped Chlorophyll data product, version 2022, with the same spatial and temporal resolution as the SST data. This Chl dataset was not used for detecting chlorophyll fronts.

The MODIS retrieval algorithm tends to flag high-gradient regions as cloud-contaminated, potentially reducing the number of fronts found where the gradients are high. While this could give rise to a slightly lower frontal probability in some places, this is unlikely to impact trends derived from this study. Conversely, the decrease in strength of persistent fronts, e.g., the shoreward edge of the Gulf Stream, may increase the probability of finding fronts.

A global fisheries landings dataset (Global Fisheries Landings V4.0) published on the Institute for Marine and Antarctic Studies (IMAS) Metadata Catalogue and the Gridded Population of the World, Version 4 (GPWv4) dataset from NASA Socioeconomic Data and Applications Center (NASA SEDAC) was used to illustrate potential socio-economic impacts of frontal changes on fisheries.

Ocean warming hotspots definition and classification

We used ocean warming hotspots defined by Hobday and Pecl29 as our study regions. Their full definition of hotspots is that areas that experience warming rates over 1.48 °C per 100 years, equivalent to the highest 10% of warming values29. Hobday and Pecl29 defined hotspots based on 50 years (1950–1999) of historical sea surface temperature data and evaluated the persistence of these hotspots into the future (2000–2050) by using global climate models. Here, we presumed that these hotspots would continue to exist. We used the hotspot boundaries (shapefile) defined and provided by Hobday and Pecl29 instead of re-defining hotspots with MODIS data. We calculated SST trends for each hotspot in case the results were biased by some hotspots that have changed. Our calculation of SST trends verifies whether the defined hotspots are still warming, not whether they remain areas with the most rapid warming in the context of the global oceans. Three hotspots – Southeast Canada b, South Africa, and Kerguelen Island – were cooling regions during our study period and were not included in our discussion but only listed in Supplementary Table 1. These cooling regions are marked by gray dashed lines in Fig. 1. Some hotspots of Hobday and Pecl29 near the Arctic and elsewhere were not included in this study due to insufficient satellite data to identify frontal positions and frontal activity trends (gray solid lines in Fig. 1).

The classification of hotspots into regions was based on Martínez-Moreno, et al.30. They are “equatorial”, “subtropical gyre”, and “boundary current”. Ocean warming hotspots that are located north of 40°N and that are not part of other regions were classified as “high-latitude”. The specific classification is shown in Supplementary Table 1. As shown in Fig. 1, some hotspots were segmented according to the boundaries of several regions that they overlap. Some segments were eliminated from this study (marked by gray solid lines in Fig. 1), such as the coastal region of the South California-Baja and the southern part of the Galapagos, according to two rules. First, if the number of rows or columns of gridded SST data was smaller than 96, we removed the segments to make sure that the 32 × 32 overlapping window of the SIED frontal detection algorithm could be successfully applied. Second, if the proportion of zeroes in the Fdens (see Methods section: “Statistical metrics”) time series was greater than 50%, we discarded the segments due to the lack of data precluding the calculation of robust trends. In several segments, we found repeated occurrences where the frontal detection algorithm could not find any fronts, due to cloud-contaminated data during the study period. As a result, there are only a few non-zero values in the corresponding time series that can be used for trend analysis, which would further create unreliable trends and bias the results. Therefore, the second rule ensures that the trends on which our results are based are robust and reliable.

Frontal detection

We applied the Cayula and Cornillon52 single image edge detection (SIED) algorithm to the 8-day averaged SST data product to detect SST fronts. This algorithm is a histogram-based method that has commonly been used to detect oceanic fronts24,27,53. In this study, the SIED method was applied using the module “Oceanographic Analysis” embedded in the Marine Geospatial Ecology Tools (MGET) (see Code Availability). The basic principle of the SIED method is to define a front as the line that separates two populations of pixels53. The method searches for bimodality in histograms on overlapping windows of pixels in a smoothed image27. First, the SST products were smoothed with a 3 × 3 pixel (~12 km) median filter to reduce variance between pixels and generate robust fronts. Then, we used 32 × 32 overlapping windows for the histogram analysis. The method checks the window for a bimodal distribution in the pixel values (temperature). If there is a bimodal distribution in the detecting window, the mean values of the two populations (water masses) are calculated. The difference between the mean values is then compared with an input detection threshold. If the difference is larger than the threshold, the method then determines optimal values for edges (optimal temperatures for fronts). Then, a spatial cohesion algorithm was applied to confirm whether the pixels of the two populations were sufficiently spatially separated. Any pixels with these optimal temperatures in the overlapping window are marked as frontal pixels. Subsequently, a contour-following algorithm was also included in the original version of the SIED method, to extend fronts on either end so long as a sufficient temperature gradient continues in the direction the front is pointing. Unfortunately, the MGET module we used does not include this step. We experimented with variable thresholds and found the outputs not sensitive to small threshold modifications and a detection threshold of 0.3°C was finally applied. With this fixed detection threshold, we treated all fronts identified by the SIED method equally instead of distinguishing between different types of fronts. An example of fronts detected by the SIED method is shown in Fig. 2a.

Statistical metrics

We investigated a series of metrics to characterize the spatial and temporal variability (frontal density and frontal frequency), intensity (frontal strength), and regional and along-front Chl of SST fronts in the twenty ocean warming hotspots (Supplementary Table 1).

Spatial density of fronts: By applying the SIED method, images with frontal pixels (pixel value = 1) and non-frontal pixels (pixel value = 0) were generated for each SST image. Black lines in Fig. 2a show fronts formed by frontal pixels. Pixels covered by land or cloud were labeled as invalid pixels (pixel value = NaN). For each frontal image, the spatial density of fronts (Fdens) was defined as:

$${{{{{\rm{Fdens}}}}}}=\frac{{{{{{{\rm{n}}}}}}}_{{{{{{\rm{front}}}}}}}}{{{{{{{\rm{n}}}}}}}_{{{{{{\rm{valid}}}}}}}}$$
(1)

\({{{{{{\rm{F}}}}}}}_{{{{{{\rm{dens}}}}}}}\): the spatial proportion of fronts or frontal density

\({{{{{{\rm{n}}}}}}}_{{{{{{\rm{front}}}}}}}\): the number of frontal pixels (pixels with value 1) in a frontal image

\({{{{{{\rm{n}}}}}}}_{{{{{{\rm{valid}}}}}}}\): the number of valid pixels (pixels with values 0 and 1) in a frontal image

Frontal frequency (Ffreq) indicates the frequency of frontal detections for an individual pixel over a given time interval. For an individual pixel, frontal frequency was defined as:

$${{{{{\rm{Ffreq}}}}}}=\frac{{{{{{{\rm{N}}}}}}}_{{{{{{\rm{front}}}}}}}}{{{{{{{\rm{N}}}}}}}_{{{{{{\rm{valid}}}}}}}}$$
(2)

\({{{{{{\rm{F}}}}}}}_{{{{{{\rm{freq}}}}}}}\): the frequency of frontal detections for an individual pixel over a given time interval

\({{{{{{\rm{N}}}}}}}_{{{{{{\rm{front}}}}}}}\): the sum of the counts of front detections over a given time interval

\({{{{{{\rm{N}}}}}}}_{{{{{{\rm{valid}}}}}}}\): the sum of the counts of valid data over a given time interval

For the time interval, we experimented with one month and three months (seasonal). With the 8-day averaged SST data, using one month as the time interval leads to the problem of insufficient samples since the maximum of \({{{{{{\rm{N}}}}}}}_{{{{{{\rm{valid}}}}}}}\) is only four. Therefore, we used three months as the time interval for generating intra-annual seasonal time series of Ffreq. For three months, the maximum \({{{{{{\rm{N}}}}}}}_{{{{{{\rm{valid}}}}}}}\) value for an individual pixel is twelve. If we suppose an individual pixel has a valid SST value in ten of the twelve images and is defined as a frontal pixel four times, then its Ffreq will be 0.40 or 40%. Ffreq used in this study does not reflect frontal persistence. Ffreq indicates the frequency of frontal activity occurring at a particular location (whether or not each individual count is caused by the same front), while frontal persistence describes how long a particular front line exists (from its development to its dissipation).

Frontal strength (Fstre) or cross-front SST gradient: The application of the SIED frontal detection relies on searching for the bimodality of pixels’ histograms instead of the gradient’s magnitude. Thus, it is challenging to distinguish strong and weak fronts by adjusting the SIED threshold, at least compared to fronts identified with gradient-based methods. However, the SST gradient can still be used to define the strength of fronts. Based on the binary frontal images, the locations of frontal pixels are known. By retrieving the SST gradient with the coordinates of frontal pixels, we successfully embed the SST gradient on frontal images to derive \({{{{{{\rm{F}}}}}}}_{{{{{{\rm{stre}}}}}}}\). The time series of Fstre was then generated by spatially averaging Fstre values.

Frontal chlorophyll or cross-front chlorophyll (FChl) was defined similarly, by retrieving Chl values at each frontal pixel based on the coordinates of frontal pixels. Regional chlorophyll (RChl) is the mean chlorophyll concentration within each hotpot. Chl was log-transformed prior to calculating FChl and RChl.

By spatially averaging pixel values, we created monthly time series of SST, Fdens, FChl, and RChl and seasonal (three-month) time series of Ffreq. We then obtained monthly anomalies of SST, Fdens, FChl, and RChl, as well as seasonal anomalies of Ffreq by subtracting the mean climatological monthly (SST, Fdens, FChl, and RChl) or seasonal (Ffreq) value from the corresponding time series. Finally, we took the annual average to get the annual anomalies time series of each metric.

Trends and significance

Linear trends of annual anomalies (SST, Fdens, Ffreq, FChl, and RChl) were calculated using a linear least-squares regression. All the observed trends for SST, SST fronts, and Chl were assessed using a Theil–Sen estimator. A modified Mann-Kendall test54 was used to assess the statistical significance of trends. Uncertainties of trends were calculated by dividing the standard deviation of the time series by the square root of the effective sample size from the Mann–Kendall test.

To emphasize the potential importance of the observed chlorophyll trends, we converted them to a % change dividing these chlorophyll trends by the 2003–2020 chlorophyll multi-year averages. The results are expressed as a % change per decade.