Arctic marine heatwaves forced by greenhouse gases and triggered by abrupt sea-ice melt

Since 2007

T he severe ecological and socioeconomic consequences of marine heatwaves [1][2][3] -events characterized by prolonged anomalously high sea surface temperature (SST) at a particular location 4 -have inspired numerous studies focused on comprehending their drivers and frequency [5][6][7] .However, MHWs occurring in the Arctic Ocean have received considerably less attention in comparison to those in other regions 8,9 .Arctic MHWs could have significant impacts on Arctic ecosystems.They can disrupt marine food chains, harm fish stocks, damage sensitive cold-water species, and cause loss of biodiversity 2,3 .These ecological changes can have cascading effects on indigenous communities and fisheries 1 .This raises an important question: to what extent can human emissions be held responsible for potentially high-impact MHW events in the Arctic?
The Arctic has warmed almost four times faster than the rest of the globe during the last few decades 10 , with a very pronounced seasonality as winter warming far exceeds summer warming 11,12 in conjunction with sea ice retreat 13,14 .Superimposed on the mean warming, the Arctic is experiencing increasingly many extreme temperature events 15 , including summer heat waves in the terrestrial Arctic 16 , and winter warming events over the Arctic Ocean 17 .Simultaneously, the extent, timing, and characteristics of Arctic sea ice cover are experiencing significant transformations [18][19][20] .In addition, since 2007 there has been a pronounced regime shift from thicker and deformed to thinner and more uniform ice cover 21 .
Arctic sea-ice variability is primarily driven by atmospheric temperature fluctuations 14 .The other drivers, such as surface albedo, clouds and water vapor, surface winds, and pole-ward atmospheric and oceanic energy transport together explain only 25% of sea-ice variability 14 .Increases in greenhouse gas (GHG) concentrations are driving the Arctic land surface temperature trends 22,23 .Nonetheless, the impact of GHG-induced warming, which accounts for around 60% of the warming, has been partially offset by the collective effects of other anthropogenic forcing, mainly anthropogenic aerosols 24,25 .After reaching the absolute minimum in September 2012, the seasonal variations in the sea-ice extent have settled at a new level, which is 22% lower than the average climatic norm of 1991-2020 19 .Through an extreme event attribution technique 26 , it was shown that an event of the magnitude of the 2012 record minimum in Arctic sea ice or more extreme is generally entirely attributable to the combination of anthropogenic and natural forcings, and such forcings are likely necessary for the occurrence of the event 26 .To our knowledge, an event attribution analysis of the MHW events in the Arctic has not yet been done.
In this study, our goal is to provide a quantitative assessment of whether GHG forcing, which is the predominant element of human-induced emissions, was a necessary factor for the emergence of MHWs in the Arctic.We also aim to determine whether GHG forcing alone is enough to predict the likelihood of these events to recur in the future.With these purposes, we use an extreme event-attribution technique based on causal counterfactual theory 27,28 to provide a quantitative assessment of whether GHG forcings are necessary or sufficient causes for Arctic MHWs.While this study's primary objective is to conduct an extreme event attribution analysis on Arctic MHWs, a preliminary analysis for detecting and attributing long-term changes in the SST time series is also carried out.Thus, we here perform two different statistical attribution methodologies to identify the human fingerprint in Arctic marginal seas SST changes both on multi-decadal timescale (changes of mean SST) and on extreme SST events on daily timescale (MHWs).
We make use of various data sources, including multiple satellite observations and reanalyses, along with harnessing three initial-condition large ensembles of coupled general circulation models with GHG-forcing only experiments: the Community Earth System Model with 20 members (CESM1-LE 29 ), the Max Planck Institute for Meteorology with 30 members (MPI-ESM-LR) 30 , and the Canadian Earth System Model with 50 members (CanESM5 31 ).Our findings provide clear evidence that the abrupt retreat of sea ice in the shallow marginal seas of the Arctic Ocean during the mid-summer maximum of downward radiative flux has resulted in unprecedented MHWs.Under the course of unfolding human-induced emissions, along with the expansion of first-year ice extent, moderate marine heatwaves in the Arctic will very likely persistently reoccur.

