Threat by marine heatwaves to adaptive large marine ecosystems in an eddy-resolving model.

Marine heatwaves (MHWs), episodic periods of abnormally high sea surface temperature (SST), severely affect marine ecosystems. Large Marine Ecosystems (LMEs) cover ~22% of the global ocean but account for 95% of global fisheries catches. Yet how climate change affects MHWs over LMEs remains unknown, because such LMEs are confined to the coast where low-resolution climate models are known to have biases. Here, using a high-resolution Earth system model and applying a "future threshold" that considers MHWs as anomalous warming above the long-term mean warming of SSTs, we find that future intensity and annual days of MHWs over majority of the LMEs remain higher than in the present-day climate. Better resolution of ocean mesoscale eddies enables simulation of more realistic MHWs than low-resolution models. These increases in MHWs under global warming poses a serious threat to LMEs, even if resident organisms could adapt fully to the long-term mean warming.

Introduction

The ocean has warmed significantly during the past few decades in most parts of the world[1]. With continuous ocean warming, prolonged extreme ocean warming events, known as marine heatwaves (MHWs), have occurred in many parts of the global ocean in the past decades[2-4]. Severe MHWs have caused negative impacts on marine ecosystems and fisheries[5-8], and the ecological responses to MHWs have been observed across a range of processes, scales, taxa and geographic regions[9]. MHWs have broader and more devastating ecological and socio-economic consequences than the impacts of long-term slower changes in the mean warming for which species might possibly adjust through adaptation[10]. Therefore, it is vital to investigate future changes in MHWs under global warming in order to develop potential mitigation strategies to reduce the overall ecological impact of climate change[11].

Both satellite and field observations of sea surface temperature (SST) have demonstrated that over the past few decades, MHWs have become longer-lasting, more frequent and extensive, primarily attributable to the increase in the mean warming of SST (refs. [3,11-13]). Both regional[14] and global[9,15] model simulations project MHWs to intensify and their incidence to increase under a warming climate. For instance, under the fossil fuel intensive scenario of Representative Concentration Pathway (RCP) 8.5 (ref. [16]), the majority of ocean areas is projected to experience almost permanent heatwaves with concomitantly stronger intensity by the end of the 21st century, with MHWs defined based on the conditions of the present climate[17], referred to as a mean warming-inclusive threshold. Similarly, many previous studies define MHWs relative to the mean climate over the historical period to investigate the changing characteristics of MHWs and their potential impact on marine life both in the past and in future[9,15]. In contrast, shifting the baseline temperature for future is useful to isolate the influence of mean background warming and higher moments[18] of temperature statistics on MHWs (ref. [17]), and a moving threshold is suggested to use to attribute MHW changes in the context of long term warming[19]. Therefore, estimating MHW changes using the mean warming-inclusive and future thresholds brackets scenarios are relevant to marine ecosystems with a range of capacity adapting to future warming.

Numerical models are important tools for elucidating the drivers and characteristics of MHWs, but the capability to reproduce MHWs in the historical record differs substantially among models at different resolutions. Low-resolution models, though computationally less intensive and useful for assessing the impact of climate change on MHWs at continental or global scales[17,20], do not resolve small scale physical processes including boundary currents and eddy transport processes[5] associated with MHWs (ref. [21,22]). High-resolution regional models with ~10 km ocean grid have much better fidelity in reproducing the magnitude and spatial structure of MHW events observed during the latter half of the 20th century[14,23]. In a comparison of global model simulations forced by atmospheric reanalysis at 1.0°, 0.25°, and 0.1° ocean grid spacing, the simulations at the 0.1° generally yield the realistic results in reproducing the frequency and duration of global MHWs during 1985-2017 (ref. [24]).

