Abstract

Previous studies have extensively investigated the impact of Arctic sea ice anomalies on the midlatitude circulation and associated surface climate in winter. However, there is an ongoing scientific debate regarding whether and how sea ice retreat results in the observed cold anomaly over the adjacent continents. We present a robust “cold Siberia” pattern in the winter following sea ice loss over the Barents-Kara seas in late autumn in an advanced atmospheric general circulation model, with a well-resolved stratosphere. Additional targeted experiments reveal that the stratospheric response to sea ice forcing is crucial in the development of cold conditions over Siberia, indicating the dominant role of the stratospheric pathway compared with the direct response within the troposphere. In particular, the downward influence of the stratospheric circulation anomaly significantly intensifies the ridge near the Ural Mountains and the trough over East Asia. The persistently intensified ridge and trough favor more frequent cold air outbreaks and colder winters over Siberia. This finding has important implications for improving seasonal climate prediction of midlatitude cold events. The results also suggest that the model performance in representing the stratosphere-troposphere coupling could be an important source of the discrepancy between recent studies.

INTRODUCTION

Over recent decades, adverse cold winters have occurred more frequently in Eurasia (1), concurrent with a pronounced warming over the Arctic, coined as the “warm Arctic cold Siberia (or continents)” (WACS) pattern (2). Some studies have suggested that these observed cold anomalies are attributable to rapid sea ice decline and Arctic warming, especially over the Barents-Kara seas (BKS) (3–6). Others have argued that the cold winter events could be a consequence of natural variability of the atmosphere (7, 8). The inconsistent and conflicting results of these studies indicate that the mechanism behind the colder Siberia is still unclear.

When examining the atmospheric circulation, studies have found a weakening or a shifting of the stratospheric polar vortex due to Arctic sea ice loss (9–12) and a subsequent downward anomaly from the stratosphere to the troposphere (13–15). Previous model studies (14, 16) have investigated the role of stratosphere-troposphere coupling induced by BKS sea ice retreat (the stratospheric pathway) in contributing to the extratropical tropospheric circulation responses in addition to the direct response of the tropospheric circulation itself (the tropospheric pathway). In observational analysis, a statistically significant relationship between Siberian cold winters and BKS sea ice loss is found (Fig. 1A). However, when the stratospheric polar vortex anomaly is removed through a simple linear regression analysis, the Siberian cold anomaly reduces to about half, suggesting a possible critical role of the stratospheric pathway (Fig. 1B). In the remainder of this study, we complement the above statistical analysis by using targeted comprehensive atmospheric general circulation model (AGCM) experiments to assess the explicit contribution of the stratospheric pathway to the Siberian cold anomaly and extreme cold events. Although previous studies (14, 16) have examined the role of the stratospheric pathway in the tropospheric circulation response to sea ice changes, our analysis expands on these previous efforts in a number of ways. First, previous analyses of the stratospheric pathway have been largely confined to the large-scale circulation response, such as the Arctic Oscillation, and have not emphasized its contribution to high-impact regional weather extremes, especially in the context of the recently hotly debated linkage between sea ice loss and the midlatitudes in the climate community [see (17) and references therein]. Our analysis provides robust evidence for the stratospheric pathway as an important component of the link between sea ice loss over the BKS and the Siberian surface climate and weather extremes. Second, our aim is for a quantitative assessment of the relative contributions of the stratospheric and tropospheric pathways and the linear additivity of the two. Exploring the connection between these two pathways in the formation of Siberian cold winters represents a step forward in our understanding of the physical mechanisms. Furthermore, the results suggest that the formation of a cold Siberia can be attributed to the regional sea ice loss over the BKS rather than the whole Arctic.

Fig. 1Observed linkage between cold Siberia and BKS sea ice loss with and without the stratospheric circulation activity.

(A) Regression of SAT (in Kelvin) in December-January-February (DJF) on the normalized early winter [November-December (ND)] BKS sea ice concentration (SIC). (B) Same as (A) after first regressing out the component of SAT variability related to DJF stratospheric vortex anomalies. Both the long-term trend and El Niño–Southern Oscillation (ENSO)–related component in all fields are removed before the regression to emphasize the interannual variability. See Materials and Methods for the data and calculation details. The sign of SIC is reversed here to emphasize the variation related to sea ice loss. Stippling denotes regions that are statistically significant at the 95% confidence level using a two-sided t test for the regression. Red sector in (A) highlights the Siberian cold anomaly region.

