15 MARCH 2024

DESER ET AL.

1877

Future Changes in the Intensity and Duration of Marine Heat and Cold Waves:
Insights from Coupled Model Initial-Condition Large Ensembles
CLARA DESER,a ADAM S. PHILLIPS,a MICHAEL. A. ALEXANDER,b DILLON J. AMAYA,b ANTONIETTA CAPOTONDI,b,c
MICHAEL G. JACOX,b,d,e AND JAMES D. SCOTTb,caNational Center for Atmospheric Research, Boulder, Coloradob
NOAA Physical Sciences Laboratory, Boulder, Coloradoc
CIRES, University of Colorado Boulder, Boulder, Coloradod
NOAA Southwest Fisheries Science Center, Monterey, CaliforniaeUniversity of California Santa Cruz, Santa Cruz, California(Manuscript received 11 May 2023, in ﬁnal form 22 December 2023, accepted 4 January 2024)
ABSTRACT: The future evolution of sea surface temperature (SST) extremes is of great concern, not only for the healthof marine ecosystems and sustainability of commercial ﬁsheries, but also for precipitation extremes fueled by moistureevaporated from the ocean. This study examines the projected inﬂuence of anthropogenic climate change on the intensityand duration of monthly SST extremes, hereafter termed marine heat waves (MHWs) and marine cold waves (MCWs),
based on initial-condition large ensembles with seven Earth system models. The large number of simulations (30–100) witheach model allows for robust quantiﬁcation of future changes in both the mean state and variability in each model. In general, models indicate that future changes in variability will cause MHW and MCW events to intensify in the northern extratropics and weaken in the tropics and Southern Ocean, and to shorten in duration in many areas. These changes are
generally symmetric between MHWs and MCWs, except for the longitude of duration change in the tropical Paciﬁc andsign of duration change in the Arctic. Projected changes in ENSO account for a large fraction of the variability-inducedchanges in MHW and MCW characteristics in each model and are responsible for much of the intermodel spread as a result of differences in future ENSO behavior. The variability-related changes in MHW and MCW characteristics noted
above are superimposed upon large mean-state changes. Indeed, their contribution to the total change in SST during
MHW and MCW events is generally ,10% except in polar regions where they contribute upward of 50%.
KEYWORDS: Atmosphere-ocean interaction; Extreme events; Sea surface temperature; Climate models;
Climate variability; Oceanic variability

1. Introduction
Prolonged sea surface temperature (SST) extremes, henceforth referred to as “marine heat waves” (MHWs) and “marine
cold waves” (MCWs), severely impact the health of aquatic ecosystems and have widespread economic repercussions for commercial ﬁsheries [see the recent review by Oliver et al. (2021)
and references therein; Smith et al. 2023]. Such events last fromweeks to months, and often recur in consecutive years over agiven region, compounding their physical and biogeochemical
effects. These SST extremes may also be accompanied by temperature anomalies at depth (e.g., Scannell et al. 2020; Amaya
et al. 2023a; Schaeffer et al. 2023). As climate change accelerates
due to human emissions of carbon dioxide and other greenhousegases, the oceans will continue to warm, thereby placing additional stress on temperature-dependent physiological processes in

Denotes content that is immediately available upon publication as open access.
Supplemental information related to this paper is available atthe Journals Online website: https://doi.org/10.1175/JCLI-D-230278.s1.
Corresponding author: Clara Deser, cdeser@ucar.edu

marine organisms (e.g., Alexander et al. 2018; Smale et al. 2019;
Cheung and Frölicher 2020).
Observational studies have documented the geographicalextent, intensity, duration, and frequency of MHWs based ondaily SST data from satellite records (1982–present), augmented
with monthly SSTs from ship-based archives going back to theearly twentieth century (Hobday et al. 2016; Holbrook et al.
2019; Oliver et al. 2021, and references therein). These studies
ﬁnd that historical MHW thresholds have been increasingly exceeded in recent years due to an underlying warming trend attributable to anthropogenic inﬂuences (Oliver et al. 2018; Oliver
2019; Oliver et al. 2021; Frölicher et al. 2018; Xu et al. 2022;
Frölicher and Laufkötter 2018; Laufkötter et al. 2020; Thoral
et al. 2022). Modeling studies have also highlighted the role oflong-term anthropogenic warming on recent and projected
trends in MHW historical threshold exceedance (Frölicher et al.
2018; Alexander et al. 2018; Yao et al. 2022; Guo et al. 2022).
However, anthropogenic effects are not limited to meanstate changes alone. As global warming progresses, the upper
ocean is expected to become increasingly stratiﬁed and themixed layer is anticipated to shoal (e.g., Capotondi et al. 2012;
Kwiatkowski et al. 2020). A shallower mixed layer will promote increased variability and decreased persistence of SST
anomalies due to its lower thermal inertia, all other factors being equal. According to this heuristic paradigm, anthropogenic
effects on long-term trends in SST extremes (e.g., MHWs and

DOI: 10.1175/JCLI-D-23-0278.1
Ó 2024 American Meteorological Society. This published article is licensed under the terms of the default AMS reuse license. For information regarding
reuse of this content and general copyright information, consult the AMS Copyright Policy (www.ametsoc.org/PUBSReuseLicenses).
Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1878

JOURNAL OF CLIMATE

MCWs) can result from both shifts in the mean (e.g., “a risingtide lifts all ships”) and from changes in the shape of the full
SST distribution as measured by the variance, skewness, andkurtosis. In addition, atmospheric circulation, air–sea exchangeof heat and momentum, oceanic thermal advection and turbulent
mixing, and remote impacts from changing El Niño–Southern
Oscillation (ENSO) teleconnections are all likely to be alteredunder anthropogenic climate change, further complicating thedrivers of SST extremes including MHWs and MCWs
(Alexander et al. 2018; Oliver et al. 2018; Sen Gupta et al.
2020; Amaya et al. 2021; Li and Thompson 2021; Oliver et al.
2021; Xu et al. 2022; Shi et al. 2022; Vogt et al. 2022; O’Brienand Deser 2023). Thus, the physical mechanisms underlying
observed and simulated changes in MHWs and MCWs arechallenging to quantify due to their complexity, varied geographical and seasonal expression, and lack of adequate data.
Only very recently have attempts been made to quantitatively
assess the various factors contributing to changes in historicaland future MHW characteristics (e.g., Oliver et al. 2021;
Amaya et al. 2021; Vogt et al. 2022).
Distinguishing between changes in the mean state andchanges in variability on the evolving characteristics of marine
temperature extremes is important not only from the standpointof physical understanding, but also for prediction and resource
management purposes (Oliver 2019; Amaya et al. 2023b). However, attributing changes in variability to anthropogenic effects
requires large amounts of data to refute the null hypothesis thatthey are due to random sampling ﬂuctuations. An innovative
study by Xu et al. (2022) employed an empirically based dynamical linear inverse model to generate thousands of synthetic realizations of SST variability consistent with the spatiotemporal
statistics of the actual observed record over the past 60 years.
With so many realizations, they were able to circumvent sampling issues and successfully isolate forced changes in variability
from forced changes in mean state on the historical evolution ofcold and warm SST extremes. They found that forced changes
in the mean state dominated those in variability; a similar conclusion was reached by Oliver (2019) based on an autoregressivemodel ﬁt to observations.
Studies of future changes in MHW characteristics have generally relied on the collection of CMIP5 and CMIP6 models,
using one simulation per model to construct a multimodel average (e.g., Frölicher et al. 2018; Plecha and Soares 2020; Qiu et al.
2021; Oliver et al. 2021; Yao et al. 2022). While this approachprovides an overall assessment, it confounds model structural
uncertainty (e.g., model-dependent forced responses in both themean state and variability) and sampling uncertainty. In contrast, “initial-condition large ensemble” simulations (hereafterreferred to simply as “large ensembles” or LEs) allow for a robust determination of the evolving forced responses in any single model at each location and time (e.g., Deser et al. 2020).
This is made possible by conducting a large number of simulations (generally 30 or more) with a given model under a given
radiative forcing protocol, where each simulation starts from adifferent initial condition. As the memory of the initial condition
fades, the simulations diverge from each other due to unpredictable internally generated variability. Thus, each ensemble
member has its own unique sequence of internal ﬂuctuations

VOLUME 37

superimposed upon an evolving forced response that is common
to all of the ensemble members. With a sufﬁcient ensemble size,
the evolving forced response can be isolated by averaging acrossthe ensemble members at each time step and at each location.
Importantly, the forced response contains not only evolving
changes in the mean background climate, but also any evolving
changes in the statistical characteristics of internal variability,
for example amplitude, autocorrelation and spatial structure.
The power of LEs has yet to be fully exploited for the studyof marine temperature extremes and their projected changes.
Frölicher et al. (2018) documented future changes in MHW
frequency, duration, spatial extent, and intensity based on a
10-member LE, while Alexander et al. (2018) used a 30-member
LE to study changes in SST extremes for large marine ecosystem
regions within the North Paciﬁc and Atlantic. More recently,
Burger et al. (2022) quantiﬁed compound MHW and ocean acidity extremes under global warming using a 30-member LE.
Here, we make use of seven different model LEs, each ofwhich contains between 30 and 100 simulations spanning theperiod 1950–2100, to isolate future changes in MHW and MCW
intensity and duration due to forced changes in variability versus
mean state in each model separately. We further examine theevolving inﬂuence of forced changes in ENSO variability on
forced changes in MHW and MCW characteristics by subsampling the LEs during ENSO-neutral states. While this subsampling procedure necessarily reduces the number of MHW and
MCW events for analysis, the large number of realizations in any
given model LE ensures that we maintain an adequate numberof events for robust results. While previous studies have shown aconnection between ENSO events and MHW events (e.g., Oliver
et al. 2018; Holbrook et al. 2019; Xu et al. 2022; Capotondi et al.
2022), none has explicitly assessed the inﬂuence of ENSO in thisway using LEs.
Finally, we analyze cold and warm extremes separately so
as not to build in any assumptions regarding linearity with respect to sign. Previous studies have generally focused exclusively on MHWs (e.g., Oliver et al. 2018; Frölicher et al. 2018;
Frölicher and Laufkötter 2018; Laufkötter et al. 2020; Guoet al. 2022) or on the symmetric component of SST variability
in general (Xu et al. 2022; Li and Thompson 2021; Shi et al.
2022). Schlegel et al. (2021) provide a global synthesis ofobserved MCW events (termed “marine cold spells” in theirstudy) based on daily SSTs over the satellite era (since 1982),
highlighting their physical characteristics, ecological impacts andrecent trends. Metrics of observed MCW and MHW events
based on daily satellite SSTs are compared in Wang et al. (2022)
and Yao et al. (2022). The latter study also compared futurechanges in MCW and MHW properties using daily SSTs fromthe CMIP6 archive; however, their multimodel mean assessment is based on averaging one simulation per model and thus
does not separate the effects of variability-induced changes
from mean-state changes.
Our study is based on monthly SST data rather than 5-day
running means as in the original MHW deﬁnition of Hobday
et al. (2016). The use of monthly data can be justiﬁed basedon the relatively slow decorrelation time scale of SST anomalies and the outsized ecological impacts from long-lasting
marine temperature extremes (Jacox et al. 2020; Amaya et al.

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

1879

DESER ET AL.

TABLE 1. Model large ensembles used in this study and their salient characteristics. Global SST change refers to the difference
between 2070–2100 and 1970–2000. MHW and MCW results for CESM2 are not sensitive to differences in biomass burning aerosol
emission protocols used for the ﬁrst and second 50 members.

Model name (CMIP
generation)

Duration; futureradiative forcing
scenario; global SST
changeCESM2 (CMIP6)

1850–2100; SSP3-7.0;
2.58C
CanESM5 (CMIP6)
1850–2100; SSP5-8.5;
4.08C
GFDL-SPEAR (CMIP6) 1921–2100; SSP5-8.5;
2.78C
MIROC6 (CMIP6)
1850–2100; SSP5-8.5;
2.08C
CESM1 (CMIP5)
1920–2100; RCP8.5;
2.68C
CanESM2 (CMIP5)
1950–2100; RCP8.5;
3.18C
MPI-ESM-LR (CMIP5) 1850–2099; RCP8.5;
1.98C

Ensemblesize

Resolution foratmosphere (ocean)

MHW/MCW
sample size (pergrid box per 31
years) for all years(ENSO-neutral
years)

Reference paper

100

1.38 3 0.98 (;1.08)

50

2.88 3 2.88 (1.48 3 0.98)

30

50 km (;1.08)

50

1.48 3 1.48 (;1.08)

699/710 (218/210) Tatebe et al. (2019)40

1.38 3 0.98 (;1.08)

619/608 (171/164) Kay et al. (2015)50

2.88 3 2.88 (1.48 3 0.98)

100

1.98 3 1.98 (;1.58)

2023b). However, shorter-lived MHW and MCW events will
be missed in our analyses. We encourage future studies tomake use of daily SST (and subsurface temperature) output
from model LEs for added insight [see Le Grix et al. (2022)
for a recent example].
The rest of this study is organized as follows. The model LEs
and our methods for deﬁning MHW and MCW events, separating forced changes due to mean state versus variability, and classifying ENSO-neutral states are described in section 2. Results
are presented in section 3, beginning with an evaluation ofmodel biases in simulating present-day MHW and MCW intensity and duration (section 3a), followed by an analysis ofmodels’ future changes in MHW and MCW characteristics(section 3b) and the role of ENSO (section 3c), and ending
with an assessment of the relative contributions from meanstate changes vs. changes in variability (section 3d). A summary and discussion are provided in section 4.

