Introduction

Understanding and quantifying the changes in climate that will occur in the next decades is of vital importance for regional adaptation and mitigation strategies and local adaptation plans. This is because this timescale is of key importance for planning different sectorial activities, such as ensuring water and food security and uninterrupted energy production1,2,3. One important factor shaping the impacts of anthropogenic global warming, particularly in the near future, is the loss of Arctic sea-ice cover. As the disappearance of sea-ice changes the surface albedo, removes the insulation between the atmosphere and surface ocean and affects the salinity profiles, these local changes in turn drive a multitude of atmospheric and oceanic teleconnections that propagate far away from the Arctic4,5. Investigations into the remote impacts of Arctic sea-ice loss have, however, often focused on long term, century-scale changes6,7,8,9,10,11,12,13,14,15,16,17 while far less attention has been given to the near-term, decadal-scale impacts16,17,18,19,20,21,22. Physically, “decadal” responses are dominated by faster adjustments of the climate system to sea-ice changes (e.g., atmospheric reorganizations and fast oceanic responses), while the “centennial” responses allow for slower deep ocean feedbacks to develop too. The distinction between “decadal” and “centennial” time-scales in this context is based on the time the sea-ice forcing has been applied for in the model before studying the resulting climate response.

In addition to the relative knowledge gap for decadal-scale impacts, another major issue is the contradictory findings regarding the remote impacts of Arctic sea-ice loss. For example, sea ice decline in the Barents-Kara seas has been linked to both negative23,24,25,26,27 and positive phases of the North Atlantic Oscillation (NAO)28,29, as well as suggested to have a time-varying relationship with the NAO30. Similarly, over the Pacific sector, some studies have suggested a strengthening of the Aleutian pressure low in response to the loss of Arctic sea-ice cover6,7,8,9, while others found the opposite response, with the loss leading to the development of the North Pacific geopotential ridge and weakening of the Aleutian Low31. These opposing circulation responses can result in a very different outcomes for future water availability over the western United States. If the sea ice loss favors the development of the North Pacific geopotential ridge, this could lead to a major drying in the American southwest, similar to the conditions experienced during the 2012–2016 California drought, affecting the water and energy supply and increasing the risk of wildfires. The strengthening of the Aleutian Low, on the other hand, would favor wetter conditions over California. Understanding the factors that could lead to such different responses to sea-ice changes, is crucial for assessing the likelihood of societally critical impacts over the North American and European continents32,33,34.

In this study, we seek to resolve some of the contradictory findings regarding the anatomy of the climate response to Arctic sea-ice changes. We present the first multi-model assessment of the decadal climate response to the loss of Arctic sea-ice cover obtained strictly using the energy conserving methodologies for isolating the impacts of sea-ice changes. We provide an evaluation of the most important climate impacts of Arctic sea-ice loss expected to occur on these timescales and detail on a robust atmospheric response found across the Pacific sector that is not dependent on the ocean model complexity.

Simulating sea-ice loss

To assess the impacts of changing sea-ice cover, we compare two sets of simulations (control and ‘low ice’) that differ only in the amounts of Arctic sea-ice cover, with all other forcings being the same. Low ice simulations are produced using energy conserving methodologies that do not involve any addition of artificial energy fluxes to melt the sea ice, but instead apply small perturbations to the sea-ice model physics parameter values19,20. We employ three different global climate models: CCSM4-SO, IPSL-CM5A2 and EC-Earth3.2. One of the models used (CCSM4-SO) employs a thermodynamic mixed-layer ocean model (so called ‘slab ocean’) while the other two (IPSL-CM5A2 and EC-Earth3.2) models utilize a full-depth dynamic-thermodynamic ocean model component. In this manner we are also able to assess whether the complexity of the ocean model influences the simulated impacts of sea-ice loss on decadal timescales. Detailed description of models and methodologies used to isolate the impacts of sea-ice loss is provided in the Methods section.

