Greenhouse Gas Forcing a Necessary Causation for Marine Heatwaves Over the Northeast Paci�c Warming Pool

8 Over the last decade, the northeast Paciﬁc (NP) experienced strong marine heatwaves (MHWs) that produced devastating marine ecological impacts and received major societal concerns. Here, we assess the link between the well-mixed greenhouse gas (GHG) forcing and the occurrence probabilities of the duration and intensity of the NP MHWs. To begin with, we apply attribution technique on the SST time series, and detect a region of systematically and externally-forced SST increase – the long-term warming pool – co-located with the past notably Blob-like SST anomalies. The anthropogenic signal has recently emerged from the natural variability of SST over the warming pool, which we attribute primarily to increased GHG concentrations, with anthropogenic aerosols playing a secondary role. With extreme event attribution technique, we further show that GHG forcing is a necessary, but not a sufﬁcient, causation for the multi-year persistent MHW events in the current climate, such as that happened in 2019/2020 over the warming pool. However, the occurrence of the 2019/2020 MHW was extremely unlikely in the absence of GHG forcing. Thus, as GHG emissions continue to ﬁrmly rise, it is very likely that GHG forcings will become a sufﬁcient cause for events of the magnitude of the 2019/2020 record event.


Introduction
On the global scale, the frequency of marine heatwaves (MHWs) -periods of anomalously high sea surface temperature (SST) at a particular location -is projected to increase further in the 21 st Century [1][2][3][4][5] .Attribution results by Laufkootter et al 6 show that there is a 20-fold increase in occurrence probability of high-impact MHWs in the historical period (1982-2017) in comparison to the pre-industrial climate.However, the link between the well-mixed greenhouse gas (GHG) forcing and the occurrence probability of MHWs is still elusive.The aim of this study is to assess to what extend GHG forcing has contributed to the intensity and duration of MHWs, which is important not only scientifically but also for decision-making regarding ocean ecosystems, marine life and the economics of regional fisheries [7][8][9] as the GHGs emissions continue to rise.
Over the last decade the northeast Pacific ocean experienced rapid resurgence of the Blob-like SST anomalies that produced devastating marine ecological impacts [7][8][9][10][11] , with extensive socioeconomic consequences 12 .The 2014/15 predominant MHW [13][14][15] , the so-called first warm "blob" 16,17 , 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 16 .The prolonged nature 14 of the 2014/15 MHW is suggested 13 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 PDO 18 , and the North Pacific Gyre Oscillation (NPGO 19 ).A stronger NPGO-PDO coupling is predicted under anthropogenic forcing 13 .The second Blob (2.0) peaked in summer 2019 and resulted from a weakened North Pacic High, which reduced the strength of the surface winds, resulting in reduced evaporative cooling in the Northeast Pacic 20 .The study by Holbrook et al 11 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,11,13,14,16,20,21 , we here perform two different statistical attribution methodologies in order to identify the human fingerprint in Northeast Pacific SST changes both on multi-decadal 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) would increase confidence in the attribution of MHWs 22 , which can be viewed as departures from this altered base state 23 .
Here, for the first time, 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 couple general circulation models (CESM1-LE 24 with 35 members and MPI-GE 25 with 100 members).These large ensembles can provide better estimates of an individual model's internal variability and response to external forcing 26,27 , and facilitate the explicit consideration of stochastic uncertainty in attribution result 28 .We also use multiple single forcing experiments from the Detection and Attribution Model Intercomparision Project (DAMIP 29 ) component of Coupled Model Intercomparison Project phase 6 (CMIP6 30 ).Two statistical methodologies are used: (1) with the first method, we analyze the observed long-term spatio-temporal 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 31 .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 32 .(2) The second statistical method is extreme event attribution 23,33 , which determines how anthropogenic forcings have changed the likelihood of occurrence of a particular event.Following Hannart et al 34,35 , we present extreme event (MHWs) attribution in terms of necessary and sufficient causation.
Thus, the objectives of the present paper are threefold: (1) to identify regions over the North Pacific where systemically forced long-term spatiotemporal changes of SST are detectable, (2) to identify the dominant external climate drivers of the detected changes, and (3) to assess the influence of GHG forcing on the intensity and duration of MHW events over the region.Our results, based on the two different attribution methods, provide a complete picture of anthropogenic influence on SST extremes over the North Pacific.