While previous studies have mostly focused on the global MHWs characteristics[3,9,17], there is an obvious increase in the frequency and duration of coastal MHWs from 1981 to 2016 based on four satellite datasets[12]; the largest impact on ecosystems is seen in LMEs found mainly in the coastal ocean[25]. The observed SST trends during the historical period in the LMEs are predominantly positive[26], and under a warming climate, the monthly SST warm extremes in LMEs over parts of the northern oceans depicted substantial increase as well using the Coupled Model Intercomparison Project Phase 5 (CMIP5) and Community Earth System Model large ensemble project (CESM-LENS) (ref. [27]).

However, ~1° resolution of the global models may not be able to resolve the crucial processes such as ocean eddies, coastal upwelling, stressing the need of higher resolution models with daily time scale to broadly investigate the changes in ocean variables such as SST due to climate change[27]. Further, there is a built-in increase in MHWs everywhere in the future ocean climate when using a mean warming-inclusive threshold, obscuring changes in characters that are specific to LMEs. Using climate simulations from a mesoscale-eddy-resolving ultra-high resolution Earth system model[28], we show an enhanced intensity and annual days of MHWs over most of the LMEs in future even using a “future threshold” above mean warming.

Need for high-resolution model simulation

Low-resolution (nominal ~ 1°) global models lack the capability of resolving small-scale processes such as boundary currents, coastal processes and ocean eddy fluxes[21] (the red circles in Fig. 1a), making them difficult to realistically simulate the characteristics of MHWs in LMEs (Fig. 1b) and their impact on marine species such as the Atlantic salmon (Salmo salar)[29,30]. In contrast to pelagic areas, the biodiversity of coastal areas is far more abundant[6], including many foundation species[7] as well as economically important fish[29,31,32] that are also vulnerable to MHWs. Observed impacts include coral bleaching[33,34], declining seagrass density[35], spawning reduction and distribution shift of marine fish[36,37].

Fig. 1

The physical processes driving MHWs and locations of LMEs.

a, Schematic diagram of processes represented in high-resolution models that allow the impact on biodiversity to be evaluated. The red, purple and brown circles indicate local, regional and teleconnection processes, with arrows illustrating the interactions between these processes and the ocean environment. b, The groups of LMEs (See “Large Marine Ecosystems (LMEs)” in Methods) based on continent including North America, South America, Europe, Africa, Asia and Australia, with a total of 54 LMEs used in this study.

With the ability to resolve small-scale processes (the red circles in Fig. 1a) and their connections to climate modes of variability (the purple and brown circles in Fig. 1a), high-resolution Community Earth System Model version 1.3 (CESM1.3) is used for simulation from 1850-2100 (See “Model descriptions” in Methods), providing valuable information for risk assessment and adaptation planning for coastal areas. The model is forced by historical forcings before 2005 and the RCP8.5 (ref. [16]) (a high-emission scenario) thereafter. Results from a low-resolution version of the same model and available CMIP5 models, which are also low-resolution, are used for a comparison.

Observed and simulated MHWs in the historical period

The frequency of MHWs (See “Definition of marine heatwaves (MHWs)” in Methods) is qualitatively similar among remotely sensed National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation Sea Surface Temperature (OISST; Fig. 2a), Group for High Resolution Sea Surface temperature (GHRSST) Multi-Product Ensemble (GMPE; Extended Data Fig. 1a) and the modeled SST from high (CESM-HR) and low (CESM-LR) resolution configuration of CESM1.3 (ref. [28]) and the CMIP5 ensemble (Fig. 2b,c,d; See “Satellite-based observational data set” and “Model descriptions” in Methods). During the historical period (1975-2004), there are one to three MHW events per year occurring over most of the globe. The obvious low frequency of MHWs in the eastern tropical Pacific is owing to El Niño-Southern Oscillation (ENSO) that can result in long period events that occur only every few years[17] (Fig. 2a). The spatial distributions of MHW frequency indicate CESM-HR (Fig. 2b) is closer to that of OISST (Fig. 2a) and GMPE (Extended Data Fig. 1a) than the low-resolution CESM-LR (Fig. 2c) and CMIP5 (Fig. 2d), as clearly delineated by the latitudinal zonal mean variations (Fig. 2e).