Simulated sea ice loss in ‘low ice’ simulations follows the seasonal cycle of the observed sea-ice loss, with the largest anomalies occurring during the summer and the smallest during the winter season (Fig. 1a–c). The total sea ice area in control simulations, as well as the resulting amounts of sea-ice loss, differ between the different models. For example, CCSM4-SO and IPSL-CM5A2 baseline Arctic sea-ice cover is higher than that in EC-Earth3.2 (black lines in panels a–c). Also, the ‘low ice’ CCSM4-SO and EC-Earth3.2 simulations reach “ice-free” conditions during the summer, while the IPSL-CM5A2 in general exhibits smaller changes when compared to the other two models (green lines in a–c). By covering this wide range of mean sea ice states and their changes (Fig. 1 and Supplementary Fig. 1), these multi-model differences additionally test the robustness of simulated climate responses. In the late autumn, the time suggested to have a strong impact on the following winter’s circulation patterns31, the three models qualitatively agree regarding the pattern of sea ice loss, where the major areas of sea ice decline make a ring encompassing the Barents-Kara, East Siberian, Chukchi, Beaufort and Greenland seas (Fig. 1d–f). The type of sea ice loss forcing applied in the ‘low ice’ experiments corresponds to the forcing expected to occur in nature in the next decades.

Fig. 1: Arctic sea-ice area and sea-ice concentration anomalies.
figure 1

a–c Monthly mean Arctic sea-ice area [×106 km2] in control (black) and low ice (green) simulations, lines denote individual ensemble members; d–f November sea-ice concentration anomalies (low ice minus control) [%]. Model names are shown above each column. Stippling indicates areas where the sea ice concentration changes are significant at the 95% confidence interval (using t-test). Monthly sea ice anomalies for December, January and February are provided in Supplementary Fig. 1.

Full size image

North Pacific response and impacts on the western US

While the greatest magnitude of Arctic sea-ice loss in our model simulations is observed in autumn, its remote impacts occur with some delay and are most pronounced during the winter season. A notable winter Pacific response in all simulations is a positive upper level (250 hPa) geopotential height anomaly over the North Pacific (Fig. 2a–c) that is accompanied by a positive sea level pressure (SLP) anomaly at the surface (Fig. 2d–f). The strongest geopotential height changes are found in CCSM4-SO and EC-Earth3.2, consistent with the largest overall sea-ice reductions in these models' simulations. The centers of the North Pacific geopotential highs, however, are located closer to the North American continent in CCSM4-SO and IPSL-CM5A2, but further to the northeast in the EC-Earth3.2. The exact location of the geopotential highs is critical for its societal impact as a rather small variation in the position of the anomalous ridge can have a substantial impact on the western United States rainfall, resulting in an extremely dry or a very wet winter in the American southwest. Thus, even when models agree on the dynamical response (i.e., establishment of a North Pacific ridge) their simulated rainfall response may still differ due to the variations in the exact position of the ridge. Precipitation decreases (significant at the 90% level, Fig. 2g–i) over the American southwest, with cumulative precipitation for DJF decreasing by 10% (not shown). A precipitation increase in the northwest is suggested in all simulations, although the magnitude and significance of rainfall changes is weaker in EC-Earth3.2 and IPSL-CM5A2 compared to the CCSM4-SO. The rainfall changes are the most prominent during the month of January in all models (Supplementary Fig. 2a–c). Weaker precipitation response and its significance in IPSL-CM5A2 may also be related to the smaller sea-ice reductions in this model (Fig. 1c and f).

Fig. 2: December-January-February (DJF) circulation and precipitation responses.
figure 2

a–c DJF 250 hPa geopotential height anomalies [m]; d–f DJF sea level pressure anomalies [hPa]; g–i DJF precipitation anomalies [mm/month]. All anomalies are the difference between the low ice and control ensemble simulations. Stippling indicates anomalies significant at the 95% confidence level (using t-test) for geopotential and sea level pressure, and 90% for precipitation. Model names are listed above each column.

Full size image

The streamfunction anomalies (shown in Fig. 3a–c) indicate that the North Pacific geopotential ridge is actually part of a much larger atmospheric structure—a Rossby wave train—featuring a large cyclonic cell in the tropical Pacific and an anticyclonic cell in the North Pacific. Plumb fluxes anomalies35, illustrating the ‘low sea ice minus control’ changes in Rossby wave propagation, are shown in Fig. 3d–f. These show an anomalous propagation of wave activity from the Pacific sector towards North America in response to the reduced Arctic sea-ice cover in CCSM4-SO and IPSL-CM5A2. In EC-Earth3.2, wave activity propagation is much less clear at this level.

