Introduction

Since the 1990s the Arctic has been warming more than twice as fast as the global average1, accompanied by rapid loss of sea ice2. This is consistent with polar amplification of climate change and is expected to continue in response to anthropogenic emissions of greenhouse gases3. Over the same period, winter temperatures over mid-latitude northern continents, especially Eurasia, have unexpectedly remained steady or cooled, with an apparent increase in severe winter weather4,5,6,7. The possibility that Arctic warming promotes more severe mid-latitude winters, by altering the atmospheric circulation, has been the subject of intense scientific debate4,5,8,9,10,11,12,13,14,15,16,17,18,19. Observational studies have suggested a clear link between Arctic sea ice loss and mid-latitude winter severity5,6,20,21,22,23,24,25,26, but dedicated numerical model experiments, which are essential to establish causality and to understand the physical mechanisms, are inconclusive, with some simulating mid-latitude cooling in response to Arctic sea ice loss27,28,29,30,31,32,33,34 and others not supporting this link16,18,35,36,37,38,39,40.

For Arctic warming to promote cooling over mid-latitudes would require changes in atmospheric circulation involving a weakening of mid-latitude westerly winds38, consistent with a negative phase of the North Atlantic Oscillation (NAO), and/or a strengthening of the Siberian High41,42. Hence, understanding and quantifying the mid-latitude atmospheric circulation response to Arctic sea ice loss is critical, but there is currently little consensus in modelling studies: the full spectrum of NAO responses has been reported including negative NAO17,27,29,30,32,41,43,44, positive NAO45,46,47,48,49,50, little response19,51,52,53,54 and a response that depends on the details of the forcing28,31,36,55,56,57,58 or the background state of the climate system59,60,61.

There are many potential reasons why previous modelling results are inconsistent, including the use of different magnitudes and patterns of imposed sea ice changes, treatment of oceanic feedbacks41, different models, and whether the simulated responses can be distinguished from internal variability. To overcome some of these limitations, the Polar Amplification Model Intercomparison Project62 (PAMIP) contribution to the sixth Coupled Model Intercomparison Project63 (CMIP6) proposed a set of coordinated experiments. Here, we analyse a large ensemble of PAMIP experiments consisting of more than 3000 simulations from 16 different models and find that all models simulate a weakening of mid-latitude tropospheric westerly winds in response to projected Arctic sea ice loss. We elucidate the main physical processes and show that the model spread depends on eddy feedback. This is 1.2 to 3 times too weak in the models, suggesting that the real-world weakening of westerly winds lies towards the higher end of the model simulations. We also show that observed relationships between Arctic sea ice and atmospheric circulation have weakened recently and are no longer inconsistent with those in models. However, the modelled response to Arctic sea ice loss is weak relative to inter-annual variability, though it is similar in magnitude and offsets the projected response to increased greenhouse gases.

Results

Multi-model response

The atmospheric response to future Arctic sea ice loss is diagnosed from two sets of global atmospheric model simulations (Methods). The first set (present-day) is driven by sea surface temperatures (SSTs) and sea ice concentrations (SICs) representing the present-day climate. The second set (future-Arctic) is the same except that Arctic sea ice and coincident SSTs are replaced with values expected if global temperatures rise by 2 ∘C. By construction, the difference (future-Arctic minus present-day) provides the model-simulated response to future Arctic sea ice loss. We assess 16 model simulations each with between 98 and 300 ensemble members (Table 1) and forced with the same SSTs and SICs. We focus on the boreal winter season (December, January and February, DJF) for which the imposed sea ice changes show reductions around the edges of the ice pack, especially in the Barents-Kara Seas, Sea of Okhotsk, Bering-Chukchi Seas, and Hudson Bay and Labrador Sea (Fig. 1a). In winter, sea ice insulates the atmosphere from the warmer ocean. Hence warm SSTs are imposed where sea ice is lost in future, producing local maxima of near surface warming response in these regions along with further warming spread throughout the Arctic and into lower latitudes (Fig. 1b).

Table 1 Models, ensemble sizes and data used.
Full size table
Fig. 1: Winter response to future Arctic sea ice loss.
figure 1

a Imposed winter sea ice concentration difference (%). b Near surface temperature (TAS) response (K). Note that surface temperature changes are imposed in regions of sea ice loss. c Mean sea level pressure (MSLP) response (hPa). All plots show the winter (December, January, February, DJF) mean, and responses are for the multi-model ensemble mean (calculated as the unweighted average of all ensemble members). Stippling indicates where the multi-model ensemble mean response is significant (95% confidence interval). Black (grey) stars indicate where 100% (90%) of the individual models agree on the sign of the response.

Full size image