2. Data and methods
a. Data from model LEs and observations
We analyze seven global coupled model LEs with ensemble
sizes ranging from 30 to 100 spanning the years 1950–2100 underhistorical and future radiative forcing (see Table 1 for modeland forcing details). We chose all available models that had atleast 30 ensemble members for the entirety of the analysis period 1950–2100 (a full listing of model LEs is provided at https://www.cesm.ucar.edu/community-projects/mmlea). Three of themodels (CESM1, CanESM2, and MPI-ESM-LR) are CMIP5
generation and the others (CESM2, CanESM5, GFDL-SPEAR,
and MIROC6) are CMIP6. Thus, our multimodel LE collection
encompasses a range of forcing scenarios and associated global
warming responses. In addition, it includes two versions of
CESM and two of CanESM, allowing for comparison across1553/1493 (391/377) Rodgers et al. (2021)
816/793 (183/175) Swart et al. (2019)553/546 (121/113) Delworth et al. (2020)

847/824 (189/180) Kirchmeier-Younget al. (2017)1619/1618 (412/404) Maher et al. (2019)

CMIP generations of the same model. In total, our multimodel
LE collection contains 420 individual simulations. All models
employed in this study utilize relatively coarse resolution(around 1.08–1.58) in their ocean component and thus do not explicitly resolve mesoscale processes that have been shown to beimportant for SST variability along western boundary currents,
the Antarctic Circumpolar Current, and other eddy-rich regions(Pilo et al. 2019; Hayashida et al. 2020; Guo et al. 2022).
Our analyses are based on monthly mean SST on the native
grid of each model (see Table 1 for details on model resolution). For the models, we use the temperature of the top ocean
model level (“tos”), which is the true SST, as opposed to surface temperature (TS) to avoid contamination from sea ice
cover. Speciﬁcally, for a grid cell that includes fractional sea
ice, tos is the water temperature directly beneath the sea ice,whereas TS is the area-weighted average of the surface temperature of the ice and the surface temperature of the water
within the grid cell. To benchmark the models against observations, we make use of SSTs from ERSSTv5 on a 28 latitude/
longitude grid (Huang et al. 2021); similar results (not shown)are obtained with HadISSTv1 on a 18 latitude/longitude grid
(Rayner et al. 2003). When comparing models against observations, and forming multimodel averages, we have regridded
the data to the CESM1 grid.

b. Methods
1) MODEL LARGE ENSEMBLES
As mentioned above, we separate future changes in MHWs
and MCWs into two components: 1) changes in variability
and 2) changes in the mean state. To isolate changes in variability, we deﬁne the internal component of SST (iSST) in
each ensemble member of a given model LE as iSST(x, t, e) 5
SST(x, t, e) 2 SST(x, t, em) where x is grid box location, t is

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1880

JOURNAL OF CLIMATE

time, e is ensemble member, and em is the ensemble mean.
Note that this procedure removes the monthly ensemble-mean
climatology from each ensemble member. We then computethe 10th and 90th percentiles of iSST(x) separately for each calendar month based on pooling together all ensemble membersof a given model LE. We compute these percentiles for the
71-yr historical period 1950–2020 in order to compare with observations, and also for three 31-yr periods (1970–2000, 2020–50,
and 2070–2100) in order to evaluate how climate change inﬂuences MHWs and MCWs. For example, for the 100-member
CESM2 LE, we determine these thresholds for the month ofJanuary for the period 1970–2000 from the set of 31 3 100 iSST
samples at each grid box, and repeat this procedure for each calendar month. Similarly, we determine the thresholds as a function of calendar month for the period 2020–50 and for theperiod 2070–2100. In this way, we obtain seasonally varying
thresholds for each period separately. Then, for each period, wedeﬁne a MHW event if iSST(x, t, e) reaches or exceeds the 90th
percentile; similarly, we deﬁne a MCW event if iSST(x, t, e) is
equal to or less than the 10th percentile. We deﬁne the intensityof each MHW and MCW by its corresponding monthly iSST
value (8C) and the duration of each MHW (MCW) event bydetermining the number of consecutive months with MHW(MCW) conditions at a given location. We then form MHW
and MCW composites of intensity (8C) and duration (months)
based on averaging all MHW and MCW events for each period
separately. To give an idea of the sample sizes for these composites, the 100-member CESM2 LE yields 3720 (12 months 3
31 years 3 100 members 3 0.1) monthly MHW samples (andMCW samples) in any given 31-yr period at each grid box. The
number of discrete MHW and MCW events (e.g., taking into
account the lifetime of each event) depends on the duration
characteristics, which vary spatially. For the CESM2 LE, eachgrid box contains, on average, 1553 discrete MHW events and1493 discrete MCW events in any given 31-yr period. Sample
size information for all seven model LEs is provided in Table 1
and spatial maps of the number of discrete MHW and MCW
events are shown in Fig. S1 in the online supplemental material.
The large sample sizes afforded by each model LE allow for arobust determination of MHW and MCW characteristics andany future changes thereof.
We have repeated our analyses using the 5th and 95th percentile thresholds to deﬁne MCW and MHW events, respectively, and ﬁnd qualitatively similar results to those based on
the 10th and 90th percentile thresholds. Thus, we adopt the10th and 90th percentile thresholds in order to provide a robust comparison with observations for which the 5th and 95th
percentile thresholds would be too restrictive (i.e., result in
too few samples; see also Jacox et al. 2020). In addition, the10th and 90th percentile thresholds allow us to robustly quantify future changes in MHW and MCW events during the subset of years that are in an ENSO-neutral state (see below).
To assess the inﬂuence of ENSO on future changes in MHWs
and MCWs, we have repeated our analyses for ENSO-neutral
states as follows. First, we compute the leading empirical orthogonal function (EOF1) of 3-month running means of iSST(x, t)over the tropical Paciﬁc domain (108N–108S, 1208E–858W) for
each time period separately (1970–2000, 2020–50, and 2070–2100),

VOLUME 37

using the pooled set of ensemble members for a given model LE.
(Note that the use of EOF1 as opposed to a ﬁxed regional index
such as Niño-3.4 SST ensures that the pattern of ENSO-relatediSST speciﬁc to each model is accommodated.) We then identify
the 30th and 70th percentile values of the associated principal
component (PC1) for each month separately. If the PC1 value ina given ensemble member e at time t lies within the 30th–70th percentile range, and if the PC1 values in each of the preceding
5 months also fall within this range, then we deﬁne the corresponding iSST(x, t, e) as being in an “ENSO-neutral” state [wechoose 5 months as a conservative estimate of the cumulative
lagged response of extratropical SST anomalies to ENSO (seeAlexander et al. 2002); however, a choice of 3 months yields
very similar results (not shown)]. We repeat this procedure forall time steps and ensemble members. We then determine the10th and 90th percentile thresholds of iSST(x) for each monthand time period from this set of ENSO-neutral samples, andthen use these thresholds to construct MHW and MCW composites from the ENSO-neutral samples. Sample size information
for the ENSO-neutral composites in each model LE is provided
in Table 1 and Fig. S2. We have checked that the results do notchange appreciably if we consider both PC1 and PC2 of tropicalPaciﬁc iSSTs in our ENSO-neutral selection procedure; however, the use of two PCs greatly reduces the number of ENSOneutral samples available for compositing, resulting in less stablestatistics and noisier patterns (not shown). While more (or less)sophisticated methods could be used to deﬁne ENSO-neutral
states, we view our approach as a valuable step in assessing theinﬂuence of future changes in tropical Paciﬁc SST variability on
MHW and MCW characteristics.

2) OBSERVATIONS
We ﬁrst compute the climatological SST for each monthat each grid box based on the period 1950–2020. We then deﬁne monthly SST anomalies by subtracting the climatologicalmonthly SST from each corresponding month. We then quadratically detrend the data to produce the estimated observediSST(x, t), hereafter iSST_Obs(x, t). We then compute the 10th
and 90th percentiles of iSST_Obs(x, t) as a function of calendar
month, and deﬁne MHW and MCW events following the procedure outlined for the model LEs above. We used quadratic
detrending in lieu of linear detrending to accommodate the accelerating pace of global warming since 1950 (NOAA National
Centers for Environmental Information 2023).
We have tested whether quadratic detrending effectively removes the forced component by applying this procedure to eachmember of a given model LE individually and comparing
the results against the “true” forced component estimated by
the ensemble mean of the LE. We ﬁnd that for the period
1950–2020, quadratic detrending produces nearly identical MHW
and MCW composites as the method outlined in section 2b(1)for each model LE (not shown), lending conﬁdence to the observational results.

3) STATISTICAL SIGNIFICANCE TESTING
We apply the false discovery rate (FDR; Wilks 2016) to a
two-sided t test at the 95% conﬁdence level to assess whether

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

DESER ET AL.

differences in observed and simulated MHW and MCW composites are statistically signiﬁcant, and whether simulated
MHW and MCW composites in the future (2020–50 or 2070–
2100) differ signiﬁcantly from those in the historical reference
period (1970–2000). We use the number of discrete events ineach composite for the t-test sample sizes. Note that the t statistic (critical t value) is relatively insensitive to the precise
number of samples for samples sizes . 75 (Student 1917) andthat our model composites always contain at least 75 discrete
samples at every grid box.