Fig. 3: DJF streamfunction and Plumb flux responses.
figure 3

a–c DJF 250 hPa streamfunction anomalies [105 m2/s]; d–f Plumb flux anomalies [-]. All anomalies are calculated as low ice minus control. Model names are listed above each set of plots.

Full size image

Previous work has suggested this Rossby wave pattern could be forced by sea-ice induced changes in tropical convection and divergence19,36,37. Changes in upper level divergence, associated with deep tropical convection changes have been shown to act as Rossby wave sources38,39,40,41. Earlier studies have also stressed the importance of certain regions in the tropical Pacific42 as starting points for the Rossby wave trains that propagate across the North American continent and into the Atlantic. For example, the large-scale tropical convection changes during the El Nino Southern Oscillation events produce a well-established extra-tropical Rossby wave response43, with the cold phase (La Niña) displaying a Rossby wave pattern very similar to that shown in Fig. 3a–c. In our simulations, areas of increased outgoing longwave radiation (OLR) in the tropical Pacific (Fig. 4a–c) indicate a large-scale convection reorganization that is also associated with the upper level divergence changes (Supplementary Fig. 3d–f). These tropical changes are central to the establishment of the Rossby wave train seen in our simulations.

Fig. 4: DJF tropical Pacific response.
figure 4

a–c U surface wind component [m/s]; d–f Latent heat flux (positive upwards) [W/m2]; g–i Outgoing Longwave Radiation (OLR) [W/m2]. Stippled areas indicate the anomalies significant at the 90% confidence level. (Note that the background easterly flow is negative, meaning that the negative anomalies are strengthening the easterly flow).

Full size image

Links between high latitude temperature changes and tropical rainfall and convection changes have been proposed in multiple contexts of past and future climate change44,45, with different mechanisms of propagation suggested for different timescales considered. In the following we describe the decadal-scale response.

Tropical convection response and the importance of tropical surface changes

Our experiments were designed to study the decadal-scale response to sea ice changes. While the initial propagation of sea-ice loss induced perturbations has been a major focus of PAMIP style experiments46,47, the response shown in our study represents the time after the initial signal of sea-ice forcing has already arrived. This developed decadal-scale response in our simulations features some major tropical Pacific anomalies. Zonal wind anomalies shown in Fig. 4a–c illustrate weakened northeasterly zonal flow just off the equator in the eastern tropical Pacific and strengthening of the southeasterly zonal flow across the central and western tropical Pacific. The areas of stronger winds coincide with the areas of increased anomalous latent heat flux into the atmosphere (Fig. 4d–f). Aloft, substantial OLR changes are present, signaling large-scale convection changes (Fig. 4g–i). Similar to the patterns accompanying the tropical Pacific convection changes during the ENSO events41, changes in the location and intensity of tropical Pacific convection in our experiments drive the northward propagating Rossby wave train responsible for the development of the anomalous geopotential ridge in the North Pacific.

In order to better understand the importance of tropical surface changes in conveying the signal of high latitude sea-ice loss and triggering the tropical convection changes that drive subsequent Rossby wave response, we also analyze two sets of idealized simulations in which tropical (30 °S–30 °N) sea-surface temperatures (SST) are prescribed, while the rest of the ocean operates in the mixed-layer ocean mode (see Supplementary Methods for detailed description). If the tropical SST and surface heat flux changes are essential for producing the sea-ice induced tropical convection response, this suggests that the Rossby wave train (or any downstream response associated with it), cannot exist in modeling configurations that do not allow tropical sea surface temperature changes (e.g., prescribed SST setups). By comparing the simulations featuring the imposed high latitude forcing in combinations with the tropical SSTs fixed to the perturbed values and the experiments featuring imposed high latitude forcing with tropical SSTs fixed to the control values it is possible to infer the importance of tropical surface changes. Perturbed and control SST values are derived from corresponding mixed layer experiments that do not prescribe tropical SSTs.