In the multi-model mean, there are minima in mean sea level pressure (MSLP) response situated over the regions of largest sea ice loss (Fig. 1c), consistent with a thermodynamic heat low response to surface warming. However, there is also a ridge of high pressure extending from Greenland to Siberia with low pressure further south, producing a response that projects onto a negative NAO and a strengthened Siberian High. Although these features are statistically significant in the multi-model ensemble mean (stippling in Fig. 1c), there is disagreement between individual models on the sign of the response in many regions (grey stars show where 90% of models agree).

Consistencies and differences among the models in the dynamical response are further illustrated in the zonally averaged zonal wind (\(\bar{u}\), where the overbar denotes the zonal average) response as a function of latitude and height (Fig. 2). In the troposphere, there is a very robust equatorward shift, with a weakening of zonal winds around 55–65∘N and a strengthening around 30–40∘N simulated by all models. However, the response is much less coherent in the stratosphere, with some models simulating a significant weakening but the majority showing an insignificant response of either sign. This suggests that a stratospheric pathway highlighted in some studies5,30,33,64,65,66 is not essential for the sign of the tropospheric and surface response. However, it could act to modulate the magnitude of the surface response, as discussed later. Even in the troposphere where the sign of the response is robust, regions of statistical significance (stippled in Fig. 2) are not consistent across models and the strength of the response varies greatly between models. This raises the key question of what the real-world response would be, and whether the differences between models can be understood in order to derive a constrained estimate. Further progress therefore requires understanding the physical processes.

Fig. 2: Consistent tropospheric response to Arctic sea ice loss.
figure 2

Zonally averaged DJF zonal wind response (\(\bar{u}\), ms−1, where the overbar denotes the zonal average) plotted as a function of latitude (°N) and height (pressure) for the ensemble mean of each of the models. The boxes show the regions used to compute the zonal wind response index (ZWRI). Stippling indicates where the ensemble mean response is significant (95% confidence interval). Contours show the climatological zonal mean winds (contour interval 5 ms−1 with negative contours dotted).

Full size image

Physical processes

We focus on explaining the robust multi-model average zonal-mean temperature and zonal wind responses in the troposphere. The largest zonal-mean warming (\(\bar{T}\)) occurs close to the surface over the Arctic (Fig. 3a), consistent with direct heating of the atmosphere by the imposed local surface warming associated with the loss of sea ice. There is a second warming maximum over the Arctic in the lower stratosphere (centred around 100–200 hPa) which has previously been suggested to be heated directly by energy transported from below5. However, in the model simulations there is a meridional overturning circulation response with air descending at high latitudes in the lower stratosphere (arrows in Fig. 3a) indicating that this region warms adiabatically as part of the dynamical response67. Furthermore, warming in this region is not robust across the models. Hence, the robust weakening of the mid-latitude tropospheric winds is unlikely to be solely caused by a simple reduction of meridional temperature gradient by heating of the high latitude atmosphere from below.

Fig. 3: Zonal mean response.
figure 3

Latitude-height cross sections of the multi-model ensemble mean zonally averaged DJF response of a temperature (\(\bar{T}\), colours, K, where the overbar denotes the zonal average) and transformed Eulerian mean (TEM, Methods) circulation (arrows), b zonal wind (\(\bar{u}\), ms−1), and c upward Eliassen-Palm (EP) flux (Fp, colours, standard deviations) and EP flux vectors (arrows representing Fϕ and Fp), d divergence of northward EP flux (∇ϕFϕ, standard deviations). Stippling as in Fig. 1. To aid visualisation the TEM circulation and EP fluxes are standardised by dividing by the internal variability of the present-day simulations (Methods).

Full size image

A meridional overturning circulation response is also seen in the mid-latitude lower troposphere (Fig. 3a) and, as we show below, is important for understanding the physical mechanisms and explaining the spread in modelled responses. This circulation is thermally indirect, with air rising over relatively cool surface conditions between 35–50∘N and descending over the warmer surface between 55 and 70∘N. It is therefore not a direct thermally driven response to the imposed high latitude warming and, as we show below, highlights a key role for changes in wave (eddy) activity.

Climatologically, transient wave activity is mostly generated near the surface by baroclinic eddies in the storm tracks around 40–60∘N and propagates vertically and meridionally until the waves break and dissipate68. Waves flux angular momentum in the opposite direction to their propagation, accelerating zonal winds in their source regions and decelerating zonal winds in their dissipation regions. Wave propagation may be depicted by Eliassen-Palm (EP) fluxes (Methods equation 6), with regions of EP flux divergence indicating where waves act to accelerate the zonal flow and regions of EP flux convergence indicating where waves act to decelerate the zonal flow (Methods equations 2 and 5).

The multi-model mean wave activity response (Fig. 3c) mainly consists of an equatorward shift of upward EP flux (Fp) in the troposphere, with increased Fp around 35–55∘N and decreased Fp around 60–75∘N, and an increased northward EP flux (Fϕ) in the mid to upper troposphere connecting these latitudes. The resulting divergence of Fϕ accelerates the zonal wind around 30–40∘N and convergence of Fϕ decelerates the zonal wind around 55–65∘N (Fig. 3d), consistent with the robust equatorward shift of the winds noted above. We will show below that the strength of the EP flux response is important for explaining the model spread, and focus first on understanding its origins.