3. Resultsa. Model validation
We begin by comparing the simulated and observed characteristics of MHW and MCW composites over the 71-yr period
1950–2020. For conciseness, we present results from the
100-member CESM2 LE in the main paper, and the remaining
six model LEs in the online supplemental material. Figure 1ashows MHW composite intensity (8C) averaged across all
100 members of the CESM2 LE. The model simulates a maximum in intensity along the equatorial Paciﬁc (values ;3.58C)
and secondary maxima along the western boundary currents
and their eastward extensions in the North Paciﬁc and North
Atlantic (values ;28C). Similar features are evident in the observations, but with reduced amplitudes, especially in theequatorial Paciﬁc (Fig. 1b). To assess whether the model is signiﬁcantly biased in its representation of MHW composite intensity, or whether the differences between the simulated andobserved MHW composites arise from sampling uncertainty(the observational record samples only one of many possible
realizations), we evaluate whether the observed composite values lie outside the 5th–95th percentile range of the distributionof 100 composite values based on each member of the CESM2
LE. According to this measure, CESM2 shows a signiﬁcantwarm bias in composite MHW intensity over much of the tropical Indo-Paciﬁc and along Antarctica, and a signiﬁcant cold
bias in the Arctic, temperate portions of the Southern Ocean,
and along the west coasts of North America and Africa(Fig. 1c). The amplitudes of these model biases are generally,0.28C except in the western equatorial Paciﬁc and offequatorial eastern Paciﬁc and along the Siberian coast where
they reach 0.48–1.08C. Similar model biases are found for
MCW composite intensity, with an even larger overestimation
across the full width of the equatorial Paciﬁc basin than is thecase for MHW intensity [Fig. 1i; note that negative (positive)
values indicate that the model overestimates (underestimates)
MCW intensity]. The biases in MHW and MCW intensity
within the tropical Indo-Paciﬁc are likely related to the model’s overestimation of ENSO amplitude (Capotondi et al.
2020), while those in the Arctic and Antarctic may be associated with deﬁciencies in the representation of sea ice cover(Danabasoglu et al. 2020).
We also examine composite MHW and MCW intensity ineach ensemble member for a direct assessment of the sampling
uncertainty due to the ﬁnite record length (note that thethresholds are still based on those obtained from all members).

1881

Figures 1d and 1e show the maximum and minimum values of
MHW composite intensity, respectively, across the 100 members of the CESM2 LE. To construct these ﬁgures, we have selected the ensemble member with the maximum composite value
and the one with the minimum composite value at each grid box
separately (i.e., the ensemble member with the maximum composite value at one grid box may differ from the ensemble member with the maximum composite value at another grid box). The
two renditions of MHW composite intensity show striking differences in amplitude throughout the global oceans, with a typicalrange of 0.48–0.88C (Fig. 1f). Larger differences (up to 1.68C) arefound in the equatorial eastern Paciﬁc, the Kuroshio and Gulf
Stream Extensions, and along the sea ice edge between Icelandand Svalbard. Similar results are found for MCW composite intensity (Figs. 1j–l), although the ensemble spread in the equatorial Paciﬁc maximizes farther west compared to the MHW
results (Fig. 1l). Such sampling variations in both MHW andMCW composite intensities over a 71-yr period are noteworthy
and imply that the observed composites may also be subject tosimilar levels of sampling uncertainty. [Note that the patterns
shown in Figs. 1d,1e,j,k may not necessarily obtain in any individual simulation.]
Next, we compare the simulated and observed MHW
and MCW composites of duration. The CESM2 ensemblemean MHW composite duration exhibits largest values(up to 6–7 months) in the tropical eastern Paciﬁc, and secondary
maxima (up to 3–4 months) in the western tropical Atlantic and
the Paciﬁc sector of the Southern Ocean (Fig. 2a). These maxima are associated with El Niño events (see section 3c). Elsewhere, the simulated composite durations are typically between
2 and 3 months. Observations show a similar, albeit noisier, pattern of MHW composite duration as the model ensemble mean,
but with generally smaller magnitudes (1–2 months) except in thetropical northeast Paciﬁc (Fig. 2b). In particular, the observed values lie outside the model’s 5th–95th percentile ensemble spread(indicative of a signiﬁcant model bias) in the southeast and northeast tropical Paciﬁc, south tropical Atlantic, central North Paciﬁc,Arctic, and portions of the Southern Ocean (Fig. 2c).
The ensemble spread in MHW composite duration in CESM2
is striking (Figs. 2d,e). In general, MHW composite duration
varies by approximately a factor of 2 between the minimum andmaximum of the 100 ensemble members throughout the global
oceans. The absolute range is largest in the eastern tropicalPaciﬁc and the Greenland Sea (up to 7 months), with secondary
maxima in the Amundsen–Bellingshausen Sea, and parts of thetropical South Atlantic and subtropical northeast Paciﬁc (Fig. 2f).
In the eastern tropical Paciﬁc, MHW composite duration can exceed 10 months or be less than 4 months, depending on the realization (Figs. 2d,e). This sampling variability is associated with
low-frequency modulation of ENSO (Capotondi et al. 2020).
Unlike intensity, duration exhibits noticeable asymmetries between MHW and MCW composites (Fig. 2). In particular, theobserved MCW composite duration shows a boomerang-shaped
pattern in the central tropical Paciﬁc, with off-equatorial maxima in both hemispheres (Fig. 2h), whereas the observed MHW
composite duration exhibits peak values along the equator in
the eastern Paciﬁc (Fig. 2b). In CESM2, the duration asymmetry
takes the form of a westward-extended maximum in the

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1882

JOURNAL OF CLIMATEVOLUME 37

FIG. 1. Composite MHW and MCW intensity (8C) during 1950–2020 from the 100-member CESM2 large ensemble and observations:(a),(g) ensemble average; (b),(h) observations; (c),(i) ensemble average minus observations; (d),(j) ensemble maximum; (e),(k) ensemble
minimum; and (f),(l) ensemble maximum minus minimum. Gray shading in (c) and (i) indicates that observations lie within the 5th–95th
percentile range of the CESM2 large ensemble.

ensemble-mean MCW composite (Fig. 2g) compared to its
MHW counterpart (Fig. 2a). Model biases in MCW composite
duration are broadly similar to those in MHW composite duration, with additional areas of signiﬁcant bias in the far westernequatorial Paciﬁc and south of Hawai’i, where the values areoverestimated and underestimated by 2–3 months, respectively(Fig. 2i). The ensemble spread in composite duration is evenmore pronounced for MCW than MHW, especially in the tropical Paciﬁc where the range between the maximum and minimum of the individual members can exceed 8 months at many
locations (Figs. 2j–l).
Analogous comparisons between simulated and observed
MHW and MCW composites for the other six model LEs areshown in Figs. S3–S8 for intensity and Figs. S9–S14 forduration. Figures 3a–g compare the model biases in MHW
composite intensity for all seven LEs, depicted as the difference between each model’s ensemble-mean composite andthe observed composite. While the regional details of the patterns of signiﬁcant bias vary across models, there are somecommon large-scale features. For example, all models show asigniﬁcant overestimation of MHW composite intensity in thetropical Indian Ocean and the western half of the tropicalPaciﬁc. This positive bias is smallest for the two CanESM LEs
and is largest in CESM2 and MIROC6 where it extends acrossmost of the Paciﬁc basin. Another common feature across models
is the signiﬁcant underestimation of MHW composite intensity in
the Arctic, and a signiﬁcant overestimation in the Antarctic.
Many of the models show a signiﬁcant underestimation of

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

1883

DESER ET AL.

FIG. 2. As in Fig. 1, but for composite MHW and MCW duration (months).

MHW composite intensity along the west coast of North
America, and all models except the two CESM LEs show asigniﬁcant positive bias in the North Atlantic. The fractional
area of signiﬁcant model bias in MHW composite intensity
regardless of sign is similar across models, ranging from 51%
in CanESM5 to 68% in MIROC6 (numbers in the upper rightof each panel in Fig. 3; see also Fig. 4). Positive biases makeup most of this total areal coverage, except in CESM1 and
CanESM2 where positive and negative biases contribute inroughly equal proportion (Fig. 4a). To summarize model behavior, we map the locations where at least 5 out of 7 models(.two-thirds) show a signiﬁcant model bias and agree on thesign of the bias (Fig. 3h). Areas of majority agreement include the tropical northern Indian and western Paciﬁc andhigh-latitude Southern Ocean where MHW intensities are
signiﬁcantly overestimated, and the Arctic Ocean where they

are signiﬁcantly underestimated; these areas comprise 34% of
the world oceans.
Model biases in MCW composite intensity are generally
analogous to their MHW counterparts, but with slightly larger
amplitudes (Figs. 3i–o). [Note that the opposite sign between
model biases in MCW and MHW composite intensity indicatesthat the biases are the same in a relative sense: i.e., a negative
bias in MCW intensity and a positive bias in MHW intensity
both denote overestimation.] There are also some asymmetries
in regional expression of the biases. For example, the signiﬁcant
overestimation of composite intensity in the equatorial Paciﬁc
extends farther east for MCW compared to MHW in most models. The areal coverage of signiﬁcant model bias regardless ofsign is also similar for MCW as MHW composite intensity, witha comparable proportion of over- and underestimation withineach model (Fig. 4b). Areas of majority agreement in the sign

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1884

JOURNAL OF CLIMATEVOLUME 37

FIG. 3. Composite MHW and MCW intensity (8C) during 1950–2020 for the ensemble mean of each model large ensemble minus observations: (a),(i) CESM2; (b),(h) CanESM5; (c),(i) GFDL-SPEAR; (d),(j) MIROC6; (e),(k) CESM1; (f),(l) CanESM2; and (g),(o) MPIESM-LR. Gray shading indicates that observations lie within the 5th–95th percentile range of the model large ensemble. The number inthe upper right of each panel denotes the fractional area (%) of signiﬁcant model bias (e.g., nongray areas). Shading in (h) and (p) shows
locations where at least two-thirds of the models show a signiﬁcant bias (pink shading for positive bias and blue shading for negative bias);
the number in the upper right denotes the fractional area (%) of the pink and blue shading.

of signiﬁcant model bias in MCW composite intensity are similar to those for MHW (Fig. 3p).
Unlike intensity, there are considerable differences amongmodels in the pattern and amplitude of signiﬁcant bias in
MHW (Figs. 5a–g) and MCW (Figs. 5i–o) composite duration.
For example, in the tropics, four of the models (CanESM5,
GFDL-SPEAR, CanESM2, and MPI-ESM-LR) show signiﬁcant underestimation MHW and MCW duration (by 1–2 months),
while the remaining three models (CESM2, MIROC6, andCESM1) show insigniﬁcant or opposite-signed biases. The former group of models also exhibits consistent asymmetries in theirpatterns of duration bias between MHW and MCW composites,with an equatorial maximum for MHW and off-equatorial(“boomerang” shaped) maxima for MCW. There is greateragreement among the 7 LEs in the sign of their extratropicalbiases, with signiﬁcant overestimation (by 1–2 months) of composite duration for both MHW and MCW within the Arctic,
North Paciﬁc and Atlantic, and Southern Ocean, although the details within each region vary considerably. There is less consensus
among models on the signiﬁcance and sign of bias in durationcompared to intensity for both MHW and MCW (Figs. 5h,p).
Only 17% and 21% of the area of the world oceans show model
consensus on MHW and MCW duration bias, respectively. Further, the regions of model consensus are less spatially coherent
for duration biases compared to intensity biases, with the North
Paciﬁc showing the most widespread areas of model agreement.
The areal coverage of signiﬁcant model bias regardless ofsign is generally smaller for duration than intensity, ranging
from 33% for MHW and 39% for MCW in CESM2 to 53%
for MHW and 55% for MCW in CanESM5 (Fig. 4b). The proportion of positive versus negative duration bias contributing to
the total areal coverage for both MHW and MCW composites
varies by model, with CESM2, GFDL-SPEAR, MIROC6, andCESM1 showing a substantial imbalance and the remaining
models showing more equivalent values (Fig. 4b).
In summary, all model LEs examined in this study show
widespread statistically signiﬁcant biases in their historicalsimulation of composite MHW and MCW intensity and duration deﬁned on the basis of monthly SSTs and taking into account sampling uncertainty due to the limited length of the

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