Comparison of the upper level stream function responses to high latitude warming indicates that the tropical Pacific Rossby wave train (with the anticyclonic cell in the North Pacific of the US west coast) is present only in simulations in which the tropical SST contain the signal of high latitude changes (Supplementary Fig. 4a), and does not exist in simulations in which the SST’s were fixed to their control values (Supplementary Fig. 4b). Tropical OLR anomalies in these two sets of experiments are shown in Supplementary Fig. 5, panels a and b. They are strongly pronounced in the experiment with prescribed perturbed SST values (panel a) and largely absent from the experiment in which the tropical SSTs were fixed to the control values (panel b).

To further constrain the role of tropical surface perturbations alone in propagating the signal of high latitude changes globally, in Supplementary Fig. 4c we show the results of simulation that features only tropical SSTs prescribed to the perturbed values, without accompanying high latitude cooling. In this case the upper level streamfunction response over the North Pacific strongly resembles the one in Supplementary Fig. 4a, suggesting that tropical surface changes are critical for the development of tropical Pacific driven Rossby wave train. The streamfunction response over the Arctic differs between panel a and c, in line with the presence and absence of high latitude forcing in the two cases. While the sea ice loss initiates the changes, the consequent tropical Pacific SST and surface flux changes are crucial for the development of the global response described in our study. In that sense, our results as described in Figs. 2 and 3, show a combined response to Arctic sea-ice and Pacific SST and surface flux changes. Some studies have suggested that atmospheric circulation variability in the Arctic is driven by the internal SST variability in the tropical Pacific36,37. The present work supports the idea that it may not have to be one or another (high latitudes driving tropical changes or tropics driving the high latitudes), but that once the changes are initiated in the tropics or high-latitudes they could potentially be reinforcing each other.

Response to the sea-ice loss induced high-latitude energy imbalance

Besides driving the circulation changes, the tropical convection changes also play an important role in compensating for sea-ice induced energy imbalance. As the sea-ice loss alters the high-latitude energy budget these changes are compensated for outside of the high latitudes. This is illustrated in Fig. 5 showing the zonal annual mean atmospheric net column and top-of-atmosphere (TOA) energy flux changes. In the high latitudes, TOA net heat flux into the space decreases (i.e., increases downward ; Fig. 5a–c) as a result of sea-ice loss induced decrease in the net upward TOA shortwave flux that is only partly compensated by an increase in net TOA upward longwave flux (not shown).

Fig. 5: Zonal annual mean top-of-atmosphere (TOA) and atmospheric net column flux changes.
figure 5

a–c TOA anomalies (positive down [W/m2]); d–f Atmospheric column net flux anomalies [W/m2]. All anomalies have been normalized by the cosine of latitude.

Full size image

Hight latitude atmospheric net column flux is positive (Fig. 5d–f) and strongly follows the anomalous net TOA flux in all simulations, indicating small or negligible impact of net surface flux changes over this area. This is a because sea-ice induced decrease in the net upward shortwave flux at the surface is largely compensated by an increase in the latent and sensible heat fluxes due to the newly open ocean areas (not shown). The dominance of high latitude TOA flux over the surface flux anomalies in response to sea-ice changes has been demonstrated in previous studies43,44. This however may not be the case in experimental setups that impose artificial surface energy fluxes in order to remove sea-ice cover, resulting in a substantial high latitude surface energy budget imbalance.

Positive anomalous high latitude flux into the atmospheric column indicates the high-latitude atmospheric energy surplus. Outside the high latitudes, negative atmospheric net column flux anomalies compensate for this surplus. This compensation however appears to be influenced by the ocean model complexity. In the CCSM4-SO, atmospheric net column flux changes continue to be dominated by the TOA net flux changes outside the high latitudes too, and the compensation for the high latitude energy surplus happens mainly in the tropics, through an increased TOA radiation to space (increased OLR due to tropical Pacific convection reorganization, see Fig. 4). In the IPSL-CM5A2 and EC-Earth3.2 experiments, outside of the high latitudes atmospheric net column flux does no longer follow TOA net changes as closely, indicating that the net surface flux changes also compensate for the high latitude energy surplus. As a result, less energy compensation is needed in the tropics. This in part explains weaker tropical Pacific convection changes (and in turn Rossby waves response) in the EC-Earth3.2 and IPSL-CM5A2 models compared to CCSM4-SO slab ocean simulations: less tropical TOA (and OLR) changes are needed for compensation since some of the energy is compensated for at the surface level too.