Previous modeling studies have suggested that low-top models or models with an insufficiently resolved stratosphere exhibit weaker stratosphere-troposphere coupling than that in observations or stratosphere-resolving models (18–20). An unrealistic simulation of stratospheric processes may limit a model’s ability to represent the complete response to sea ice forcing (13). Therefore, our approach is to use an AGCM with a good representation of stratospheric dynamics to explore the atmospheric circulation responses to BKS sea ice decline and its impact on the WACS pattern.

RESULTS

Colder Siberia in response to BKS sea ice forcing

The AGCM SC-WACCM4 (Specified Chemistry Whole Atmosphere Community Climate Model version 4), which includes a good representation of the stratosphere, is used to conduct three perturbation experiments (BKS_FL, BKS_TP, and BKS_SP; see Materials and Methods and table S1). These experiments aim to examine the linkage between anomalous surface forcing due to BKS sea ice anomalies and Siberian surface climate and the role of the stratospheric and tropospheric pathways in this linkage. The BKS SIC anomalies are obtained from the differences between projected future sea ice and sea ice in historical runs following the previous work of Sun et al. (13) (see Materials and Methods and Fig. 2). The BKS sea ice decline reaches its maximum in November (about 22%; see Fig. 2A), during which the ice edge retreats from the Barents Sea to the Kara Sea (Fig. 2B). The magnitude of BKS sea ice forcing is comparable to that in previous modeling studies, which used the observed sea ice anomaly during the past three decades (see the Supplementary Materials for more discussion on the forcing) (6).

Fig. 2SIC anomalies (%) in the BKS_FL and BKS_TP runs compared to those in the CTRL run.

(A) Monthly SIC changes averaged over the BKS region [highlighted by the black sector in (B)]. (B) Spatial distribution of BKS SIC changes in November. The orange and red lines in (B) denote the ice edge (15% contour of SIC) in the CTRL and BKS_FL (BKS_TP) experiments, respectively.

In response to the prescribed SIC anomalies shown in Fig. 2, the surface air temperature (SAT) anomaly in the BKS_FL (full response) run is characterized by a local surface warm anomaly and a cold anomaly in northern Eurasia (more statistically significant over Siberia than over Europe) (Fig. 3A), as detected in observations (compare with Fig. 1A). The Siberian cold anomaly persists throughout the whole winter, achieves its maximum in December and January, and disappears at the end of February (fig. S1).

(A) SAT (color shadings; in Kelvin) and surface wind (vectors; in m s−1) anomalies in the BKS_FL experiments during DJF. (B) Geopotential height anomalies (500 hPa) [color shadings; in geopotential meters (gpm)] in the BKS_FL experiments during DJF. (C) Same as Fig. 1A but for regression of geopotential height on the BKS SIC index. The contours in (B) and (C) show the climatological DJF 500-hPa geopotential height in the CTRL and the observations, respectively. The contour interval in (B) and (C) is 120 gpm. Stippling denotes regions that are significantly different from zero at the 95% confidence level using a two-sided t test for the color-shading variables. The surface wind anomaly that is significant at the 95% confidence level is shown. Red sector in (A) highlights the Siberian cold anomaly region similar to that in Fig. 1A.

Associated with the colder Siberia are anomalous surface winds in the sense of a local-scale cyclonic anomaly over the Arctic Ocean, a cyclonic anomaly over the midlatitude North Atlantic Ocean, and a continental-scale anticyclonic anomaly over northern Eurasia, indicating an enhanced Siberian High (vectors in Fig. 3A). Correspondingly, there is an increase in geopotential height over the Arctic and a decrease over midlatitude Asia in the mid-troposphere (Fig. 3B), which is in line with the observed geopotential height anomaly (Fig. 3C). The Siberia surface cold response predicted by our modeling experiments and the anomalous circulation pattern (Fig. 3, A and B) closely resemble that seen in the observations (Figs. 1A and 3C), enhancing our confidence that the observed anomalies are actually associated with BKS sea ice loss. The Arctic positive geopotential height anomaly extends to northern Europe to favor a ridge anomaly near the Ural Mountains, while the East Asia trough becomes much deeper compared to the CTRL state, favoring southward cold advection from the Arctic. The anomalous ridge and trough pattern and associated near-surface winds would enhance the advection of cold air from the Arctic over the Eurasian midlatitudes contributing to the anomalous Siberian cold anomaly.

