Recent marine heatwaves in the North Pacific warming pool can be attributed to rising atmospheric levels of greenhouse gases

Abstract

Over the last decade, the northeast Pacific experienced marine heatwaves that caused devastating marine ecological impacts with socioeconomic implications. Here we use two different attribution methods and show that forcing by elevated greenhouse gases levels has virtually certainly caused the multi-year persistent 2019–2021 marine heatwave. There is less than 1% chance that the 2019–2021 event with ~3 years duration and 1.6 ^∘C intensity could have happened in the absence of greenhouse gases forcing. We further discover that the recent marine heatwaves are co-located with a systematically-forced outstanding warming pool, which we attribute to forcing by elevated greenhouse gases levels and the recent industrial aerosol-load decrease. The here-detected Pacific long-term warming pool is associated with a strengthening ridge of high-pressure system, which has recently emerged from the natural variability of climate system, indicating that they will provide favorable conditions over the northeast Pacific for even more severe marine heatwave events in the future.

Introduction

On the global scale, the frequency of marine heatwaves (MHWs)—periods of anomalously high sea surface temperature (SST) at a particular location that lasts for at least 5 consecutive days—is projected to increase further in the twenty-first Century^1,2,3,4,5. The probability of MHWs exceeding the pre-industrial 99th-%tile will increase under future global warming¹. Attribution results by Laufkoetter et al.⁶ further show that the occurrence probabilities of the duration and intensity of high-impact MHWs have increased 20-fold in the historical period (1982–2017) in comparison to the pre-industrial climate. MHWs are extreme events with wide-ranging impacts on marine ecosystems, including geographical shifts of species⁷, and mortality of marine mammals and sea birds⁸. Beyond their impacts on marine ecosystems, ecological responses to MHWs can have socioeconomic implications⁹, such as loss of fisheries income, food availability, erosion of essential ecosystem services (e.g., carbon capture, water quality), and mass mortality of iconic species⁹. The aim of this study is to identify the fraction of the likelihood of a MHW event’s magnitude that is attributable to GHG forcing, which is important not only scientifically, but also for decision-making regarding ocean ecosystems, the economics of regional fisheries, and marine life^8,9,10,11 as the GHGs emissions continue to rise.

Over the last decade, the northeast Pacific ocean experienced a rapid resurgence of the Blob-like SST anomalies that caused devastating marine ecological impacts^11,12, such as very low-ocean productivity¹³; dramatic mortality events in seabird species¹⁴; and major outbreaks of harmful algal blooms that produce extremely dangerous toxins¹⁵. Furthermore, they contributed to severe drought conditions across the US West Coast¹⁶. The 2014–2015 MHW^17,18,19, the so-called first warm “blob”^20,21, which appeared off the coast of Alaska, originated in winter and was the result of a resilient atmospheric ridge in the Northeast Pacific, which weakened the climatological Aleutian Low and related surface winds²⁰. The prolonged nature¹⁸ of the 201–2015 MHW is suggested ¹⁷ to be linked to the dynamics and coupling between the two dominant modes of winter (January–March) SST variability in the North Pacific, namely the Pacific Decadal Oscillation (PDO²²), and the North Pacific Gyre Oscillation (NPGO²³). A stronger NPGO-PDO coupling is predicted under anthropogenic forcing¹⁷. The second Blob (2.0) peaked in summer 2019 and resulted from a weakened North Pacific High, which reduced the strength of the surface winds, resulting in reduced evaporative cooling in the Northeast Pacific^24,25. The study by Holbrook et al.¹² identifies multiple drivers that can enhance the occurrence of MHWs over the North Pacific, including El Niño and the positive PDO phases.

Unlike previous studies which have focused on linking the SST patterns in the North Pacific to changes in the oceanic circulation and the extratropical/tropical teleconnections^{2,12,17,18,20,24,26}, we here perform two different statistical attribution methodologies in order to identify the human fingerprint in Northeast Pacific SST changes both on multidecadal timescale (changes of mean SST) and on extreme SST events on daily timescale (Marine Heatwaves). Evidence that anthropogenic forcing has altered the base state (long-term changes of mean SST) over the northeast Pacific, which is characterized by strong low-frequency SST fluctuations, would increase confidence in the attribution of MHWs²⁷, since rising mean SST is the dominant driver of increasing MHW frequency and intensity, outweighing changes due to temperature variability^1,2.

In this study, we provide a quantitative assessment of whether GHG forcing, the main component of anthropogenic forcings, was necessary for the North Pacific high-impact MHWs (the Blob-like SST anomalies) to occur, and whether it is a sufficient cause for such events to continue to repeatedly occur in the future. With these purposes, we use two high-resolution observed SST datasets, along with harnessing two initial-condition large ensembles of coupled general circulation models (CESM1-LE^28,29 with 35 members, and MPI-GE³⁰ with 100 members). These large ensembles can provide better estimates of an individual model’s internal variability and response to external forcing^31,32, and facilitate the explicit consideration of stochastic uncertainty in attribution results³³. We also use multiple single-forcing experiments from the Detection and Attribution Model Intercomparision Project (DAMIP³⁴) component of Coupled Model Intercomparison Project phase 6 (CMIP6³⁵).

Two statistical methodologies are used: (1) with the first method, we analyze the observed long-term spatiotemporal changes of SST to (a) detect the presence of a signal beyond changes solely due to natural (internal) variability, and to (b) attribute the detected changes in long-term SST to external climate drivers. The climate over North Pacific is potentially influenced by two external climate drivers: well-mixed greenhouse gases (GHGs) and anthropogenic aerosols (AER), which have opposing effects on the radiation energy balance³⁶. Our attribution analysis is based on assessing the amplitude of the response of SST to each external forcing from the observations via the estimation of scaling factors^37,38. (2) The second statistical method is extreme event attribution^27,39, which determines how anthropogenic forcings have changed the likelihood of occurrence of a particular event. Following Hannart et al.^40,41, we present extreme event attribution in terms of necessary and sufficient causation. Our results, based on the two different attribution methods, provide a complete picture of anthropogenic influence on SST extremes over the North Pacific.