Several previous studies30,31,32,36,69,70,71 have highlighted an increase in Fp in response to Arctic sea ice loss. However, the cause is unclear given that the reduced meridional temperature gradient between the equator and the pole (Fig. 3a) would weaken the baroclinic generation of eddies in the storm tracks72,73 and hence reduce Fp, as seen in other studies43,61,67,74. There is a small increase in Fp north of 80∘N (Fig. 3c) that is consistent with zonal asymmetries in near surface temperature and sea level pressure response over regions of sea ice loss (Fig. 1), but the main signal is an equatorward shift resulting in an increase in Fp around 35–55∘N that also extends into the stratosphere. Hence, understanding the main physical processes involved in this equatorward shift is key for understanding the mid-latitude atmospheric response to Arctic sea ice loss.

Although causality cannot be unequivocally determined, much insight can be gained by considering the evolution from autumn to winter, focussing on processes that are robustly simulated across the models. In October there is a weakening of both \(\bar{u}\) and Fp around 50–75∘N (Fig. 4b, c) consistent with the imposed reduction in meridional surface temperature gradient, whereas an equatorward shift is first seen in November (Fig. 4f, g) and develops into the DJF pattern (Fig. 3). To explain this evolution we highlight the following processes which appear to be simulated by the majority of models, though we note that they do not explain every aspect of the circulation response and other processes75 are also likely operating:

Fig. 4: Physical mechanisms of zonal mean response.
figure 4

Latitude-height cross sections of the multi-model ensemble mean zonally averaged October response of a temperature (\(\bar{T}\), colours, K) and TEM circulation (arrows), b zonal wind (\(\bar{u}\), ms−1), c upward EP flux (Fp, colours) and EP flux vectors (arrows representing Fϕ and Fp), d divergence of upward EP flux (∇pFp). (e,f,g,h) As (a, b, c, d) but for November. Stippling as in Fig. 1. To aid visualisation the TEM circulation and EP fluxes are standardised by dividing by the internal variability of the present-day simulations (Methods).

Full size image
  1. 1.

    Reduced zonal wind shear and eddy formation. Arctic warming decreases the surface meridional temperature gradient in October (Fig. 4a), reducing the wind shear on the poleward side of the jet (around 60–70∘N, Fig. 4b) via the thermal wind relation (Methods equation 1). Reduced wind shear reduces baroclinic eddy formation, weakening the storm track and reducing Fp at the surface around 50–75∘N (Fig. 4c).

  2. 2.

    Meridional overturning circulation. A lower tropospheric mid-latitude meridional overturning circulation anomaly (Fig. 4a, e, highlighted above) develops between October and December such that the resulting flow maintains thermal wind balance and is consistent with changes in eddy activity. Some aspects of this circulation can be understood by considering that the reduced Fp at the surface around 50–75∘N (Fig. 4c, g) results in a positive ∇pFp immediately above (Fig. 4d, h) because the flux into this region reduces more than the flux out of this region. An increase in ∇pFp must be balanced (Methods equation 2) by an increase in zonal wind and/or an equatorward flow (negative \({\bar{v}}^{* }\)). However, zonal wind tends to be reduced in response to the imposed weakening of the surface temperature gradient (Fig. 4b) and ∇pFp is at least partly balanced by an equatorward flow near the surface around 45–55∘N (Fig. 4a, e). To maintain mass continuity, a meridional circulation develops with ascent further equatorward (35–50∘N), poleward flow in the mid to upper troposphere (50–70∘N) and descent around 65–75∘N. Adiabatic cooling of the ascending air around 35–50∘N would tend to enhance the meridional temperature gradient and, via thermal wind, strengthen the wind shear further to the south, and reduce the meridional temperature gradient and weaken the wind shear further to the north. Hence, this meridional circulation tends to promote an equatorward shift of the storm track and source of Fp, further reducing Fp on the poleward side of the jet and producing an increase in Fp on the equatorward side as seen in the DJF pattern of wave activity response (Fig. 3c).

  3. 3.

    Positive eddy feedback. The wave activity response (Fig. 3c) results in a divergence of Fϕ on the equatorward side of the jet and convergence of Fϕ on the poleward side of the jet (Fig. 3d) that reinforces the equatorward shift of the storm track and hence enhances the response. Increased Fp can also propagate into the stratosphere, weakening the polar vortex and subsequently further enhancing the equatorward shift by weakening the tropospheric winds on the poleward side of the jet, especially in late winter.

Sensitivity across models