Results
Arctic MHWs characteristics.Satellite observations from the NOAA OISSTv2 32 dataset unveil that the summer of 2007 marked the beginning of a shift towards a new era of marine heatwaves over the Arctic.We detect eleven MHW events from 2007 to 2021 (Fig. 1) with an average duration of 37 days, peak temperature anomaly of 3 °C, intensity of 2.2 °C, and an average cumulative heat intensity of 95 °C days.The detected Arctic MHWs have each been accompanied by a record decline in Arctic sea ice, in particular in the years 2007, 2012, and 2020 (Fig. 1c).The 2007 MHW with 91 days duration and 3.4 °C max.intensity (188 °C days cumulative intensity; Fig. 1a) is spatially compound with an extreme minimum Arctic sea ice extent (SIE), which reached a minimum of 4.28 × 10 6 km 2 in September 2007 33 .The next minimum was reached in 2012, with the Arctic experiencing a September average SIE of only 3.62 × 10 6 km 2 34 , which is 22% below the 1991-2020 mean, concurrent with an MHW with 30 days duration and 2.1 °C intensity.In 2020 the average Arctic SIE was 3.9 × 10 6 km 2 34 , which is 19% below the 1991-2020 mean for September (Fig. 1c).This extreme minimum SIE event is a compound event with the 2020 MHW that lasts 103 days with 4 °C max.intensity, and a large 300 °C day cumulative heat intensity.During the 2020 MHW, SST over the Kara and the Laptev seas reached 6 °C above climatology.The 2020 MHW is the most severe event both in terms of intensity and duration that has ever been detected since the year 1982 (for the evolution of the 2020 MHW over the Kara see Supplementary Fig. 1c).
As illustrated in Fig. 1a, MHW events mainly happened over the Arctic marginal seas, including the Kara, Laptev, East Siberian, Chukchi, and part of the Beaufort seas.Thus, the region of interest in this study is the Arctic marginal seas, here defined as 45°E-140°W and 68°N-80°N, excluding the largely icefree Barents, and Norwegian seas.The majority of the Arctic marginal seas with a shallow mixed-layer depth (10-11 m in July-August; Supplementary Fig. 1a) are predominantly covered by first-year ice (Supplementary Fig. 1b 8 ).The recent prevalence of first-year ice phenomena 21,35,36 creates extensive areas where MHW events can occur and develop.Extreme MHW events in the Arctic are superimposed on a long-term SST increase with 1.2 °C decade −1 warming over 1996-2021 in JAS (July-September; Fig. 1b), pronounced over the Kara and the Laptev seas (Fig. 1d).In addition to the long-term warming, over the last two decades the eastern Arctic marginal seas have experienced more +25 days summer −1 with extreme SST (>95thpercentile of climatology; 1983-2012) in comparison with previous decades (Fig. 1e).
Results reveal that 82% of the SST variability over the marginal seas, where MHWs are prone to occur, can be explained by variability in the net atmospheric surface fluxes (Q net ; Fig. 1b).This is estimated with regression of the normalized values of SST (i.e., minus mean and divided by the standard deviation) as predictand against Q net as a predictor.Q net is driven by changes in the shortwave radiation (R SW ), long-wave radiation (R LW ), and turbulent heat fluxes (sensible heat (SH) + latent heat (LH)).The Q net has exhibited positive anomalies since 2006, with pronounced peaks in 2007 and 2020, indicating that since 2006 the increase in radiative heat gain by the ocean is larger than the increase in turbulent heat loss to the atmosphere (Fig. 1b).In other words the net seasonal surface heat flux in summer is being stored in the ocean.In MHW years (2007, 2012, 2016, and 2020)  the positive anomalies of Q net , which leads to positive anomalies in SST, point to the dominant effect of air-sea heat fluxes in the MHWs development.(see the section "The underlying mechanisms of Arctic marine heatwaves" for further discussions on the Arctic MHW's underlying mechanisms).In the following analysis, we utilize an event-attribution technique to identify the fraction of the likelihood of detected MHW events magnitude that is attributed to GHG forcing.
Marine heatwaves attribution.We use an extreme event attribution technique 27 to identify the fraction of the likelihood of MHWs intensity (°C), duration (days), and cumulative heat intensity (°C days) that is attributable to GHG forcing 37 .We estimate the probabilities of MHW events with specific characteristics occurring in the presence and absence of GHG forcing.These probabilities are calculated for both actual (ALL-forcing includes anthropogenic and natural external forcing) and counterfactual (fixed GHG forcing) scenarios using observations and model simulations.The estimated probabilities are used to calculate two event-attribution metrics, namely the probability of necessary causation (PN, Eqs.(2,3)), and the probability of sufficient causation (PS, Eqs.(2,3)) 38,39 (for details see the "Methods" section).A summary of the PN and the PS values for the duration, intensity, and cumulative heat intensity of the three selected MHWs is presented in Table 1.
The probability of necessary causation curve (PN), which describes the probability that GHG forcing is a necessary cause of the event, saturates to 1.0 by the time MHW intensity reaches 1.5 °C, indicating that any MHW event with intensity larger than this value could not occur without GHG forcing with ≥99% probability (Fig. 2a).The observed 2007, 2012, 2020 MHW events with 3.4, 2.1, and 4 °C intensity, respectively, are more intense  than the intensity rate at which the PN curve reaches 1.0.Thus, we conclude that for these events GHG forcing is virtually certainly a necessary cause.For instance, a PN value of approximately 0.99 ([0.98-1.0])for the intensity of the 2020 MHW means that there is a 99% chance that GHG forcing is required for this event to occur (Fig. 2a).
In terms of cumulative intensity (Fig. 2b, d), the PN curve saturates to 1.0 by the time MHW cumulative intensity reached 205 °C days, indicating that any MHW event with a cumulative intensity larger than this value could not occur without GHG forcing with ≥99% probability.For instance, the PN value of the 2007 MHW with 188 °C days is equal to 0.93 ([0.90-0.99]),suggesting that an MHW with 188 °C days cumulative intensity has <7% occurrence probability under the no-greenhouse gas effect.The PN value of the 2012 MHW with a 30-day duration is equal to 0.58 ([0.57-0.78]),indicating that an MHW with a 30day duration would still have a 42% chance to occur in a world without GHG forcing.
There is still an important question: Is GHG forcing sufficient causation for the occurrence of these MHWs?To address this question we calculate the probability of sufficient causation (PS, Eq. ( 3)), which describes the probability that the inclusion of GHG forcing is sufficient for the event's occurrence.Comparing the PS curves (Fig. 2c) with those for PN (Fig. 2a), the PN is enhanced by an event being rare in the counterfactual world (without GHG forcing), whereas PS is enhanced further by the event being frequent in the real world.The PS curve saturates to zero for intensities larger than 2.3 °C, as such events are rare even with ALL forcing scenarios.However, the PS values equal to 0.66 for events with 1 °C intensity indicate that it is likely (with ≥66% probability) that GHG forcing is a sufficient cause for these events' occurrence.In other words, there is more than a 66% chance that the events will recur when GHG forcing is present.
In summary, while for extreme events in the current climate, such as those in 2007, 2012, and 2020, the presence of GHG forcing is necessary (PN = 1.0), it is not sufficient (PS is small).This means that GHG forcing must be present (with ≥99% probability) for such events to happen, but the inclusion of GHG forcing alone is not enough to guarantee the event's occurrence.However, for moderate events, with an intensity in the range of 0.5-1 °C, GHG forcing emerges as a sufficient cause (with 66-99% probability).This implies that if GHG forcing continues to rise, events with moderate intensity will consistently recur.
The underlying mechanisms of Arctic marine heatwaves.We explore how the amount of solar energy absorbed in areas of open water has varied spatially and temporally over the Arctic MHWs era (F rw = F r (1−α)(1−SIC), See Methods for details).Results show that there is a persistent increase in the amount of solar energy deposited in the upper Arctic Ocean and surrounding seas at a maximum rate of 70% decade −1 over 1996-2021 (Fig. 3e), despite considerable inter-annual variability (Fig. 3a).Specifically, during the 2020 (2007) MHW the percent anomalies of annual cumulative solar energy absorbed was 120% (110%) higher than the average rate observed during 1983-2012 (Fig. 3a).As expected, the largest increases in solar heat input occurred over the areas with largest sea ice loss (25% decade −1 ; not-shown), and longer open water season (40 days decade −1 ; Fig. 3f).The response of SST to the absorbed solar energy is strongest when the mixed layer is shallow and thus the ocean can rapidly respond to the thermal forcing 40 .The majority of shallow marginal seas (with a summer-time mixed layer in the order of 10-11 m (Supplementary Fig. 1a)) are predominantly covered by first-year ice (Supplementary Fig. 1b).The stratification within the firstyear ice region hinders the thorough mixing and downward dispersion of energy from solar radiation 8,41 , leading to unusually high SST in this area.This extensive coverage offers ample opportunity for the occurrence of MHW events.
Our findings show that MHWs are primarily triggered by an abrupt sea ice retreat, which coincides with the midsummer (July) maximum downward radiative fluxes.The rate of sea ice melt in June/July has increased from 18 × 10 2 km 2 day −1 in 1996 to 25 × 10 3 km 2 day −1 in 2021 at a speed of 38% in 25 years (Fig. 3b).In regions with first-year ice, achieving an early sea ice retreat is more feasible.An abrupt early summer ice melt allows more ocean surface warming, particularly when it coincides with the midsummer (July) maximum surface net flux [41][42][43] .The observed net surface radiative fluxes in June/July strongly covaries with MHWs intensities (Fig. 3d).The anomalously warm SST that leads to MHWs onset are related to abnormally high radiative fluxes into the ocean, while non-MHW years such as 2010 and 2014 experienced a negative anomaly in downward radiative fluxes (Fig. 3d).However, why was SST higher in 2020 compared to 2012, which leads to a more intense MHW in 2020, despite 2012 experiencing more total ice loss by the end of summer (Fig. 1c)?Indeed, the extent of upper Arctic Ocean warming during each summer is significantly influenced by the interaction between two key factors: the timing of sea ice retreat and the atmospheric heat input (downward radiative fluxes) 41 .The 2020 MHW was stronger than the 2012 MHW owing to earlier sea ice loss in 2020 over a wide area near the peak of atmospheric heat flux (July) (Fig. 1f).The anomalously warm SST in summer 2007 can likely be attributable to northward warm Pacific ocean currents which exert a more significant influence on the SST compared to the timing of sea ice retreat.
In autumn (with a peak in October) the extra heat stored in the ocean is released back into the atmosphere through surface upward LW radiation (LW_up) and turbulent fluxes (LH + SH) 42 .According to the observed record, the LW_up and SH + LH fluxes have increased substantially over the marginal seas with a trend in the order of +5.6 W/m 2 /decade from 1996 to 2021.Specifically, during the 2020 MHW the outgoing surface fluxes were +15 W/m 2 larger than the climatology (1983-2012), which indicates an anomalously warm SST and extra heating to the air (Fig. 3c).This result is further supported by the close spatial proximity of the most substantial SST increase (Fig. 1d) with turbulent flux, and longwave radiation increases (Fig. 3g).

Is the effect of GHG-forcing detectable in Arctic SST warming?
The robustness of single MHW attribution to GHG forcing is corroborated with evidence that the basal state has also been altered by GHG forcing.Extreme MHW events in the Arctic are superimposed on a long-term SST increase with 1.2 °C decade −1 warming over 1996-2021 in August (Fig. 4a).In this section we focus on the question of whether the signal of GHG-forcing is detectable in the recently observed long-term SST increase.Here we focus on August mean SSTs since it provides the most appropriate representation of Arctic Ocean summer SSTs.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 response to GHG forcing by using uni-variate total least square (TLS) regression analysis (Eq.(1) in the "Methods" section).
The response of SST to GHG response is derived from two ensembles of model simulations: the MPI-ESM-LR 30-member 30 and the CanESM5 31 50-member ensembles forced with historical well-mixed greenhouse gases only.The trend patterns of SST in response to increasing GHG forcing from 1995 to 2020 are presented in Fig. 4b, c.The estimate of "time-evolving" internal variability is derived from the ALL forcing simulations (see the "Methods" section).Since the long-term variations of the coupled atmosphere, sea ice, and ocean in the Arctic climate system are subject to the influence of natural variability 44 such as the Pacific Decadal Oscillation 45 , and Arctic Oscillation 46 , the large size of the ensembles is a crucial requirement to robustly sample historical (transient) internal variability of SST over the Arctic.We follow the same approach used in previous studies by Barkhordarian et al. 37 and estimate the amplitude of the response of SST to external forcing from the observations via the estimation of scaling factors 47 .The resulting scaling factors, and their 95% confidence intervals conducted separately for MPI-ESM-LR (30-member), and CanESM5 (50-member) are shown in Fig. 4d, e, respectively.
The SST time series attribution results clearly illustrate the emergence of a detectable GHG signal in the SST trends ending in 2011 (2009) and later on with < 5% risk of error (Fig. 4).We reach this conclusion, as the 95% uncertainty range of scaling Detection of GHG 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).
factors, derived from fits of the regression model (Eq.( 1)) to the estimate of internal variability onto the GHG signal, does not include the zero line but is consistent with unity (a i ≠ 0 ⋂ a i = 1) for 25-year trends ending in 2011 ( 2009) and later on.
The models show striking similarities in their anthropogenically forced emergence from internal variability in SST.The clear detection of GHG signal in the observed SST trends is in line with the studies that have established a linear correlation between the decrease in summer Arctic sea-ice coverage and the rise in globalmean temperature [48][49][50] , atmospheric CO 2 concentration 51 , and cumulative anthropogenic CO 2 emissions 52 .Thus, we reinforce the event attribution findings of MHWs by showcasing that GHG forcing has significantly altered the underlying state in which extreme events take place.