The Northeast Pacific long-term "warming pool"
We identify an outstanding warming pool over the Gulf of Alaska in the high-resolution NOAA OISST 36 dataset (Figure 1 and Supplementary Figure S1).The detected warming is more prominent in SON months (September-November) when the SST warms by ∼ 2 • C over the 1996-2020 period (Figure 1a), which could have major ecological consequences.Hereafter, we refer to this region as the Pacific long-term warming pool.
We identify significant changes in daily mean SST and seasons over the warming pool.SSTs are now higher on average for every day of the year than they were in the late 20 st Century (1982-1999).Over the last two decades, summers are on average 1 • C warmer and 37 days longer.SSTs that marked the start of summer (June 1 st ) now come 11 days earlier, and SSTs that marked the end of summer (September 1 st ) now come around 27 days later.Spring and Autumn have shifted and winters are on average 0.5 • C warmer and 11 days shorter (Supplementary Figure S2).
The continuing changes in seasons will have increasingly profound implications over the region.
We further identify concurrent forced long-term changes in SST and large-scale atmospheric patterns in the warming pool region.We detect a reversal from cooling trends in winter, spring and autumn in an earlier period (1983-2007) to significant warming trends over the 1996-2020 period (Supplementary Figure S3).This detected reversal of trend signal is more pronounced in cold months (September-to-February), likely due to enhanced atmospheric stability.In Sep-to-Feb months, the observed trends in 500-hPa geopotential heights (500 Gph) show a tendency towards a ridge with a magnitude on the order of ∼+30 m/decade increase, located over the warming pool area (Figure 2a).The strengthening ridge is associated with a ∼+150 hPa/decade increase in SLP (Figure 2c).The pronounced increase in 500 Gph and SLP enhances atmospheric stability, suppresses the loss of heat from the ocean to the atmosphere, and thus lead to a lack of the usual cold advection in the upper ocean 16 .Therefore, the observed reversal of a cooling to a significant warming over the last decades, resulting in the configuration of the here-detected long-term Pacific warming pool, is not only due to warming, but also due to lack of usual cold advection over the region (less winter/autumn cooling).However, it remains unclear whether the here-detected Pacific long-term warming pool, and the associated trends in large-scale circulation patterns, 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.