Atlantic atmospheric circulation changes and Mediterranean rainfall

Propagation of the Rossby wave train from the Pacific Sector can also be traced by considering the Plumb flux anomalies35,48 shown in Fig. 3d–f. From the tropical Pacific, the wave energy is propagated across the Northern American continent and into the tropical Atlantic. Anomalous Plumb fluxes indicate a strong westward propagation in the North Atlantic, just south of Greenland, in the EC-Earth3.2 and IPSL-CM5A2 models, but not in CCSM4-SO. Negative SLP anomalies northwest of the Mediterranean Basin are more pronounced in the two models coupled to a full dynamic-thermodynamic ocean (Fig. 2d–f) and are centered in the Atlantic (as opposed to being centered in the Mediterranean in CCSM4-SO). The respective locations of SLP anomalies are consistent with the corresponding surface temperature anomalies that are advected further south in the slab model simulations compared to the coupled simulations (Supplementary Fig. 6a–c). Similarly, a patch of high latitude SLP increase, encompassing the Greenland seas and Scandinavia in the EC-Earth3.2 and ISPL-CM5A2 models, is also not present in CCSM4-SO. In accordance with this, drier winter conditions over Scandinavia are seen only in EC-Earth3.2 and IPSL-CM5A2 simulations. It is possible that the atmospheric response over the Atlantic region is additionally shaped by emerging ocean feedbacks that are not simulated in a CCSM4-SO setup, such as decadal changes in the Atlantic Meridional Overturning Circulation (AMOC). For example, weakened AMOC means less heat is transported to the high latitudes partly compensating for the excess of heat induced by the loss of Arctic sea-ice. In this case, the atmospheric response in the fully coupled simulations would be expected to be weaker compared to those obtained with the slab ocean, similar to what is found in this study.

Despite this, due to the weakening of the Azorean high across all the models, a winter response that resembles negative NAO-like conditions is present in all simulations. Over the European continent, the only region where statistically significant (at the 90% level) precipitation changes are consistently found in all three model ensembles is the western Mediterranean, where wetter winter conditions are suggested (Fig. 2g–i). The anomalous Plumb fluxes over this region are westward, suggesting a possible decrease in Rossby wave propagation in response to Arctic-sea ice loss. The dipole temperature signature over Europe and Mediterranean, typically associated with the negative NAO winters49 is, however, not present: mean DJF temperature anomalies are positive over northern Europe and show no change over southern Europe (Supplementary Fig. 6).

Comparison of Atlantic vs. Pacific atmospheric responses suggests that model complexity, i.e., employing the slab ocean instead of a fully coupled dynamic-thermodynamic ocean, does not substantially affect the simulated decadal-scale response to sea-ice changes in the North Pacific, but it makes a difference in the Atlantic sector. This differs from the previous findings comparing the slab ocean and fully coupled model responses in a non-energy conserving setups for isolating the impacts of sea ice loss22.

The recent decade

Over the last forty years, Arctic sea-ice cover has been rapidly disappearing, with statistically significant downward trends now emerging in all months50,51. However, detection of the sea-ice induced teleconnections outside the model world remains challenging. The ongoing climate change is not only influenced by the loss of Arctic sea-ice cover, but is also shaped by multiple other factors (greenhouse gas forcing, natural and anthropogenic aerosol emissions, different cycles of natural variability, etc.). All these different factors interfere with the detection of the climate signals associated exclusively with sea-ice loss. In Fig. 6, we use the ERA5 reanalysis data52 to compare the most dominant atmospheric features over the last decade with the ones from the 25-year period at the start of the satellite observations of Arctic sea-ice cover. We consider the preceding autumn’s sea-ice cover and the following winter’s SLP, upper level geopotential and streamfunction anomalies over the recent winters (2009–2019) relative to their corresponding 25-year means starting in 1979. The seasonal cycle of sea-ice anomalies over the chosen periods are shown in Supplementary Fig. 7.

Fig. 6: Sea-ice concentrations and circulation anomalies from ERA5.
figure 6