Conclusions
The summer of 2007 marked the beginning of a shift towards a new era of marine heatwaves over the shallow marginal seas of the Arctic Ocean.Severe marine heatwave events are predominantly concentrated in the first-year-ice area along the edges, where stratification in summer inhibits the downward dispersion of energy from solar radiation 8 resulting in unusual fluctuations in SST.Here we show that 82% of the sea surface temperature variability over the shallow Arctic marginal seas can be attributed to changes in the net atmospheric heat fluxes.Marine heatwaves in the Arctic are primarily triggered by an early and abrupt sea-ice retreat, which coincides with the midsummer (July) maximum downward radiative fluxes.The 2020 marine heatwave was stronger than the 2012 event, even though 2012 had experienced more total sea-ice loss, owing to an abrupt sea-ice retreat in 2020 across a broad area near the peak of downward radiative fluxes in July.
Event attribution analysis reveals that any marine heatwave event with an intensity larger than 1.5 °C has <1% occurrence probability under no-greenhouse gas effect.Thus, the occurrence of the recently observed extreme marine heatwave, such as that occurred in 2007 (with 3.4 °C intensity) and in 2020 (with 4 °C intensity), would have been exceptionally unlikely in the absence of GHG forcing.Thus, GHG forcing is necessary for the occurrence of these events, though not sufficient.This means that GHG forcing must be present for such events to happen, but the inclusion of GHG forcing alone is not enough to guarantee the event's occurrence.However, for moderate events, with an intensity in the range of 0.5-1 °C, GHG forcing emerges as a sufficient cause (with 66-99% probability).This suggests that if GHG forcing continues to increase, along with the expansion of first-year ice extent, events with moderate intensity will very likely consistently reoccur.
The long-term SST detection results clearly show the emergence of a GHG signal in the observed 3 °C ocean warming (1996-2021) over the shallow Arctic marginal seas where marine heatwaves are prone to occur.The fact that marine heatwaves are superimposed on a GHG-induced systematic SST warming confirms the marine heatwave attribution results and implies that the Arctic region will face increasingly frequent and intense extreme SST events, which will exacerbate climate change impacts in the Arctic, and cause Arctic sea ice extent to shrink even faster in the near future, especially given the continued increase in GHG emissions.