Detection of systematic changes in the long-term Pacific warming pool
We demonstrate that the SST increase over the warming pool is stronger than could be just due to natural (internal) variability, and systematically forced changes of SST are detectable at the 95% significant level (Figure 1b-d).This result is robust across comparisons of observed SST trends with (1) unfocred trends derived from Pre-industrial control simulations, which provide up to 280 pseudo-realizations of how SST might have changed in the absence of external climate drivers (CMIP6 30 , Figure 1b), (2) naturally-forced trends derived from Paleo simulations over the 850-1850 period (CCSM4 37 , Figure 1c), and (3) trends obtained from time-evolving internal variability over 1996-2020, derived from a single-model 100-member ensemble (the MPI-ESM Grand Ensemble 25 , Figure 1d) (See

Method).
On the one hand, a substantial portion of SST trends can be explained by the natural variability in winter (DJF) and spring (MAM) (Supplementary Figure S1).On the other hand, in summer (JJA) and 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).
The pronounced increases in 500Gph (+30 m/decade) and SLP (+150 hPa/decdae) are also found to be externally forced at the 95% significant level.Such changes serve to systematically suppress the loss of heat from the ocean to the atmosphere, supporting the configuration of the warming trends.The Northeast Pacific (Golf of Alaska) stands out as a region with trends in 500Gph (Figure 2b) and SLP (Figure 2d) that are beyond the range of trends due to natural variability of the climate system, as obtained in 280 pseudo-realizations of unforced trends in 500 Gph and SLP derived from CMIP6 pre-industrial control simulations.
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 SON months where the most significant and pronounced reversal of cooling to warming trends in SST is detectable over the last decades (Supplementary Figure S3).
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) by using uni-variant total least square (TLS) regression analysis (Eq. 1 in Method).The forced response (ALL signal) is derived from two ensembles of model simulations: the CESM1 35-member Large Ensemble (CESM-LE 24 ) and the MPI-ESM 100-member Grand Ensemble (MPI-GE 25 ).
The externally forced changes in SST over the NP warming pool are attributable to anthropogenic forcing in the 21 st century.The 1-year moving scaling factors (a i ) over 1983-2020 time period demonstrates that the 25-year SST trends in the warming pool (the blue box in Figure 3a) emerge from natural variability starting in 1994-2018 onward (Figure 3b).We reach this conclusion, as the 95%-tile range of the internal variability-generated uncertainty of scaling factors (gray shaded area in Figure 3b, assessed from fits of the regression model 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 (1994-2018, 1995-2019 and 1996-2020; Fig. 3b).This indicates the emergence of a detectable ALL signal in SST changes (anthropogenic forcing being the dominant) in the 21 st century (with < 5% risk of error) over the warming pool.The similarity of the forced response in CESM-LE (Figure 3b) and MPI-GE (Figure 3c) 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, as the effect of anthropogenic forcing intensifies (Figure 3b).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 sensitive and responsive to anthropogenic forcing.
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 31 .In the following, we focus on 1996-2020 time period and assess the response of SST to the GHGs-only forcing and the AER-only forcing in isolation.

Uni-variant trend detection
Over the last two decades, both GHG-induced radiative warming and AER-induced local warming have contributed, individually, to the SST increase observed in the Northeast Pacific .The SST trend pattern driven by AER-forcing only shows that AER cools most of the global oceans over the 1980's and offsets GHG-induced radiative warming, with the most pronounced SST declines over the North Pacific ocean (Supplementary Figure S4).In the 21 st Century, however, the aerosols-induced warming acts as an amplifier for GHG-induced warming in the Northeast Pacific (Figure 4a and Figure S4).
Interestingly, the AER-forcing only simulation over 1996-2020 shows a horseshoe-like pattern of SST warming that is reminiscent of the PDO (Figure 4a).The prevailing view is that the variability of PDO occurred through internal climate variability 38 .In the 21 st Century warming climate, however, PDO trends are clearly distinguishable from internal climate variability 39,40 , with the amplitude reduced and the decadal cycle shifted toward a higherfrequency band 41 .In addition, regional trends in anthropogenic aerosols over North Pacific with aerosols-cloud interaction being the dominant aerosols forcing 42 , can influence the PDO through modulation of the Aleutian low [42][43][44][45] .Thus, it is plausible that the external forcings (GHG and AER forcing) have influenced both North Pacific SSTs and the PDO index itself.
Uni-variant detection analysis (See Method) reveals that irrespective of the models used, the elevated GHG concentration has a robust detectable influence on the observed SST increase over the region (with <5% risk of error).As the internal variability-generated uncertainty of scaling factors (a i in Eq. 1) does not include the zero In summary, we detect a signal from both GHG and AER in the observed long-term trend of SST and demonstrate that AER-forcing can act in favour of GHG-induced warming, amplifying this warming.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 bivariant (two-dimensional) attribution analysis is required, where the observed trend of SST is projected onto two hypothetical signals (GHG and AER) simultaneously 32,46 .

Bi-variant trend attribution
In order to consider the combined influence of GHG and AER signals, we carry out a two-dimensional attribution analysis 32,46 when observed SST data are projected onto two response patterns of GHG and AER simultaneously.
We demonstrate that forcing by elevated GHG levels and AER are attributed as key causes for the observed SST increase.The two-dimensional (bivariate) uncertainty contour for the GHG and AER is shown with an ellipse (Figure 4c).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 (Figure 4c).However, the uni-variant scaling factor for AER signal doesn't include unity, thus, we attribute the pattern of increased SST over the warming pool primarily to increased GHG concentrations, with anthropogenic aerosols playing a secondary role.
Providing evidence that the GHG forcing has a demonstrable influence on the base climate state in the region, it is reasonable to pursue an event attribution analysis for MHWs.In the following, we use "Extreme Event Attribution 23,33 " framework to quantify the influence of GHG forcing on the intensity, and duration of each individual detected MHWs.

Marine heatwaves characteristics
We detect 51 MHW events over the Northeast Pacific from January 1982 to April 2021 in NOAA OISSTv2 36 daily SST data (Figure 5).The detected MHWs over the 21 st century (2000-2020) are 2-fold more frequent, 2.4-fold longer lasting and 1.8-fold more intense in comparison with those occurred in the previous decades.This raises the question of whether their occurrence, intensity and duration have been influenced by external climate drivers.
The maximum SST anomaly (above the climatology in the peak date of the event) pattern of the three major MHWs, 2014/2015, 2016 and 2019/2020 is displayed in Figure 5a.The 2014/2015 MHW, which received major societal concerns, lasts for 350 days with 1.5 • C intensity (average SST anomaly during the event).The 2019/2020 MHW, which lasts for more than a year (480 days) with 1.6 • C intensity, is the most severe MHW both in terms of intensity and duration that has ever detected since the year 1982 and will have major ecological consequences.
The co-location of the Blob-like 17,20 anomalies (MHWs, Fig. 5a) over the Northeast Pacific with the here-detected long-term warming pool points towards a possible positive feedback between the mean SST warming and the changing in MHWs properties 2 .
Up to 60% of the events detected over the last decade (2010-2020) are either more intense (Figure 5b) and/or longer lasting (Figure 5c) 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 pseudo-realizations of MHWs characteristics in the absence of external climate drivers.In this manner, we can identify events that are significantly distinguishable from the distribution of events detected in Pre-industrial climate (Ho: Obs. ± P 95% = 0).
Here, we conclude, therefore, that external forcing has significant contribution (P < 0.05) to the intensity and duration of the detected MHWs, including the three major 2014/2015, 2016 and 2019/2020 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.

Marine Heatwaves Attribution (Extreme event attribution)
Extreme event attribution 33 is used to assess the influence of GHG forcing on the duration and intensity of the three major observed MHWs.Following Hannart et al 34,35 , we present event attribution in terms of necessary and sufficient causation.According to their definitions, requiring the presence of a particular forcing scenario for an event to occur would be necessary causation.In contrast, if the particular forcing scenario always produces the event in question, this would be sufficient causation 34,35 .In current study the event is MHWs and that particular forcing scenario is GHG forcing.A summary of the PS,PN and PR values (See Merthod) for duration and intensity of the three selected MHWs is presented in Table 1.
The probability ratio values (PR, Eq.3), 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 6 for the both 6/23 2014/2015 and 2019/2020 MHWs implies that the occurrence of such events is 10 6 times more likely under the influence of GHG forcing.The probability of necessary causality (PN, Eq.3), 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/2020 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.3), which describes the probability that inclusion of GHG forcing is sufficient for the event's occurrence, is small (< 0.01) for the more extreme (Table 1), as such events are rare even with ALL forcing scenarios.
In summary, event attribution results provide evidence that the occurrence of a multi-year persistent MHWs, such as that in 2014/2015 (with 350 days duration and 1.5 • C intensity) and in 2019/2020 (with 480 days duration and 1.6 • C intensity), are entirely attributable to the combination of anthropogenic and natural forcings (ALL forcing), and that 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 a sufficient, causation (PS < 0.01).In other words, a MHW of this magnitude requires GHG forcing to occur, but the inclusion of GHG forcing alone is not enough to guarantee the event's occurrence.