Distinct role of the stratospheric pathway

The experiments BKS_TP and BKS_SP are designed to distinguish the roles that the stratospheric and tropospheric pathways play in the formation of the colder Siberia. In BKS_TP, the same sea ice anomalies as in BKS_FL are imposed (Fig. 2), but the zonal mean state of the stratospheric circulation is restored to that in the CTRL, thus numerically largely suppressing the stratospheric response to the forcing and the stratosphere-troposphere interaction (see Materials and Methods for details). When shutting down the stratospheric pathway and leaving the tropospheric pathway acting in isolation in this manner, the Siberian cold anomaly becomes much weaker than that in BKS_FL, although the same magnitude of the local BKS warming still exists (Fig. 4A). The smaller magnitude of the Siberian cold anomaly is due to much weaker ridge and trough over Eurasian midlatitudes and corresponding weaker surface northerly wind and thus much weaker cold air advection from the Arctic region (Fig. 4, A and B). Therefore, the full response in the troposphere is significantly reduced in the absence of the coupling with the zonal mean stratospheric response, and the tropospheric pathway alone cannot produce the Siberian cold anomaly as seen in the observations and BKS_FL run.

(A, C, and E) Same as Fig. 3A but for the SAT and surface wind anomalies in the BKS_TP (A) and BKS_SP (C) and the sum of the BKS_TP and BKS_SP (E) runs. (B, D, and F) Same as Fig. 3B but for the geopotential height anomaly in the BKS_TP (B) and BKS_SP (D) and the sum of the BKS_TP and BKS_SP (F) runs. Stippling denotes regions that are significantly different from zero at the 95% confidence level using a two-sided t test for the color-shading variables. Red sector in (E) highlights the Siberian cold anomaly region similar to that in Fig. 1A.

To explore the explicit contribution of the stratospheric pathway to the colder Siberia, we conducted a stratospheric pathway experiment, namely, the BKS_SP run, in which the zonal mean stratosphere is nudged to a perturbed state, derived from the BKS_FL run as a result of BKS sea ice loss, but the CTRL sea ice is prescribed. A surface cold anomaly and anomalous anticyclonic winds are found over northern Eurasia in BKS_SP compared with the CTRL (Fig. 4C). Northern Europe actually ends up much colder in the BKS_SP run. The 500-hPa geopotential height anomaly pattern over Asia also resembles that in the BKS_FL run and is statistically significant, although slightly weaker than the full response (Fig. 4D). The deepening of the East Asian trough favors the northerly wind over northern Asia advecting cold air masses from the Arctic to Siberia.

Comparison between BKS_TP and BKS_SP provides further insights into the relative contributions from the stratospheric and tropospheric pathways. Without the stratospheric response, the response mainly reflects local SAT changes as a direct response to the retreated sea ice (BKS_TP; Fig. 4A). A full representation of the stratospheric response to sea ice loss enables the model to simulate the Siberian cold anomaly (BKS_SP; Fig. 4C). When summing the results from both BKS_TP and BKS_SP, we find that the SAT, surface winds, and 500-hPa geopotential height patterns are comparable with those simulated by BKS_FL in both the spatial structures and magnitudes, suggesting the linear additivity of the stratospheric and tropospheric pathways (compare Fig. 3, A and B, with Fig. 4, E and F). The results, thus, reinforce the important role of the stratospheric pathway in driving the formation of the teleconnection pattern between the BKS sea ice loss and Siberian midlatitudes.