Difference in November sea-ice concentrations (a) and DJF sea level pressure (b), DJF 250 hPa geopotential height (c) and DJF 250 hPa streamfunction (d) between the recent decade and the 25-year period since the start of sea-ice satellite observations from ERA5 re-analysis. (November means are calculated for a 10-year period 2009–2018 relative to a 25-year period 1979–2003, DJF means that are calculated for a 10-year period starting in December 2009 and ending in February 2019 relative to a 25-year period starting in December 1979 and ending in February 2004). Crosses highlight the locations of the centers of anomalous geopotential highs from CCSM4 slab ocean (black), EC-Earth3.2 (green) and IPSL-CMPA2 (white) ensemble means, stippled areas indicate the anomalies significant at the 90% confidence level. JRA version of the same figure is provided in Supplementary Fig. 8.

Full size image

The recent decade shows a statistically significant increase in the SLP over the North Pacific and central Siberia and a decrease over the Mediterranean Basin. Aloft, geopotential height increases in the North Pacific and Barents-Kara-Greenland seas regions and decreases over north-east United States and western Europe (Fig. 6c). This pattern is reminiscent of the “North American Winter Dipole” that has increased in frequency and intensity during recent decades53,54. The stream function anomalies (Fig. 5d) show a Rossby wave train propagating from the tropical Pacific. A similar picture emerges when considering the JRA data55 (Supplementary Fig. 8).

While the atmospheric features dominating the last decade (Fig. 6) bare some interesting similarities with the features of the atmospheric response to Arctic sea-ice loss shown in this study (Figs. 2 and 3), their respective magnitudes differ substantially, and are 2 to 3 times larger in the reanalysis data. The “real” world, full-forced with GHG changes, in addition to the Arctic sea-ice loss features multiple other processes and feedbacks that can contribute to Arctic warming. Arctic amplification arises as a consequence of sea ice loss, but also atmospheric and oceanic poleward heat transport changes, cloud and water vapor feedbacks56,57. It also occurs in the absence of sea-ice changes - Arctic amplification is present even in the aquaplanet simulations that have no sea ice when forced with greenhouse gas increase58,59. However, in the “real” world, there are also multiple other competing factors and influences. For example, El Niño-Southern Oscillation (ENSO) is known to have strong remote effect on both Pacific and North Atlantic circulation and could potentially affect the Arctic sea-ice loss induced Rossby wave train. In that sense La Niña conditions could strengthen the sea-ice impacts and El Niño weaken it. In fact, our simulations show that during a strong La Niña event the response in the low ice simulations is very similar to the mean response shown in Fig. 3, but substantially stronger (Supplementary Fig. 9a). However, during a strong El Niño event the response in the North Pacific and North Atlantic is the opposite to the mean response shown in Fig. 3 (Supplementary Fig. 9b) in agreement with the typical teleconnection pattern during warm ENSO events60. Recent studies61,62,63 suggested dominance of a La Niña-like mean response to global warming in those CMIP6 model projections that correctly account for effects linked to the ocean-thermostat and aerosol forcing. This response is suggested to be gradually replaced by the El Niño-like one in the longer term. These effects could over the forthcoming decades first strengthen the La Niña-like response to Arctic sea-ice loss we report here, and then later weaken it.

With further increase in global temperatures, other competing factors may become important too. Antarctic sea-ice loss has also been suggested to result in the opposite patterns to the Arctic one, which could lead to partial cancellation of their mutual effects if the Antarctic sea-ice loss intensifies19. In the upper troposphere, tropical warming is projected to dominate over the polar warming, meaning that high to mid-latitude temperature gradient will increase at upper levels. This too could exert opposite influences on midlatitude circulation patterns compared to the ones arising from the Arctic warming alone64. Finally, we note that Arctic sea-ice loss has been suggested to affect the planetary-scale wave propagation during the weak polar vortex events through a stratospheric pathway27. Further examination is needed to clarify the possible interplay between the stratospheric and the tropical surface processes described in this study.

Given the complexity of processes involved, potential similarities between the observed and simulated patterns discussed here (geopotential ridge in the North Pacific as a part of tropical Pacific driven Rossby wave train) are no short of surprising. This calls for a detection and attribution tailored experiments that would help untangle the causes behind recent persistent geopotential highs in the North Pacific responsible for California’s droughts, including the potential role of Arctic sea-ice loss.

Discussion and conclusions