We further assess the processes described above by investigating whether they explain the sensitivity of the response across the models. To quantify the response for each model, we define a zonal wind response index (ZWRI, Methods, the difference between the zonally averaged zonal wind responses in the boxes in Figs. 2 and 3b). The correlation between ZWRI and zonal mean temperature response across the models (Fig. 5a) is consistent with a sensitivity to the strength of the meridional circulation (process 2, Fig. 3a): models with a larger ZWRI show enhanced lower tropospheric adiabatic cooling by ascent around 35–50∘N and warming by descent around 65–75∘N compared to models with a smaller ZWRI. In addition, models with a larger ZWRI show greater mid to upper tropospheric warming due to enhanced southerly advection around 50–70∘N compared to models with a smaller ZWRI.

Fig. 5: Sensitivity of zonal mean response.
figure 5

Latitude-height cross sections of the correlation across models between ZWRI and response in (a)\(\bar{T}\), b ∇ϕFϕ and c refractive index (n2). All data are for DJF.

Full size image

ZWRI is also correlated with the divergence of Fϕ around 30–40∘N and convergence of Fϕ around 45–65∘N (Fig. 5b) consistent with a positive feedback from wave driving in determining the magnitude of the simulated response (process 3). The tropospheric response may also be enhanced by eddy feedback involving the stratosphere (process 3): increased Fp enters the stratosphere around 35–55∘N (Fig. 3c) and ZWRI is correlated with increases in the refractive index (n2, Methods equation 7) in these latitudes (Fig. 5c), allowing more waves to weaken the stratosphere and affect the troposphere via downward propagation, strengthening the negative NAO and equatorward jet shift through further eddy feedback76.

Constrained response

The model spread can potentially be exploited to obtain an estimate of the true response using the concept of emergent constraints (Methods). In this approach, the key physical processes that explain the differences between modelled responses must be understood and related to quantities that can be observed. If a robust relationship exists, then the true response may be inferred by comparing observations of these quantities with those in the models. Results above suggest an important role for eddies: the DJF zonal wind response is initiated by thermal wind balance but is sensitive to eddy feedback which alters ∇ϕFϕ in the troposphere (Figs. 3d and 5b) and potentially involves the stratosphere via changes in Fp and n2 (Figs. 3c and 5c). Eddies are generated by the mean flow but can also feedback onto the mean flow in such a way that increases or reduces their generation (i.e. a positive or negative feedback77,78). Hence, we hypothesise that the response to sea ice loss may be related to the strength of the eddy feedback simulated by the different models.

We estimate eddy feedback by examining the role of eddies in driving internal variability (Methods). For each model, we compute the local correlation across the ensemble members between DJF zonal mean zonal wind (\(\bar{u}\)) and the divergence of northward EP flux (∇ϕFϕ) in the present-day simulations. Averaged across the models, this correlation is positive throughout most of the troposphere (Fig. 6a) consistent with acceleration of the mean flow by a convergence of eddy momentum flux. The square of this correlation represents the fraction of the total variability of \(\bar{u}\) that is explained by variations in ∇ϕFϕ and it would be expected to increase as the eddy feedback becomes more strongly positive. The squared correlation varies greatly across the models (Fig. 6b), suggesting a wide range of eddy feedbacks simulated by the different models.

Fig. 6: Eddy feedback.
figure 6

a Latitude-height cross section of the multi-model mean local correlation between DJF \(\bar{u}\) and ∇ϕFϕ. Correlations are computed across ensemble members for the present day simulation for each model separately, and then averaged to make the multi-model mean. b Fraction of mid-upper troposphere zonal wind variance explained by ∇ϕFϕ as a function of latitude (the square of the local correlations in a averaged over 600 to 200 hPa) for each of the models (coloured curves) and in the reanalyses (black curves). Reanalyses values are computed from interannual time series over the period 1979–2016 inclusive. The eddy feedback parameter (M) is computed as the average over latitudes 25–72∘N (shown by the box in a).

Full size image

We define a measure of eddy momentum feedback strength (M) as the variance of DJF \(\bar{u}\) explained by ∇ϕFϕ (i.e. the squared correlation) averaged over the mid to upper extratropical troposphere in the region shown by the box in Fig. 6a (Methods). We find that M is positively correlated with ZWRI across models (r = 0.49, p = 0.03, Fig. 7a), supporting our hypothesis that the response to sea ice loss is strengthened by eddy feedbacks and allowing a constrained estimate of the zonal wind response to be obtained via the ensemble regression approach (ER, Methods). We diagnose the observed eddy feedback using three reanalyses and computing the correlations in time rather than across ensemble members (Methods). The observed eddy feedback is between 1.2 and 3 times larger than that in any of the climate models. Using the observed eddy feedback strength to scale the zonal wind responses in the models, we find that the multi-model ensemble mean ZWRI increases from 0.7 (± 0.1) ms−1 to 0.9 (± 0.2) ms−1 (95% confidence intervals) (Fig. 7a, green shading), suggesting that the real world response would lie towards the higher end of the model simulations.

Fig. 7: Emergent constraints.
figure 7