DESER ET AL.

FIG. 4. Fractional area (%) of signiﬁcant positive (red bars) andnegative (blue bars) model bias in composite MHW (open bars)
and MCW (ﬁlled bars) (a) intensity and (b) duration during 1950–
2020. Note that a positive (negative) bias in MHW (MCW) intensity indicates model overestimation.

observational record. These biases cover ;50%–70% of thearea of the global oceans for intensity and ;35%–55% for duration in every model. The majority of model LEs signiﬁcantly overestimate MHW and MCW intensity in the tropicalwestern Paciﬁc and Indian Oceans and across portions of theSouthern Ocean, and signiﬁcantly underestimate it in the
Arctic. The models show less consistency in their durationbiases both in terms of sign and location, although the North
Paciﬁc is a region of general overestimation.

b. Future changes
Next, we examine future changes in MHW and MCW composite intensity and duration due to changes in variability(mean state changes are considered in section 3d). As before,
we present results from the 100-member CESM2 LE in the main
paper and the remaining six model LEs in the supplemental
material. Note that all model LEs use relatively high-emissions
scenarios for their projections (Table 1). Figure 6 shows composite intensity maps for three periods (1970–2000, 2020–50, and
2070–2100) and the differences between each future period andpresent-day (e.g., 2020–50 minus 1970–2000, and 2070–2100 minus

1885

1970–2000) for the CESM2 LE. All three periods show similar
spatial patterns of MHW composite intensity, but their amplitudes
evolve over time (Figs. 6a–c). Relative to the present day, theMHW composite intensity at midcentury is projected to increasesigniﬁcantly in the Arctic and along the Antarctic coastline, andwithin portions of the North Paciﬁc, North Atlantic and thetropics, and decrease signiﬁcantly within the Southern Ocean andparts of the tropical Indian, western Paciﬁc, and eastern Atlanticbasins, with amplitudes generally ,0.28C (Fig. 6d). The increasein MHW intensity within the Arctic and along the Antarctic coastline is likely related to the melting of sea ice and subsequent exposure of the ocean to the atmosphere, allowing for variability inheat exchange across the air–sea interface. We speculate that thepatterns of MHW intensity change within the North Paciﬁc andNorth Atlantic are potentially related to altered patterns of atmospheric circulation variability (O’Brien and Deser 2023), while thedecreases in MHW intensity in the Southern Ocean may be due
in part to a poleward shift of the storm track and associated reduction in air–sea heat and momentum ﬂux variability (Wu et al.
2012; Guo et al. 2022). A similar pattern of future change is
projected for late century as for midcentury, with an approximately twofold increase in amplitude, except in the tropicalPaciﬁc where MHW composite intensity is projected to decreasesigniﬁcantly (up to 0.68C), opposite to the increase seen at midcentury (Fig. 6e). The diverging mid- and late-century projections in the tropical Paciﬁc are related to changes in the
behavior of ENSO as will be shown in section 3c. In CESM2,ENSO is projected to intensify until about 2050 and weaken
thereafter (Maher et al. 2023), consistent with the changes in
MHW intensity in the tropical Paciﬁc. Future changes in MCW
composite intensity are largely similar to those in MHW intensity, both in pattern and amplitude (Figs. 6i,j). Analogous comparisons between present-day and future MHW and MCW
composite intensity for each of the six other model LEs areshown in Figs. S15–S20.
Future changes in MHW and MCW composite duration in
CESM2 are shown in Fig. 7. Statistically signiﬁcant changes begin to emerge midcentury and become widespread by late century for both MHW and MCW (Figs. 7d,e,i,j). The late-century
changes in composite duration, which are largely similar between MHW and MCW, include reductions over the tropicalIndo-Paciﬁc of up to 2–3 months and smaller decreases over theSouthern Ocean and eastern North Paciﬁc, accompanied by increases over the North Atlantic and near Antarctica (Figs. 7e,j).
Notable asymmetries occur over the tropical Paciﬁc, where latecentury duration changes maximize in the east for MHW and in
the west for MCW, likely related to differences in their presentday climatological mean states (Fig. 7a vs. Fig. 7f). Anotherasymmetry is apparent in the western tropical Atlantic, where
signiﬁcant late-century increases in duration occur for MCW
but not for MHW. We speculate that this difference may be related to the signiﬁcant reduction in MHW duration over theeastern tropical Paciﬁc, which may act to offset (and mask) the
increase that would otherwise occur over the western tropicalAtlantic via Paciﬁc-to-Atlantic interbasin teleconnection mechanisms (similar to those which occur for ENSO events; see, e.g.,
Alexander et al. 2002). Analogous maps of MHW and MCW

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1886

JOURNAL OF CLIMATEVOLUME 37

FIG. 5. As in Fig. 3, but for composite MHW and MCW duration (months).

composite duration in present-day and future periods for eachof the six other model LEs are shown in Figs. S21–S26.
Next, we compare the projected late-century changes incomposite intensity across the seven model LEs for both
MHW and MCW (Fig. 8); midcentury changes are shown in
Fig. S27. The models show considerable diversity in their patterns of projected change, especially at low latitudes. Overthe tropical Paciﬁc, for example, CanESM2 and CESM2 and
MPI-ESM-LR show relatively pronounced reductions in intensity for both MHW and MCW, while the other models show
more muted changes of mixed sign (Fig. 8). There is some correspondence between the patterns of future change in MHW/
MCW intensity over the tropical Paciﬁc and the projected behavior of ENSO across models, but it is not exact. For example,
CanESM2 and CESM2 show pronounced late-century reductions in ENSO variance (Maher et al. 2023), consistent with decreasing MHW/MCW intensity, but CanESM5 and MIROC6
show increases in ENSO variance and MPI-ESM-LR shows little change in ENSO variance (Maher et al. 2023), yet thesemodels also simulate decreasing MHW/MCW intensity (albeit
less than in CanESM2 and CESM2). Models agree that MHW
and MCW composite intensity will increase in the Arctic andalong Antarctica, and diminish in the Indo-Paciﬁc sector of theSouthern Ocean by late century (Figs. 8h,p). There is also model

consensus that composite intensity will increase over portions of
the subtropics and that it will decrease in the tropics on eitherside of the equator, especially for MHWs (Figs. 8h,p). The fractional area of the World Ocean that is projected to experience asigniﬁcant late-century change in MHW and MCW composite
intensity (regardless of sign) varies from a low of 64%–65% in
GFDL-SPEAR to a high of 84%–85% in CESM2 (Fig. 8). Of
this fractional area, the proportion of signiﬁcant increase outweighs that of signiﬁcant decrease in four of the models, andvice versa in two of the models (Figs. 9a,b). Model consensus on
the signiﬁcance and sign of projected changes in MHW and
MCW intensity is found over 36%–38% of the area of theWorld Ocean (Figs. 8h,p).
Like intensity, models show a variety of patterns and amplitudes of projected late-century change in MHW and MCW
composite duration (Fig. 10; midcentury changes are shown in
Fig. S28). For example, CESM2 features an interbasin contrast
between the North Atlantic and the tropical Paciﬁc, while CanESM5 shows out of phase changes between the tropics and extratropics (excluding the Antarctic coastline). Models disagree
on the sign of duration change over the Arctic, where signiﬁcant
increases are simulated by CESM1, CanESM2, and MPI-ESM-LR,
and signiﬁcant decreases by CanESM5, MIROC6, and GFDLSPEAR. The reasons for this are unclear. There is little agreement

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

1887

DESER ET AL.

FIG. 6. Composite MHW and MCW intensity (8C) from the 100-member CESM2 large ensemble during (a),(f) 1970–2000, (b),(g) 2020–2050,
and (c),(h) 2070–2100, and differences (d),(i) 2020–50 minus 1970–2000, and (e),(j) 2070–2100 minus 1970–2000. Note that the color bar range istwice as large in (a)–(c) and (f)–(h) compared to (d), (e) and (i),(j). Gray shading in (d), (e), (i), and (j) indicates that the differences are not statistically signiﬁcant according to the false discovery rate applied to a two-sided t test at the 95% conﬁdence level. Boxes in (e) and (j) outline regions used for Figs. 17 and 18. Note that positive (negative) values in (d) and (e) and negative (positive) values in (i) and (j) indicate an increase(decrease) in intensity.

among models on the spatial pattern of duration changes within
the tropical Paciﬁc, although modest decreases in MHW duration
are found in the majority of models in the far west (Figs. 10h,p).
Likewise, there is model consensus that MHW and MCW duration will decrease over much of the North Paciﬁc and at scatteredlocations within the Southern Ocean except along Antarctica. Fiveof the seven models indicate that the fraction of the world
oceans projected to experience a signiﬁcant decrease in MHW
and MCW duration will outweigh the fraction with a signiﬁcant
increase (Figs. 9c,d). Finally, model consensus on the sign of signiﬁcant projected changes in MHW and MCW duration isfound for 27%–28% of the area of the World Ocean, somewhat
less than for intensity (Figs. 10h,p).

c. Role of ENSO
It is clear from the preceding results that the tropical Paciﬁc
is a region of high-amplitude MHW and MCW intensity andduration, and that the projected changes in this region are
subject to large intermodel spread. Model LEs also differ substantially in their ENSO projections, including the sign, spatial
pattern, and time dependence of future change, for reasons

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1888

JOURNAL OF CLIMATEVOLUME 37

FIG. 7. As in Fig. 6, but for composite MHW and MCW duration (months).

that are not well understood (Maher et al. 2023). Given that
SST variability in the tropical Paciﬁc associated with ENSO is
known to impact remote regions via atmospheric and oceanic
teleconnections, this raises the question: How much do projected changes in ENSO inﬂuence projected changes in MHW
and MCW characteristics worldwide? To address this question, we
recompute the MHW and MCW composites for ENSO-neutral
conditions following the procedures outlined in section 2b(1). As
before, we begin by showing results for the CESM2 LE, followed by an intercomparison of late-century changes in ENSOneutral MHW and MCW composite intensity and duration
across all seven model LEs.
As expected, the tropical Paciﬁc maxima in MHW andMCW composite intensity amplitudes in CESM2 are greatly
reduced in the ENSO-neutral samples compared to all

samples in all three time periods (Figs. 11a–c,f–h comparedwith Figs. 6a–c,f–h). Other regions, most notably the North
Paciﬁc, also show a slight decrease in amplitude in each time
period. It is clear that changes in ENSO variability impact future changes in MHW and MCW composite intensity. For example, at midcentury, the areal coverage of signiﬁcant change
is reduced from 70% in the all-sample composites (Figs. 6d,i)to 53% in the ENSO-neutral composites (Figs. 11d,i), alongwith some alterations in the sign and pattern of change in the
Paciﬁc sector. By late century, the differences between the
all-sample (Figs. 6e,j) and ENSO-neutral (Figs. 11e,j) composite changes over the North Paciﬁc are particularly striking.
The ENSO-neutral composites show a pattern reminiscent of
the positive phase of the Paciﬁc decadal oscillation (PDO;
Newman et al. 2016), with negative values over the western

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

DESER ET AL.

1889