In this study, we show that forcing by elevated well-mixed GHG levels, has virtually certainly caused the multiyear persistent 2019–2021 marine heatwave, in a necessary causation sense. There is less than 1% chance that the 2019–2021 event with ~3 years duration and 1.6 °C intensity could have happened in the absence of GHG forcing. We further discover an outstanding warming pool over the Northeast Pacific, co-located with the multi-year persistent MHW events. The warming pool is marked by concurrent and pronounced increase in annual mean, and variance of SSTs, decrease in cold-season low-cloud’s cooling effect, and strengthening cold-season atmospheric ridge of high-pressure system. There is less than 5% chance that natural (internal) variability is responsible for the observed cold-season strengthening high-pressure system. Our long-term SST detection and attribution results further reveal that forcing by elevated GHG levels and the recent industrial aerosol-load decrease are the key causes for the configuration of the here-detected long-term warming pool (P <  0.05).

Results

Marine heatwaves characteristics

We detect 40 MHW events over the Northeast Pacific from January 1982 to March 2022 in NOAA OISSTv2⁴² daily SST data (Fig. 1). The detected MHWs over the twenty-first century (2000–03/2022) are 4.5-fold more frequent, ninefold longer-lasting, and threefold more intense in comparison with those occurred in the previous decades (Fig. 1c). During the three major MHWs over the northeast Pacific in 2014–2015 (Fig. 2b), 2016 (Fig. 2c), and 2019–2021 (Fig. 2a) the maximum SST anomalies (above the climatology in the peak date of the event) reached 5.6 °C, 5 °C, and 5.8 °C, respectively. The 2014–2015 MHW, which received major societal concerns, lasts for 600 days with 1.4 °C intensity (average SST anomaly during the event). The 2019–2021 MHW, which lasts for almost 3 years (1000 days) with 1.6 °C intensity, is the most severe MHW both in terms of intensity and duration that has ever been detected since the year 1982. This raises the question of whether their occurrence, intensity, and duration have been influenced by external climate drivers.

Fig. 1: External forcing had a significant contribution to the observed marine heatwaves characteristics.

a Intensity, b duration of the observed 40 marine heatwave events (MHWs) detected between January 1982 to March 2022 over the northeast Pacific (160°W–100°E, 10°N–60°N). Each black (red) bar represents a MHW event and the gray (blue) whiskers show the 98%-tile distribution of MHW characteristics in an undisturbed Pre-industrial climate. Detection is claimed in cases where the whiskers do not include the zero line (Ho: Obs. ± P_98% ≠ 0). Events that are statistically distinguishable (P < 0.02) from the distribution of events in Pre-industrial climate are denoted by red bars with blue whiskers. Up to 60% of the events detected over the last decade (2010–2021) are either more intense and/or longer-lasting than could solely be attributed to climate variability in the absence of external climate drivers. The two major MHWs with high impact occurred in 2014–2015, and 2019–2021 are denoted with stars. c Changes in MHWs characteristics (duration, frequency, and intensity) over 2000–03/2022 in comparison with those detected over 1982–1999.

Fig. 2: Evolution of the MHWs, and the long-term trends of the background state.

a SST anomaly patterns (above the 1983–2012 climatology) during the evolution of the MHW from 2019 to 2021, b the MHW from 2014 to 2015, and c the MHW from 05/2016 to 11/2016. Observed pattern of trends over 1996–2021 in (d) annual mean SST (OISST; Units: ^∘C decade⁻¹), e SST annual variance (OISST; units: °C² decade⁻¹), f cold-season (November-to-February) geopotential height at 500-hPa (ERA5 reanalysis; units: m decade⁻¹).

Our results reveal that up to 60% of the MHW events detected over the last decade (2010–03/2022) are either more intense (Fig. 1a) and/or longer-lasting (Fig. 1b) than could solely be attributed to climate variability in the absence of external climate drivers. We arrive at this conclusion by comparing the detected MHWs with events occurring solely in response to the natural variability of the climate system. To obtain these estimates, we use CMIP6 pre-industrial control simulations of daily SST, which provide pseudorealizations of MHWs characteristics in the absence of external climate drivers^43,44. In this manner, we can identify events that are significantly distinguishable from the distribution of events detected in Pre-industrial climate (Ho: Obs. ± P_98% ≠ 0).

Here, we conclude, therefore, that external forcing has a significant contribution (P < 0.02) to the intensity and duration of the detected MHWs, including the three major 2014–2015, 2016, and 2019–2021 MHWs. However, it remains unclear whether the well-mixed GHG forcing is a sufficient causation for the occurrence of these MHWs and to what extend GHG forcing has changed the likelihood of the observed events. In the following, we use “Extreme Event Attribution^27,39” framework to identify the fraction of the likelihood of these events duration and intensity that is attributable to GHG forcing.

Marine heatwaves attribution (extreme event attribution)

We use extreme event-attribution technique³⁹, and assess the influence of GHG forcing on the duration and intensity of the observed MHWs. In extreme events attribution approach^27,45,46, Y is defined as an observed extreme situation based on exceedance over a threshold u of a relevant climate index Z. We evaluate the extent to which an external climate forcing f, for instance greenhouse gas (GHG) forcing, has changed the probability of occurrence of the event Y. The variable X_f is introduced to indicate whether or not the forcing f is present. The probability p1 = P(Y = 1∣X_f = 1) of the event occurring in the real world, with f present, is referred to as factual, while p₀ = P(Y = 1∣X_f = 0) is referred to as counterfactual (with forcing f being absent)^45,46. Following Hannart et al.^40,41, we present event attribution in terms of necessary and sufficient causation. According to their definitions, requiring the presence of a particular forcing scenario (f) for an event to occur would be necessary causation. In contrast, if the particular forcing scenario (f) always produces the event in question, this would be sufficient causation^40,41. In this study, the event is MHW’s intensity and duration (Y) and that particular forcing scenario (f) is GHG forcing.

For each detected MHW, we estimate the probabilities that a MHW has occurred that equals or exceeds the duration and intensity of the observed MHW in ALL forcing (actual world) and fixGHG forcing (counterfactual world) scenarios. That is, we calculate the probability, ${P}_{fixGHG}^{duration}$, of the threshold (duration of the observed MHW) being exceeded without GHG forcing, and the probability, ${P}_{ALL}^{duration}$, of exceeding the threshold with GHG forcing. Similarly, we calculate the probability, ${P}_{fixGHG}^{intensity}$, of the threshold (intensity of the observed MHW) being exceeded without GHG forcing, and the probability, ${P}_{ALL}^{intensity}$, of exceeding the threshold with GHG forcing. Choosing the 1982–2021 time period, the MHW characteristics (duration and intensity) from each year and each 20 realizations of CESM1 model are pooled together (780 years total) and probabilities are estimated.