In this study we have shown that the multi-model decadal response to Arctic sea-ice loss results in development of a North Pacific geopotential ridge and negative NAO-like conditions in the Atlantic Basin during winters. This circulation response is characteristic of predominantly dryer conditions over the American southwest (and the Scandinavian peninsula in fully coupled simulations) and wetter conditions over the Mediterranean Basin. This is opposite to the negative surface pressure and/or upper level geopotential anomalies in the Northern Pacific found in a number of previous studies6,7,8,9,10,11,12,13,14,15,16,17 that considered centennial responses to sea ice loss and applied non-energy conserving experimental setups for isolating its impacts.

Other than the PAMIP initiative that provides a common framework for testing the responses to sea-ice loss in a multi-model setup46, previous studies feature a number of methodological differences making it difficult to discern the exact causes of the discrepancies in simulated responses to sea-ice loss. These study design differences include: different time-scales of the responses considered (i.e., centennial6,7,8,9,10,11,12,13,14,15,16,17 vs. decadal18,19,20), the type of method used for isolating the impacts of sea-ice loss (energy conserving8,19,20 versus non-energy conserving6,7,9,18), consideration of Arctic changes alone6,7,9,18,19,20 as opposed to considering both Arctic and Antarctic changes simultaneously8,65, model complexity (atmosphere-only65 vs. slab ocean19 vs. fully coupled model simulations6,7,18,20,31), consideration of sea-ice changes over only part of the Arctic as opposed to the whole Arctic sector31,47, emphasizing different seasons21, different background ocean states as well as any intermodel differences. Arctic sea-ice response can also be affected by decadal SST variability5,66,67 and by internal variability68. In particular, ENSO is known to have strong remote effect on both Pacific and North Atlantic circulation and could potentially interfere with the Arctic-midlatitude linkages, modulating the SST-induced Rossby wave train (see Supplementary Fig. 9).

PAMIP style experiments represent a very important step in overcoming these issues. However, current PAMIP protocol focuses on the initial arrival of sea ice induced perturbations by using a large ensemble of very short (up to 1-year long) simulations. Physically, the response shown in this study represents the next phase, after the initial signal of sea-ice forcing has already arrived and got established. Thus, by providing a robust multi-model responses obtained using the same type of energy conserving methodologies for isolating the impacts of sea ice loss, this study makes a major step towards understanding the impacts of Arctic sea-ice loss on decadal time-scales.

It is interesting to note that the majority of the studies showing the strengthening of the Aleutian pressure low investigated the responses happening a few centuries after the sea-ice loss commenced, while the studies indicating the geopotential ridge strengthening mostly focus on the response occurring within couple of decades after the sea-ice loss. The substantial difference between the reported short and long term responses indicates that climate adaptation efforts aimed at the impacts of a changing Arctic may need to consider these time-scales separately, since there could be a need to adjust differently to the short term and long term impacts.

The Mediterranean Basin and California, two regions whose rainfall shows a striking connection to the Arctic sea-ice loss in our simulations, also represent some of the top ‘uncertainty hotspots’ in terms of future projections of rainfall changes69,70. Given the rather large inter-model differences in the future projections of the Arctic sea-ice loss, this raises the question as to whether the large spread in projections of future sea-ice loss could be linked to the large spread in the future projections of other climate variables far away from the Arctic, like Californian and Mediterranean rainfall. Understanding these links could provide important insights for the development of new emergent constraints71 that would help decrease these uncertainties. The Mediterranean region has already been named ‘a climate change hot spot’ as a highly populous region that could experience greater aridification than any other area in the world72,73, while California’s precipitation, influenced by changing midlatitude and tropical circulation patterns, has been suggested to change anywhere between −30% and +75% by the end of century (2060–2089 relative to 1960–1989)70. Finding the ways to decrease the uncertainties in the future projections of Mediterranean and California rainfall is a burning issue while one, possibly overlooked solution, may be linked to better constraining the impacts of Arctic sea-ice loss.

Methods

Perturbed sea-ice physics simulations

In order to isolate the atmospheric impacts of sea-ice loss, we compare two sets of experiments that differ only in the amount of sea-ice cover. All other factors - the greenhouse gas concentrations, solar insulation and land use are kept constant between the ensemble members and representative of present-day climate. Greenhouse gases are prescribed values from the year 2000.