FIG. 8. Composite MHW and MCW intensity (8C) differences between 2070–2100 and 1970–2000 for the ensemble mean of each model
large ensemble: (a),(h) CESM2; (b),(i) CanESM5; (c),(j) GFDL-SPEAR; (d),(k) MIROC6; (e),(l) CESM1; (f),(m) CanESM2; and (g),(n)MPI-ESM-LR. Gray shading indicates that the differences are not statistically signiﬁcant according to the false discovery rate applied to a
two-sided t test at the 95% conﬁdence level. The number in the upper right of each panel denotes the fractional area (%) of signiﬁcant differences (e.g., nongray areas). Note that positive (negative) values in (d) and (e) and negative (positive) values in (i) and (j) indicate an increase (decrease) in intensity. Shading in (h) and (p) shows locations where at least two-thirds of the models show statistically signiﬁcantvalues (pink for positive and blue for negative); the number in the upper right denotes the fractional area (%) of the pink and blue
shading.

and central North Paciﬁc and positive values in the eastern
North Paciﬁc, for both MHW and MCW and opposite in sign
to the all-sample composites. The fact that the ENSO-neutral
late-century changes in intensity are of the same sign for
MHW and MCW indicates that their responses are strongly
asymmetric (e.g., opposite directions of change) in this region.
The similarity of future changes in intensity over the Arctic,
North Atlantic, and Southern Ocean between the ENSO-neutral
and all-sample composites indicates that changes in ENSO variability play only a minor role in these areas. It is interesting tonote that the late-century decrease in ENSO-neutral composite
intensity over the tropical Paciﬁc and Atlantic is larger for MCW
than MHW (cf. yellow shading in Fig. 11j with blue shading in
Fig. 11e), a result which is insensitive to whether one or two tropical Paciﬁc PCs are used in the deﬁnition of ENSO-neutral samples (not shown).
MHW and MCW durations in CESM2 are considerably reduced in the ENSO-neutral composites compared to their“all-sample” counterparts in each time period (1970–2000, 2020–50,
and 2070–2100), with values ranging from 1 to 3 months(Figs. 12a–c,e–f; compare with Figs. 7a–c,f–h). As expected, thelargest changes occur in the tropical Paciﬁc, where the prominent
duration maximum virtually disappears in the ENSO-neutralcomposites. Secondary duration maxima in regions with known
teleconnections to the tropical Paciﬁc such as the tropicalAtlantic, northeast Paciﬁc, and Amundsen–Bellinghausen Sea(e.g., west of the Antarctic Peninsula) also show pronounced
decreases in the ENSO-neutral composites compared to the
all-sample composites. The resulting patterns of ENSO-neutralcomposite duration are very similar between MHW and MCW
for each time period (cf. Figs. 12a–c,e,f), more so than their
all-sample counterparts (cf. Figs. 7a–c,e,f). While these ENSOneutral patterns do not change appreciably over time, theiramplitudes diminish slightly (note, for example, the reductionin areal extent of values between 2 and 3 months shown in thedarker green shading in Fig. 12).

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1890

JOURNAL OF CLIMATEVOLUME 37

FIG. 9. Fractional area (%) of signiﬁcant positive (red bars) and negative (blue bars) differences between
2070–2100 and 1970–2000 for each model large ensemble for composite (a) MHW intensity, (b) MCW intensity,(c) MHW duration, and (d) MCW duration. Open (ﬁlled) bars are for all (ENSO-neutral) MHW and MCW samples.

Signiﬁcant future changes in MHW and MCW composite
duration during ENSO-neutral states (Figs. 12d,e,i,j) are considerably weaker in magnitude and less widespread comparedto their all-sample counterparts (Figs. 7d,e,i,j). Indeed, signiﬁcant late-century changes in ENSO-neutral MHW and MCW
duration are mainly conﬁned to the Southern Hemisphere,
Arctic Ocean, and portions of the North Atlantic. It is notablethat the pattern of late-century duration change is much more
symmetric between MHW and MCW events in the ENSOneutral composites compared to the all-sample composites,
consistent with our earlier conjecture that changes in ENSO
variability are responsible for asymmetric changes in duration
over the tropical Atlantic. The ENSO-neutral MHW andMCW composites for the six other model LEs are shown in
Figs. S29–S34 for intensity and Figs. S35–S40 for duration.
Next, we summarize the late-century changes in ENSO-neutral
MHW and MCW composite intensity and duration across allseven model LEs. Like CESM2, the other model LEs all show
reduced areal coverage of signiﬁcant future changes in MHW
and MCW composite intensity and duration during ENSOneutral states (Figs. 13 and 14) compared to their all-sample
counterparts (Figs. 8 and 10). For intensity, the areal fractionof signiﬁcant change is reduced by 10%–40% depending on themodel (cf. open and shaded bars in Figs. 9a,b). Duration shows
even more pronounced decreases in areal fraction of signiﬁcantchange, with reductions of 42%–84% depending on the model(cf. open and shaded bars in Figs. 9c,d). Notably, the areal fraction of signiﬁcant future change in duration for ENSO-neutralcomposites falls below 25% for all models except CESM2(shaded bars in Figs. 9c,d). The relative proportion of positive vs.
negative contributions to the areal fraction of signiﬁcant futurechange remains similar for the ENSO-neutral composites compared to the all-sample composites for a given model (Fig. 9).
Finally, the areal fraction of the world oceans with at least twothirds model consensus on the sign of signiﬁcant future change in
MHW and MCW intensity (duration) decreases from 36% to38% (27%–28%) in the all-sample composites to 20%–21%(1%) in the ENSO-neutral composites (Figs. 13h,p and 14h,p). In
terms of amplitude, modest reductions are found between the
ENSO-neutral and all-sample MHW and MCW composites foreach model, both for intensity (cf. Figs. 13 and 8) and duration
outside of the tropical Paciﬁc (cf. Figs. 14 and 10). In summary,future changes in ENSO-related variability account for a considerable fraction of the projected late-century MHW and MCW

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

1891

DESER ET AL.

FIG. 10. As in Fig. 8, but for composite MHW and MCW duration (months).

composite changes in intensity, and nearly all of the projectedchange in duration. Similar conclusions are found for midcentury
changes (Figs. S41 and S42).

d. Mean state versus variability
Up to now, we have focused on the future evolution of MHW
and MCW characteristics due to changes in variability. However,this evolution will be superimposed upon mean state shifts (“arising tide lifts all ships”). Figure 15 shows the late-century meanstate SST change in each model LE computed as the ensemblemean difference between 2070 and 2100 minus 1970 and 2000.
The large-scale patterns of mean state change are generally similar across the models [see Manabe et al. (1991), Xie et al. (2010),
and Armour et al. (2016) for a discussion of pattern formation
mechanisms], but their amplitudes vary considerably due to differences in climate sensitivity and/or radiative forcing scenario(recall Table 1). Of the three CMIP6 models with the same
(SSP5-8.5) radiative forcing scenario, CanESM5 shows the largest change in global-mean SST (4.048C) and MIROC6 shows
the smallest (2.038C). Similarly, of the three CMIP5 models
with the same (RCP8.5) radiative forcing scenario, CanESM2
shows the largest change in global-mean SST and MPI the smallest(3.078 and 1.948C, respectively; Table 1 and Fig. 15). All models
simulate intensiﬁed warming in the equatorial Paciﬁc and greater

SST increases in the Northern Hemisphere compared to theSouthern Hemisphere. However, they differ in their representation of the North Atlantic “warming hole” response and sign of
the SST response in the Weddell and Amundsen–Bellinghausen
Seas (Fig. 15). The amplitude of global-mean SST warming does
not appear to be directly correlated with the magnitude or spatial
extent of signiﬁcant future changes in MHW and MCW intensityand duration across models (cf. Fig. 15 with Figs. 8 and 10).
How large are the relative contributions of changes in variability versus changes in mean state to projected changes in
MHW and MCW characteristics? Figure 16 shows maps ofthe ratio of late-century changes in MHW and MCW composite intensity due to changes in variability divided by those due
to changes in variability-plus-mean state for each model LE(e.g., Fig. 8 divided by the sum of Figs. 8 and 15). The number
at the upper right of each map indicates the fractional area ofthe World Ocean containing ratios within 60.1. As can beseen, the models range from 82% to 91% in this metric for
MHW and from 83% to 90% for MCW. Thus, mean statechanges dominate future changes in MHW and MCW composite
intensity over most of the World Ocean. The Arctic and Antarctic
are exceptions, where ratios reach up to 60.5 for MHW and evenhigher for MCW in many models, indicating that contributions
from changes in variability are comparable to those from changes

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1892

JOURNAL OF CLIMATEVOLUME 37

FIG. 11. As in Fig. 6, but for ENSO-neutral samples.

in the mean state for projected late-century changes in MHW andMCW composite intensity. This is likely due to melting of sea ice
and subsequent exposure of the ocean to atmosphere, which thenallows for variability in heat exchange across the air–sea interface.
Other regions with relatively large contributions from changes
in variability include the subpolar North Atlantic, an area inﬂuenced by sea ice–related changes in density and deep ocean
convection. More modest contributions from changes in variability are found over portions of the tropical Paciﬁc and along
the Kuroshio–Oyashio front in some models, with values up to
about 60.25. The relative contribution of changes in variability tends to be larger for midcentury changes in MHW andMCW composite intensity (Fig. S43) than late-century ones,
likely a consequence of the high-emissions scenarios used toforce the models. In some models, this contribution can reachup to 60.4 in parts of the equatorial Paciﬁc and up to 60.8 in
the Arctic (Fig. S43). The fractional area of the World Ocean
containing ratios within 60.1 is also correspondingly lower for
midcentury changes, ranging from 72% to 86% for MHW and
69%–86% for MCW across models.
To further illustrate the relative contributions from changes
in the mean state versus variability, we show histograms ofarea-averaged iSST from the CESM2 LE for three regions ofinterest: the Arctic (poleward of 678N), the western tropicalPaciﬁc (88S–68N, 1558E–1758W), and the northeast Atlantic(358–628N, 308W–08) in Figs. 17a, 17c, and 17e, respectively.
For each regional average, we consider only those grid boxes
whose late-century MHW composite intensity changes
are statistically signiﬁcant and of the same sign (positive in thecase of the Arctic and northeast Atlantic, and negative in the

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

1893

DESER ET AL.

FIG. 12. As in Fig. 7, but for ENSO-neutral samples.

case of the western tropical Paciﬁc; regional boundaries are outlined in Figs. 6e,j). The distributions shown in Figs. 17a, 17c,
and 17e are based on all monthly iSST samples (not just MHW
and MCW samples) from all 100 ensemble members during
1970–2000 (gray bars) and 2070–2100 (blue bars); the purple bars
are based on the same set of samples as the blue bars, but withthe ensemble-average mean state change (2070–2100 minus
1970–2000) added back in. The 10th and 90th percentiles of each
distribution are shown as vertical solid lines, and the 50th percentile as a vertical dashed line.
In the Arctic, there is a notable widening of the iSST distribution in late century compared to present day (blue vs graybars), and this widening is approximately symmetric withrespect to sign, as measured by the difference between the10th percentiles of the future and present-day distributions(0.158C) and between the 90th percentiles of the two distributions (0.168C) (see Fig. 17a). Adding in the mean state change
shifts the blue distribution by 1.278C, and results in a complete
separation of the present-day (gray bars) and late-century(purple bars) histograms (Fig. 17a). Unlike the Arctic, thepresent-day western tropical Paciﬁc iSST distribution is distinctly skewed, with a longer tail of negative values comparedto positive ones (gray bars in Fig. 17c). We speculate that thisnegative skewness is a reﬂection of the fact that La Niña
events extend farther west than El Niño events, both in natureand the model (Capotondi et al. 2020). This distribution narrows in late century, with a larger change at the lower end
than the upper end (the 10th percentile increases by 0.418C
while the 90th percentile decreases by 0.258C): however, thenegative skewness is maintained (blue bars in Fig. 17c).

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1894