The estimated probabilities are used to calculate three event-attribution metrics^40,47 (see “Methods”). The probability ratio (PR), which describes how many times as likely the event occurrence is with ALL forcing than ALL minus GHG forcing (fixGHG). The probability of necessary causation (PN) describes the probability that GHG forcing is a necessary cause of the particular event; that is, that GHG forcing is required for the event’s occurrence, and finally the probability of sufficient causation (PS) describes the probability that GHG forcing is sufficient for the event, such that a scenario with GHG forcing will see the occurrence of this event every time. A summary of the PS, PN, and PR values for the duration and intensity of the three selected MHWs is presented in Table 1.

Table 1 Attribution of MHWs duration and intensity to GHG forcing.

The probability ratio values (PR, Eq. (1), Eq. (2)), which quantifies the additional probability of an event’s intensity (duration) due to GHG forcings, increases for the more extreme MHWs (Table 1). A PR value of 10⁶ for both 2014–2015 and 2019–2021 MHWs implies that the occurrence of such events is 10⁶ times more likely under the influence of GHG forcing. The probability of necessary causality (PN, Eq. (1), Eq. (2)), which describes the probability that GHG forcing is a necessary cause of the detected MHWs, increases with increasing the severity of MHWs in terms of both duration and intensity. For instance, a PN value of approximately 0.99 ([0.98–1.0]) for the 2019–2021 MHW means that 99% of the probability of such an event is due to GHG forcing. In other words, there is a 99% chance that GHG forcing is required for this event to occur. The probability of sufficient causation (PS, Eq. (1), Eq. (2)), which describes the probability that the inclusion of GHG forcing is sufficient for the event’s occurrence, is small (<0.01) for the more extreme events (Table 1), as such events are rare even with ALL-forcing scenarios.

In summary, event-attribution results provide evidence that the occurrence of a multiyear persistent MHW, such as that occurred in 2014–2015 (with 600 days duration and 1.4 °C intensity) and in 2019–2021 (with 1000 days duration and 1.6 °C intensity), are entirely attributable to the combination of anthropogenic and natural forcings (ALL forcing). The occurrence of the MHWs was extremely unlikely in the absence of GHG forcing (with less than 1% occurrence probability under no-GHG effect), and GHG forcing is necessary for the occurrence of these events (PN = 1.0 [0.98–1.0]; Table 1). Those MHWs are extreme in the current climate, so the inclusion of GHG forcing is a necessary, but not just merely sufficient, causation (PS < 0.01; Table 1). In other words, a MHW of this magnitude requires GHG forcing to occur (with >99% probability), but the inclusion of GHG forcing alone is not enough to guarantee the event’s occurrence.

Discovering an outstanding long-term warming pool

We discover an outstanding warming pool over the Gulf of Alaska in the high-resolution NOAA OISST⁴² dataset, co-located with the MHWs induced SST anomalies. The region is marked by a concurrent and pronounced increase in annual mean (0.4 °C decade⁻¹; Fig. 2d), and annual variance of SSTs (0.6 °C² decade⁻¹; Fig. 2e) over 1996–2021. Hereafter, we refer to this region as the Pacific long-term warming pool.

The co-location of the major MHWs over the Northeast Pacific (Fig. 2a–c) with the here-detected long-term warming pool points towards possible positive feedback, and suggests that the prominent MHWs occur where mean warming (Fig. 2d), and higher variance overlap (Fig. 2e). We also detect significant changes in daily mean SST and seasons over the warming pool. SSTs are now (2000–2021) higher on average for every day of the year than they were in the late twentieth Century (1982–1999). Over the last two decades, summers are on average 1.5 °C warmer and 37 days longer. SSTs that marked the start of summer (June 1) now come 11 days earlier, and SSTs that marked the end of summer (September 1) now come around 27 days later. Winters are on average 0.5 °C warmer and 11 days shorter (Fig. 3f).

Fig. 3: Detection of externally forced trends in large-scale circulation patterns.

a Observed trend pattern in cold-season (November-to-February) over 1996–2021 of 500-hPa geopotential height (500 Gph; units: m decade⁻¹), and b of mean sea-level pressure (SLP; units: Pa decade⁻¹). Regions where systematically and externally forced changes are detectable at 95% confidence level in (c) 500 Gph, and (d) SLP, in comparison with 280 pseudorealizations of unforced trends in 500 Gph and SLP due to pre-industrial variability derived from CMIP6 control runs. e Observed pattern of trends over 1995–2018 in total cloud-cover fraction (EUMETSAT; units: % decade⁻¹). f Changes in daily mean SST and seasons based on NOAA OISST over the Northeast Pacific warming pool (indicated with a black box on Fig. 4c). Expansion of summers by more than one month and shrinking of winters over the warming pool.

Large-scale atmospheric circulations play an important role in the occurrence of extremes⁴⁸. In the cold season (November-to-February), the observed trends in geopotential height at 500-hPa (500 Gph) according to ERA5⁴⁹ reanalysis show a tendency towards a ridge with a magnitude on the order of ~ +30 m decade⁻¹ increase, located over the warming pool area (Fig. 3a), resembling the ridiculously resilient ridge^50,51, which we show that is getting intensified over the 1996–2021 time period. The cold-season enormous strengthening ridge is associated with a ~+250 Pa decade⁻¹ increase in mean sea-level pressure (SLP, Fig. 3b), resulting in an amplification of the background state of the MHWs. The pronounced increase in 500 Gph and SLP, which enforce heat-trapping systems, are also associated with a decreasing trend in cloud cover starting in 1995/1996. The EUMETSAT⁵² satellite data shows a 5% decade⁻¹ decreasing trend in cold-season cloud cover during 1995–2018 (Fig. 3e). Low-cloud cover reduction is the major contribution to the observed decline in total cloud fraction (not shown), resulting in decreases of winter-time low-cloud’s cooling effect.