The underlying mechanisms lie in the response of the midlatitude circulation (Fig. 4B versus Fig. 4D and Fig. 4A versus Fig. 4C), which is attributed to the downward effect of the stratospheric zonal mean change on the tropospheric circulation. In the BKS_FL run, first, the stratospheric polar vortex is weakened by the BKS sea ice retreat mostly via enhanced planetary-scale upward wave propagation into the stratosphere, as can be seen in the circulation and Eliassen-Palm (EP) flux response in fig. S2 (A and B). The increased wave flux due to the sea ice loss is achieved by the linear constructive interference of planetary-scale waves (fig. S2, C and D), especially via the contribution of the zonal wave 1 component (fig. S2C). The wave anomalies are largely in phase with the climatological waves so that more wave energy is transported into the stratosphere (fig. S2, B to D). This stratospheric circulation anomaly persists in the stratosphere throughout the winter and is accompanied by downward migration back to the troposphere and the surface (fig. S2A). Unsurprisingly, this downward migration disappears when the stratospheric zonal mean flow is nudged back to the CTRL and the stratosphere-troposphere dynamical interaction is decoupled in the BKS_TP run (fig. S3, A and B), but largely recurs in the BKS_SP run (fig. S3C). The downward influence of the stratospheric processes significantly amplifies the weakening of the tropospheric zonal mean flow via anomalous poleward wave propagation (fig. S3D). This downward effect is mainly manifested over the North Atlantic–Eurasia sector as an intensified ridge and trough pattern, as shown in Fig. 4D. Therefore, the results support the argument that there is a dynamical linkage between sea ice loss and a cold Siberia and further suggest that the coupling between the stratosphere and the troposphere is key to producing this connection. Ultimately, the changes in the stratospheric zonal mean circulation owe their existence to upward-propagating planetary wave anomalies excited by the sea ice loss below, through the troposphere, and these results highlight the importance of capturing the stratospheric zonal mean response and its subsequent downward feedback for the full tropospheric circulation response and the cooling over Siberia.

Responses in cold air outbreaks over Siberia

The above analysis has demonstrated that BKS sea ice loss results in a colder winter climatology over Siberia and that the altered stratospheric zonal mean state is a key component of this cold anomaly. Given a weaker stratospheric polar vortex favoring cold events over Eurasia (21), what are the implications of BKS sea ice loss and the role of the stratosphere-troposphere coupling for cold extremes on shorter, synoptic time scales? Figure 5 (A and B) shows the probability density function (PDF) of daily SAT over Siberia in all experiments and observations, respectively. There is a shift of the PDF toward colder SAT with BKS sea ice loss in BKS_FL and BKS_SP and in observations. Specifically, in observational analysis, the probability of 1 standard deviation (σ) below the DJF mean SAT of the normal BKS SIC years increases from the normal years value of 17.8 to 27.4% in the lower BKS SIC years. Similarly, in the modeling experiments, the probability of 1σ below the CTRL DJF mean increases from the CTRL value of 15.7 to 22.9% in BKS_FL and 21.4% in BKS_SP, while BKS_TP is similar to CTRL (16.8%). To address whether the change in the cold tail of the distribution in Siberian daily SAT results from the climatological shift of the distribution or a change in the shape of the distribution, we used the two-sample Kolmogorov-Smirnov (KS) test. The results (P > 0.05) suggest that the shape of the PDF in the perturbation runs is about the same as that in the CTRL run. Thus, changes in the PDF of Siberian daily SAT are largely due to the shift of the climatological mean in BKS_FL and BKS_SP in response to BKS sea ice loss. The KS test result for the observations is the same. This suggests that BKS sea ice loss forces a seasonal mean climatological cold anomaly over the Siberian surface through stratosphere-troposphere coupling, superimposed on which the occurrence of daily-scale cold events follows the same distribution as in the CTRL climatology but with a shift toward colder SAT.