JOURNAL OF CLIMATEVOLUME 37

FIG. 13. As in Fig. 8, but for ENSO-neutral samples.

Adding in the mean state change shifts the blue distribution
by 2.28C (purple bars in Fig. 17c). Unlike the Arctic, the western tropical Paciﬁc distributions show some overlap between
present-day (gray bars) and late-century (purple bars) results
due to the long negative tail in the future histogram. Like theArctic, the northeast Atlantic iSST distributions widen in the
future, and this widening is approximately symmetric withrespect to sign (0.188C for the 10th percentiles and 0.198C
for the 90th percentiles; Fig. 17e). However, due to the pronounced “warming hole” response in this region (recall
Fig. 15a), the mean state change shifts the distribution only
slightly (0.188C), resulting in a large overlap between presentday and late-century (Fig. 17e, gray vs. purple bars). In summary, the three regions show distinctive signatures of futurechange. In the Arctic, projected changes in variability signiﬁcantly widen the SST distribution, while those in the mean
shift the distribution to an entirely new climatic state (e.g.,values that are warmer than any during present day). In thewestern tropical Paciﬁc, projected changes in variability narrow the negatively skewed distribution (and slightly reduce
the skewness), while changes in the mean state shift the distribution to warmer but not completely unprecedented values.
In the northeast Atlantic, the projected widening of the distribution is equal to the magnitude of the mean state shift.

As stated in section 2b, our thresholds for deﬁning MHW
and MCW events (based on the 10th and 90th percentiles of
the iSST distribution in a given time period, respectively) arecomputed for each month separately. Here, we brieﬂy show
the seasonal variation in these thresholds for the three regions
selected above (Figs. 17b,d,f). In the Arctic, the 10th and 90th
percentiles of the iSST distribution during 1970–2000 exhibit
pronounced seasonal cycles, with minimum amplitudes(60.028C) during January–April and largest magnitudes
during June–November (maximum of 60.228C in August);
Fig. 17b (gray bars). The model’s seasonal variation is similar tothat in observations, although the minimum values are underestimated (Fig. 17b, orange bars). This seasonal variation is likely
associated with the amount of open water available for exchanging heat with the atmosphere (i.e., during the months ofJanuary–April when sea ice cover is typically around 100%, thewater directly beneath the ice is buffered from air–sea interaction and hence is able to maintain a relatively constant temperature). A similar seasonal variation in the 10th and 90th
percentiles is found during 2070–2100, but with considerablylarger amplitudes compared to the present day (Fig. 17b, bluebars). For example, winter values increase to around 60.18C
while late summer values reach to nearly 60.58C. There is nooverlap in any month between the present-day percentiles andUnauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

1895

DESER ET AL.

FIG. 14. As in Fig. 10, but for ENSO-neutral samples.

the late-century percentiles that include mean-state changes(Fig. 17b, gray and purple bars).
The present-day MHW and MCW thresholds in the western tropical Paciﬁc also vary seasonally: the 90th percentile
reaches maximum values (1.08–1.18C) from September toMarch and minimum values from May to July (0.68–0.78C); asimilar timing is seen for the 10th percentile but the amplitudes are larger by about 0.48C, reﬂecting the negative skewness in the full iSST distribution noted earlier (Fig. 17d, graybars). While negative skewness is also found in observations,the model overestimates the 10th and 90th percentiles in all
months, consistent with the model biases shown in Figs. 1cand 1i. Unlike the Arctic, the seasonal cycle of the projected
decrease in threshold magnitudes is inversely related to thepresent-day seasonal cycle, with maximum decreases from June
to August (0.78 and 0.458C for the 10th and 90th percentiles, respectively; Fig. 17d, gray and blue bars). Like the Arctic, the present-day percentiles and the late-century percentiles that include
mean-state changes are completely separate in every month(Fig. 17d, gray and purple bars). Compared to the other regions,the northeast Atlantic shows a muted seasonal cycle of futurechange in MHW and MCW threshold amplitudes, ranging from
0.28 to 0.38C (cf. gray and blue bars in Fig. 17f), with considerable overlap between the present-day and future distributions inevery month (cf. gray and purple bars in Fig. 17f). The model realistically simulates present-day thresholds in almost every
month (cf. orange and gray bars in Fig. 17f).
Our MHW and MCW composites convey the average intensity and duration across all samples in a given time period.
In Fig. 18, we show the range of intensity values across all individual MHW and MCW events that make up these composite values for each of the three selected regions. The 5th–95th
percentile range of intensity values for individual MHW
events in the Arctic increases from 0.028 to 0.358C during
1970–2000 to 0.108–0.798C during 2070–2100 (Fig. 18a, grayand blue bars, respectively). The positive skewness in the future intensity distribution of individual MHW events is particularly pronounced. Individual MCW events are largely the
mirror image of MHW events, with notable negative skewness
in the future compared to present day (Fig. 18b, gray and blue
bars, respectively). The present-day and future distributionsof individual MHW (Fig. 18a) and MCW (Fig. 18b) event intensity remain clearly separated even when mean statechanges are included (cf. gray and purple bars).
In the western tropical Paciﬁc, the distribution of individual
MHW event intensity shifts toward lower values in the futurecompared to the present day, even without mean state changes
included (Fig. 18c, gray and blue bars). While this result is

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1896

JOURNAL OF CLIMATEVOLUME 37

FIG. 15. Ensemble-mean SST differences (8C) differences between 2070–2100 and 1970–2000 for each model large ensemble: (a) CESM2,(b) CanESM5, (c) GFDL-SPEAR, (d) MIROC6, (e) CESM1, (f) CanESM2, (g) MPI-ESM-LR, and (h) multimodel ensemble (MME)average. The text in parentheses after each model name indicates the radiative forcing scenario applied, and the number in the upper right ofeach panel denotes the global-mean SST difference.

somewhat counterintuitive, it is consistent with the behavior of
the upper half of the full iSST distribution shown in Fig. 17b. A
similar future shift (toward less negative values) is found for thecorresponding distributions of individual MCW events, with ahint of a bimodal shape in the future distribution (Fig. 18d, grayand blue bars). The present-day and future distributions of individual event intensity remain distinct when mean state changes
are included (gray and purple bars in Figs. 18c,d).
In the northeast Atlantic, the individual MHW event intensity distributions have a similar degree of positive skewness inpresent day and late century (Fig. 18e, gray vs. blue bars).
While the MHW distribution is shifted toward warmer values
in the future, there remains some overlap with the present-day
distribution (Fig. 18e, gray vs. purple bars). Analogous behavior
is seen for the individual MCW event distributions, with an evengreater overlap between future and present-day and negativeskewness (Fig. 18f). In summary, the Arctic, western tropicalPaciﬁc, and northeast Atlantic show diverse behavior in their
MHW and MCW intensity distributions, including shape (degreeof normality), seasonal dependence, and future change as a resultof changes in variability in addition to changes in the mean state.

4. Summary and discussion
We have investigated future changes in the intensity andduration of composite MHWs and MCWs based on monthly
SSTs in seven coupled model initial-condition large ensembles(LEs) over the period 1970–2100 under a range of radiative
forcing scenarios. The large number of simulations (30–100)
available for each LE provides voluminous samples of MHW
and MCW events, allowing for robust quantiﬁcation of theirprojected evolution as a result of changes in variability andcomparison with changes in the mean state in each model. We
ﬁnd that mean state changes are much larger than variabilityrelated changes in MHW and MCW characteristics, consistent

with previous work based on daily SSTs from the CMIP archives (e.g., Frölicher et al. 2018; Oliver 2019; Yao et al. 2022)
and a high-resolution coupled model simulation (Guo et al.
2022). However, we also identify widespread statistically signiﬁcant contributions from projected changes in variability inevery model. These changes in variability generally accountfor ,10% of the total change in intensity by late-century, except over the Arctic and Antarctic where they contribute up
to 50%–75%. By midcentury in some models, changes in variability contribute up to 40% in parts of the equatorial Paciﬁcand up to 80% in the Arctic. For ecosystems that have some
capacity to adapt to the slowly evolving SST mean state,variability-induced changes in MHW/MCW properties may
pose additional challenges in the coming decades (Guo et al.
2022; Oliver 2019; Oliver et al. 2021).
The patterns and magnitudes of projected forced changes
in MHW and MCW intensity and duration due to variability
are complex and model dependent. While it is difﬁcult to generalize, the majority of LEs show that future changes in variability will cause MHW and MCW events to intensify in theextratropics (excluding the Southern Ocean) and weaken in
the tropics and Southern Ocean, although the regional details
and spatial extent of these changes differ across models. The
majority of LEs also show that future changes in variability
will generally cause MHW and MCW events to signiﬁcantly
shorten in duration, although these changes are less geographically widespread than those in intensity. However, there isconsiderable model disagreement in the sign of signiﬁcant duration change over the Arctic, the North and South Atlantic,
and the south Indian Ocean. Future changes in intensity andduration of composite MHWs and MCWs are generally symmetric in all models, except for the longitude of maximum duration change in the tropical Paciﬁc. This asymmetry is likely
related to the different spatial and temporal characteristics of
El Niño and La Niña events (e.g., Okumura and Deser 2010;

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

DESER ET AL.

1897

FIG. 16. Ratio of late-century changes (2070–2100 minus 1970–2000) in MHW and MCW composite intensity due to changes in variability divided by that due to changes in variability-plus-mean state for (a),(h) CESM2; (b),(i) CanESM5; (c),(j) GFDL-SPEAR; (d),(k)MIROC6; (e),(l) CESM1; (f),(m) CanESM2; and (g),(n) MPI-ESM-LR. (h),(p) The multimodel ensemble (MME) average. The numberin the upper right of each panel denotes the fractional area (%) of values within the range from 20.1 to 10.1.

Capotondi et al. 2020). There does not appear to be a direct
correspondence between models’ amplitude of global-mean
SST warming and their magnitude and spatial extent of forcedchanges in MHW and MCW characteristics due to variability.
Projected changes in ENSO variability have a large impact
on future changes in MHWs and MCWs worldwide in every
model. All LEs show reduced areal coverage of signiﬁcant future changes in MHW and MCW intensity and duration during ENSO-neutral states compared to all states (late-century
reductions of 10%–40% for intensity and 42%–84% for duration, depending on the model). Relatedly, the areal fractionof the world oceans with at least two-thirds model consensus
on the sign of signiﬁcant late-century change in MHW and
MCW intensity (duration) decreases from 37% (28%) in the
all-sample composites to 20% (1%) in the ENSO-neutral composites. Thus, future changes in ENSO exert a major control on the
future evolution of MHW and MCW characteristics, underscoring the importance of reducing model spread in ENSO projections, which is considerable and poorly understood (Maher et al.
2023). We note that the collection of models used hereencompasses a range of projected ENSO behavior, from littlechange in the MPI-ESM-LR model to decreasing ENSO variance in CanESM2 to nonmonotonically increasing ENSO variance in the remaining models (Maher et al. 2023). The physicalmechanisms underlying the diversity of projected ENSO behavior across models remain to be elucidated, although we note that
ENSO results from a sensitive balance of competing effects andits response to anthropogenic emissions may likewise depend sensitively on the relative magnitudes and rates of change of different
factors including the Bjerknes feedback, atmospheric damping,
upper ocean stratiﬁcation, and thermocline depth [see Cai et al.
(2021) for a recent review]. Without a sound understanding of
how and why ENSO may change under anthropogenic emissions,this uncertainty is unlikely to be narrowed.
A potential caveat of our ﬁndings is the fact that every
model LE shows statistically signiﬁcant biases in composite
MHW and MCW intensity and duration deﬁned on the basis ofmonthly SSTs over the historical period 1950–2020. These biases
are widespread in every model, covering ;50%–70% of the areaof the global oceans for intensity and ;35%–55% for duration.

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1898

