Changing ocean seasonal cycle escalates destructive marine heatwaves in a warming climate

Marine heatwaves (MHWs) can cause various adverse effects on marine ecosystems associated with complicated social ramifications. It has been well established that the gradually rising sea surface temperature (SST) due to anthropogenic carbon emission will cause an increase of the MHW duration and intensity. However, for species with strong adaptation capacity or mobility, MHW changes due to the altered SST variability under greenhouse warming are more crucial but so far remain poorly assessed. Under the high carbon emission scenario, we show that the cumulative duration (intensity) of MHWs, with the effect of secular SST increase excluded, is projected to be 60% (100%) higher by the end of this century than in the 1990s due to an amplified SST seasonal cycle. This increase becomes more evident for stronger MHWs, reaching up to 8 (30) folds for the extreme MHW category. The amplified SST seasonal cycle also causes pronounced seasonality of MHWs, making them more active in summer-autumn than winter-spring. Our results suggest that MHWs are likely to have increasingly devastating impacts on a wide range of marine species in the future without taking effective steps for carbon emission reduction.


Introduction
Both satellite and in-situ observations have consistently revealed a significant mean sea surface temperature (SST) increase over the past three decades in many parts of the global ocean due to the rising greenhouse gases (Rogelj et al 2012, Cheng et al 2019. If the baseline for defining marine heatwaves (MHWs) (Hobday et al 2016, Oliver et al 2020 is chosen as a fixed historical period as many studies did (Frölicher et al 2018, Oliver et al 2018a, Pilo et al 2019, Hayashida et al 2020, this rising mean SST exerts cascading impacts on the likelihood, duration and intensity of MHWs. In particular, existing literature indicates that MHWs have intensified and become longer in the past three decades mainly ascribed to the mean SST increase (Oliver et al 2018a). Such trends in the MHW statistics are projected to become more prominent in the future under the high carbon emission scenario because of the accelerated sea surface warming (Frölicher et al 2018, Hayashida et al 2020. By the end of this century, many parts of the global ocean are projected to reach a permanent MHW state (Frölicher et al 2018, Hayashida et al 2020.
Although the mean SST increase plays a dominant role in the MHW changes under greenhouse warming, understanding MHW changes caused by the changing SST variability in the future is important from dynamical and biological concerns. Dynamically, an MHW, as an anomalous event, should be objectively determined relative to a contemporaneous equilibrium position that also shifts as the mean SST rises (Jacox 2019, Jacox et al 2020. Furthermore, MHW changes caused by the mean SST increase and changing SST variability could have different biological impacts (Pinsky et al 2013, Habary et al 2017, Jin and Agustí 2018, given the great diversity of adaptation capacity or mobility among marine species. For instance, mobile species can migrate gradually in response to the increasing mean SST at a centennial scale and find their suitable alternative habitats in the future with mean SST similar to the present condition, but may be harder to do so for SST changes at much shorter time scales (Jacox et al 2020, Payne 2020. Therefore, for species that adapt or move fast, MHW changes caused by changing SST variability are more relevant , 2020, Jacox et al 2020. So far, it remains poorly assessed how MHWs will respond to the anthropogenic changes of SST variability. Here we address this important issue based on a suite of coupled global climate model (CGCM) simulations (Taylor et al 2012, Eyring et al 2016.

Observation products
An observation-based SST product, National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation Sea Surface Temperature V2.0 (OIS-STv2) (Reynolds et al 2007, Banzon et al 2016, is used as a benchmark to assess the performance of CGCMs in simulating MHWs. This daily SST dataset has a spatial resolution of 0.25 • and covers a period from 1 September 1981 to the present. In addition, to diagnose the CGCMs' performance in simulating the historical SST trend, Hadley Centre Sea Ice and Sea Surface Temperature dataset version1.1 (HadISST v1.1) (Rayner et al 2003) is used, which has a spatial resolution of 1 • and covers 1870-present.

CMIP5/6 data
A suite of CGCM simulations during 1850-2100 in CMIP5 (Taylor et al 2012) and CMIP6 (Eyring et al 2016) are used to assess the MHW changes under the high carbon emission scenario (RCP8.5 scenario in CMIP5 and SSP5-8.5 scenario in CMIP6) caused by changing SST variability. Two selection criteria are applied. First, only CGCMs with daily SST available are retained, guaranteeing an accurate calculation of MHWs. Second, CGCMs showing poor performance in simulating the historical SST trend during 1900-2020 or MHW statistics during 1982-2020 are discarded, as these CGCMs may be problematic in representing MHWs and their response to anthropogenic climate changes. Two widely used statistics for MHWs are considered here. The first one is the annual cumulative duration (denoted as CD henceforth) or equivalently the annual MHW days. It is the product of the annual mean duration and frequency of MHWs. We do not consider the latter two quantities separately because they would be strongly constrained by each other when a substantial fraction of a year is occupied by MHWs. The second statistic is the annual cumulative intensity of MHWs (denoted as CI henceforth). We use CI instead of annual mean intensity as the former is a better measurement of environmental stress of MHWs on marine ecosystems (Frölicher et al 2018, Oliver et al 2018b, 2020, Filbee-Dexter et al 2020, Gupta et al 2020, Hayashida et al 2020. On the one hand, the CGCM's performance in simulating the historical SST trend is measured as the regression coefficient β of the simulated global mean SST time series onto the observed one derived from the HadISST. On the other hand, the CGCM's performance in simulating MHWs is measured as the difference of the global mean MHW statistics between the simulations and OISSTv2. CGCMs with β < 0.6 or difference of MHW statistics larger than its interannual variation in the observation (figure S1 available online at stacks.iop.org/ERL/17/054024/mmedia) are discarded. Based on the above two criteria, 15 CGCMs in CMIP5 and 19 CGCMs in CMIP6 are selected for MHW analysis (table S1). The MHW statistics for each CGCMs are calculated on its native grid and interpolated to a regular 1 • by 1 • grid when calculating the multi-model average.

CESM simulation
The limited model output in CMIP5/6 simulations makes it difficult to uncover the dynamical causes of changing SST variability under greenhouse warming. For this sake, we run a long-term CGCM simulation that outputs the diagnostic terms in the temperature governing equation to facilitate an SST budget analysis. This simulation uses the Community Earth System Model (CESM) version 1.3, of which the atmosphere component is Community Atmosphere Model version 5 (CAM5), the ocean component is Parallel Ocean Program version 2 (POP2), the ice component is the Community Ice Code version 4, and the land component is Community Land Model version 4. The nominal horizontal resolutions of CAM5 and POP2 are 1 • . The Gent-McWilliams scheme (Gent and Mcwilliams 1990), FK08 scheme (Fox-Kemper et al 2008) and KPP scheme (Large et al 1994) are adopted to parameterize the heat transport by mesoscale eddies, submesoscale eddies and microscale turbulence, respectively. The CESM simulation consists of a 500 year pre-industry control experiment and a 250 year historical and future transient climate experiment, following the design protocol of the CMIP5 experiments (Taylor et al 2012). The historical experiment uses time-varying historical forcing from 1850 to 2005. The future experiment uses external forcing following the RCP8.5 future scenario for the period 2006-2100. More details can be found in Chang et al (2020).

Definition of MHWs
In this study, MHWs are identified based on the approach proposed by Hobday et al (2016) with some modifications to exclude the effects of mean SST increase on MHWs. The CGCM simulations in CMIP5 and CMIP6 cover 1850-2100 during which period the SST trend is highly nonlinear. For this reason, rather than applying a linear detrending to the entire SST time series, we adopt a 31 year moving linear detrending based on the premise that the SST trend can be approximated as being linear within any 31 year period. In order to compute the seasonally varying 90th-percentile threshold θ, a linear trend is fitted and subtracted from the SST series in a prescribed 31 year baseline period  at each grid cell. By detrending, we also remove the climatological mean SST within the baseline period. Then, for any year in question Y c , a linear trend is fitted based on the SST series within the 31 year period centered on Y c . Finally, MHWs within Y c are identified based on the criterion (Hobday et al 2016) that the detrended SST anomaly (SSTA) exceeds the value of seasonally varying θ for at least five consecutive days or a longer period. Once individual MHWs within Y c are identified, the CD is computed as the total number of MHW days and the CI is computed as the integral of MHW intensity over all the MHW days.
To assess projected changes of MHWs in the 21st century, we compute the values of CD and CI during 2000-2085 with the last 15 years of the 21st century not available due to the boundary effect of the moving detrending. It should be noted that choosing a different baseline period may alter the absolute values of CD and CI but will not affect their long-term trends, as long as the baseline period does not overlap with 2000-2085 (Zhang et al 2005).
As to the validation of CGCM-simulated MHWs against observations, the above moving linear detrending procedure is not applicable given the short observational SST records. Instead, we linearly detrend the SST time series in CGCM simulations and observations during 1982-2020. Then MHWs are computed based on the SSTA with the baseline period chosen as 1982-2020.

SST budget
Following previous studies (Lee et al 2004), the vertical mean temperature in the upper 20 m is used as a proxy for SST. Its controlling equation is derived as: (1) where ⟨. . .⟩ denotes a vertical average in the upper 20 m, ADV = − ⟨∇ · (uT)⟩ is the advection by model-resolved flows with T the temperature, u the three-dimensional velocity vector and ∇ = (∂/∂x, is the effect of air-sea heat exchange with Q net (0) the net downward heat flux at the sea surface, Q sw (−h) the shortwave radiation flux penetrating out of the layer base, ρ 0 the ocean reference density, C p the ocean specific heat capacity and h = 20 m, EDDY is the parameterized heat flux convergence by ocean mesoscale eddies (Gent and Mcwilliams 1990) and submesoscale eddies (Fox-Kemper et al 2008), and MIX is the parameterized vertical mixing by microscale turbulence (Large et al 1994). It should be noted that the change of SST due to snow and sea ice melting effects is also parameterized in SHF. But such effects are generally negligible equatorward of 60 • . The anomalous SST budget for each month under greenhouse warming is obtained by subtracting the mean SST budget during 1990-2020 in this month from its counterpart during 2070-2100. Then the annual mean value over the 12 months is further subtracted from each term in the anomalous SST budget to remove the effect of mean SST increase. Finally, integrating the resultant anomalous SST budget with time converts the effects of ADV, SHF, EDDY, MIX on SSTA tendency to their effects on SSTA itself (denoted as ADV-T, SHF-T, EDDY-T, and MIX-T, respectively). In this study, the winter-spring and summer-autumn in the northern hemisphere correspond to February-April and July-September, respectively. The opposite is true for the southern hemisphere.
Existing studies suggest that the seasonally varying mixed layer depth may account for the seasonal difference in SST trends in some parts of the ocean (Alexander et al 2018). Although the mixed layer depth does not explicitly appear in equation (1), its effects on the SST changes are included in the MIX term. A deeper mixed layer generally corresponds to a stronger turbulent vertical mixing and vice versa.

Decomposition of SSTA into seasonal and nonseasonal components
The seasonal cycle of SSTA (T ′ SE ) in any year in question Y c is first calculated for each day of the year by averaging the corresponding SSTA records during the 31 year period centered on Y c , and then smoothed by applying a 31 day moving average. Similarly, the climatological mean seasonal cycle of SSTA during the baseline period 1950-1980 is derived by averaging SSTA records during 1950-1980 for each day of the year, followed by a 31 day moving average. Once T ′ SE is obtained, the nonseasonal component T ′ NS is computed by subtracting T ′ SE from SSTA.

Projected MHW changes under the greenhouse warming
To quantify changes of CD and CI caused by changing SST variability under greenhouse warming, we introduce the absolute and relative increasing rates. The absolute increasing rate is derived from the linear fit of CD or CI time series during 2000-2085. The relative increasing rate is defined as the absolute increasing rate divided by the mean CD or CI during 1990-1999. The global mean absolute/relative increasing rates are derived from the globally averaged time series of MHW statistics instead of averaging the increasing rates at individual grids. The simulated mean CD and CI during 1982-2020 averaged over all the CGCM simulations show general agreement with the observed ones (figures S1 and S2), although there are quantitative discrepancies, especially in the equatorial regions partially due to the model deficiencies in simulating the El Niño-Southern Oscillation (Oliver et al 2018a. In terms of the global average, the simulated MHWs are generally weaker and longer than the observations (figure S1), as already reported by previous studies (Pilo et al 2019). But the difference between the simulated and observed CD (CI) values is less than 12% (18%) and much smaller than the interannual variation of CD (CI) in the observations, suggesting that the simulated MHWs are qualitatively reliable in the sense of their global mean features.
The projected CD exhibits significant increases over 54.7% area of the global ocean (figures 1(a) and (c)). The absolute and relative increasing rates are spatially inhomogeneous and generally more evident in the mid-and high-latitude than tropical regions and in the northern than the southern hemisphere. Especially, the largest increasing rates occur in the Arctic Ocean, where MHWs are projected to occupy more than half a year by the end of this century. Averaged over the global ocean, the value of CD at the end of this century is projected to be 60% (15.6 d) larger than its mean value during 1990-1999 ( figure 2(a)). We remark that the increasing rates of globally averaged CD are qualitatively similar among individual CGCMs with a small inter-model spread. Therefore, these signals are unlikely to be an artifact due to model deficiencies.
The spatial pattern of projected CI changes resembles that of CD, with 66.7% area of the global ocean exhibiting significant increases (figures 1(b) and (d)). But its relative increasing rate is systematically larger than that of CD, suggesting that the annual mean intensity of MHWs is also projected to become higher in response to the anthropogenic climate changes even after removing the effect of mean SST increase. By the end of this century, the globally averaged CI is projected to increase by 100% (27.3 • C·day) compared to its mean value during 1990-1999 (figure 2(e)). Again, different CGCMs consistently project the increasing rates of globally averaged CI. Therefore, we conclude that MHWs will occupy a longer period and exert larger environmental stress under greenhouse warming even if the effect of mean SST increase is excluded.
Not all MHWs exert devastating effects on marine ecosystems. Stronger MHWs are more likely to induce acute biological and economic impacts. For this reason, we classify MHWs into different intensity categories (Hobday et al 2018), that is, moderate, strong, severe and extreme events, and evaluate their projected changes separately. As illustrated in figure 2, the relative increasing rates of CD and CI become progressively larger as the intensity category rises. For MHWs categorized as strong or higher-level events, the relative increasing rate of CD (CI) averaged over the global ocean is 240% (410%) per century, about four folds of the value with all-category MHWs counted. The value is further increased to 800% (3000%) for the extreme MHWs. In addition, with ongoing anthropogenic climate changes, the fraction of the ocean experiencing an extreme MHW at least once per two decades is projected to increase (figure S3). The fraction value increases from 13.6% during 1990-2020 to 27.0% during 2055-2085. This expansion of the extreme-MHW affected regions would cause a severe threat to the large marine ecosystems (LMEs) (Alexander et al 2018, Nielsen et al 2021. During 2055-2085, the region experiencing 5% or higher annual extreme-MHW occurrence frequency is projected to cover a substantial fraction of LMEs around the subarctic gyres in the North Pacific Ocean and North Atlantic Ocean, and the major fishing grounds of the world (figure S3). In contrast, the extreme MHWs are typically centennial or rarer events in these regions during 1990-2020.

Projected MHW changes caused by amplified SST seasonal cycle
Our analysis has revealed that the changing SST variability under greenhouse warming causes significant increases of CD and CI over a large fraction of the global ocean. The changing SST variability can be attributed to changes of either standard deviation (σ SSTA ) or skewness (b SSTA ) of SSTA (Oliver et al 2018a(Oliver et al , 2020. We find that it is the former that accounts primarily for the projected increases of CD and CI (figure 3). The projected σ SSTA exhibits a significant positive trend over many parts of the global ocean (figures 3(a) and (b)). The spatial distribution of the relative increasing rate of σ SSTA agrees remarkably well with those of CD and CI, with a pattern correlation coefficient being 0.94 and 0.93, respectively. In contrast, the value of b SSTA decreases with time ( figure 3(c)), acting to weaken MHWs and partially offsetting the effect of increasing σ SSTA on MHW changes. The negative trend of b SSTA is most   1990-1999. pronounced in the Arctic Ocean and marginal ice zone along the Antarctic continent. Nevertheless, its effect on MHWs in these regions is still outweighed by that of increasing σ SSTA (figure 1).
In terms of the global average, the dominant role of increasing σ SSTA on MHW changes becomes more evident. The projected global mean CD and CI increase monotonically with the increasing σ SSTA (figure S4). We remark that CD and CI increase at a higher rate than σ SSTA . For instance, a 12% increase of σ SSTA corresponds to a 50% and 98% increase of CD and CI, respectively. These values are further increased to 14 and 40 folds when only the extreme MHWs are considered. It thus suggests that even a small-to-moderate increase of σ SSTA under the anthropogenic climate changes would cause a pronounced increase of most destructive MHWs.
To uncover the causes for the increasing σ SSTA under greenhouse warming, we decompose the SSTA into seasonal and nonseasonal components (T ′ SE and T ′ NE ), respectively. Over most parts of the global ocean, the projected σ SSTA increase is mainly ascribed to T ′ SE , whereas T ′ NE makes a negligible contribution (figures 3(a) and (d)). The only exception occurs in the Arctic Ocean where T ′ NE accounts for about one-third of the increase of σ SSTA under greenhouse warming.
To examine whether the fundamental contribution of the changing T ′ SE to the increasing σ SSTA under greenhouse warming can be extended to the MHW changes, we recalculate MHWs by referencing MHWs at the year in question to the contemporaneous SST seasonal cycle instead of the climatological mean seasonal cycle during the baseline period. In this case, the projected MHW changes only depend on the changing T ′ NE . Once the effect of changing T ′ SE is eliminated, the projected increases of CD and CI are substantially reduced almost everywhere (figures 4(a) and (b)). In terms of the global average, the relative increasing rate of CD is reduced from 60% to 20% per century without the effects of changing T ′ SE (figures 2(a) and 4(c)). This reduction becomes even more evident for the relative increasing rate of CI, that is, from 100% to 15% per century (figures 2(e) and 4(d)). Again, such features are consistently simulated by individual CGCMs. Therefore, we conclude that the increased standard deviation of T ′ SE plays a dominant role in the projected increases of CD and CI in the global ocean under the anthropogenic climate changes.
The increasing standard deviation of T ′ SE under greenhouse warming suggests an amplified SST seasonal cycle. Consistent with the existing literature (Gupta et al 2015, Alexander et al 2018, the latter is found to be due to the faster SST increase in summer-autumn than winter-spring (figure S5). Especially in the mid-and high-latitude North Pacific and Atlantic regions collocated with the pronounced increase of CD and CI, the change of SST seasonal cycle amplitude is even comparable to the mean SST increase (figures 1 and S5). The seasonally differed SST increasing rates could originate from either airsea heat exchange or heat redistribution through internal ocean processes. To shed light on its underlying dynamics, an SST budget analysis is performed using the diagnostic output from the CESM simulation. The CESM-simulated SST trends in summerautumn and winter-spring as well as their difference agree reasonably well with the ensemble mean of CGCMs (figure S5), lending supports to the representativeness of CESM simulation.
The underlying dynamics for the faster sea surface warming in summer-autumn than winter-spring are region-dependent (figure 5). In the mid-and highlatitude North Atlantic, the amplified SST seasonal cycle is mainly due to the seasonally differed effects of turbulent vertical mixing (i.e. the MIX term). The latter reflects the role of seasonally varying mixed layer depth in SST changes: as the mixed layer is generally shallower in summer-autumn than winter-spring, an equal amount of heat uptake from the atmosphere would result in a larger SST increase in summerautumn than winter-spring (Alexander et al 2018). In most of the remaining ocean, including the mid-and high-latitude North Pacific, the faster SST increase in summer-autumn is primarily caused by the larger heat surplus injected by the air-sea heat flux (i.e. the SHF term), implying a driver from the changes of seasonality in the atmosphere (Santer et al 2018, Chen et al 2019.

Discussion
The timing of MHWs is an important factor for the biological and economic impacts of MHWs, as MHWs occurring in a particular season may exert more adverse influences on marine species than those in other seasons (Mills et al 2017, Jones et al 2018, Laurel and Rogers 2020, Santora et al 2020. Although the rising mean SST may cause a uniform increase of MHW statistics over different seasons, the effect of amplified SST seasonal cycle on MHW changes is likely to be season-dependent. To test this conjecture, we compare the projected changes of MHW statistics in winter-spring and summer-autumn. Unlike section 3, the MHWs are identified using a fixed baseline following Frölicher et al (2018) to include the MHW changes caused by both the rising mean SST and changing SST variability. Only CI is analyzed here because CD becomes saturated due to the permanent MHW state by the end of 21st century in this case (Frölicher et al 2018, Hayashida et al 2020. Figure 6 displays the projected increases of CI averaged over all the CGCM simulations in the midand high-latitude North Pacific and Atlantic with pronounced amplification of SST seasonal cycle. By the end of this century, the spatially averaged CI in the mid-and high-latitude North Pacific (Atlantic) is projected to increase to 450 (360) • C·day in summerautumn, significantly larger than 280 (180) • C·day in winter-spring (figures 6(a) and (b)). This suggests an evident seasonal difference in MHW statistics in the future under the high carbon emission scenario. Consistent with our conjecture, such seasonality of MHW statistics is primarily attributed to the amplified SST seasonal cycle. Once the effect of changing SST seasonal cycle on MHW changes is eliminated, the projected increases of CI in winter-spring and summer-autumn are almost equal (figures 6(c) and (d)). Therefore, the seasonality of MHW statistics under greenhouse warming is strongly affected by the amplified SST seasonal cycle but less so by the rising mean SST.

Conclusions
This study reveals that the amplified SST seasonal cycle in response to the anthropogenic climate changes will cause a substantial increase of MHW statistics over many parts of the global ocean, including a large fraction of LMEs. Although the impact of amplified SST seasonal cycle on MHW statistics is quantitively overwhelmed by that of the mean SST increase (Frölicher et al 2018, Oliver 2019, Hayashida et al 2020, it does not mean that the former effect is negligible as they may have different biological implications. Adaptation of marine species to the changing climate takes time. It would be easier for marine species to find suitable habitats in the case of a slowly rising mean SST at a centennial scale than a change of SST with comparable magnitude at a seasonal scale ( figure S5). Furthermore, the amplified SST seasonal cycle causes a pronounced seasonal difference in the MHW statistics in the future, whereas the seasonality of MHWs is nearly unaffected by the rising mean SST. Finally, although MHWs in this study are defined based on a seasonally varying percentile threshold, the amplified SST seasonal cycle would also probably exert strong influences on the MHW changes if a fixed threshold (e.g. an absolute temperature value) were used for definition. In particular, the mean SST increase and amplified SST seasonal cycle will work in concert to increase the extreme high SST in summer-autumn, escalating MHWs defined with a fixed threshold under the greenhouse warming. Currently, a comprehensive understanding of the biological impacts of MHW changes caused by the amplified SST seasonal cycle is still lacking but will be necessary to evaluate the influence of greenhouse warming on marine ecosystems.
All data that support the findings of this study are included within the article (and any supplementary files).