Fig. 2

Observed and simulated frequency of MHWs during the historical period.

Spatial distribution of annual MHW frequency based on OISST data (a), as well as the simulation outputs in CESM-HR (b), CESM-LR (c) and CMIP5 (d). e, Zonal mean MHW frequency for OISST (solid black), GMPE (dashed gray), CESM-HR (red), CESM-LR (blue) and CMIP5 (orange). For CMIP5, the yellow shading is added to represent one standard deviation calculated based on the 20 models (Supplementary Table 1). f, The box-and-whisker plot of MHW frequency grouped by continent for the LMEs, with the minimum and maximum (line end points), 25th and 75th percentile (boxes), medians (horizontal lines), and average (black points). Note that due to data availability, the OISST and GMPE data used in this study spans from 1982 to 2011, slightly different from the historical period of 1975-2004 used in climate modeling. The overlapping period of 1982-2011 that encompasses the model simulations and OISST or GMPE was also used in model evaluation, yielding similar results. Results show that CESM-HR is more realistic in reproducing the frequency of MHWs compared to the CESM-LR and CMIP5.

Extended Data Fig. 1

Observed frequency and mean intensity of MHWs.

Spatial distribution of annual MHW frequency (a) and mean intensity (b) of MHWs during 1982–2011 based on Group for High Resolution SST Multi–Product Ensemble (GMPE) (refs. [38,39]). Results show that one to three MHW events per year occurring over most of the globe, and large spatial heterogeneity in mean intensity with higher values over areas such as the western boundary current regions.

For instance, biases of MHW frequency in CESM-HR, CESM-LR and CMIP5 relative to OISST are -0.31, -0.55 and -0.60 (15%, 27% and 30%) times per year (Fig. 2e), respectively, although negative bias still exists across majority of the latitude bands in CESM-HR. The zonal mean of GMPE only includes the latitudes within 55 degrees of the equator due to the diminished consensus among the multiple data sets that are included in GMPE at high latitudes[38-40].

Below we classify the LMEs into 6 groups (Fig. 1b and See “Large Marine Ecosystems (LMEs)” in Methods), and show that CESM-HR (red) more closely captures the frequency of MHWs when compared to OISST/GMPE than its low-resolution counterparts: CESM-LR (blue) and the CMIP5 ensemble (orange in Fig. 2f) (P<0.05). The improved simulated global mean SST from CESM-HR relative to CESM-LR[28] is partly explained by the difference in computing eddy vertical heat transport, that is, it is explicitly computed in CESM-HR but parameterized in CESM-LR[41,42], and the better simulated mixed layer depth by CESM-HR[28].

The spatial (Fig. 3a-d; Extended Data Fig. 1b) and zonal (Fig. 3e) mean distributions of MHW mean intensity in CESM-HR are much closer to those of OISST and GMPE compared to CESM-LR and CMIP5, consistent with the comparison over the LMEs (Fig. 3f). It is noteworthy that large differences exist in the intensity over the tropics and subtropics between OISST and GMPE, with smaller model biases when benchmarked against GMPE (Fig. 3e). The mean intensity of MHWs also exhibits large spatial variations, with higher intensity occurring in the western boundary current (WBC) regions, the eastern and central equatorial Pacific boundary current regions and the eastern boundary current regions (Fig. 3a; Extended Data Fig. 1b). The intensity over the WBC regions in CESM-LR and CMIP5 is underestimated, while in CESM-HR it is overestimated (Fig. 3b-d), leading to higher peak of the intensity in LMEs (that is, over South America and Asia in Fig. 3f) compared to OISST and GMPE. Nevertheless, the comparable maximum intensity exhibited in GMPE and CESM-HR over the LMEs in North America strongly support improvements in CESM-HR relative to CESM-LR and CMIP5.