(A) PDF of the daily SAT over Siberia in the CTRL, BKS_FL, BKS_TP, and BKS_SP experiments. (B) Same as (A) but for the observations in which the red line is that in the 10 lowest BKS SIC years (referred to as LowSIC) and the gray line is that in all other years (regarded as normal years and referred to as NorSIC) (see Materials and Methods for details). (C) Frequency of CAOs per winter in each experiment. (D) Same as (C) but for the duration of the CAO events (days). The vertical line in (A) and (B) shows the 1σ of daily Siberian SAT below the DJF mean in the CTRL and NorSIC years. In (B), the DJF SAT component related to long-time trend and ENSO (calculated as linear regression on Niño 3.4 index) are first removed in the daily time series of Siberian SAT. The bars in (C) and (D) show the 2.5th to 97.5th percentile interval of the mean values based on 10,000 bootstrap samples. The asterisks in (C) and (D) denote that the difference between the perturbation runs and the CTRL are statistically significant at 90% (*) and 95% (**) confidence levels based on 10,000 times bootstrap samples.

To further address the change in extreme weather events over Siberia in winter, we investigated changes in the characteristics of cold air outbreaks (CAOs) (see Materials and Methods for the CAO definition). A Siberian CAO event is characterized by a consecutive decrease in local SAT and a strengthening of the Siberian High over a period of about 6 days (fig. S4, A and B). The CAO life cycle is such that when the Siberian High reaches its maximum, Siberian SAT exceeds the CAO criterion and obtains its minimum value about 2 days later (fig. S4, A and B). Here, the frequency and duration of CAO events are presented (Fig. 5, C and D). Compared to the CTRL run (about 2 CAO events per winter; gray bar), there is a significant increase in CAO events (about 2.5 per winter; red bar) in the BKS_FL run as a result of BKS sea ice loss. When the stratospheric pathway is numerically shut down in BKS_TP, the frequency of CAO events remains comparable to that in CTRL (about 2 per winter; orange bar). In contrast, there are more frequent CAO events (about 2.5 per winter; blue bar) in the BKS_SP run, similar to that in the BKS_FL run. As discussed above, this change in occurrence of CAOs, as defined by the CTRL threshold, results from a shift of the whole distribution to colder temperatures when BKS sea ice is reduced and the stratospheric pathway is present. The key driver of CAOs is the presence of a blocking high anticyclone, and observational analysis has indicated that the increase in frequency of winter cold events over Siberia is accompanied by intensified anticyclonic activities (22). The BKS sea ice loss and stratospheric pathway favor anticyclonic conditions over Siberia through the response of the ridge-trough pattern. As for the duration of CAO events, an extension can be found in BKS_FL, which can be attributed to the combined contribution by the stratospheric and the tropospheric pathways (Fig. 5D). In addition, the CAO severity (that is, the minimum SAT) is roughly the same in each experiment (fig. S4A). Therefore, the increase in the frequency of CAO events is primarily due to a stratospheric pathway–dominated shift in the climatological mean state of Siberian SAT such that the CAO threshold as defined from the CTRL climate is more frequently crossed.

DISCUSSION

Through specifically targeted experiments using the advanced, stratosphere-resolving AGCM SC-WACCM4, this study reveals the predominant role of the stratospheric pathway in linking the BKS sea ice loss to cold Siberia. The results demonstrate that inclusion of the stratospheric pathway produces a coherent teleconnection pattern ranging from the Arctic Ocean to the Eurasian midlatitudes via stratosphere-troposphere dynamical coupling, characterized by an enhanced Ural Ridge and eastern Asian trough, in agreement with the response to BKS sea ice loss inferred from observations. Correspondingly, the Siberian High intensifies, leading to southward cold air advection and cold anomalies over Siberian midlatitudes, also favoring an increase in the frequency of CAOs.

Our modeling experiments have exhibited a robust “cold Siberia” pattern as forced by the BKS sea ice loss, but the presence of the stratospheric pathway is the key for this. As discussed in the beginning of this study, inconsistencies have been continually reported among various modeling studies. A detailed discussion of the discrepancies in the atmospheric response to Arctic sea ice loss across climate models can be found in (17). Here, through our study, we suggest that those discrepancies could, in part, be attributed to the following differences across studies. First, the representation of the stratosphere in the models used in previous studies may be an important source of uncertainty. We note that the modeling studies that cannot reproduce the Eurasian cold anomaly used low-top or coarse–vertical resolution models (7, 8). Most of the studies that simulated the Eurasian cold anomaly, including the present one, used high-top or high–vertical resolution models (6, 14). In particular, an additional experiment with Community Atmosphere Model version 4 (CAM4) (a low-top version of SC-WACCM4 with a similar physics package), forced with identical BKS sea ice decline, fails to capture the Siberian cold anomaly (fig. S5A). This may stem from more than simply model vertical resolution, as details of how the modeled stratosphere is tuned, in the presence of parameterized gravity waves, to achieve a realistic climate state may matter (23). Therefore, the model’s ability to capture the Eurasian cold anomaly could be significantly affected by its performance in simulating stratospheric processes.