The detected long-term pacific warming pool could have major ecological consequences. For instance, marine species shift their distributions in response to warming trends^7,11, and the size of various fishes decline due to increasing water temperature and reduced oxygen levels^53,54, causing changes to the entire ecosystem.

Detection of systematic changes in the Pacific warming pool

However, it remains unclear whether the here-detected Pacific long-term warming pool, and the associated trends in large-scale circulation patterns (SLP and 500 Gph), are a simple manifestation of internal climate variability or are externally and systematically forced. To address this question, we begin with comparing the observed trends with estimates corresponding to the natural (internal) variability of the climate system, which ENSO, PDO, NPGO etc. are part of. To obtain estimates of the natural (internal) variability we use three data sources (see “Methods”). We test the null hypothesis that the observed trends in SST, SLP, and 500Gph over 1996–2021 are within the 2.5–97.5%-tile distribution of unforced trends of those variables (as derived either from the pre-industrial climate or historical variability) or naturally forced trends (as derived from the 850–1850 millennium simulation)^38,55.

There is less than 5% chance that internal variability is responsible for the observed cold-season strengthening high-pressure system in the region. In spite of the large interannual variability over the north Pacific region, our detection analyses (displayed in Fig. 3c, d) show that the pronounced increases in 500 Gph (+30 m decade⁻¹) and SLP (+250 hPa decade⁻¹), respectively, over 1996–2021 are systematically forced and larger than could be due to natural (internal) variability alone (with <0.05 risk of error). Such changes serve to systematically enhance atmospheric stability, supporting the configuration of the warming pool.

Detection analyses on the SST trends (displayed in Fig. 3c–e) reveal that the SST increase over the warming pool is as well stronger than could be due to natural (internal) variability alone, and systematically forced changes of SST are detectable at the 95% confidence level in summer and autumn (JJASON; June–November). This result is robust across comparisons of observed SST trends with (1) unforced trends derived from CMIP6³⁵ pre-industrial control simulations, which provide up to 280 pseudorealizations of how SST might have changed in the absence of external climate drivers (Fig. 4b), (2) trends obtained from time-evolving historical variability over 1996–2021, derived from the MPI-ESM GE³⁰ 100-member ensemble (Fig. 4c), and (3) naturally forced trends derived from CCSM4⁵⁶ Paleo simulations over 850–1850 (Fig. 4d).

Fig. 4: Detection of systematically forced trends in observed sea surface temperature (SST) record.

a Observed trend pattern of SST over 1996–2021 in JJASON (June-to-November) based on OISSTv2. Regions where systematically forced changes of SST are detectable at 95% significant level (b) in comparison with 280 pseudorealizations of unforced trends of SST due to pre-industrial variability derived from CMIP6 control runs, (c) in comparison with trends due to historical variability over 1996–2021 estimated from the 100 realizations of MPI-GE and, (d) in comparison with 39 pseudorealizations of naturally (internal + external) forced trends derived from CCSM4 850–1850 Paleo simulation (P < 0.05). e One-year moving scaling factors of observed 25-year SST trends onto the ALL signal (anthropogenic + natural) derived from the ensemble mean of CESM1-LE 35 members over 1983 to 2021 in JJASON over the warming pool (black box in c). The gray shaded area displays the 95%-tile range of the internal variability-generated uncertainty of scaling factors (a) assessed from fits of the regression model to 280 control run segments for the raw (dark gray) and double (light gray) the model variance. Detection of ALL signal is claimed when the gray shaded area does not include the zero line but is consistent with unity (a_i ≠ 0 ⋂ a_i = 1, with <5% risk of error).

On the one hand, in winter (JDF), and in spring (MAM) a substantial portion of SST trends over the North Pacific can be explained by the natural variability alone, over 80% and 75% of the grids, respectively (Supplementary Fig. 1). On the other hand, annually, in summer (JJA), and in autumn (SON), the configuration of the Pacific warming pool emerges from the natural variability and is found to be externally and systematically forced (P < 0.05).

Having demonstrated that the here-detected warming pool significantly deviates from natural variability, that is, externally forced changes are detectable over the region (P < 0.05), we proceed to an attribution step by checking whether the detected changes are consistent with what climate models describe as the expected response to anthropogenic forcing. The following attribution analysis focuses on JJASON months, where the most significant and pronounced warming trends in SST are detectable over the last decades (Supplementary Fig. 1).

As a first step to determining if the observed SST trends constitute the system’s forced response to external climate drivers, we project the observed SST changes on the model simulated forced response, ALL signal (anthropogenic + natural) by using univariate total least square (TLS) regression analysis (Eq. (3) in “Methods”). The forced response (ALL signal) is derived from two ensembles of model simulations: the CESM1 35-member Large Ensemble (CESM-LE²⁸) and the MPI-ESM 100-member Grand Ensemble (MPI-GE³⁰).

The 1-year moving scaling factors (a_i) during 1983 to 2021 time period demonstrates that the 25-year SST trends in the warming pool emerge from natural variability in 2018 and later on (Fig. 4e). We reach this conclusion, as the 95%-tile range of the internal variability-generated uncertainty of scaling factors for the row (dark gray), and double (light gray) the model variance (assessed from fits of the regression model (Eq. (3)) to control run segments) does not include the zero line but is consistent with unity (a_i ≠ 0 ⋂ a_i = 1) for 25-year trends ending in 2018 and later on (1993–2018, 1994–2019, 1995–2020, and 1996–2021; Fig. 4e). This indicates the emergence of a detectable ALL signal (anthropogenic forcing being dominant) in SST changes in the twenty-first century (with <5% risk of error) over the warming pool. The similarity of the forced response in CESM-LE (Fig. 4e) and MPI-GE (Supplementary Fig. 2) suggests that the smaller ensemble size of the CESM-LE does not affect the forced signal.

It is interesting to note that, the range of internal variability-generated uncertainty in the scaling factors (the gray shaded area) decreases when more recent years are included in the analysis (Fig. 4e). This could imply that the signal of SST increase is becoming dominated by the external forcing, which makes it easier to separate the signal’s large contribution from the internal climate variability (noise). In other words, the region might become more responsive to anthropogenic forcing as a result of, for instance, decreasing mixed layer depth that could cause a stronger SST response for the same heat flux⁵⁷.