Fig. 3

Observed and simulated mean intensity of MHWs during the historical period.

Spatial distribution of annual MHW mean intensity based on OISST data (a), as well as the simulation outputs in CESM-HR (b), CESM-LR (c) and CMIP5 (d). e, Zonal mean MHW mean intensity for OISST (solid black), GMPE (dashed gray), CESM-HR (red), CESM-LR (blue) and CMIP5 (orange). For CMIP5, the yellow shading is added to represent one standard deviation calculated based on the 20 models (Supplementary Table 1). f, The box-and-whisker plot of MHW mean intensity grouped by continent for the LMEs, with the minimum and maximum (line end points), 25th and 75th percentile (boxes), medians (horizontal lines), and average (black points). Note that due to data availability, the OISST and GMPE data used in this study spans from 1982 to 2011, slightly different from the historical period of 1975-2004 used in climate modeling. The overlapping period of 1982-2011 that encompasses the model simulations and OISST or GMPE was also used in model evaluation, yielding similar results. Results show that simulated MHW intensity by CESM-HR is in general closer to that in OISST and GMPE than simulated by CESM-LR and CMIP5.

The apparently high MHW intensity in the eddy-rich WBC regions, that is, the Scotian Shelf and Patagonian Shelf close to the North and South America, respectively, has been reported previously for a 0.1° high resolution ocean only model[43], or coupled ocean and sea ice model[24], possibly linked to the stronger internal variability of SST exhibited in high resolution models[43]. However, the relatively coarser resolution (0.25°) of the satellite SST data (see “Satellite-based observational data set” in Methods) might lead to underestimation of SST variability at mesoscale eddy-resolving resolution[21] (that is, 0.1°).

In general, both in the LME regions and on a global scale, CESM-HR is more skillful in reproducing the frequency and mean intensity of MHWs compared with the coarse resolution models (CESM-LR and CMIP5), lending some confidence for the following analysis.

Changes in MHWs under a warming climate

We define future MHWs using the “mean warming-inclusive threshold” and “future threshold” determined from simulations of the future (See “Definition of marine heatwaves (MHWs)” in Methods) to isolate the effect of mean background warming and higher statistical moments of SST on MHWs. The number of annual MHW days is more highly correlated with changes of important basic biological groups in the world than the mean or maximum SST (ref. [6]). On the other hand, the intensity or the temperature anomaly during a MHW event can represent the level of acute heat stress for marine ecosystems and is closely linked to mortality of organisms such as intertidal barnacles[17].

Based on the mean warming-inclusive threshold, under RCP8.5, CESM-HR projects strong increases in annual MHW days (Fig. 4c) and average MHW intensity (Fig. 4d). Compared to the historical period, the mean annual MHW days between 70°N/S are projected to increase by 287.2 days during 2071-2100 (Fig. 4c), and a permanent MHW state will be reached in many areas of the equatorial and subtropical regions. Many marine species in the equatorial region live near their high temperature ceiling and are highly sensitive to MHWs (refs. [44,45]). In contrast, the increase in annual MHW days in the North Atlantic and WBC region is much smaller. The mean intensity shows a mean increase of 1.2 °C based on CESM-HR (2071-2100; Fig. 4d), Moreover, the increase of intensity is much larger over the northern hemisphere with faster mean SST warming[46] (Extended Data Fig. 2a).

Fig. 4

Projected changes in annual MHW days and mean MHW intensity.

Shown are from CESM-HR based on mean warming-inclusive threshold and future threshold. Projected changes in annual days (left) and mean intensity (right) of MHW in 2071-2100, based on future threshold (a,b) and mean warming-inclusive threshold (c,d). The areas surrounded by the black solid line and coastline represent the LMEs. Note that the color bar range is different in each plot. Results show mild increases based on future threshold in the annual MHW days and mean MHW intensity, while the changes are much stronger based on mean warming-inclusive threshold.