Methods
Defining marine heatwaves.We identify marine heatwaves (MHWs) from daily NOAA OISST time series available from January 1982 to February 2022 and follow the standardized and widely used [53][54][55][56] MHWs definition developed in Hobday et al. 4 .
Arctic MHWs are detected when the following three criteria are satisfied: (a) SSTs exceed a seasonally varying threshold, defined as the 95th-%tile of SST variations based on a 30-year climatological period (1983-2012), (b) the extreme SSTs are sustained for at least five consecutive days with gaps of less than 3 days, and (c) the SSTs are warmer than long-term mean summer temperature 9 .Note that the region with SST standard deviation <0.25 °C has been masked out before the analysis since the MHW activities are less reliable due to the high concentration of sea ice and low SST standard deviation 9 .At each location and for each MHW, we calculated the event duration (time between start and end dates), intensity (SST anomaly above the threshold average over the event duration), and cumulative intensity (the integral of the SST anomalies over time for the duration of the event).The MHW definition used in this study is available as software modules in R (heatwaveR 57 ).
Open-water period calculation.As in several past studies [58][59][60][61] , the open-water period is defined here as the time period between the last day of observed sea ice concentration (SIC; NOAA/ NSIDC 62 ) above the 15% threshold prior to the day of annual minimum, and the first day with SIC above the 15% after the annual minimum.The annual minimum day is defined as the median of all days August-October that equals the minimum SIC for the year.A 5-day moving average is applied to the daily SIC time series at each grid cell prior to detection of the open-water period to reduce the impact of short-term SIC fluctuations 58,60 .
Cumulative heating of the upper ocean.As defined in the study by Perovich et al. 43 , the flux of solar heat input directly to the ocean (F rw ) depends on the incident solar irradiance (F r ), the ice concentration (C), and the albedo of the ocean (a) and can be expressed as F rw = F r (1−α)(1−C).The ocean albedo is set to 0.07 43,63 .Pegau and Paulson (2001) determined that while the albedo of open water in Arctic pack ice had modest variations due to solar zenith angle and cloud conditions, a value of 0.07 was typical and representative.F rw is evaluated at every grid cell each day from 1982 to 2021.Annual cumulative amounts of solar heat are estimated by summing the daily values for each year.Mean values of Q w were subtracted from each annual value to produce annual anomalies.
Long-term trend detection and attribution.We follow the same approach used in previous studies by Barkhordarian et al. [64][65][66] and estimate the amplitude of the response of SST to external forcing from the observations via the estimation of scaling factors 47,67 , which is a linear regression model as follows: where y obs represent the observations and each x i is the modeled response to one of the m forcings that are 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.Uncertainties in a i are estimated by accounting for the effect of internal climate variability on y obs , using timeevolving (historical) internal variability, which is estimated at each time step by removing the ensemble mean from each of the ensemble members of ALL-forcing simulations.The response of SST to GHG response is derived from two ensembles of model simulations: the MPI-ESM-LR 30-member ensemble 30 and the CanESM5 50-member ensemble 31 forced with historical wellmixed greenhouse gases only.MPI-ESM-LR and CanESM5 were chosen as they provide large ensembles for GHG-forcing only simulations with SST output, and represent opposite sides of the equilibrium climate sensitivity with 2.8 K 68 for MPI-ESM-LR and 5.6 K for CanESM5 31,69 .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.If the null hypothesis H DE is rejected, it indicates that the observed change y obs cannot be explained by internal variability u obs alone.Once detection has been established, attribution is assessed by testing the null hypothesis H AT : a = 1.When there is insufficient evidence to reject H AT , the attribution of changes to the respective forcing is claimed 70,71 .
Time-evolving (historical) internal variability.To have an estimate of the natural internal variability of SST, of which the Pacific Decadal Oscillation 72,73 , and Arctic Oscillation 46 are part, we utilize the 30 members of MPI-ESM-LR 30 , and the 50 members of the CanESM5 31 model's response to ALL forcing.The "evolving internal variability" is estimated at each time step by removing the ensemble mean from each 30 (50) ensemble members.
Extreme event attribution.In the extreme event attribution approach 28,74,75 based on the Causal counterfactual theory 38,76 , the variable Y represents an observed extreme event that exceeds a threshold u for a relevant climate index Z.We use this method to assess the degree to which an external climate forcing f, such as greenhouse gases (GHG) forcing, has altered the likelihood of the occurrence of the event Y.Following Hannart et al. 38,76 , we present event attribution in terms of necessary and sufficient causation.The probability of the threshold (duration of the observed MHW) being exceeded without GHG forcing is denoted as P duration fixGHG , and the probability of exceeding the threshold with GHG forcing is denoted as P duration ALL .Similarly, the probability of the threshold (intensity of the observed MHW) being exceeded without GHG forcing is denoted as P intensity fixGHG , and the probability of exceeding the threshold with GHG forcing is denoted as P intensity ALL .The estimated probabilities are used to calculate two event-attribution metrics 38,39 .
These event attribution metrics are as follows: • Probability of necessary causation (PN): This metric represents the probability that GHG forcing is a necessary cause of the event, meaning that the event would not have occurred in the absence of GHG forcing.
• Probability of sufficient causation (PS): This metric represents the probability that GHG forcing is a sufficient cause of the event, meaning that the event always occurs when GHG forcing is present.
To obtain reliable estimates of the probabilities, we use daily SST output from the Community Earth System Model Large Ensemble (CESM1-LE) 29 with a 20-member ensemble with ALL forcing and a 20-member ensemble with excluded time evolution of GHG forcing (LE-fixGHG).The differences between the two 20-member ensembles are due to internal variability and the 20 simulations can be considered as 20 plausible realizations of the real world 77 .The CESM1-LE simulations are running from 1920 to 2100.From 2006 to 2100 the Representative Concentration Pathway 8.5 forcing (RCP8.5 78 ) is used.To the best of our knowledge, the CESM1-LE is the only comprehensive model available with complementary historical single-forcing large ensembles with SST output on a daily time-scale.This version of CESM has been widely used for Arctic sea-ice studies 26,79,80 , and generally performs well.In addition, Arctic sea ice extent 81 and sea ice thickness 82 in the CESM1-LE have been shown to be realistic when compared to satellite observations post-1978.Given that the magnitude of daily SST variability in CESM1-LE compares well with OISSTv2 observations (based on detrended data during 1983-2021), this model ensemble is appropriate for MHWs attribution analysis.
We calculate the PN and, PS separately for MHW intensity as and similarly for MHWs cumulative intensity: It should be noted that the equations presented here only apply as PN or PS if the resulting values are >0; if negative, the PN or PS is assigned a probability of 0. To calculate the uncertainty on the PN (PS) a resampling 27 method is used.
Observations.The daily sea surface temperature (SST) data is from the National Oceanic and Atmospheric Administration (NOAA) daily optimum interpolation sea surface temperature (DOISST 32,83 ) Version 2.1 bias-corrected and improved 84 product with 0.25°resolution available from September 1981.In icecovered regions, the in situ and satellite-advanced very highresolution radiometer (AVHRR) SSTs were blended with proxy SSTs from ice concentrations 85 .Arctic sea surface temperatures (SSTs) are estimated mostly from satellite sea ice concentration (SIC) estimates 85 .We use the NOAA/NSIDC Climate Data Record of Passive Microwave Sea Ice Concentration in 25 km gridded resolution from the NASA Team sea ice algorithm 62,86,87 , distributed by the National Snow and Ice Data Center (NSIDC).The surface heat fluxes are derived from ERA5 88 reanalysis.
ERA5 simulates observed atmospheric profiles more accurately than other reanalysis data sets, such as ERA-Interim, JRA-55, CFSv2, and MERRA-2 89,90 .), and arguably represent one of the best data sets available for the Arctic region.Estimates of ocean heat budgets in the ERA5 data set are good, and the improved measurements of air temperatures by radiosonde and other sounding techniques have proved that the data set has significant improvements on its former predecessors 88,90,92 .
Climate models.Three large ensembles of ocean-atmosphere coupled models with single forcing experiments are used; the Community Earth System Model with 20 members (CESM1-LE 29 ), the Max Planck Institute-Earth System Model with 30 members (MPI-ESM-LR 30 ), and the Canadian Earth System Model with 50 members (CanESM5 31 ).Climate model datasets consisted of historical simulations (1950-2014) and future projections forced with the emission scenario SSP2-4.5.The summary of single-forcing experiments used in this study is as follows: • ALL signal: Ensemble of 20 simulations from CESM1-LE, ensemble of 30 simulations from MPI-ESM-LR, and ensemble of 50 simulations from CanESM5 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 simulations from CESM1-LE, 30 simulations from MPI-ESM-LR, and 50 simulations from CanESM5 forced with historical changes in wellmixed greenhouse gases only.
To determine the effect of the forcing factor that was held fixed in the "fix" ensembles, we subtract the ensemble mean of the "fix" ensemble from the ensemble mean of ALL and term these residuals GHG.Thus, the combined effects of internal variability and forced response in any individual member of the GHG ensembles can be calculated as 93