a Emergent constraint based on the ensemble regression (ER) between eddy momentum feedback and the zonal wind response index (ZWRI). Black line shows the regression with hatching showing the 95% confidence interval. Horizontal green line shows the constrained ensemble mean response, with the shading showing its 95% confidence interval (Methods). Vertical black line and grey shading shows the mean and range of eddy feedback from the reanalyses. Ellipses show the 95% uncertainties obtained by bootstrapping with replacement the ensemble members. b As a but for the stratospheric polar vortex (SPV) response. A one-sided test is used to calculate p values since we expect the response to increase as eddy feedback strengthens. All data are for DJF.

Full size image

We now apply the ER method in a more general sense to assess the response to sea ice loss in different quantities and regions: for any variable and location, we regress our estimate of eddy momentum feedback strength against the response in that variable and use the observed eddy feedback to diagnose the constrained response. Note, that ER will have little impact where the regression is small. We find that ER enhances the zonal wind response throughout the atmosphere, including the stratospheric polar vortex (Fig. 7b, r = −0.44, p = 0.04) which increases in magnitude from 0.1 to −0.6 ms−1 in the simple ensemble mean (EM) to −0.1 to −1.7 ms−1 in ER, accompanied by an increase in regions showing a significant zonal wind response (compare Fig. 8c, d). This is consistent with an increased refractive index response in ER particularly around 35–40∘N (Fig. 5c) in the upper troposphere, which allows more waves to propagate into the stratosphere and weaken the SPV (Fig. 8a, b). It is well established that changes in the stratosphere propagate downwards and affect the troposphere and surface winds76, and this stratospheric pathway likely enhances the near surface wind response. Overall, ER results in a greater reduction of the SPV (Fig. 8c, d), stronger equatorward shift of the storm tracks (Fig. 8e, f), and a stronger negative NAO response (increased from −0.4 to −1.2 hPa in EM to −0.3 to −2.1  hPa in ER). However, ER does not constrain the Eurasian temperature response (Eurasia T, Methods), which ranges from −0.2 to +0.4 ∘C.

Fig. 8: Impact of emergent constraint.
figure 8

Refractive index (n2) response (standardised) for a the multi-model ensemble mean (EM) and b the ensemble regression (ER). c, d As a, b but for stratospheric zonal wind (ms−1) at 10 hPa. e, f As a, b but for near-surface zonal wind (ms−1) at 925 hPa. Arrows in a and b show the EP flux response (standardised). Colour bars represent each column. Stippling shows where the EM or ER response is significant (95% confidence interval, Methods). All data are for DJF.

Full size image

We assess the robustness of our proposed constraint to several sources of uncertainty. The simulated response to Arctic sea ice loss is small relative to internal variability, leading to substantial uncertainties in individual models (ellipses in Fig. 7). We therefore further assess the statistical significance of the regressions by bootstrapping the individual members for each model with replacement. This is repeated 1000 times and results in p values that are slightly less significant (p = 0.08 and 0.06 for ZWRI and SPV respectively), highlighting the need for very large ensembles to obtain robust results79. We also tested the sensitivity of ER to outliers80 by removing each model in turn and repeating the regression. This is most sensitive to removing E3SMv1 and CanESM5, which increases the p values to 0.16 and 0.07 respectively for ZWRI, and to 0.06 and 0.10 respectively for SPV. The calculation of EP fluxes can be sensitive to the frequency of data used, which ranges from 20 min to daily means depending on the model (Table 1) and 6 h in the reanalyses81. We tested this sensitivity by recalculating the regressions (Fig. 7) for subsets of models with similar sampling frequencies (20–30 min, 6 hourly, and daily) and found the sign to be the same for each subset as for the full set (though the regressions were no longer significant at the 95% confidence level). The eddy feedback parameters computed from the Reanalyses are based on 38 years compared to more than 100 members for the models. However, the three reanalyses are in good agreement and together provide a similar sample size to the models. Finally, the eddy feedback parameters for reanalyses are based on time series that include the effects of coupled modes of internal variability such as El Niño, in addition to changes in natural and anthropogenic radiative forcings that are not present in the model experiments. We recalculated the eddy feedback using time series from AMIP simulations that also include these factors for five models (CanESM5, HadGEM3-LL, HadGEM3-MM, IPSL-CM6A-LR, MIROC6) and found a small increase (0.09 on average, compared to a multi-model mean of 0.30), but that all model values remain lower than the reanalyses.

Our proposed emergent constraint is therefore reasonably robust and is strongly linked to the physical processes. However, we note that other emergent constraints are possible. For example, ZWRI is correlated with the background SPV (r = 0.50, p = 0.02, Supplementary Fig. 1) but in this case the observations are near the middle of the model range so that ER is close to the simple ensemble mean. Had we only considered this relationship we would have erroneously constrained the response to be near the EM, highlighting the importance of basing constraints on the physical processes61,82. It is also possible that the response could be further enhanced by coupling with the ocean17,41,44, albeit with greater variability79, and/or by the pattern of sea ice loss36,57,58,83, and could depend on the phase of the quasi-biennial oscillation84. However, initial results from a subset of models for PAMIP experiments that are designed to investigate the effects of coupling and the pattern of sea ice loss (Methods) do not show consistent differences in ZWRI across models compared to the standard experiments. This suggests that these effects may not be large, but results from more models are needed for further assessment, and could provide an important out of sample test82 of the emergent constraint proposed here.