Extended Data Fig. 2

Future changes in SST and SST standard deviation.

Spatial distribution of future changes (2071–2100 minus 1975–2004) in SST (a) and detrended daily SST standard deviation (b) based on CESM–HR. Results show larger increase of SST over the northern hemisphere compared to the southern hemisphere, and dipole features over the WBC regions in the future changes of SST standard deviation, likely attributable to the shifts in the frontal position in WBC regions.

The mean warming dominates the changes in MHWs and explains 94% or more of the simulated changes, as inferred by the similarity between the changes directly estimated from the simulations (Fig. 4c, d) and the changes calculated based on the pseudo scenario with mean warming alone (Extended Data Fig. 3). The pseudo scenario is designed by adding a perturbation calculated based on the thirty-year mean SST differences between future (2071–2100) and historical period (1975–2004) to the historical daily SST. The dominant role of mean warming in enhancing the changes in MHWs is consistent with previous results showing that the changes in the mean SST was the primary driver of the changes in MHW globally[11] or in coastal regions[12] during the historical period.

Extended Data Fig. 3

Projected changes in MHW annual days and intensity due to mean warming.

The calculation is based on a pseudo scenario by adding a perturbation to the historical daily SST, with the perturbation equivalent to thirty–year mean SST differences between future (2071–2100) and historical (1975–2004) periods.

In contrast, by utilizing the future threshold, much smaller increases in annual MHW days and average MHW intensity are projected (Fig. 4a, b) by CESM-HR. The mean annual MHW days over 70°N/S are projected to increase by only 2.8 days during 2071-2100 compared to 1975-2004 (Fig. 4c), which is much lower than the result obtained using a mean warming-inclusive threshold. Besides annual days of MHWs, the mean intensity shows consistently mild increases over the oceans worldwide, with a mean increase of 0.2 °C (2071-2100) (Fig. 4b). Likewise, the application of similar methods yields comparably small changes in the annual MHW days over northeast Pacific Blob and MHW intensity over North Atlantic Ocean, respectively by the end of this century in RCP8.5 using coarse resolution simulations[13,22]. Moreover, the dipole feature of higher increase over northern hemisphere and lower increase over southern hemisphere exhibited from analysis using the mean warming-inclusive threshold disappears, resulting in more uniform increases in MHW intensity.

Analysis using CESM-LR and the CMIP5 multi-model ensemble in general supports the findings discussed above, further illustrating the dependence of MHW characteristics relative to the baseline climate (Extended Data Figs. 4, 5). However, CMIP5 and CESM-LR do not show WBC regions having distinct behavior relative to other mid-latitude regions (Extended Data Figs. 4b, 5b), while CESM-HR projects much larger changes over the major WBC areas including the Kuroshio Extension, the Gulf Stream, the Zapiola Anticyclone, the Agulhas Return Current, the East Australian Current and the South Pacific storm-track (Fig. 4b). Moreover, almost all these major WBC regions except the East Australian Current delineate a distinctive meridional dipole intensity changes with increase over the poleward flank and decrease over the equatorward flank. The dipole feature can be explained by changes in the detrended SST variance (Extended Data Fig. 2b), consistent with a previous study that demonstrated the relationship between SST variance and MHW intensity[13], with changes in SST variance and MHWs intensity likely attributable to the shifts in the frontal position in WBC regions[47]. Detrending SST prior to calculation of SST variance excludes the influence of greenhouse-induced long-term trends on SST variability[27]. The 30-year SST trends in the historical and future climate over LMEs are listed in Supplementary Table 2.

Extended Data Fig. 4

Projected changes in annual MHW days and mean MHW intensity.

Shown are from CESM–LR based on mean warming–inclusive threshold and future threshold. Projected changes in annual days (left) and mean intensity (right) of MHW in 2071–2100, based on future threshold (a,b) and mean warming–inclusive threshold (c,d). The areas surrounded by the black solid line and coastline represent the LMEs. Note that the color bar range is different in each plot. Results show much smaller increases in the annual MHW days and mean MHW intensity based on future threshold in comparison to those obtained from mean warming–inclusive threshold.