The sea-ice cover is altered without imposing artificial energy fluxes, thus ensuring that the observed atmospheric response originates solely from the sea-ice changes and is not additionally affected by imposing energy flux perturbations. The sea-ice loss is achieved through sea-ice physics parameter perturbations. The selected sea-ice physics parameters for each model are listed in Table 1. Parameter perturbations are applied in the Northern Hemisphere only (“low Arctic ice” simulations). The anomalies represent the ensemble mean differences between “low Arctic ice” and control simulations.

Table 1 Sea-ice physics parameter values in control and perturbed (low ice) simulations
Full size table

Statistical significance is estimated using Welch’s t-test in order to account for different sample size and variances between the control and ‘low Arctic ice’ simulations. Prior to t-test calculations, the equivalent sample size is determined at each grid point (accounting for temporal autocorrelation).

Models

Community Climate System Model version 4 (CCSM4) in slab ocean mode – CCSM4-SO

The CCSM474 configuration employed incorporates the Community Atmosphere Model version 4 (CAM4)75, Community Land Model version 476 and Los Alamos Sea Ice Model version 4 (Community Ice CodeE 4—CICE4)77 coupled to a mixed layer (slab) ocean. The atmosphere and land models are run on a 1.9° × 2.5° (latitude × longitude) finite volume grid with 26 atmospheric levels in the vertical dimension. The ice and ocean models use a 1° displaced pole grid. The CCSM4-SO simulations were run for 40 years with the last 20 years considered in the analysis. We consider 2 sets of low ice simulations (using two sets of perturbed sea-ice physics parameter values) and one set of control simulations, with each set consisting of 6 ensemble members. Individual ensemble members were produced through different initializations of atmospheric model components. More details on these simulations are provided in Cvijanovic et al.19.

European community Earth-System Model version 3.2 – EC-Earth3.2

The EC-Earth3.2 model78 incorporates the Integrated Forecast Systems (IFS) model (6)—simulating atmospheric processes, the land surface exchange and hydrology model (H-TESSEL)79, coupled to a full-depth ocean model (Nucleus for European Modelling of the Ocean version 3.6—NEMO3.6)80 and a multilayer dynamic-thermodynamic sea-ice model (Louvain-la-Neuve Sea Ice Model version 3—LIM3)81. These components are coupled using the OASIS coupler system82. The IFS model was run on a finite volume grid with spectral resolution T255L91, that is approximately equal to an 80 km grid size, with 91 atmospheric levels in the vertical dimension (up to 0.1 hPa). The ocean model employs ORCA1L75 grid (with nominal resolution of 1 degree and 75 vertically discretized levels). Simulations were run for 30 years, with the last 20 used in the analysis. The ocean initial conditions originate from years 50 and 90 of the present-day control run (branched from a year 2000 of the historical simulation). In total, there are 10 control and 10 ‘low ice’ ensemble members produced by 5 different initializations of atmospheric model components for each of the two ocean initializations.

Institut Pierre-Simon-Laplace Model CM5A version 2 – IPSL-CM5A2

The IPSL-CM5A283 configuration employed utilizes the atmospheric general circulation model LMDZ5A84, the land surface module ORCHIDEE (ORganizing Carbon and Hydrology In Dynamic EcosystEms)85, the Nucleus for European Modelling of the Ocean version 3.6 (NEMO3.6) ocean model80 and the dynamic and thermodynamic sea-ice model The Louvain-la-Neuve Sea Ice Model version 2 (LIM2)86. The components are coupled using the OASIS coupler system (10). The atmospheric and land models are run on a ~ 1.9° × 3.7° (latitude × longitude) finite volume grid, with 39 atmospheric levels in the vertical dimension. The ocean model is run on a “tripolar” ORCA2.3 grid with 31 vertically discretized levels. The simulations were run for 30 years, with the last 20 used in the analysis. Oceanic initial conditions correspond to year 90 of a present-day control simulation (branched from a 500-yr preindustrial spin-up simulation). Ten initial atmospheric conditions are chosen randomly from a stabilized present-day control run and used to produce 10 control and 10 ‘low ice’ ensemble members. More details on these simulations can be found in Simon et al.20.