GHG
where the subscript i refers to an individual ensemble member, and the subscript em refers to the ensemble mean.

Fig. 1
Fig. 1 Arctic marine heatwave characteristics.a Cumulative heat intensity map of the major marine heatwaves (MHWs) in 2007, 2012, 2019, and 2020.b MHWs maximum intensity (gray bars; second Y axes), actual SST in July-September (OISSTv2; JAS; red curve), and reconstructed SST in JAS with net sea surface heat fluxes (Q net ; blue curve).Q net is driven by changes in shortwave radiation, long-wave radiation, and turbulent heat fluxes.The percent variance of SST variability explained by Q net is 82%.c Relative Arctic sea ice extent anomalies (blue bars) as a percentage of the 1991-2020 mean, and the detected 11 MHWs maximum intensity over 1982-2021 (red bars; second Y axes).d Linear trend in SST based on OISSTv2 over 1996-2021 in JAS (July-September) in °C/decade, e Changes in the number of days with SST > 95th-%tile of 1983-2012 climatology, over 2001-2021 minus 1982-2000.f Areas where the sea-ice melt onset in 2012 and 2020 coincides with maximum downward radiative fluxes (June/July).

Fig. 2
Fig. 2 Marine heatwaves attribution.The probability of necessary causation (PN) of greenhouse gases forcing for marine heatwave event's a intensity, b cumulative intensity.The probability of sufficient causation (PS) of greenhouse gases forcing for marine heatwave event's c intensity, d cumulative intensity.Uncertainties are estimated by bootstrap-re-sampling with replacement.Likelihood scales of ≥66% (likely), ≥90% (very likely), and ≥99% (virtually certain) are represented with horizontal lines.

Fig. 4
Fig. 4 Detection of GHG signal in observed changes of SST. a Observed trend in the month of August SST over 1995-2020, according to OISSTv2.Response of SST to GHG signal over 1995-2020 in August, derived from ensemble means of b 50 CanESM5 GHG-forcing only realizations, c 30 MPI-ESM-LR GHG-forcing only realizations.1-year moving scaling factor of observed SST trends (derived from OISSTv2) onto the GHG signal derived from d MPI-ESM-LR (30 members), e CanESM5 (50 members) over 1982-2021.The gray shaded area displays the 95%-tile range of historical (transient) internal variability-generated uncertainty of scaling factors, derived from ALL-forcing experiments of 30 members of MPI-ESM-LR and 50 members of CanESM5.Detection of GHG 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).

Table 1
Attribution of marine heatwaves duration, intensity, and cumulative intensity to greenhouse gas forcing.We present the event attribution results as the probability of necessary (PN) and sufficient (PS) causation for three MHWs detected over the Arctic in 2007, 2012, and 2020.