Extended Data Fig. 5

Projected changes in annual MHW days and mean MHW intensity.

Shown are from CMIP5 based on mean warming–inclusive threshold and future threshold. Projected changes in annual days (left) and mean intensity (right) of MHW in 2071–2100, based on future threshold (a,b) and mean warming-inclusive threshold (c,d). The areas surrounded by the black solid line and coastline represent the LMEs. Note that the color bar range is different in each plot.

Future changes of MHWs in LMEs with CESM-HR

Fisheries catch varies by two orders of magnitude among the LME regions, so it is useful to present the MHW days and intensity for the present and future for the different LME regions separately (Fig. 5). The mean and standard deviation of MHWs are shown in Supplementary Table 2. Historically, over the LME regions between 70°N/S, annual MHW days range from 27.4 to 40.6, with an average of 33.2 (x axis in Fig. 5a, c). With the mean warming-inclusive threshold, the annual MHW days over LMEs soar to 351.4 (y axis in Fig. 5c). The results using the future threshold largely suppresses the dominant effect of mean SST changes on the future MHW days. A total of 98% (except one) of LMEs show more MHWs days, with a mean annual increase of 2.8 days by the end of this century compared to the historical period (Fig. 5a). The increase in the mean annual MHW days is contributed by the increase in the persistence of MHWs, as indicated by the increase in the autocorrelation of SST despite of decreases in the frequency of MHWs[13].

Fig. 5

Comparison of MHW days and intensity between future and historical period.

Shown are for the LMEs based on mean warming-inclusive threshold and future threshold. The mean annual days (a,c) and intensity (b,d) among the LMEs for MHWs defined based on threshold for each period (1975-2004 and 2071-2100; a,b), respectively, and mean warming-inclusive threshold (c,d). The smaller squares and larger circles represent the small and big fishery catch, respectively, with color indicating the index of the LMEs. Results show an enhanced MHW intensity and annual days over most of the LMEs in future even with a “future threshold” above mean warming.

Consistently, using mean warming-inclusive thresholds, we find that the mean intensity in LMEs increases by more than 100%, from 1.2 °C during the historical period to 2.9 °C by the end of the century. As expected, the mean intensity over the LMEs based on future thresholds yields a small increase of 0.2 °C. Despite this, 93% of LMEs display an increase in intensity (Fig. 5b). The increase is primarily contributed by the changes in the SST variance, reflected by a statistically significant correlation (P<0.05) between MHW intensity and SST variance in both historical and future periods (Extended Data Fig. 6). Given the vast diversity of geographical locations of the LMEs, forcing of the increased SST variance is equally diverse. Under greenhouse warming, dominant modes of climate variability, which strongly influence MHWs across the global ocean, are generally projected to increase in their variance. As such, increased SST variance, hence an increased intensity of MHWs over majority of the LMEs, at least in part attributable to strengthened ENSO variability[48], enhanced SST variability over north tropical Atlantic[49], increased frequency of stronger positive Indian Ocean Dipole[50] and stronger nonlinear relationship between evaporation and SST over the North Pacific[51] under greenhouse warming.

Extended Data Fig. 6

Relationships between SST standard deviation and MHW mean intensity.

Both SST standard deviation and MHW mean intensity based on future threshold are averaged in each LME during historical period (1975–2004, a) and future period (2071–2100, b). The asterisk on the top left of the correlation coefficient R indicates statistical significance (P<0.05). Results show strong correlations between the MHW mean intensity and SST standard deviation.