JOURNAL OF CLIMATEVOLUME 37

FIG. 17. Histograms of area-averaged SST (8C) from the CESM2 large ensemble for (a) Arctic (poleward of 678N),(c) western tropical Paciﬁc (88S–68N, 1558E–1758W), and (e) northeast Atlantic (358–628N, 308–08W) based on all
months from all 100 ensemble members during 1970–2000 (gray) and 2070–2100 (blue) after removing the ensemblemean climatological seasonal cycle for each period (see text for details). Purple histograms are the same as the blue
histograms but with the mean state change (2070–2100 minus 1970–2000) added back in. The 10th and 90th percentiles of each distribution are shown as vertical solid lines, and the 50th percentile is shown as a vertical dashed line.
(b),(d),(f) Seasonal variation of the 10th and 90th percentiles of the distributions shown in (a), (c), and (e). Observed
values for the period 1970–2000 are denoted in orange.

The majority of model LEs signiﬁcantly overestimate intensity
in the tropical western Paciﬁc and Indian Oceans and across portions of the Southern Ocean, and signiﬁcantly underestimate it in
the Arctic. The models show less consistency in their durationbiases both in terms of sign and location, although the North
Paciﬁc is a region of general overestimation. The physical mechanisms underlying these model biases remain to be uncovered, but
deﬁciencies in the representation of tropical and extratropicalmodes of variability including ENSO and its teleconnections(Fasullo et al. 2020; Maher et al. 2023) and inadequacies in thesimulation of mean mixed layer depth (Guo et al. 2022) and sea
ice properties (Notz and SIMIP Community 2020) are likely toplay a role. An important aspect of our assessment of model bias
is the consideration of sampling ﬂuctuations inherent in the short(71 year) historical record. All model LEs exhibit substantial ensemble spread in historical MHW/MCW composites, underscoring

the importance of having LEs for robust statistics and cautioning
against overinterpretation of the observational record, which may
be subject to similar levels of uncertainty due to limited sampling.
Another potential caveat of our study is the relatively coarse
spatial resolution (approximately 18–28) of the model LEs, which
limits their ability to represent mesoscale ocean processes andsmall-scale air–sea interactions (Pilo et al. 2019; Hayashida et al.
2020; Guo et al. 2022). In particular, the dynamically downscaled
ocean model experiments of Pilo et al. (2019) and Hayashida et al.
(2020) and the fully coupled high-resolution modeling experimentof Guo et al. (2022) demonstrate that explicit representation ofocean mesoscale processes improves the simulation of historical
MHW characteristics based on daily SSTs along the westernboundary currents, the Antarctic circumpolar current and othereddy-rich regions of the world oceans. Some of the improvements
are also attributable to the better simulation of mixed layer depth

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

1899

DESER ET AL.

FIG. 18. As in Fig. 17, but for (a),(c),(e) MHW and (b),(d),(f) MCW intensity (8C) based on all monthly samples.
The 5th and 95th percentiles of each distribution are shown as vertical solid lines, and the 50th percentile is shown asa vertical dashed line.

(atial resolution affects the simulation of MHW/MCW characteristics based on monthly SSTs remains to be assessed. The enhanced ﬁdelity of the high-resolution models in turn lends
conﬁdence to their future projections, although it does not guarantee that modes of variability such as ENSO and their evolution
under global warming are also better simulated.
Our study leaves many important questions for future investigation as follows. What is the seasonal dependence of future changes in MHW and MCW intensity and duration
across model LEs? Are future changes larger and more robustin summer than winter due to the shallower depth of the upper ocean mixed layer? What role does ENSO play in setting
the seasonality of future changes in MHW and MCW characteristics, and are there any asymmetric inﬂuences between El
Niño and La Niña events? Do more sophisticated approaches
to removing ENSO inﬂuences such as multivariate linear inverse modeling (Zhao et al. 2021) reveal additional insights
into the role of future changes in ENSO on MHWs and
MCWs? Are there additional insights to be gained on the future behavior of MHW and MCW events by using daily data?

What is the subsurface structure of the projected changes inintensity and duration, and does it differ between MHWs and
MCWs? What are the physical mechanisms governing projected changes in MHW and MCW variability in different locations, and do the dominant processes depend on whether
ENSO is playing an active role? Do future changes in largescale modes of atmospheric circulation variability such as
those documented in O’Brien and Deser (2023) affect the future evolution of MHW and MCW characteristics? How doprojected changes in MHW and MCW characteristics found here
relate to future changes in variance and persistence of daily andmonthly SST anomalies documented in Li and Thompson (2021)
and Shi et al. (2022)? We hope to address many of these outstanding questions in future studies by leveraging the unique setof information contained in Earth system model initial-condition
large ensemble archives.
Acknowledgments. We thank the three anonymous reviewers
for their constructive comments and suggestions. We acknowledge
the efforts of all those who contributed to producing the model

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1900

JOURNAL OF CLIMATE

simulations and observational datasets used in this study. A. C.
and A. P were supported by the NOAA Climate Program Ofﬁce
Modeling, Analysis, Predictions and Projections (MAPP) program. NCAR is sponsored by the National Science Foundation
under Cooperative Agreement 1852977.
Data availability statement. The ERSSTv5 dataset is publicly available from www.ncei.noaa.gov/products/extendedreconstructed-sst. All model simulations used in this study are
publicly available from https://www.cesm.ucar.edu/communityprojects/mmlea and https://esgf-node.llnl.gov/search/cmip6/. Allanalyses were performed using NCL; analysis codes are available from Adam S. Phillips upon request.

REFERENCES
Alexander, M. A., I. Bladé, M. Newman, J. R. Lanzante, N.-C.
Lau, and J. D. Scott, 2002: The atmospheric bridge: The inﬂuence of ENSO teleconnections on air–sea interaction over the
global oceans. J. Climate, 15, 2205–2231, https://doi.org/10.1175/
1520-0442(2002)015,2205:TABTIO.2.0.CO;2.
}}, J. D. Scott, K. D. Friedland, K. E. Mills, J. A. Nye, A. J. Pershing, and A. C. Thomas, 2018: Projected sea surface temperatures over the 21st century: Changes in the mean, variability
and extremes for large marine ecosystem regions of northernoceans. Elementa, 6, 9, https://doi.org/10.1525/elementa.191.
Amaya, D. J., M. A. Alexander, A. Capotondi, C. Deser, K. B.
Karnauskas, A. J. Miller, and N. J. Mantua, 2021: Are longterm changes in mixed layer depth inﬂuencing North Paciﬁcmarine heatwaves? [in “Explaining Extremes of 2019 from a
Climate Perspective”] Bull. Amer. Meteor. Soc., 102 (1), S59–
S66, https://doi.org/10.1175/BAMS-D-20-0144.1.
}}, M. G. Jacox, M. A. Alexander, J. D. Scott, C. Deser, A. Capotondi, and A. S. Phillips, 2023a: Bottom marine heatwaves
along the continental shelves of North America. Nat. Commun., 14, 1038, https://doi.org/10.1038/s41467-023-36567-0.
}}, and Coauthors, 2023b: Marine heatwaves need clear deﬁnitions so coastal communities can adapt. Nature, 616, 29–32,
https://doi.org/10.1038/d41586-023-00924-2.
Armour, K., Marshall, J., Scott, J., A. Donohoe, and E. Newsom,
2016: Southern Ocean warming delayed by circumpolar upwelling and equatorward transport. Nat. Geosci., 9, 549–554,
https://doi.org/10.1038/ngeo2731.
Burger, F. A., J. Terhaar, and T. L. Frölicher, 2022: Compound
marine heatwaves and ocean acidity extremes. Nat. Commun., 13, 4722, https://doi.org/10.1038/s41467-022-32120-7.
Cai, W., and Coauthors, 2021: Changing El Niño–Southern Oscillation in a warming climate. Nat. Rev. Earth Environ., 2, 628–
644, https://doi.org/10.1038/s43017-021-00199-z.
Capotondi, A., M. A. Alexander, N. A. Bond, E. N. Curchitser,
and J. D. Scott, 2012: Enhanced upper ocean stratiﬁcation
with climate change in the CMIP3 models. J. Geophys. Res.,
117, C04031, https://doi.org/10.1029/2011JC007409.
}}, C. Deser, A. S. Phillips, Y. Okumura, and S. M. Larson,
2020: ENSO and Paciﬁc decadal variability in the Community
Earth System Model version 2. J. Adv. Model. Earth Syst.,
12, e2019MS002022, https://doi.org/10.1029/2019MS002022.
}}, M. Newman, T. Xu, and E. Di Lorenzo, 2022: An optimal
precursor of northeast Paciﬁc marine heatwaves and central Paciﬁc El Niño events. Geophys. Res. Lett., 49,
e2021GL097350, https://doi.org/10.1029/2021GL097350.

VOLUME 37

Cheung, W. W. L., and T. L. Frölicher, 2020: Marine heatwaves
exacerbate climate change impacts for ﬁsheries in the northeast Paciﬁc. Sci. Rep., 10, 6678, https://doi.org/10.1038/s41598020-63650-z.
Danabasoglu, G., and Coauthors, 2020: The Community Earth System Model version 2 (CESM2). J. Adv. Model. Earth Syst., 12,e2019MS001916, https://doi.org/10.1029/2019MS001916.
Delworth, T. L., and, Coauthors, 2020: SPEAR: The next generation
GFDL modeling system for seasonal to multidecadal prediction
and projection. J. Adv. Model. Earth Syst., 12, e2019MS001895,
https://doi.org/10.1029/2019MS001895.
Deser, C., and Coauthors, 2020: Insights from Earth system model
initial-condition large ensembles and future prospects. Nat.
Climate Change, 10, 277–286, https://doi.org/10.1038/s41558020-0731-2.
Fasullo, J., A. S. Phillips, and C. Deser, 2020: Evaluation of leading modes of climate variability in the CMIP archives. J. Climate, 33, 5527–5545, https://doi.org/10.1175/JCLI-D-19-1024.1.
Frölicher, T. L., and C. Laufkötter, 2018: Emerging risks from marine heat waves. Nat. Commun., 9, 650, https://doi.org/10.
1038/s41467-018-03163-6.
}}, E. M. Fischer, and N. Gruber, 2018: Marine heatwaves under global warming. Nature, 560, 360–364, https://doi.org/10.
1038/s41586-018-0383-9.
Guo, X., and Coauthors, 2022: Threat by marine heatwaves toadaptive large marine ecosystems in an eddy-resolving model.
Nat. Climate Change, 12, 179–186, https://doi.org/10.1038/
s41558-021-01266-5.
Hayashida, H., R. J. Matear, P. G. Strutton, and X. Zhang, 2020:
Insights into projected changes in marine heatwaves from ahigh-resolution ocean circulation model. Nat. Commun., 11,4352, https://doi.org/10.1038/s41467-020-18241-x.
Hobday, A. J., and Coauthors, 2016: A hierarchical approach todeﬁning marine heatwaves. Prog. Oceanogr., 141, 227–238,
https://doi.org/10.1016/j.pocean.2015.12.014.
Holbrook, N. J., and Coauthors, 2019: A global assessment of marine heatwaves and their drivers. Nat. Commun., 10, 2624,
https://doi.org/10.1038/s41467-019-10206-z.
Huang, B., C. Lsiu, V. Banzon, E. Freeman, G. Graham, B. Hankins, T. Smith, and H.-M. Zhang, 2021: Improvements of the
Daily Optimum Interpolation Sea Surface Temperature(DOISST) version 2.1. J. Climate, 34, 2923–2939, https://doi.
org/10.1175/JCLI-D-20-0166.1.
Jacox, M. G., M. A. Alexander, S. J. Bograd, and J. D. Scott,2020: Thermal displacement by marine heatwaves. Nature,584, 82–86, https://doi.org/10.1038/s41586-020-2534-z.
Kay, J. E., and Coauthors, 2015: The Community Earth System
Model (CESM) large ensemble project: A community resource for studying climate change in the presence of internal
climate variability. Bull. Amer. Meteor. Soc., 96, 1333–1349,
https://doi.org/10.1175/BAMS-D-13-00255.1.
Kirchmeier-Young, M. C., F. W. Zwiers, and N. P. Gillett, 2017:
Attribution of extreme events in Arctic sea ice extent. J. Climate, 30, 553–571, https://doi.org/10.1175/JCLI-D-16-0412.1.
Kwiatkowski, L., and Coauthors, 2020: Twenty-ﬁrst century ocean
warming, acidiﬁcation, deoxygenation, and upper-ocean nutrient and primary production decline from CMIP6 model
projections. Biogeosciences, 17, 3439–3470, https://doi.org/10.
5194/bg-17-3439-2020.
Laufkötter, C., J. Zscheischler, and T. L. Frölicher, 2020: High-impact
marine heatwaves attributable to human-induced global warming.
Science, 369, 1621–1625, https://doi.org/10.1126/science.aba0690.

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

