Introduction

The severe ecological and socioeconomic consequences of marine heatwaves1,2,3—events characterized by prolonged anomalously high sea surface temperature (SST) at a particular location4—have inspired numerous studies focused on comprehending their drivers and frequency5,6,7. However, MHWs occurring in the Arctic Ocean have received considerably less attention in comparison to those in other regions8,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 biodiversity2,3. These ecological changes can have cascading effects on indigenous communities and fisheries1. 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 decades10, with a very pronounced seasonality as winter warming far exceeds summer warming11,12 in conjunction with sea ice retreat13,14. Superimposed on the mean warming, the Arctic is experiencing increasingly many extreme temperature events15, including summer heat waves in the terrestrial Arctic16, and winter warming events over the Arctic Ocean17. Simultaneously, the extent, timing, and characteristics of Arctic sea ice cover are experiencing significant transformations18,19,20. In addition, since 2007 there has been a pronounced regime shift from thicker and deformed to thinner and more uniform ice cover21.

Arctic sea-ice variability is primarily driven by atmospheric temperature fluctuations14. 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 variability14. Increases in greenhouse gas (GHG) concentrations are driving the Arctic land surface temperature trends22,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 aerosols24,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–202019. Through an extreme event attribution technique26, 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 event26. 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 theory27,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-LE29), the Max Planck Institute for Meteorology with 30 members (MPI-ESM-LR)30, and the Canadian Earth System Model with 50 members (CanESM531). 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 OISSTv232 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 × 106 km2 in September 200733. The next minimum was reached in 2012, with the Arctic experiencing a September average SIE of only 3.62 × 106 km2 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 × 106 km2 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).

Fig. 1: Arctic marine heatwave characteristics.
figure 1

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 (Qnet; blue curve). Qnet is driven by changes in shortwave radiation, long-wave radiation, and turbulent heat fluxes. The percent variance of SST variability explained by Qnet 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).