Importantly, there is a significant correlation of 0.9 (P<0.05) between the future and historical mean intensity over all LMEs, compared to otherwise 0.6 based on mean warming-inclusive threshold (Fig. 5b, d), emphasizing comparable severity of MHWs at present and future for the majority of LMEs. In other words, LMEs that are under stress now will continue to be so in the future but in addition to the stress due to the mean warming and the increased intensity. The change from a more scattered distribution (Fig. 5d) to the alignment almost in a straight line (Fig. 5b) is to a large extent because of the built-in increase in MHW intensity due to the mean warming in Fig. 5d that dominates the response when the mean warming-inclusive threshold is used.

Our result based on the future threshold show that marine species in most of the LMEs would still experience an increase in the threat of MHWs, if they were able to adapt to the slowly increasing mean warming. To highlight this point, we compare result for LMEs of the category I, which are the 15 LMEs with the largest fishing capacity (Extended Data Fig. 7) and for all other LMEs defined as category II (Fig. 1A in ref. [52]). There is a generally larger changes in the mean annual MHW days and intensity in category I, compared to that of category II, indicative of a potentially more intensified impact of MHWs on LMEs with higher catches.

Extended Data Fig. 7

Comparison of MHW days and intensity between future and historical period.

Shown are for the top 15 catch LMEs based on mean warming–inclusive threshold and future threshold. The mean annual days (a,c) and intensity (b,d) for the top 15 catch LMEs (yellow dots, refer to category I) for MHWs defined based on threshold for each period (1975–2004 and 2071–2100; a,b), respectively, and mean warming–inclusive threshold (c,d). The yellow and blue triangles represent the average annual days or intensity of MHWs over the high (category I) and low catch (category II) areas, respectively. Results show larger changes in the mean annual MHW days and intensity in category I compared to that of category II, implicative of a more intensified impact of MHWs on LMEs with higher catches.

Organisms might adapt to climate change to a certain extent[44], but the rate of adaptation can vary widely among species[53,54], and spatial heterogeneity in the changes of SST might lead to differences in the extent to which marine species must adapt. Changes in MHWs defined using the future threshold are relevant if species can adapt fully to the future mean warming, which might not be possible[13] due to the rate at which SST is changing relative to what ecosystems have experienced in the past[55,56]. The increase in MHWs under the future threshold can be considered the “most optimistic” scenario for establishing the lower bound of climate change impact on marine ecosystems.

Conclusions

We find an increased intensity and annual days of MHWs over the majority of the LMEs in the future climate by applying a “future threshold” definition of MHWs. Our result of a widespread increase of MHWs over LMEs implies that even if we assume that organisms in LMEs were able to adapt fully to the impact of the long-term mean warming, the LMEs would still face serious threats under global warming. Our result is based on a high fidelity simulation using a high-resolution model that provides improved simulation of MHWs in the LME regions. As computational power continues to improve, we expect that a multi-model ensemble of high-resolution model simulations will soon be possible to project future MHW changes under multiple climate forcing scenarios, to assess the associated uncertainty, and to provide early warning of the likely changes. Importantly, our initial result indicates that even under the most optimistic assumption, risks to LMEs are substantial. The result therefore has far-reaching ecological, social, and economic implications and calls for a response strategy from the impacted communities and policy makers.

Methods

Model descriptions

Here, we use a high-resolution Earth System Model simulation, spanning 250 years from 1850 to 2100 (ref. [28,57]) that uses CMIP5 historical forcings until 2005 and the RCP8.5 (ref. [16]) (high-emission scenario) thereafter. The models in the high-resolution configuration of the CESM1.3 were used for the simulation. The atmospheric and land models have a nominal horizontal resolution of 0.25°, while a nominal horizontal resolution of 0.1° is used for the ocean and sea ice components. This high-resolution configuration of CESM1.3 is referred to as CESM-HR. At the resolutions of the individual components, the model allows for mesoscale eddies in the ocean to better delineate the interactions between the mesoscale phenomenon and large-scale circulation[28]. For comparison, simulations with CESM at a coarser spatial resolution of 1° in both atmosphere and ocean, referred to as CESM-LR, as well as the multi-model ensemble of 20 models participating in CMIP5 (ref. [16]) are also used in this study (Supplementary Table 1). All simulations during the historical period (1975-2004) for CMIP5, CESM-LR and CESM-HR, as well as the OISST/GMPE (See “Satellite-based observational data set” in Methods) for 1982-2011 are interpolated to 1°. Please note that the start time of OISST and GMPE is September 1981, thus January 1982 is selected as the start and December 2011 as the end for the thirty-year comparison. Comparison of the future and historical MHW based on CESM-HR is performed at the spatial resolution of 0.25°.