Discussion

We have analysed the boreal winter atmospheric circulation response to future Arctic sea ice loss using a very large multi-model ensemble, comprising 16 different models and more than 3000 ensemble members with the same experimental protocol. We find a robust response in the troposphere, with a weakening of mid-latitude tropospheric westerly winds and an equatorward shift of the storm tracks, consistent with a negative phase of the NAO. However, the strength of the modelled response varies between models and is proportional to the strength of the eddy momentum feedback, enabling a constrained estimate of the real-world response to be obtained. Since all the models underestimate the observed eddy feedback strength, the real-world response is likely to be at the higher end of the model range. The stratosphere response is not consistent across the models and therefore not essential for determining the sign of the tropospheric circulation response to sea ice loss. However, the eddy feedback constraint indicates a robust weakening of the stratospheric polar vortex, suggesting that in the real world the stratosphere may play an active role that amplifies the surface response, consistent with other studies30,33,36,64,66,74.

Our emergent constraint suggests that models may underestimate the response to sea ice loss and is consistent with other evidence that models underestimate the predictable fraction of NAO variability in seasonal85,86,87, interannual88, and decadal forecasts89,90,91, and in historical climate simulations92,93,94. This has been referred to as the signal-to-noise paradox95 since models are unexpectedly able to predict the real world better than they can predict one of their own ensemble members. We find that the eddy momentum feedback is underestimated in all of the models, suggesting that this could be a potential cause of this model error as previously suggested96. Possible reasons for underestimated eddy feedback include errors in wave propagation, imperfect representation of interactions with the mean flow, and unresolved waves, but further analysis is left for future work. We also note that our emergent constraint could potentially explain some of the model spread in jet shifts forced by a variety of factors, including changes in greenhouse gases, but this is beyond the scope of the present study.

Diagnosing the response to sea ice loss is not possible from observations alone since causality cannot be established61. Nevertheless, there is a perception that models and observations do not agree5. Our results suggest that model errors could be responsible for some of the apparent inconsistencies and that taking these into account yields a robust modelled response in the troposphere and stratosphere. However, it is important to note that the short observational record contains considerable sampling uncertainty97,98. For example, the period 1979 to 2012 suggests a strong link between autumn sea ice and winter atmospheric circulation and Eurasian temperature (Fig. 9a, c, e) consistent with previous studies6,20,21,22,23,24,25. If we extend the analysis period by adding the most recent 8 years of data (2013–2020), we find that the observed relationships between autumn sea ice extent and DJF NAO, SPV and Eurasian T indices are weaker in magnitude (Fig. 9b, d, f). This is consistent with recent evidence that the observed relationships are modest7 and intermittent66 and may not reflect a causal link with Arctic sea ice99,100,101.

Fig. 9: Weakened observed relationships.
figure 9

a Observed winter (DJF) NAO anomaly time-series for the period 1979 to 2012 (black) along with the variability that is linearly related to autumn (September-November) sea ice extent in the Arctic (blue) and Barents-Kara (BK) Seas (red). Pearson correlation (r) and regression coefficients (reg, 95% confidence interval, per standard deviation of sea ice extent) are indicated. b As a but for the period 1979 to 2020. c, d As a, b but for SPV. e, f As a, b but for Eurasia T. Indices are defined in Methods.

Full size image

Our model estimates are within, and therefore consistent with, the magnitudes of the observed relationships based on the extended period of record (1979–2020). However, the modelled response to sea ice loss remains relatively weak, amounting to 30% or less of the observed interannual standard deviation (σ) of the NAO and SPV. Assuming linearity and given that the imposed reduction in DJF Arctic sea ice extent in our experiments is around 4σ, this implies that interannual Arctic sea ice variations account for less than 10% of the interannual variations in NAO and SPV, and are thus unlikely to drive large seasonal mean impacts in individual winters. Nevertheless, these values are similar in magnitude to changes in the NAO and SPV expected by the end of the century in response to increases in greenhouse gases65,102,103, and will therefore impact long term projections.

Methods

Model experiments