Summary and Conclusions
The aim of this study is to assess whether greenhouse gas (GHG) forcing was necessary for the North Pacific MHWs to occur and whether it is a sufficient cause for such events to continue to repeatedly occur in the 21st Century.
Before addressing this question and applying Extreme event attribution on MHWs, a detection and attribution analysis is performed for the SST time series.The detection and attribution of long-term change in SST strengthens event attribution findings 22 by demonstrating that external forcing has discernibly changed the background state against which MHW events occur.
We detect the Northeast Pacific as a hotspot region for marine heatwaves (MHWs) -the long-term warming pool -with an SST increase on the order of +2 • C over 1996-2020.Seasons also show significant changes over the warming pool, including the expansion of summer by more than one month (37 days) and shrinking of winters by 11 days.We further identify concurrent pronounced increases in 500 Gph and SLP over the warming pool region, which are larger than expected from natural (internal) variability alone, and consequently suppresses the loss of heat from the ocean to the atmosphere.Therefore, the observed reversal of a cooling to a significant warming over the last decades, resulting in the configuration of the here-detected long-term Pacific warming pool, is not only due to warming, but also due to lack of usual cold advection over the region (less winter/autumn cooling).
Consequently, greater exposure to heat over the region leads to 2-fold more frequent, 2.4-fold longer lasting, and 1.8-fold more intense MHWs in the 21 st century, in comparison with those detected in the 20 st century.The co-location of the Blob-like 17,20 anomalies (MHWs) over the Northeast Pacific with the here-detected long-term warming pool points towards a possible positive feedback between the mean SST warming and the changing in MHWs properties 2 .We show that up to 60% of the MHW events detected over the last decade (2010-2020) are more intense and longer lasting than could solely be attributed to climate variability in the absence of external climate drivers.
Further results indicate that, over the past decades, the SST in the warming pool has undergone a change that 7/23 is well beyond the estimated range of changes solely due to natural variability defined in both the undisturbed pre-industrial climate and the climate over 850-1850 perturbed with natural external forcing.We attribute the pattern of increased SST over the warming pool primarily to increased GHG concentrations, with anthropogenic aerosols playing a secondary role.In addition, regional trends in anthropogenic aerosols over North Pacific with aerosols-cloud interaction being the dominant aerosols forcing 42 , can influence the PDO through modulation of the Aleutian low [42][43][44][45] .Thus, it is plausible that the external forcings (GHG and aerosols forcing) have influenced both North Pacific SSTs and the PDO index itself.
After we demonstrate that the GHG forcing has a dominant influence on the base climate state in the region, we pursue an event attribution analysis for MHWs on a more localized region.Extreme event attribution analysis revleas that GHG forcing is a necessary but not sufficient causation for the multi-year persistent MHW events in the current climate, such as that happened in 2014/2015 and 2019/2020.That is, a MHW of this magnitude requires GHG forcing to occur, but the inclusion of GHG forcing alone is not enough to guarantee the occurrence of the event.
However, given that the occurrence of the 2019/2020 MHW was extremely unlikely in the absence of GHG forcing (under no-GHG forcings), combined with increasing trends in MHWs duration and intensity, it is very likely that future MHW events more extreme than that in 2019/2020 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/2020 record event.These results are important for communicating climate change, developing adaptation plans, and understanding the causes of extreme events.Table 1.Attribution of MHWs duration and intensity to GHG forcing.We present the event attribution results as probability of necessary (PN) and sufficient (PS) causation as well as probability ratio (PR) for three high-impact MHWs detected over the North Pacific, 2014/2015, 2016 and 2019/2020.