Second, the forcing used in previous studies is diverse. Sea ice perturbations in different geographical areas of the Arctic Ocean may result in different midlatitude responses (24, 25). The literature suggests that, for AGCM simulations, BKS sea ice forcing, which exerts the largest impact on the midlatitude circulation (25), causes a Siberian cold anomaly (as shown in Fig. 3A), while pan-Arctic sea ice loss leads to a weaker cold anomaly (compare Fig. 3A and fig. S5B; the pan-Arctic SIC forcing can be seen in fig. S5C) or insignificant response over Siberia (7, 26). The impact of sea ice loss over other regions, such as the Sea of Okhotsk and the Labrador Sea, could possibly nonlinearly interfere with the impact of BKS sea ice decline and result in different midlatitude responses (26, 27). For coupled ocean-atmosphere GCM experiments with imposed sea ice loss, studies have found an important role of the oceanic processes in the overall response, but the induced Siberian cooling is about the same (28) regardless of whether the model is coupled to a dynamic ocean or not. However, in future warming scenarios, the Siberian cooling due to Arctic sea ice loss is likely to be overwhelmed by the increase of radiative forcing, as seen in Coupled Model Intercomparison Project phase 5 (CMIP5) future projections (6).

This study investigates the cold Siberia pattern and demonstrates that, with a stratosphere-resolving atmospheric model and imposed sea ice loss over the BKS, the downstream region including Siberia becomes colder with an increased occurrence of CAO events. Our work demonstrates the critical role of stratosphere-troposphere coupling in creating the colder Siberia pattern, and the results point to improving the representation of the stratosphere and its response to sea ice loss as an important step toward more reliable projections of regional climate and extreme weather changes in the context of future polar climate change. In addition, because the largest sea ice loss usually occurs in late autumn–early winter and its impact persists throughout the whole winter due to the long time scale of stratospheric processes, this study has important implications for improving the seasonal climate prediction of midlatitude cold events. Our results also aid in reconciling the discrepancies in previous studies and advance the scientific understanding behind the ongoing debate on the linkage between the Arctic and midlatitude climate and extreme weather.

MATERIALS AND METHODS

Experimental design

Models. The stratosphere-resolving (that is, “high-top”) AGCM used in this study is the SC-WACCM4, a component of the National Center for Atmospheric Research (NCAR) Community Earth System Model version 1.2 (CESM1), with specified chemistry rather than the full interactive chemistry version WACCM4 (29). The removal of interactive chemistry does not change the climatology and variability of the atmospheric circulation below 55 km but much reduces the computational cost, which makes SC-WACCM4 suitable for studies of stratosphere-troposphere dynamical coupling (29). Like WACCM4, SC-WACCM4 has 66 vertical levels and a horizontal resolution of 1.9° latitude by 2.5° longitude, and the model lid is at 5.1 × 10−6 hPa (approximately 140 km). For comparison to a “low-top” model, the CAM4, which shares the same physical parameterizations as WACCM4 (SC-WACCM4), was used. CAM4 has 26 vertical levels, and the model lid is at about 3.5 hPa (30, 31).