We assess coordinated experiments from the Polar Amplification Model Intercomparison Project62 (PAMIP). PAMIP experiment 1.1 simulates the present day climate using global atmosphere models constrained at the surface by present day estimates of sea surface temperature (SST) and sea ice concentration (SIC). PAMIP experiment 1.6 is the same as 1.1 except that Arctic SIC is replaced with values expected with global temperatures 2 degrees Celsius warmer than pre-industrial conditions. Where sea ice is lost in 1.6 relative to 1.1, future SSTs are used; elsewhere the SSTs are the same in 1.6 and 1.1. By construction, the difference between 1.6 and 1.1 provides the simulated response to future Arctic sea ice loss and associated local changes in SST. Each experiment starts on 1st April 2000 and runs for 14 months. We analyse results from 16 models, each with at least 98 ensemble members (Table 1) and forced with the same SSTs and SICs. We also analyse a small subset of models to assess the effects of ocean-atmosphere coupling (PAMIP experiments 2.1 and 2.3, EC-EARTH3, HadGEM3-MM, NorESM2-LM) and regional sea ice changes (PAMIP experiments 3.1 and 3.2, CNRM-CM6-1, ECHAM6.3, EC-EARTH3, HadGEM3-MM). All data were re-gridded to the resolution of the coarsest model (3∘ latitude by 3∘ longitude) before comparison.

Observations

We use Eliassen-Palm (EP) fluxes from three reanalyses (ERA-Interim, NCEP-NCAR and JRA-55) available from the Centre for Environmental Data Analysis81. Sea ice observations are taken from HadISST1.1104. Sea level pressure, near-surface temperature and stratospheric wind observations are taken from the ERA5 reanalysis105.

Indices

The North Atlantic Oscillation (NAO) index is calculated as the difference in mean sea level pressure between two small boxes located around the Azores (28–20∘W, 36–40∘N) and Iceland (25–16∘W, 63–70∘N) with the average over the whole time series removed to create anomalies88. To quantify the response to Arctic sea ice loss, we define a Zonal Wind Response Index (ZWRI) that is calculated as the difference in zonally averaged zonal wind response between 30–39∘N and 54–63∘N averaged over 600 to 150 hPa, which correspond to the regions with the largest zonal wind responses in the multi-model mean. The strength of the stratospheric polar vortex (SPV) is computed as the zonal-mean zonal wind averaged over 54–66∘N at 10 hPa. Eurasian temperatures (Eurasia T) are averaged over the region34 60–120∘E, 40–60∘N. The Barents–Kara (BK) Seas region58 is taken as 10–100∘E, 65–85∘N.

Thermal wind

The geostrophic balance between pressure gradient and Coriolis forces, combined with the hydrostatic equation leads to the thermal wind relationship, showing that a reduction in meridional temperature gradient as the Arctic warms is accompanied by a decrease in vertical wind shear68:

$$f\frac{\partial u}{\partial p}=\frac{R}{ap}\frac{\partial T}{\partial \phi }$$
(1)

where u, T, p, ϕ are zonal wind, temperature, pressure and latitude, f is the Coriolis parameter, a is the radius of the Earth and R is the gas constant.

Transformed Eulerian Mean momentum equation and Eliassen-Palm fluxes

Much insight into extratropical atmospheric circulation can be gained from the Transformed Eulerian Mean (TEM) form of the zonal mean quasi-geostrophic (QG) momentum equation68,81:

$$\frac{\partial \bar{u}}{\partial t}=f{\bar{v}}^{* }+\frac{1}{a\cos \phi }{{{{{{{\boldsymbol{\nabla }}}}}}}}{{{{\cdot }}}}{{{{{{{\boldsymbol{F}}}}}}}}+\bar{\epsilon }$$
(2)

where the overbar represents the zonal mean and the * denotes the TEM circulation:

$${\bar{v}}^{* }=\bar{v}-\frac{\partial }{\partial p}\left[\frac{\overline{{v}^{\prime}{\theta }^{\prime}}}{\partial \bar{\theta }/\partial p}\right]$$
(3)
$${\bar{w}}^{* }=\bar{w}+\frac{1}{a\cos \phi }\frac{\partial }{\partial \phi }\left[\frac{\overline{{v}^{\prime}{\theta }^{\prime}}\cos \phi }{\partial \bar{\theta }/\partial p}\right]$$
(4)

where the prime represents departures from the zonal mean, u, v, w are the zonal, meridional and vertical velocities, θ is potential temperature, f is the Coriolis parameter, t is time, \(\bar{\epsilon }\) represents the effects of friction and parameterised processes, and the Eliassen–Palm (EP) flux divergence is given by

$${{{{{{{\boldsymbol{\nabla }}}}}}}}{{{{\cdot }}}}{{{{{{{\boldsymbol{F}}}}}}}}={\nabla }_{\phi }{F}_{\phi }+{\nabla }_{p}{F}_{p}=\frac{1}{a\cos \phi }\frac{\partial {F}_{\phi }\cos \phi }{\partial \phi }+\frac{\partial {F}_{p}}{\partial p}$$
(5)

where the QG northward and vertical EP fluxes are given by

$$\big\{{F}_{\phi },{F}_{p}\big\}=a\cos \phi \left\{-\overline{{u}^{\prime}{v}^{\prime}},\frac{\overline{{v}^{\prime}{\theta }^{\prime}}}{\partial \bar{\theta }/\partial p}f\right\}$$
(6)