15 MARCH 2024

DESER ET AL.

Le Grix, N., J. Zscheischler, K. B. Rodgers, R. Yamaguchi, and T. L.
Frölicher, 2022: Hotspots and drivers of compound marine heatwaves and low net primary production extremes. Biogeosciences,19, 5807–5835, https://doi.org/10.5194/bg-19-5807-2022.
Li, J., and D. W. J. Thompson, 2021: Widespread changes in surface temperature persistence under climate change. Nature,599, 425–430, https://doi.org/10.1038/s41586-021-03943-z.
Maher, N., and Coauthors, 2019: The Max Planck institute grandensemble: Enabling the exploration of climate system variability. J. Adv. Model. Earth Syst., 11, 2050–2069, https://doi.
org/10.1029/2019MS001639.
}}, and Coauthors, 2023: The future of the El Niño–Southern
Oscillation: Using large ensembles to illuminate time-varying
responses and inter-model differences. Earth Syst. Dyn., 14,
413–431, https://doi.org/10.5194/esd-14-413-2023.
Manabe, S., R. J. Stouffer, M. J. Spelman, and K. Bryan, 1991:
Transient responses of a coupled ocean–atmosphere model
to gradual changes of atmospheric CO2. Part I. Annual mean
response. J. Climate, 4, 785–818, https://doi.org/10.1175/15200442(1991)004,0785:TROACO.2.0.CO;2.
Newman, M., and Coauthors, 2016: The Paciﬁc decadal oscillation, revisited. J. Climate, 29, 4399–4427, https://doi.org/10.
1175/JCLI-D-15-0508.1.
NOAA National Centers for Environmental Information, 2023:
State of the climate: Global climate report for 2022. Accessed
21 December 2023, https://www.ncei.noaa.gov/access/monitoring/monthly-report/global/202213.
Notz, D., and SIMIP Community, 2020: Arctic sea ice in CMIP6.
Geophys. Res. Lett., 47, e2019GL086749, https://doi.org/10.
1029/2019GL086749.
O’Brien, J. P. and C. Deser, 2023: Quantifying and understanding
forced changes to unforced modes of atmospheric circulation
variability over the North Paciﬁc in a coupled model large
ensemble. J. Climate, 36, 17–35, https://doi.org/10.1175/JCLID-22-0101.1.
Okumura, Y. M., and C. Deser, 2010: Asymmetry in the durationof El Niño and La Niña. J. Climate, 23, 5826–5843, https://doi.
org/10.1175/2010JCLI3592.1.
Oliver, E. C. J., 2019: Mean warming not variability drives marine
heatwave trends. Climate Dyn., 53, 1653–1659, https://doi.org/
10.1007/s00382-019-04707-2.
}}, and Coauthors, 2018: Longer and more frequent marine
heatwaves over the past century. Nat. Commun., 9, 1324,
https://doi.org/10.1038/s41467-018-03732-9.
}}, J. A. Benthuysen, S. Darmaraki, M. G. Donat, A. J. Hobday, N. J. Holbrook, R. W. Schlegel, and A. Sen Gupta,
2021: Marine heatwaves. Annu. Rev. Mar. Sci., 13, 313–342,
https://doi.org/10.1146/annurev-marine-032720-095144.
Pilo, G. S., N. J. Holbrook, A. E. Kiss, and A. M. Hogg, 2019:
Sensitivity of marine heatwave metrics to ocean model resolution. Geophys. Res. Lett., 46, 14 604–14 612, https://doi.org/
10.1029/2019GL084928.
Plecha, S. M., and P. M. M. Soares, 2020: Global marine heatwave
events using the new CMIP6 multi-model ensemble: From
shortcomings in present climate to future projections. Environ. Res. Lett., 15, 124058, https://doi.org/10.1088/1748-9326/
abc847.
Qiu, Z., F. Qiao, C. Jang, L. Zhang and Z. Song, 2021: Evaluationand projection of global marine heatwaves based on CMIP6
models. Deep-Sea Res. II, 194, 104998, https://doi.org/10.1016/
j.dsr2.2021.104998.
Rayner, N., D. E. Parker, E. B. Horton, C. K. Folland, L. V. Alexander, D. P. Rowell, E. C. Kent, and A. Kaplan, 2003: Global

1901

analyses of sea surface temperature, sea ice and night marine
air temperatures since the late nineteenth century. J. Geophys.
Res., 108, 4407, https://doi.org/10.1029/2002JD002670.
Rodgers, K. B., and Coauthors, 2021: Ubiquity of human-inducedchanges in climate variability. Earth Syst. Dyn., 12, 1393–
1411, https://doi.org/10.5194/esd-12-1393-2021.
Scannell, H. A., G. C. Johnson, L. Thompson, J. M. Lyman, andS. C. Riser, 2020: Subsurface evolution and persistence of marine heatwaves in the northeast Paciﬁc. Geophys. Res. Lett.,
47, e2020GL090548, https://doi.org/10.1029/2020GL090548.
Schaeffer, A., A. Sen Gupta, and M. Roughan, 2023: Seasonal stratiﬁcation and complex local dynamics control the sub-surface
structure of marine heatwaves in eastern Australian coastal waters. Commun. Earth Environ., 4, 304, https://doi.org/10.1038/
s43247-023-00966-4.
Schlegel, R. W., S. Darmaraki, J. A. Benthuysen, K. Filbee-Dexter,and E. C. J. Oliver, 2021: Marine cold-spells. Prog. Oceanogr.,
198, 102684, https://doi.org/10.1016/j.pocean.2021.102684.
Sen Gupta, A., and Coauthors, 2020: Drivers and impacts of themost extreme marine heatwaves events. Sci. Rep., 10, 19359,
https://doi.org/10.1038/s41598-020-75445-3.
Shi, H., and Coauthors, 2022: Global decline in ocean memory
over the 21st century. Sci. Adv., 8, eabm3468, https://doi.org/
10.1126/sciadv.abm3468.
Smale, D. A., and Coauthors, 2019: Marine heatwaves threaten
global biodiversity and the provision of ecosystem services.
Nat. Climate Change, 9, 306–312, https://doi.org/10.1038/s41558019-0412-1.
Smith, K. E., and Coauthors, 2023: Biological impacts of marine
heatwaves. Annu. Rev. Mar. Sci., 15, 119–145, https://doi.org/
10.1146/annurev-marine-032122-121437.
Student, 1917: Tables for estimating the probability that the meanof a unique sample of observations lies between 2‘ and any
given distance of the mean of the population from which thesample is drawn. Biometrika, 11, 414–417, https://doi.org/10.
1093/biomet/11.4.414.
Swart, N. C., and Coauthors, 2019: The Canadian Earth System
Model version 5 (CanESM5.0.3). Geosci. Model Dev., 12,4823–4873, https://doi.org/10.5194/gmd-12-4823-2019.
Tatebe, H., and Coauthors, 2019: Description and basic evaluationof simulated mean state, internal variability, and climate sensitivity in MIROC6. Geosci. Model Dev., 12, 2727–2765,
https://doi.org/10.5194/gmd-12-2727-2019.
Thoral, F., S. Montie, M. S. Thomsen, L. W. Tait, M. H. Pinkerton,and D. R. Schiel, 2022: Unravelling seasonal trends in coastal
marine heatwave metrics across global biogeographicalrealms. Sci. Rep., 12, 7740, https://doi.org/10.1038/s41598-02211908-z.
Vogt, L., F. A. Burger, S. M. Grifﬁes, and T. L. Frölicher, 2022:
Local drivers of marine heatwaves: A global analysis with an
Earth system model. Front. Climate, 4, 847995, https://doi.org/
10.3389/fclim.2022.847995.
Wang, Y., J. B. Kajtar, L. V. Alexander, G. S. Pilo, and N. J.
Holbrook, 2022: Understanding the changing nature ofmarine cold-spells. Geophys. Res. Lett., 49, e2021GL097002,
https://doi.org/10.1029/2021GL097002.
Wilks, D. S., 2016: “The stippling shows statistically signiﬁcantgrid points”: How research results are routinely overstatedand overinterpreted, and what to do about it. Bull. Amer.
Meteor. Soc., 97, 2263–2273, https://doi.org/10.1175/BAMS-D15-00267.1.

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC

1902

JOURNAL OF CLIMATE

Wu, L., and Coauthors, 2012: Enhanced warming over the global
subtropical western boundary currents. Nat. Climate Change,
2, 161–166, https://doi.org/10.1038/nclimate1353.
Xie, S.-P., C. Deser, G. A. Vecchi, J. Ma, H. Teng, and A. T. Wittenberg, 2010: Global warming pattern formation: Sea surface
temperature and rainfall. J. Climate, 23, 966–986, https://doi.
org/10.1175/2009JCLI3329.1.
Xu, T., M. Newman, A. Capotondi, S. Stevenson, E. Di Lorenzo,
and M. A. Alexander, 2022: An increase in marine heatwaves
without signiﬁcant changes in surface ocean temperature

VOLUME 37

variability. Nat. Commun., 13, 7396, https://doi.org/10.1038/
s41467-022-34934-x.
Yao, Y., C. Wang, and Y. Fu, 2022: Global marine heatwaves
and cold-spells in present climate to future projections.
Earth’s Future, 10, e2022EF002787, https://doi.org/10.
1029/2022EF002787.
Zhao, Y., M. Newman, A. Capotondi, E. Di Lorenzo, and D.
Sun, 2021: Removing the effects of tropical dynamics from
North Paciﬁc climate variability. J. Climate, 34, 9249–9265,
https://doi.org/10.1175/JCLI-D-21-0344.1.

Unauthenticated | Downloaded 02/18/24 08:19 PM UTC