The climate over North Pacific is potentially influenced by two external climate drivers: well-mixed greenhouse gases (GHGs) and anthropogenic aerosols (AER), which have evolving contribution in driving SST over the Gulf of Alaska. In the following, we focus on 1996–2021 time period and assess the response of SST to the GHGs-only forcing and the AER-only forcing in isolation.

Evolution of forced SST trends

Here we are isolating the evolving contributions of industrial aerosols, AER_indust (Note that biomass burning aerosols are not considered), and greenhouse gases (GHGs) in driving annual SST trends during the 1920–2021 time period (Fig. 5a). We use two ensembles that are identical to the CESM-LE^28,32, but each ensemble excludes the time evolution of one forcing agent: greenhouse gases (LE-fixGHG, include 20 members) and anthropogenic aerosols (LE-fixAER_indust, include 20 members). We define the response to the withheld forcing as i.e., ALL_em − fixGHG_em²⁹. The time series of GHG-induced SST anomalies shows a steep positive trend, with an increase in magnitude over 1920 to 2025 time period (Fig. 5a). The SST trend pattern driven by AER_indust-forcing only shows that AER_indust cools North Pacific over the 1980’s (Fig. 5a), and offsets GHG-induced radiative warming (Fig. 5b), with the most pronounced SST declines over the Northeast Pacific during 1956–1980 (Fig. 5c). Noteworthy, aerosols contribute negatively to the radiative forcing response (cooling effect), and there is indeed an observed increase in aerosols in this period³⁶. However, in the most recent period, the cooling due to AER-forcing is replaced by warming (Fig. 5a, g). The reversal of the cooling trends due to AER-forcing to warming after around 1995/1996 is very likely a result of concurrent aerosols optical depth (AOD) declines over North America and Europe^58,59. Consequently, the ALL-forcing trend pattern over 1996–2021 is characterized by amplified warming, when both GHG-induced radiative warming and AER_indust-induced local warming have contributed, individually, to the SST increase (Fig. 5e).

Fig. 5: Evolving roles of the greenhouse gases forcing and industrial aerosols in driving annual SST over the North Pacific.

a Area mean annual SST anomalies over the Gulf of Alaska (black dotted box in e), derived from the ensemble mean of 20 realizations of CESM1-LE model, ALL_em minus fixGHG_em in red, and ALL_em minus fixAER_em in blue. Response of annual SST to the ALL forcing (anthropogenic + natural), greenhouse gases (GHG) forcing and industrial aerosols (AER) forcing (b–d) over 1956–1980, and (e–g) over 1996–2021. h Response of mean sea-level pressure (SLP) to AER-forcing in DJF (December–February). The Aleutian low region defined by 160–220°E, 30–65°N latitude–longitude is denoted by a red box in h.

Interestingly, the AER_indust-forcing only simulation over 1996–2021 shows a horseshoe-like pattern of SST warming that is reminiscent of the PDO (Fig. 5g). Despite the PDO being commonly recognized as an internally generated mode of the climate system, recent studies have suggested that PDO phase transition is also significantly affected by external climate drivers^60,61,62. Regional trends in anthropogenic aerosols over North Pacific, with aerosols-cloud interaction being the dominant aerosols forcing⁶³, can influence the PDO through modulation of the Aleutian low^63,64,65,66. The CESM1-LE model simulates a strengthening of the Aleutian Low in DJF (decrease in mean sea-level pressure in the north Pacific index (NPI) region defined by the 160–220°E, 30–65°N latitude–longitude domain) in response to AER-forcing over 1996–2021 (Fig. 5h). The winter-time strengthening of the Aleutian Low (−50 Pa decade⁻¹ decreasing trend in SLP), induce SST trends in the northeast Pacific that resemble the PDO (Fig. 5g). This makes it difficult to distinguish PDO-related anomalies from climate change signals. Thus, it is plausible that the anthropogenic aerosols have influenced both on North Pacific SSTs and the PDO index itself.

Univariate trend detection and bivariate trend attribution

The univariate detection analysis, which is based on estimating the amplitude of the response of SST to each external forcing from the observations via the estimation of scaling factors^37,38,67 (see “Methods”), reveals that irrespective of the models used, the elevated GHG concentration has a robust detectable influence on the observed SST increase over the region (P < 0.05, Fig. 6a). As the internal variability-generated uncertainty of scaling factors (a_i in Eq. (3)) does not include the zero line for the GHG signal derived either from CESM1-LE (ensemble mean of 20 simulations, green bar, Fig. 6a) or CMIP6-DAMIP models (ensemble mean of 34 simulations conducted by 5 DAMIP models, cyan bar, Fig. 6a). Neither the scaling factors of the AER signal derived from CESM1-LE (red bar) nor from CMIP6-DAMIP models (purple bar) include the zero line, suggesting that the AER-induced local warming over the Gulf of Alaska also contributes significantly to the observed increase in SST over the region in JJASON.

Fig. 6: Detection and attribution of sea surface temperature trends over the warming pool.

a Univariate attribution: Scaling factors (a) of observed SST changes over the warming pool (black box in b), and its 95-%tile uncertainty range, assessed from fits of the regression model to control run segments, in JJASON (June-to-November) over 1996–2021 against the ALL signal (anthropogenic + natural) derived from MPI-GE (100 runs; black bar), and CESM1-LE (35 runs; red bar). Scaling factors (a) of observed SST changes and its 95-%tile uncertainty range against the greenhouse gases signal (GHG) and the anthropogenic aerosols signal (AER) derived from CESM1-LE (20 runs; green and blue bars, respectively), and derived from DAMIP-CMIP6 (34 runs conducted by 5 models; cyan and purple bars, respectively). b Observed trend pattern of SST over 1996–2021 time period. c Bivariate attribution: The ellipse displays the 90% of the estimated joint distribution of scaling factors for the GHG & AER signals when observed data are regressed onto two signals simultaneously during 1996–2021. The one-dimensional uncertainty intervals for the univariate and bivariate analysis for two signals are shown as red and black whiskers, respectively. Bivariate attribution is claimed in cases where the ellipse excludes the origin (0,0) but the point (1,1) lies inside the ellipse.