Observations
The area of research is the North Pacific ocean defined by the 10 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 24 , but each ensemble excludes the time evolution of one forcing agent: greenhouse gases (LE-fixGHG, include 20 members) and anthropogenic aerosols (LE-fixAER, include 20 members).

CMIP6-DAMIP
We also analyze five models that participate in the Detection and Attribution Model Intercomparison Project (DAMIP 29 ) component of the Coupled Model Intercomparison Project Phase 6 (CMIP6 30 ).From the DAMIP 29 single forcing runs we consider 2 groups of simulations.One group (GHG) includes 34 simulations conducted by 5 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 = d 2 ∑ d i=1 1 l i , where d is the number of models and l is the ensemble size 50 .The final internal variance is then just 1/n the internal variance.Thus, in the multi-model ensemble mean of 34 simulations conducted by 5 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 summary of model data used in this study is presented in Table 2.
The summary of single-forcing experiments used in this study is as follows: • ALL signal: Ensemble of 35 runs from CESM-LE, and ensemble of 100 runs from MPI-GE forced with ALL forcing, which includes anthropogenic factors such as human emissions of greenhouse gases, atmospheric 9/23 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 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 ensemble of 34 runs conducted by 5 CMIP6-DAMIP models forced with historical changes in anthropogenic aerosols.