Experiments. The modeling experiments in this study are summarized in table S1. Sea surface temperature (SST) and SIC were prescribed to define the model’s surface boundary condition. In the control run (hereafter referred to as CTRL), the repeating climatological seasonal cycle of SST and SIC, from the CESM1-WACCM4 historical outputs (average of seven ensembles) from the CMIP5 and averaged during 1980–1999, were prescribed. In the first perturbation run, referred to as BKS_FL (full response), the settings are identical to those of the CTRL, but the SIC in BKS is replaced by that in the CMIP5 Representative Concentration Pathway (RCP) 8.5 outputs (average of four ensembles) averaged during 2080–2099 (see Fig. 2). The SST in BKS was also replaced with RCP8.5 SST in the open water areas that used to be covered by sea ice in CTRL run, similar to previous studies (13). The second perturbation experiment, BKS_TP, has the same SST and SIC forcings as those in BKS_FL run, but a nudging method was used to largely remove the contribution of the stratospheric pathway and to isolate the tropospheric adjustment by itself (tropospheric pathway). In this nudged experiment, the zonal means of temperature, zonal wind, meridional wind, and specific humidity above 90 hPa were nudged toward those in the CTRL with a time scale of 6 hours. The fields were fully nudged above 54 hPa (the nudging coefficient is 1), and the nudging strength was linearly decreased to zero between 54 and 90 hPa, with no nudging applied below 90 hPa. The tapering region in the vertical, that is, the lower stratosphere (54 to 90 hPa), was chosen to relax the strength of nudging over some vertical scale to prevent issues that would arise at the sharp boundary between no nudging and nudging (32). The nudging was performed at every time step of the model integration, but the target states from the CTRL run were read in every 6 hours, and the model fields were nudged toward the linear interpolation between consecutive target states, which, in this case, were the full time-evolving zonal means of the CTRL simulation. On the basis of Hitchcock and Haynes (32), both the vertical levels of the nudging and the tapering that we used here and the nudging time scale of 6 hours are reasonable choices that achieve what we aim for. In the BKS_TP run, we, therefore, numerically suppressed stratosphere-troposphere dynamical coupling that occurred via a change in the zonal mean stratospheric circulation, allowing the tropospheric pathway to be viewed in isolation. To explicitly quantify the contribution of the stratospheric pathway, another nudging experiment, referred to as BKS_SP (stratospheric pathway), was conducted. The BKS_SP run was the same as the CTRL run, except that the zonal mean state in the stratosphere was nudged toward that in the BKS_FL run with the same nudging coefficients as in the BKS_TP run. The BKS_SP run provides an opportunity to examine the impacts of the perturbed stratospheric circulation by the BKS sea ice anomaly on tropospheric circulation and surface climate when direct surface sea ice forcing is absent.

In addition, to examine the impact of the nudging method, we conducted an additional experiment in which all the settings were the same as the CTRL, except that nudging was applied in the stratosphere and above to nudge the zonal mean state to that of the CTRL itself, referred to as CTRL_NUDG. The results of CTRL_NUDG (zonal wind in December for example) showed that nudging neither changed the zonal mean states nor markedly altered the stationary planetary-scale wave structure (fig. S6). A theoretical discussion of the nudging impact can be found in the study of Hitchcock and Haynes (32).

We also conducted the ARC run, in which the pan-Arctic SIC and SST changes were prescribed. For comparison, two additional experiments with the low-top model version CAM4, that is, CTRL_CAM4 and BKS_CAM4, were performed. These two experiments were the same as CTRL and BKS_FL, except that CAM4 was used.

The anomaly (or response) in the perturbation run was defined as the difference between each perturbation run and the CTRL run. Because the most significant midlatitude response to early winter sea ice changes occurred during the following winter (25, 33), we focused on the responses in DJF. All of the experiments were integrated for 60 model years, with the first 10 years discarded as spin-up. There are thus 49 whole winters in the last 50 model years analyzed in each run in this study. We checked the robustness of the results with two different statistical methods (see the “Statistical analysis” section). The results showed that the conclusions in this study were robust. In addition, we tested whether the same response was obtained when dividing our 49 winters data set in two halves, and it was.

Cold air outbreak. A CAO event is defined as two or more consecutive days during which area-averaged Siberian daily SAT (60°E to 140°E and 45°N to 65°N; highlighted by the red box in Fig. 3) is at least 1 standard deviation (σ) below the DJF mean SAT in the CTRL run, where σ is the average of the daily standard deviation of Siberian daily SAT in 49 winters in the CTRL run (34). The Siberian High index is defined as the weighted area-averaged sea-level pressure (SLP) over 65°E to 100°E and 60°N to 75°N, the region where the maximum in the leading empirical orthogonal function mode of winter SLP in the CTRL run is located (fig. S4C).