Furthermore, in order to separate the driver’s contributions to the response, a combined influence of GHG, and AER should be considered. In this case, a bivariate (two-dimensional) attribution analysis is required, where the observed trend of SST is projected onto two hypothetical signals (GHG and AER) simultaneously^37,38. The two-dimensional (bivariate) uncertainty contour for the GHG and AER is shown with an ellipse (Fig. 6c). The ellipse containing 90% of the estimated joint distribution of scaling factors for the two signals (GHG and AER) excludes the origin (0,0), indicating that the effects of GHG and AER signals are detectable simultaneously. In addition, the scaling factors are consistent with unit amplitude since the point (1,1) lies within the ellipse (Fig. 6c).

In summary, we detect a signal from both GHG and AER in the observed long-term trend of SST in JJASON, and demonstrate that AER-forcing can act in favor of GHG-induced warming, amplifying this warming. Therefore, we conclude that forcing by elevated GHG levels and the recent industrial aerosol-load decrease are identified as key causes for the here-detected long-term warming pool (P < 0.05).

Conclusions

Over the last decade, the North Pacific experienced strong marine heatwaves (MHWs) that produced devastating marine ecological as well as socioeconomic impacts, and received major societal concerns. We use the  extreme event-attribution technique based on causal counterfactual theory^40,41,47, and provide a quantitative assessment of whether GHG forcing, the main component of anthropogenic forcings, was necessary for the North Pacific high-impact MHWs (the Blob-like SST anomalies) to occur, and whether it is a sufficient cause for such events to continue to repeatedly occur in the future.

In this study, we show that forcing by elevated well-mixed GHG levels, has virtually certainly caused the multiyear persistent 2019–2021 marine heatwave. In other words, there is less than 1% chance that the 2019–2021 event with ~3 years duration and 1.6 °C intensity could have happened in the absence of GHG forcing. Our extreme event-attribution analysis further reveals that the GHG forcing is a necessary, but not sufficient, causation for the multiyear persistent MHW events in the current climate, such as that happened in 2014–2015, and 2019–2021. That is, a MHW of this magnitude requires GHG forcing to occur (with 99% probability), but the inclusion of GHG forcing alone is not enough to guarantee the occurrence of the events. However, given that the occurrence of the 2019–2021 (2014–2015) MHW was extremely unlikely in the absence of GHG forcing (with <1% occurrence probability under no-GHG forcings), combined with increasing trends in MHWs duration and intensity, it is likely that future MHW events will be attributable to GHG forcing, and that the inclusion of GHG forcings will become a sufficient cause for events of the magnitude of the 2019–2021 (2014–2015) record event.

We further discover a systematically-forced outstanding warming pool in the high-resolution NOAA OISST⁴² dataset, co-located with the MHWs maximum SST anomalies. The region is marked by concurrent and pronounced increase in annual mean (0.4 °C decade⁻¹), and variance of SSTs (0.6 °C² decade⁻¹), decreases in winter-time low-cloud’s cooling effect, and increases in atmospheric stability with a strengthening ridge. The region is further identified with 0.5 ^∘C warmer and shorter winters with 37 days longer summers. Consequently, greater exposure to heat, and lack of usual winter-time cooling leads to 4.5-fold more frequent, ninefold longer-lasting, and threefold more intense MHWs in the twenty-first century, in comparison with those detected in the twentieth century. We further show that up to 60% of the MHW events detected over the last decade (2010–2021) are either more intense and/or longer-lasting than could solely be attributed to climate variability in the absence of external climate drivers (P < 0.05). We refer to this region as the Pacific long-term warming pool.

The univariate detection analysis indicate that anthropogenic signal has recently emerged from the natural variability of SST over the warming pool in JJASON (June-to-November). The pronounced increases in 500 Gph (+30 m decade⁻¹) and SLP (+250 Pa decade⁻¹) in the cold season are also found to be externally forced, and larger than could be as a result of natural (internal) variability alone (P < 0.05). Such changes enforce heat-trapping systems, and serve to systematically enhance atmospheric stability in cold season together with decreases in winter-time low-cloud’s cooling effect, supporting the configuration of the warming pool.

The evolution of anthropogenic aerosols plays an important role in the northeast Pacific SST trends subsequently (post-1996) when East Asian aerosol emissions decreases. The cooling to warming SST transition after around 1995/1996 is reproduced in the aerosol-forced historical simulation, when the ALL-forcing trend pattern is characterized by amplified warming, with the contribution of both GHG-induced radiative warming and AER_indust-induced local warming to the SST increase over the region. Our bivariate attribution analysis over 1996–2021, based on multiple model data sources, demonstrates that forcing by elevated GHG levels and the recent industrial aerosol-load decrease are key causes for the here-detected warming pool (P < 0.05). These results further strengthen the MHWs attribution findings, that GHG forcing has virtually certainly caused the 2019–2021 MHW, by demonstrating that GHG forcing has also discernibly changed the background state against which MHW events occur. This indicates that the region will very likely experience more intense MHW events as the GHGs emissions continue to rise.

The here-detected long-term Pacific warming pool together with the strengthening high-pressure system on the background state of the high-impact MHWs, which we show that significantly deviate from natural variability, and that will continue and intensify in the course of unfolding anthropogenic climate change, will provide favorable conditions over the northeast Pacific for even more sever marine heatwave events in the future. Such change could have profound societal and marine ecosystem impacts over the region.

Methods

Observations

The area of research is the North Pacific ocean defined by the 10°N–60°N, 90°E–90°W latitude–longitude domain. Two datasets of sea surface temperature are used: (1) NOAA OISSTv2⁴², which is daily remotely sensed National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation (OI) SST V2 high-resolution (0.25^∘) gridded SST data available over January 1982 to April 2021. (2) The Hadley Centre HadISST⁶⁸, which is produced monthly with 1° grid resolution since 1870. The mean sea-level pressure (SLP) and 500 hPa geopotential height (500Gph) are derived from ERA5⁴⁹ reanalysis. Cloud-cover data is from EUMETSAT’s Satellite Application Facility on Climate Monitoring (CM SAF) available over 1982–2018⁵². The summary of observation and model data used in this study is presented in Table 2 and Supplementary Table 1.

Table 2 Observation and model data used in this study.

Large ensembles