Estimating natural (internal) variability
We used three sources for the estimation of "time-invariant" and "time-evolving" internal variability 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 51 .The 7,000-year pre-industrial control (PIC) runs, which are the concatenated PIC runs of 7 models, provide up to 280 pseudo-realizations (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").For all piControl simulations, a linear trend is subtracted, to reduce a possible tiny influence of model drift.
2: Time-evolving internal variability.We further utilize MPI-GE 25 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 "time-evolving" internal variability and response to external forcing 26,27 and facilitate the explicit consideration of stochastic uncertainty in attribution results 28 .The transient forced response (F t ) is quantified by taking the ensemble mean of 100 members at each time step, F t = ∑ e=100 e=1 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 25 .
3: Natural (internal+external) variability.We use the Paleo simulations over the 0850-1850 millennium derived from CCSM4 37 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 non-overlapping 25 year segments, which provides 39 pseudo-realizations of how SST might have changed in the absence of anthropogenic influences.
With the estimated internal variability we test the null-hypothesis that the observed trend in SST over 1996-2020 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.We note here the adoption of a risk of false rejection (P < 0.05) of the null hypothesis of "no external forcing".
threshold average over the event duration).In the analysis two events with a peak of less than 2 days are considered as a single event.The MHW definition used in this study is available as software modules in R (heatwaveR 59 ).

Extreme event attribution
Extreme event attribution 33 is used to assess the influence of GHG forcing on the duration and intensity of the observed MHWs.In extreme events attribution approach, 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 of the event occurring in the real world, with f present, is referred to as factual, while p 0 = P(Y = 1|X f = 0) is referred to as counterfactual (with forcing f being absent) 60,61 .
Following Hannart et al 34,35 , 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.In current study the event is MHW's intensity and duration (Y ) and that particular forcing scenario ( f ) is GHG forcing.
In order to obtain reliable estimates of the probabilities, we use daily SST output from CESM1-LE with 20-member ensemble with ALL forcing and 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 26 .It is important to first assess a model's ability to reproduce the events in question before performing an event attribution analysis 22,62 .To the best of our knowledge, the CESM1-LE is the only comprehensive model available with complementary historical single-forcing large ensembles in daily time-scale.In addition, the results of detection and attribution analysis of log-term SST presented in Sections 2.2 and 2.3, indicate that the CESM-LE is suitable for an analysis of the attribution of Marine Heatwaves in the region because of its strong detection of the anthropogenic signal (Scaling factor a i is very close to unity; Figure 3b).
Choosing the 21 st century, 2000-2020, the MHW characteristics (duration and intensity) from each year and each 20 realizations of CESM1 model are pooled together (400 years total) and probabilities are estimated.We are pooling the data to have a larger sample size to get a more robust estimate of the probabilities.
For each detected MHW, we estimate the probabilities that an 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 duration f ixGHG , of the threshold (duration of the observed MHW) being exceeded without GHG forcing, and the probability, P duration ALL , of exceeding the threshold with GHG forcing.
Similarly, we calculate the probability, P intensity f ixGHG , of the threshold (intensity of the observed MHW) being exceeded without GHG forcing, and the probability, P intensity ALL , of exceeding the threshold with GHG forcing.
The estimated probabilities are used to calculate three event attribution metrics 34,63 .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,      (35) and MPI-GE 25 (100) Historical simulations in daily time-scale CESM1-LE 64 (20 members; ALL and fixGHG) GHG-forcing only, AER-forcing only runs 34 runs conducted by 5 models of CMIP6 30 -DAMIP 29 20 runs conducted by CESM1-LE 64 Pre-industrial control simulations (7,000 year) CMIP6 30 Paleo simulations (0850-1850 millennium) CCSM4 37 Paleo simulations line for the GHG signal derived either from CESM1-LE (ensemble mean of 20 simulations, green bar, Fig 4b) or CMIP6-DAMIP models (ensemble mean of 34 simulations conducted by 5 DAMIP models, cyan bar).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 Golf of Alaska also contribute significantly to the observed increase in SST over the region.