Observation data

ERA-Interim. The European Centre for Medium-Range Weather Forecasts (ECMWF) reanalysis ERA-Interim was used in Fig. 1 on a resolution of 1.5° longitude × 1.5° latitude grid (35).

Optimum Interpolation Sea Surface Temperature. The SST data set used to calculate the Niño 3.4 index was from the National Oceanic and Atmospheric Administration (NOAA) Optimum Interpolation (OI) SST version 2 on a resolution of 1° longitude × 1° latitude (37). Monthly variables during 1982–2016 were used in Fig. 1. We defined the early winter BKS SIC index using the area-averaged SIC in the region of 15°E to 100°E and 70°N to 82°N in ND. Here, the sign of SIC was reversed to emphasize the changes associated with sea ice loss. For all the fields in Fig. 1, both the long-term trend and ENSO-related component (calculated as linear regression on Niño 3.4 index) were first removed. The contemporaneous (DJF) averaged stratospheric contribution was removed via linear regression before calculating the regression onto BKS SIC in Fig. 1B. We simply used the polar averaged (65°N to 90°N) zonal mean zonal wind at 50 hPa to define the stratospheric vortex variability. The result in Fig. 1B was not sensitive to the definition of stratospheric vortex variability (not shown). The indices used in Fig. 1 were normalized before the regression. For the calculation of the PDF of daily SAT over Siberia in the observations (Fig. 5B), the DJF mean SAT components related to long-term trend and ENSO (calculated as linear regression on Niño 3.4 index) were first removed from the daily time series in each winter. Then, we chose the daily SAT in the 10 lowest BKS ND SIC years as the lower than normal SIC samples, similar to (6), referred to as LowSIC. The daily SAT in all other years was regarded as the normal BKS SIC years, referred to as NorSIC.

Statistical analysis

Statistical significance was assessed by a two-sided Student’s t test. A bootstrap resampling with replacement was used to calculate the confidence interval of the CAO frequency and duration in each experiment. The confidence interval of the mean value was calculated as the 2.5th to 97.5th percentile range of 10,000 bootstrap samples (Fig. 5, C and D). The bootstrap method was also applied to test the statistical significance of the differences between the perturbation runs and the CTRL run in Fig. 5 (C and D). We also assessed the significance of the results in Figs. 3 and 4 using the bootstrap method and found similar results (not shown). The two-sample KS test was used to compare the PDF distribution shape of Siberian daily SAT in each perturbation run to that in the CTRL run. The mean was removed before conducting the KS test to compare the shape of the distribution only. The null hypothesis, that the distributions of two samples were the same, was rejected when P ≤ 0.05.

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.

Acknowledgments: We thank M. Ting, J. Lu, P. Kushner, and the three reviewers for their helpful suggestions. We also thank L. Sun for the suggestions for the preparation of the sea ice forcing. We would like to acknowledge high-performance computing support from Cheyenne (doi:10.5065/D6RX99HX) provided by NCAR’s Computational and Information Systems Laboratory, sponsored by the NSF. Funding: P.Z., Y.W., and B.D. were funded by the NSF (AGS-1406962). I.R.S. and P.C. were supported by NSF funding to the NCAR. K.L.S. was funded by the NSF Office of Polar Programs, Arctic Research Opportunities, PLR-1603350. X.Z. was supported by the NSF (ARC-1023592). Author contributions: P.Z. and Y.W. conceived the article and experiments. P.Z. performed experiments and analyzed data. I.R.S. helped P.Z. to verify the nudging method. Y.W., I.R.S., and K.L.S. provided feedback on analysis. P.Z. and Y.W. wrote the manuscript. I.R.S., K.L.S., and X.Z. contributed to constructive revision. I.R.S. and P.C. provided guidance on the nudging experiments. B.D. contributed to a part of data preparation. Competing interests: The authors declare that they have no competing interests. Data and materials availability: The CESM-WACCM4 sea ice and SST outputs were accessed at NCAR CMIP Analysis Platform (doi:10.5065/D60R9MSP). All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials. Additional data related to this paper may be requested from the authors.