Two ensembles of ocean-atmosphere coupled model simulations are used; the Community Earth System Model 35-members Large Ensemble (CESM-LE)²⁸, and the Max Planck Institute for Meteorology 100-members Grand Ensemble (MPI-GE)³⁰. The CESM-LE includes 35 simulations (members) running from 1920 to 2100. From 2006 to 2100, the Representative Concentration Pathway 8.5 forcing (RCP8.5⁶⁹) is used. The MPI-GE is conducted by using the Max Planck Institute Earth System Model (MPI-ESM1.1) and includes 100 simulations (members) running from 1850 to 2100 under the RCP8.5 scenarios. For disentangling the forced response and the internal climate variability each member in the ensembles is initialized with different initial conditions. The large size of the ensemble is a crucial requirement to robustly sample internal variability. The mean of each ensemble averages out the internal variability, thus represents the forced response of the system.

In addition, for attributing the detectable signal to a specific forcing agent we make use of two ensembles that are identical to the CESM-LE²⁸, but each ensemble excludes the time evolution of one forcing agent: greenhouse gases (LE-fixGHG, include 20 members) and anthropogenic aerosols (LE-fixAER_indust, include 20 members).

CMIP6-DAMIP

We also analyze five models that participate in the Detection and Attribution Model Intercomparison Project (DAMIP³⁴) component of the Coupled Model Intercomparison Project Phase 6 (CMIP6³⁵). From the DAMIP³⁴ single-forcing runs, we consider two groups of simulations. One group (GHG) includes 34 simulations conducted by five models forced with historical well-mixed greenhouse gases only. A second group (AER) includes 34 simulations conducted by 5 models forced with anthropogenic aerosols only. The DAMIP models (ensemble size) are CNRM-CM6-1 (6), CanESM5 (10), GISS-E2-1-G (4), HadGEM3-GC3-LL (4), and IPSL-CM6A-LR (10). In these simulations, aerosol and GHG emissions (or concentrations) are allowed to vary in time whereas all other forcing variables are set to pre-industrial values. The formula to give equal weighting of the individual models is, $n=\frac{{d}^{2}}{\mathop{\sum }\nolimits_{i = 1}^{d}\frac{1}{{l}_{i}}}$, where d is the number of models and l is the ensemble size⁷⁰. The final internal variance is then just 1/n the internal variance. Thus, in the multimodel ensemble mean of 34 simulations conducted by five models (GHG-only and AER-only simulations), the internal variability is reduced by about 70%, which leads to an enhanced signal-to-noise ratio in estimated signal patterns. The name of the models from DAMIP-CMIP6, the number of model ensemble members and the pre-industrial control years used in the study is presented in Table 2 and Supplementary Table 1.

The summary of single-forcing experiments used in this study is as follows:

• ALL signal: Ensemble of 35 runs from CESM-LE, and the ensemble of 100 runs from MPI-GE forced with ALL forcing, which includes anthropogenic factors such as human emissions of greenhouse gases, atmospheric aerosols, ozone, land-use changes, and natural external factors such as stratospheric aerosols due to the large volcanic eruptions and solar forcing.

• GHG signal: Ensemble of 20 runs from CESM-LE, and the ensemble of 34 runs conducted by 5 CMIP6-DAMIP models forced with historical changes in well-mixed greenhouse gases.

• AER signal: Ensemble of 20 runs from CESM-LE, and the ensemble of 34 runs conducted by 5 CMIP6-DAMIP models forced with historical changes in anthropogenic aerosols.

Defining marine heatwaves

We identify marine heatwaves (MHWs) from daily NOAA OISST time series available from January 1982 to March 2022 and follow the standardized and widely used^1,2,12,71 MHWs definition developed in Hobday et al.¹⁹. MHWs occur when SSTs exceed a seasonally varying threshold, defined as the 95th-%tile of SST variations based on a 30-year climatological period (1983–2012), for at least five consecutive days. At each location and for each MHW, we calculated the event duration (time between start and end dates), frequency and intensity (SST anomaly above the threshold average over the event duration). In the analysis, two events with a peak of less than 3 days are considered as a single event. The MHW definition used in this study is available as software modules in R (heatwaveR⁷²).

Extreme event attribution

Extreme event attribution³⁹ is used to assess the influence of GHG forcing on the duration and intensity of the observed MHWs. Following Hannart et al.^40,41, we present event attribution in terms of necessary and sufficient causation. In order to obtain reliable estimates of the probabilities, we use daily SST output from CESM1-LE with a 20-member ensemble with ALL forcing and a 20-member ensemble with excluded time evolution of GHG forcing (LE-fixGHG). The differences among 20 ensemble members are due to internal variability and the 20 simulations can be considered as 20 plausible realizations of the real world³¹. To the best of our knowledge, the CESM1-LE is the only comprehensive model available with complementary historical single-forcing large ensembles in daily timescale. Given that the magnitude of daily SST variability in CESM1-LE compares well with observations (standard deviation of 1.91 vs 1.94 °C, respectively, based on detrended data during 1982–2021), this model ensemble can be used for MHWs attribution analysis. In addition, the results of detection and attribution analysis of long-term SST presented in “Detection of systematic changes in the Pacific warming pool” and “Univariate trend detection and bivariate trend attribution”, indicate that the CESM-LE is suitable for an analysis of the attribution of MHWs in the region because of its strong detection of the anthropogenic signal (scaling factor a_i is very close to unity; Figs. 4e and  6a).

Choosing the 1982–2021 time period, the MHW characteristics (duration and intensity) from each year and each 20 realizations of CESM1 model are pooled together (780 years total) and probabilities are estimated. For each detected MHW, we estimate the probabilities that a MHW has occurred that equals or exceeds the duration and intensity of the observed MHW in ALL forcing (actual world) and fixGHG forcing (counterfactual world) scenarios.

We calculate the probability of necessary causation (PN), sufficient causation (PS) and probability ratio (PR) separately for MHWs duration as:

$$PR=\frac{{P}_{ALL}^{duration}}{{P}_{fixGHG}^{duration}} ;\quad PN=1-\frac{{P}_{fixGHG}^{duration}}{{P}_{ALL}^{duration}};\quad PS=1-\frac{1-{P}_{ALL}^{duration}}{1-{P}_{fixGHG}^{duration}},$$

(1)

and similarly for MHWs intensity:

$$PR=\frac{{P}_{ALL}^{intensity}}{{P}_{fixGHG}^{intensity}} ;\quad PN=1-\frac{{P}_{fixGHG}^{intensity}}{{P}_{ALL}^{intensity}};\quad PS=1-\frac{1-{P}_{ALL}^{intensity}}{1-{P}_{fixGHG}^{intensity}}$$

(2)

It should be noted that the equations presented here only apply as PN or PS if the resulting values are greater than 0; if negative, the PN or PS is assigned a probability of 0.

Calculating uncertainties on PN

To calculate the uncertainty on the probability of necessary causation (PN) a resampling³⁹ method is used. The pool of data from each year in 1982–2021 and from each 20 ensemble member in CESM1-LE is re-sampled to reproduce a set of data to use for the calculation of the probabilities. With 1000 estimates of the density curves, a nonparametric 90% confidence interval for each of the metrics can be determined.

Estimating natural (internal) variability

We used three sources for the estimation of time-invariant and time-evolving (historical) internal variability of climate system, of which PDO, ESNO, NPGO, etc. are part, as well as natural (internal + external) variability.

1: Time-invariant internal variability

The long pre-industrial control simulations from global climate models (GCMs) participating in the CMIP6 project are performed under control conditions (i.e., with constant atmospheric composition, no episodic volcanic influences, and no variation in solar output). The models (number of years used from control integration) are CESM1-LE (1000), MPI-GE (2000), CNRM-CM6-1 (500), CanESM5 (1000), GISS-E2-1-G (1000), HadGEM3-GC3-LL (500) and IPSL-CM6A-LR (1000). Since climate models may underestimate variability, we double the simulated variance prior to the attribution analysis⁷³. The 7000-year pre-industrial control (PIC) runs, which are the concatenated PIC runs of 7 models, provide up to 280 pseudorealizations (the series is split into 280 nonoverlapping 25-year segments) of how the climate might have changed in the absence of external influences ("Pre-industrial variability"). The pre-industrial variability are calculated using concatenated control runs in order to capture the very low-frequency 10–20 years fluctuation of SST. For all piControl simulations, a linear trend is subtracted, to reduce a possible tiny influence of model drift.

2: Time-evolving (historical) internal variability

We further utilize MPI-GE³⁰ large initial-condition ensemble, which produces 100 realizations of single model’s response to ALL forcing. This ensemble can provide better estimates of an individual model’s historical internal variability and response to external forcing^31,32. The transient forced response (F_t) is quantified by taking the ensemble mean of 100 members at each time step, ${F}_{t}=\frac{\mathop{\sum }\nolimits_{e = 1}^{e = 100}{f}_{et}}{100}$, where f_et is a single ensemble member at time step t. The “evolving internal variability” is estimated at each time step by removing the ensemble mean (F_t) from each 100 ensemble members (f_et) and calculate the standard deviation across the ensemble³⁰.

3: Natural (internal + external) variability

We use the Paleo simulations over the 0850–1850 millennium derived from CCSM4⁵⁶ model to obtain an estimate of natural external variability of SST, associated with solar variability and stratospheric aerosols due to large volcanic eruptions. This millennium simulation is split into 39 nonoverlapping 25-year segments, which provides 39 pseudorealizations of how SST might have changed in the absence of anthropogenic influences⁵⁵.

With the estimated internal variability, derived from the three sources mentioned above, we test the null hypothesis that the observed trend over 1996–2021 is within the 2.5–97.5%-tile distribution of unforced trends (as derived either from the pre-industrial control simulations or MPI-GE 100-member historical runs) or naturally forced trends (as derived from the 850–1850 millennium simulation). If the null hypothesis is rejected, it indicates that the observed SST changes deviates significantly from internal variability, that is, the observed change in SST cannot be explained by internal variability alone, and externally (systematically) forced changes are detectable in observed SST record at 5% significant level. We note here the adoption of risk of false rejection (P <  0.05) of the null hypothesis of “no external forcing”. Since when the regional null hypothesis is valid, on an average n = 0.05 m local alternatives will falsely be rejected (m number of grid points)⁷⁴.

Long-term trend detection and attribution

We follow the same detection and attribution approach used in previous studies by Barkhordarian et al.^37,38,55,67 that is a two-step process. In the first step—detection—consists of assessing whether systematically (externally) forced changes are detectable in the observed SST change. This is achieved by testing the null hypothesis that the observed change in SST is drawn from the population of an undisturbed climate. In the second step—attribution—we examine the null hypothesis that the observed change in SST is drawn from a hypothetical population of a climate disturbed by a specific external influence.

The attribution analysis here is based on estimating the amplitude of the response of SST to each external forcing from the observations via the estimation of scaling factors^75,76, which is a linear regression model as follows:

$${y}_{obs}=\mathop{\sum }\limits_{i=1}^{m}({x}_{i}-{u}_{i}){a}_{i}+{u}_{obs},$$

(3)

where y_obs represents the observations and each x_i the modeled response to one of m forcings that is anticipated by climate models. a_i is an unknown scaling factor. The noise on y_obs, denoted by u_obs, is assumed to represent internal climate variability, while the noise on x_i, denoted by u_i, is a result of both internal variability and the finite ensemble used to estimate the model response. In order to account for the noise in response patterns Total Least Squares (TLS) methodology⁷⁷, which minimizes the perpendicular distance between the scatter points and the best-fit line is used.

Uncertainties in a_i are estimated by accounting for the effect of internal climate variability on y_obs, using samples from climate model control simulations. Therefore, the distribution of the scaling factor a_i is assessed from fits of the regression model (Eq. (3)) to 280 nonoverlapping control run segments derived from 7000-year control simulations (conducted by 7 global climate models).

Detection of a climate change signal occurs if the uncertainty range around a scaling factor a_i is shown to be significantly different from zero. This is handled by testing the null hypothesis H_DE: a = 0 (where 0 is a vector of zeros). If the null hypothesis H_DE is rejected, it indicates that the observed representation of climate change deviates significantly from internal variability, that is, the observed change y_obs cannot be explained by internal variability u_obs alone. Once detection has been established, attribution (consistency of observed changes with a combination of external forcing) is assessed by testing the null hypothesis H_AT: a = 1 (where 1 denotes a vector of unit). When there is insufficient evidence to reject H_AT, the attribution of changes to the respective forcing is claimed^37,38,67.