Figure 1 .
Figure 1.Detection of systematically forced trends in sea surface temperature (a) Observed trend pattern of SST over 1996-2020 in SON (September-November) based on NOAA OISSTv2.(b) Regions where systematically forced changes of SST are detectable (in comparison with 280 pseudo-realizations of unforced trends derived from 7,000-year CMIP6 Pre-industrial control runs, P < 0.05).(c) Regions where anthropogenically forced changes of SST are detectable (in comparison with 39 pseudo-realizations of naturally forced trends derived from CCSM4 850-1850 Paleo simulation, P < 0.05).(d) Regions where systematically forced changes of SST are detectable (in comparison with unforced trends due to transient variability over 1996-2020 estimated from the 100 realizations of MPI-GE, P < 0.05).Units are C/decade..

Figure 2 .
Figure 2. Detection of systematically forced trends in large-scale circulation patterns.Observed pattern of changes in (a) 500 hPa geopotential height (500 Gph; units in m/decade) and (c) mean sea level pressure (SLP; units hPa/decade)) according to the ERA5 reanalysis over the period 1996-2020 in Autumn/winter (September to February).Regions where systematically forced changes of 500Gph (b) and SLP (d) are detectable with less than 5% risk of error (in comparison with 280 pseudo-realizations of unforced trends derived from CMIP6 Pre-industrial control runs).

Figure 3 .Figure 4 .
Figure 3. Detection of climate change signal over the warming pool.(a) Observed trend pattern of SST over 1996-2020.((b) 1-year moving scaling factors of 25-year trends of observed SST onto ALL signal (anthropogenic + natural) derived from ensemble mean of CESM1-LE 35 members over 1983 to 2020 in SON over the warming pool (blue box in Figure 3a).The gray shaded area displays the 95%-tile range of the internal variability-generated uncertainty of scaling factors in a stationary climate assessed from fits of the regression model to 280 control run segments for the raw (dark gray) and double (light gray) the model variance.(c) 1-year moving scaling factors of 25-year trends of observed SST onto ALL signal derived from ensemble mean of MPI-GE 100 members.The gray shaded area displays the 95%-tile range of the internal variability-generated uncertainty of scaling factors in a stationary climate assessed from MPI-GE transient variability over 1983-2020.Years on the horizontal axis indicate the end year of the moving 25-year time window.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)

Figure 5 .
Figure 5. Characteristics of marine heatwaves over the North Pacific.(a) The maximum SST anomaly patterns (above the climatology in the peak date of the event) for 2014/2015, 2016 and 2019/2020 MHWs.(b) Intensity, (c) duration of the detected 51 MHWs between January 1982 to April 2021.Each black (red) bar represents a MHW and the gray (blue) whiskers show the 95%-tile distribution of MHW characteristics in an undisturbed Pre-industrial climate.Detection is claimed in cases where the whiskers does not include the zero line (Ho: Obs. ± P 95% = 0).Events that are statistically distinguishable (P < 0.05) from the distribution of events in Pre-industrial climate are denoted by red bars with blue whiskers.The three major MHWs with high-impacts occurred in 2014/2015, 2016 and 2019/2020, are denoted with green stars.

Table 2 .
49N -60 • N, 90 • E -90 • W latitude-longitude domain.Two data sets of sea surface temperature are used: (1) NOAA OISSTv236, which is daily remotely sensed NationalLarge ensemblesTwo ensembles of ocean-atmosphere coupled model simulations are used; the Community Earth System Model 35-members Large Ensemble (CESM-LE)24, and the Max Planck Institute for Meteorology 100-members Grand Ensemble (MPI-GE)25.The CESM-LE includes 35 simulations (members) running from 1920 to 2100.From 2006 to 2100 the Representative Concentration Pathway 8.5 forcing (RCP8.549) is used.The MPI-GE is conducted by using the Max Planck Institute Earth System Model (MPI-ESM1.1)and includes 100 simulations (members)

Table 2 .
Observation and model data used in this study