Large Marine Ecosystems (LMEs)

The LMEs (http://www.lme.noaa.gov/) mainly refer to the coastal areas and the outer edge of coastal currents, including river basins and estuaries up to the seaward boundary of the continental shelfs or well-defined systems of currents without continental shelfs[25]. A LME usually include an area of 200,000 square kilometers or more[25]. The LMEs are rich in biodiversity, including 95% of the global fish catch, though they cover only 22% of the total ocean area[52], providing goods and services to billions of people worth more than US$12.6 trillion annually[58]. In this study, the LMEs within 70°N/S are divided into 6 groups according to the adjoined continents (see Fig. 1b), North America, South America, Europe, Africa, Asia and Australia. Some LMEs might be located between two continents. For example, the Mediterranean Sea lies between Europe and Africa, but to simplify our analysis, it is considered part of Europe.

Definition of marine heatwaves (MHWs)

A MHW is a prolonged, discrete, anomalously warm water event[59]. Specifically, for each grid cell, a threshold for each day of a year is first determined based on the 90th percentile using daily-mean SSTs in the 11-day moving window centered on the specific day over a long, 30-year segment to ensure a sufficient sample size. This is followed by a 31-day moving average of the daily threshold. Five consecutive days or more with SST above the threshold is identified as a MHW event, and two events separated by an interval of two or fewer days are considered as one event. Note that a MHW defined above might not only occur in the warmer months, MHWs in colder months are also fatal for some creatures[60-62]. The number of days per event, denoted as duration, the number of annual MHW events (frequency) and the average intensity representing the mean deviation of SST from the climatological mean within the event are calculated first, and then the total number of annual MHW days, as well as the mean intensity are derived. Note that the threshold over the high latitudes is affected by the melting rate of ice and snow which is not taken into account; therefore, the analysis in this study primarily focuses on regions within 70°N/S which are less likely affected by ice and snow cover.

Two thresholds are used in this study to calculate the future (2071-2100) MHWs: (1) the 90th percentile over the historical period[24,43] (1975-2004), called mean warming-inclusive threshold, as defined above, (2) the future (2071-2100) 90th percentile[22], called future threshold, provides a delineation of the extent to which MHW changes are associated with the mean warming or non-seasonal temperature changes[19].

Satellite-based observational data set

To compare the MHW index calculated from the simulations, the satellite-based NOAA OISST V2.1, referred to as OISST (https://www.ncdc.noaa.gov/oisst/)[63,64], is used to calculate the MHW index globally in the study. The OISST product has been widely used in MHW studies[5,15,59]. This data set was derived from remotely sensed SSTs by the Advanced Very High-resolution Radiometer (AVHRR) infrared satellite data and in-situ measurements. With a spatial resolution of 0.25° on daily scale globally, this product represents the water temperature in the top 0.5 m of the ocean. In addition, another global daily data set named Group for High Resolution Sea Surface temperature (GHRSST) multi-product ensemble (GMPE) (refs.[38,39]) V2.0 (https://dap.ceda.ac.uk/neodc/esacci/sst/data/gmpe/; last access: September 20, 2021) including multiple SST data such as MyOcean OSTIA reanalysis, CMC 0.2 degree, AVHRR ONLY Daily 1/4 degree OISST, and MGDSST (ref. [40]) with the spatial resolution at 0.25° is also used.