Full size image

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 ice-free 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. 1b8). The recent prevalence of first-year ice phenomena21,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 (>95th-percentile 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 (Qnet; 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 Qnet as a predictor. Qnet is driven by changes in the shortwave radiation (RSW), long-wave radiation (RLW), and turbulent heat fluxes (sensible heat (SH) + latent heat (LH)). The Qnet 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 Qnet, 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 technique27 to identify the fraction of the likelihood of MHWs intensity (°C), duration (days), and cumulative heat intensity (°C days) that is attributable to GHG forcing37. 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.

Table 1 Attribution of marine heatwaves duration, intensity, and cumulative intensity to greenhouse gas forcing.
Full size table

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).

Fig. 2: Marine heatwaves attribution.
figure 2

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.

Full size image

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 30-day 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 (Frw = Fr(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 forcing40. 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 first-year ice region hinders the thorough mixing and downward dispersion of energy from solar radiation8,41, leading to unusually high SST in this area. This extensive coverage offers ample opportunity for the occurrence of MHW events.

Fig. 3: The underlying mechanisms of marine heatwaves.
figure 3

a marine heatwaves (MHWs) maximum intensity (gray bars; second Y-axis), and percent anomalies in cumulative annual solar heat input ([Frw = Fr(1−α)(1−SIC)]; blue curve). b MHWs maximum intensity (gray bars; second Y-axis), and the rate of early summer (June–July) sea ice melt in km2 day−1 (blue curve). c Net outgoing fluxes (longwave upward (LWup) + sensible heat (SH) + Latent heat (LH)) over 1982–2021 averaged over the red box. d MHWs maximum intensity (gray bars; second Y axes), and anomalies in net downward radiative fluxes in June/July (blue curve). Linear trend over 1996–2021 in e the annual cumulative solar heat input in % decade−1 (Frw = Fr(1−α)(1−SIC)). f The number of open water days (sea ice < 15%; days decade−1). g the outgoing fluxes in October (longwave upward (LWup) + sensible heat (SH) + latent heat (LH)) in W m−2 decade−1. All the time-series are average over the red box in panel c (Arctic marginal Seas).

Full size image

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 × 102 km2 day−1 in 1996 to 25 × 103 km2 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 flux41,42,43. The observed net surface radiative fluxes in June/July strongly co-varies 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/m2/decade from 1996 to 2021. Specifically, during the 2020 MHW the outgoing surface fluxes were +15 W/m2 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).

Fig. 4: Detection of GHG signal in observed changes of SST.
figure 4

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 (ai ≠ 0 ⋂ ai = 1, with < 5% risk of error).

Full size image

The response of SST to GHG response is derived from two ensembles of model simulations: the MPI-ESM-LR 30-member30 and the CanESM531 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 variability44 such as the Pacific Decadal Oscillation45, and Arctic Oscillation46, 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 factors47. 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 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 (ai ≠ 0 ⋂ ai = 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 global-mean temperature48,49,50, atmospheric CO2 concentration51, and cumulative anthropogenic CO2 emissions52. 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 radiation8 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 used53,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 temperature9. 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 deviation9. 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 (heatwaveR57).

Open-water period calculation

As in several past studies58,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/NSIDC62) 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 fluctuations58,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 (Frw) depends on the incident solar irradiance (Fr), the ice concentration (C), and the albedo of the ocean (a) and can be expressed as Frw = Fr(1−α)(1−C). The ocean albedo is set to 0.0743,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. Frw 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 Qw 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 factors47,67, which is a linear regression model as follows:

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

where yobs represent the observations and each xi is the modeled response to one of the m forcings that are anticipated by climate models. ai is an unknown scaling factor. The noise on yobs, denoted by uobs, is assumed to represent internal climate variability, while the noise on xi, denoted by ui, is a result of both internal variability and the finite ensemble used to estimate the model response. Uncertainties in ai are estimated by accounting for the effect of internal climate variability on yobs, using time-evolving (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 ensemble30 and the CanESM5 50-member ensemble31 forced with historical well-mixed 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 K68 for MPI-ESM-LR and 5.6 K for CanESM531,69.

Detection of a climate change signal occurs if the uncertainty range around a scaling factor ai is shown to be significantly different from zero. This is handled by testing the null hypothesis HDE: a = 0. If the null hypothesis HDE is rejected, it indicates that the observed change yobs cannot be explained by internal variability uobs alone. Once detection has been established, attribution is assessed by testing the null hypothesis HAT: a = 1. When there is insufficient evidence to reject HAT, the attribution of changes to the respective forcing is claimed70,71.

Time-evolving (historical) internal variability

To have an estimate of the natural internal variability of SST, of which the Pacific Decadal Oscillation72,73, and Arctic Oscillation46 are part, we utilize the 30 members of MPI-ESM-LR30, and the 50 members of the CanESM531 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 approach28,74,75 based on the Causal counterfactual theory38,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}_{{\rm {fixGHG}}}^{{\rm {duration}}}\), and the probability of exceeding the threshold with GHG forcing is denoted as \({P}_{{\rm {ALL}}}^{{\rm {duration}}}\). Similarly, the probability of the threshold (intensity of the observed MHW) being exceeded without GHG forcing is denoted as \({P}_{{\rm {fixGHG}}}^{{\rm {intensity}}}\), and the probability of exceeding the threshold with GHG forcing is denoted as \({P}_{{\rm {ALL}}}^{{\rm {intensity}}}\). The estimated probabilities are used to calculate two event-attribution metrics38,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 world77. The CESM1-LE simulations are running from 1920 to 2100. From 2006 to 2100 the Representative Concentration Pathway 8.5 forcing (RCP8.578) 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 studies26,79,80, and generally performs well. In addition, Arctic sea ice extent81 and sea ice thickness82 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

$${\rm {PN}}=1-\frac{{P}_{{\rm {fixGHG}}}^{{\rm {intensity}}}}{{P}_{{\rm {ALL}}}^{{\rm {intensity}}}};\quad {\rm {PS}}=1-\frac{1-{P}_{{\rm {ALL}}}^{{\rm {intensity}}}}{1-{P}_{{\rm {fixGHG}}}^{{\rm {intensity}}}}$$
(2)

and similarly for MHWs cumulative intensity:

$${\rm {PN}}=1-\frac{{P}_{{\rm {fixGHG}}}^{{\rm {Cum-intens}}}}{{P}_{{\rm {ALL}}}^{{\rm {Cum-intens}}}};\quad {\rm {PS}}=1-\frac{1-{P}_{{\rm {ALL}}}^{{\rm {Cum-intens}}}}{1-{P}_{{\rm {fixGHG}}}^{{\rm {Cum-intens}}}}$$
(3)

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 resampling27 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 (DOISST32,83) Version 2.1 bias-corrected and improved84 product with 0.25° resolution available from September 1981. In ice-covered regions, the in situ and satellite-advanced very high-resolution radiometer (AVHRR) SSTs were blended with proxy SSTs from ice concentrations85. Arctic sea surface temperatures (SSTs) are estimated mostly from satellite sea ice concentration (SIC) estimates85. 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 algorithm62,86,87, distributed by the National Snow and Ice Data Center (NSIDC). The surface heat fluxes are derived from ERA588 reanalysis.

ERA5 simulates observed atmospheric profiles more accurately than other reanalysis data sets, such as ERA-Interim, JRA-55, CFSv2, and MERRA-289,90. The ERA5 data have been used to study Arctic climate and sea-ice concentration change in many previous studies (e.g., refs. 17,91), 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 predecessors88,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-LE29), the Max Planck Institute-Earth System Model with 30 members (MPI-ESM-LR30), and the Canadian Earth System Model with 50 members (CanESM531). 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 well-mixed 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 as93

$${\rm {GH{G}}}_{{\rm {i}}}={\rm {fixGH{G}}}_{{\rm {i}}}-2({\rm {fixGH{G}}}_{{\rm {em}}})+{\rm {AL{L}}}_{{\rm {em}}}$$
(4)

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