Note, that we use the full primitive equation EP fluxes81 in our analysis but the main effects are captured by the QG form given above.

Standardisation

To aid visualisation the TEM circulation and EP fluxes are scaled by dividing by the internal variability of the present-day simulations, computed as the variance across ensemble members, averaged over all of the models, and then square-rooted to obtain the standard deviation. This standardises the magnitude of signals across the depth of the atmosphere and indicates the general sense of the propagation (i.e. north or south, upwards or downwards) but does not necessarily indicate the precise geometric direction.

Refractive index

Insight into the propagation of wave activity can be gained by examining the refractive index106,107:

$${n}_{k}^{2}=\frac{{\bar{q}}_{\phi }}{\bar{u}}-{\left(\frac{k}{a\cos \phi }\right)}^{2}-{\left(\frac{f}{2NH}\right)}^{2}$$
(7)

where

$${\bar{q}}_{\phi }=\frac{2{{\Omega }}}{a}\cos \phi -\frac{1}{{a}^{2}}{\left[\frac{{\left(\bar{u}\cos \phi \right)}_{\phi }}{\cos \phi }\right]}_{\phi }-\frac{{f}^{2}}{{\rho }_{0}}{\left({\rho }_{0}\frac{{\bar{u}}_{z}}{{N}^{2}}\right)}_{z}$$
(8)

is the meridional gradient of the zonal mean potential vorticity, and k, N, H and Ω denote the zonal wave number, buoyancy frequency, scale height, and Earth rotation frequency respectively. Note, that waves are refracted towards high values of \({n}_{k}^{2}\).

Eddy feedback

Eddies are generated by baroclinic processes in the storm tracks and eddy activity propagates horizontally and vertically. Eddies flux angular momentum in the opposite direction to their propagation, driving (accelerating) the zonal wind in their source regions and dragging (decelerating) the zonal wind in their dissipation regions. These interactions with the mean flow can promote or reduce further eddy generation, producing a positive or negative eddy feedback77,78. In simple models, it is possible to switch off the effects of eddies allowing eddy feedback to be diagnosed by comparing to the full response108. Previous studies have assessed eddy feedback in comprehensive models using lagged relationships in daily data focussing on large scale patterns of variability109,110. However, lagged relationships potentially include persistence that is unrelated to feedbacks and may miss fast feedbacks that occur before the diagnosed lag. They also require analysing large volumes of daily data.

Here, we propose a new approach that can be used with seasonal mean data and does not require lags to be specified. We reason that the fraction of seasonal mean zonal wind variability that is related to eddies will increase as the eddy feedback becomes more positive, and measure the eddy feedback strength as the local correlation squared (i.e. the variance explained) between DJF \(\bar{u}\) and ∇ϕFϕ averaged over the mid to upper northern troposphere (25–72∘N, 600–200 hPa, the box shown in Fig. 6a). Although this measure is imperfect since it includes eddy driving, we show that it explains some of the differences in modelled responses. For the model simulations, the eddy feedback parameter is calculated across the ensemble members for the present-day simulations, whereas for the reanalyses it is calculated from time series covering the period 1979–2016. The reanalyses eddy feedback parameters are insensitive to removing a linear trend, but we find some sensitivity to the different ways of calculating eddy feedback using models for which time series data are also available (see main text). However, our main conclusions remain valid.

Constrained estimates

We use emergent constraints82,111,112 to estimate the real world response. We seek an observable quantity (x) that provides a physical explanation for differences in the model estimates of the response to Arctic sea ice loss (yi, the ensemble mean estimate of the response for model i) such that

$${y}_{i}={\bar{y}}_{i}+\beta {x}_{i}+{\varepsilon }_{i}$$
(9)

where the overbar denotes the average of the model estimates, β is the slope of the linear regression between xi and yi, and εi is an identically independently distributed random variable with zero expectation (i.e. noise).

If such a quantity exists, then the simple multi-model ensemble mean (EM i.e. with β = 0) is an inappropriate estimate of the true response because εi would not be independent of x112. Instead, the observable response (yO) may be estimated using ensemble regression112 (ER):

$${y}_{O}=\bar{{y}_{i}}+\beta ({x}_{O}-\bar{{x}_{i}})$$
(10)

where xO is the observed estimate of x. The error variance of the multi-model ensemble mean is

$${S}_{y}^{2}={\sigma }_{\varepsilon }^{2}\left(\frac{1}{n}+\frac{{({x}_{O}-\bar{{x}_{i}})}^{2}}{\mathop{\sum }\nolimits_{i = 1}^{n}({x}_{i}-\bar{{x}_{i}})}\right)$$
(11)

where n is the number of models and \({\sigma }_{\varepsilon }^{2}\) is the variance of εi. The second term in the parentheses accounts for estimation error in regression slope and grows quadratically with the error in the model mean estimate of x.