Atmospheric circulation regimes modulating Eurasian winter decadal cooling

Abstract

The Eurasian winter cooling observed in the early 21st century, which contrasted with ongoing global warming, was initially attributed to Arctic warming until atmospheric internal variability became the prevailing explanation when the cooling ceased despite unabated Arctic warming. Yet the atmospheric variability responsible for it remains elusive. Here, three comparable cooling episodes, including the early 21st-century one, are identified in reanalysis since 1901. By synthesizing them with cooling episodes in large-ensemble simulations, we demonstrate two intrinsic atmospheric regimes as modulators: the Scandinavian pattern (SCAND) and North Atlantic Oscillation (NAO). The substantial positive SCAND trend strongly favors Eurasian cooling, with the NAO trend modifying cooling magnitude. Our logistic regression model indicates that the likelihood of at least one analogous cooling episode before 2050 is approximately 30% due to internal variability, with only half potentially overriding the anthropogenic warming. These findings provide important insights for energy security strategies for midlatitude Eurasian regions.

Introduction

The intensification of extreme winter cold spells across midlatitude Eurasia during the early 21st century occurred paradoxically against a backdrop of unequivocal global warming^1,2,3, imposing significant socioeconomic burdens through cascading disruptions to transportation, energy infrastructure, and agriculture^4,5,6. These severe winters were anchored by a pronounced interdecadal cooling trend over Eurasia that emerged in the late 1980s^2,7,8,9,10 and culminated during 1998–2012^11,12,13,14. During 1998–2012, winter temperatures over a large area of central Eurasia decreased by over 4 °C per decade¹⁴—a trend magnitude comparable to the concurrent Arctic warming.

The early 21st-century Eurasian cooling was initially attributed to accelerated Arctic warming and associated sea-ice loss^{10,15,16,17,18,19,20,21,22}, as the Arctic warmed 2–4 times faster than the global mean (also known as Arctic amplification^{2,23,24,25,26}) during the same period. Proposed mechanisms include both tropospheric pathways^{27,28,29,30,31,32,33} (e.g., decelerated westerly jet, enhanced Ural blockings) and stratospheric connections^34,35,36,37 (e.g., weakened polar vortex). However, the paradigm was challenged when Eurasian cooling ceased post-2012 despite continuous Arctic warming^{38,39,40,41,42,43,44,45}, leading to alternative interpretations emphasizing atmospheric internal variability^{12,46,47,48,49,50,51,52}. Moreover, in large-ensemble atmospheric general circulation model simulations with prescribed Arctic sea-ice loss, the ensemble mean is largely masked by internal variability among different ensemble members, leading to little/no robust Eurasian cooling response to declined Arctic sea ice^48,53,54. Using three sets of large-ensemble simulations forced by seasonal sea ice conditions from preindustrial, present-day, and future periods, a recent estimate demonstrates that atmospheric internal variability explains about 70% of the observed Eurasian cooling⁵⁵.

However, the specific internal atmospheric circulation pattern most conducive to Eurasian decadal cooling remains elusive^56,57. Three key limitations have obscured the cause of decadal Eurasian cooling in the literature: (1) the conflation of Eurasian cooling across disparate timescales, from synoptic cold spells to decadal trends^{13,51,58,59,60,61}, (2) the preconceived notion of warming Arctic-cooling Eurasia (WACE) mode^49,61,62,63, and (3) the incomplete separation of influences of external forcing and internal variability using linear detrending method^64,65. In this study, we specifically focus on decadal-scale Eurasian cooling and, for the first time, identify its periodicity driven by internal variability throughout reanalysis and the Community Earth System Model (CESM) 2 Large Ensemble (LENS2)⁶⁶. The LENS2 multi-member ensemble mean (MMEM) is used to isolate internal variability in reanalysis and individual simulations (“Methods”). Our results indicate that Eurasian cooling episodes are not always accompanied by pronounced Arctic surface warming, contrary to the prevailing WACE concept; nevertheless, these cooling episodes share common characteristics regarding atmospheric circulations, which provides an opportunity to pinpoint the primary atmospheric internal variability responsible for cooling episodes and thereby improve predictions of potential Eurasian cooling.

Results

Periodical episodes of decadal Eurasian winter cooling

The wintertime (January–February mean, JF) Eurasian surface air temperature (SAT) time series (SAT_Eurasia; black solid line, Fig. 1a) is decomposed into externally forced (EX; black dashed line) and internal variability (IV; gray bars) components based on a regression method and LENS2 MMEM (“Methods”). The JF mean is adopted because these months exhibit both stronger cooling trends and higher intermonth correlation, unlike December which shows weaker cooling and low correlations with JF months (Supplementary Table 1). To isolate the atmospheric circulation regimes most conducive to Eurasian cooling, we first focus exclusively on the IV component. We applied the 15-year rolling trend analysis (“Methods”) to IV SAT_Eurasia, revealing the three most significant cooling episodes (bars in upper panel of Fig. 1b): P1 (1917–1931), P2 (1958–1972), and P3 (1998–2012). All three cooling episodes exhibit salient IV cooling, with domain-averaged rates of –1.90 (P1), –2.90 (P2), and –2.77 (P3) °C decade^–1 and extremes exceeding –4.0 °C decade^–1 east of the Caspian Sea (Fig. 1c–e). This reinforces that such decadal Eurasian cooling episodes, including the 1998–2012 one, are principally governed by atmospheric internal variability^{12,46,48,49,50,51,55,60} rather than external forcings. Notably, three past cooling episodes are consistent across observational SAT datasets (Supplementary Fig. 1) and occur with a cycle of approximately 40 years, whether or not the EX component is considered (Fig. 1a–b). LENS2 historical simulations produce 338 analogous cooling episodes across 100 ensemble members, showcasing similar spatiotemporal characteristics (Supplementary Figs. 2 and 3a), with the recurrence period of 42 ± 17 years (mean ± one standard deviation across 100 LENS2 members).

Fig. 1: Pronounced decadal cooling episodes over Eurasia in reanalysis.

a Time series of JF SAT anomalies (relative to the 1961–1990 mean) averaged over Eurasia (SAT_Eurasia; black solid line; outlined by black boxes in c–e) from 1901 to 2023, and its internal variability (IV; gray bars) and externally forced (EX; black dashed line) components calculated using the LENS2 MMEM and linear regression (“Methods”). Their 15-year trends during P1 (1917–1931), P2 (1958–1972), and P3 (1998–2012) are given, with corresponding p values in parentheses. b Rolling 15-year linear trends (marked on the eighth year) of SAT_Eurasia (upper; 40°–60°N, 45°–110°E) and SAT_BKS (lower; 70°–85°N, 30°–90°E), with gray bars, black solid, and dashed lines refer to IV, EX, and total trends, respectively. Significant IV trends exceeding the 95% confidence level are marked by asterisks. c–e Trend patterns of IV SAT (colors) during c P1, d P2, and e P3, while white hatches outlined by white contours indicate SAT trends are significant at the 95% confidence level.

During Eurasian cooling episodes, while the Barents-Kara Seas (BKS) exhibit warming trends in reanalysis (lower panel of Fig. 1b and Fig. 1c–e), SAT trends over the BKS in LENS2 simulations vary greatly, with only 20% of simulated episodes showcasing significant BKS warming (Supplementary Fig. 3b). This indicates that Eurasian cooling episodes are not always accompanied by pronounced Arctic warming and thus emphasizes the necessity to examine the driver of decadal Eurasian cooling independent of Arctic warming.

SCAND’s leading role with NAO as an amplifier

The P1–P3 Eurasian cooling is directly related to the low-level cold advection, particularly the intense meridional advection (cyan dots in Fig. 2g–i). Such cold advection results from northeasterly anomalies caused by the intensified Siberian High (dark purple colors in Fig. 2g–i) and increased Ural blocking around 50°–70°E (Fig. 2d–f). The quasi-barotropic high-pressure anomaly over Scandinavia and BKS is accompanied by two low-pressure anomalies over western Europe and eastern Russia (Fig. 2a–c), constituting a canonical Scandinavian (SCAND) pattern⁶⁷. The trend patterns of 500-hPa geopotential height anomalies (HGT500) during P1, P2 and P3 (Supplementary Fig. 4a–c, much the same as Fig. 2a–c) show strong spatial correlations (r = 0.76, 0.71, and 0.73, p < 0.01) with the positive SCAND pattern (SCAND⁺; Fig. 3a), which is obtained by rotated empirical orthogonal function analysis (REOF, see details in “Methods”). The trends of SCAND index and SAT_Eurasia—computed over successive 15-year subperiods—show strong anti-correlation in reanalysis (r = –0.82, p < 0.01), with P1, P2, and P3 reaching the largest SCAND trends and the lowest SAT_Eurasia trends (Fig. 3e). This anti-phase relationship between SCAND and SAT_Eurasia trends also holds in LENS2 simulations with significant negative correlations [r = –0.52 (p < 0.01) ± 0.12; Fig. 3g], whereby the overwhelming majority of simulated cooling episodes feature positive SCAND trends (Supplementary Fig. 5). These results suggest a crucial role of SCAND⁺ in Eurasian cooling, through surface anticyclonic anomalies over northern Eurasia^14,15,30 and enhanced Ural blocking frequency or duration^{31,68,69,70,71} (Fig. 2). The surface cold anomalies thus intensified, in turn, act to reinforce the above favorable atmospheric circulations, forming positive feedback^72,73.

Fig. 2: Atmospheric circulation trend patterns during P1, P2, and P3 in reanalysis.

a–c Trend patterns of JF streamfunction at 200 hPa (colors) and associated wave activity flux (Eqs. 4 and 5; arrows, with those greater than 1.25 m² s^–2 decade^–1 being shown) during a P1, b P2, and c P3. d–f Same as (a–c), but for blocking frequency at each longitude (“Methods”; black solid line), with gray bars showing the 95% confidence level. g–i Same as (a–c), but for geopotential height at 1000 hPa (colors). Cyan dots signify regions where 850-hPa cold advection trends equal or stronger than –0.8 °C day^–1 decade^–1. In (a–c, g–i), white hatches outlined by white contours indicate trends shown in colors are significant at the 95% confidence level.

Fig. 3: Relationship between trends of SAT_Eurasia vs. SCAND/NAO index in reanalysis and LENS2.

a Positive phase of the Scandinavian (SCAND) pattern and b corresponding JF-mean SCAND index (gray bars) during 1901–2023. c, d Same as (a, b), but for the North Atlantic Oscillation (NAO) pattern and index. See “Methods” for the definition of SCAND/NAO pattern and index. P1, P2, and P3 trends of SCAND/NAO index are given, with p values in parentheses. e Scatter plot of JF 15-year trends (black dots) between SAT_Eurasia vs. SCAND index along with the linear fit (black line). f Same as (e), but for SAT_Eurasia vs. NAO index. Correlation coefficients and corresponding p values are given, and the P1/P2/P3 trends are marked by gold/blue/red. g Same as the black line in (e), but for linear fits calculated from LENS2 100 individual members (gray lines). The ensemble mean of linear fits (black line), correlation coefficients, and corresponding p values are given. The ensemble mean results of correlation coefficients are similar between those calculated as the average of 100 simulated coefficients and calculated by scattering all 100 members’ simulations in one time.

Far away from our targeted region, the North Atlantic Oscillation (NAO) can exert an impact on temperatures over central Eurasia and Northeast Asia (dominated by horizontal advection) as strong as its impact on the European countries nearby, which is referred to as the NAO downstream teleconnection^{74,75,76,77,78}. This NAO downstream influence on Eurasian SAT is mediated by continuous advection and Rossby waves along the subtropical/polar jet, mainly during winter and on the interannual time scale^79,80,81. Besides, from the synoptic perspective, NAO could modulate the Ural blocking frequency, thereby altering Eurasian temperatures^82,83,84.

Stemming from our analysis, NAO exhibits distinct temporal evolutions in three cooling episodes (Fig. 3d): P1 shows a positive trend of NAO index (NAO⁺); P2, neutral (NAO⁰); and P3, negative (NAO^–). Consistently, 15-year SAT_Eurasia trends show no dependence on NAO trends in reanalysis (r = 0.10, p = 0.74; Fig. 3f), suggesting that NAO trends are not essential for major Eurasian cooling episodes.

While NAO trends are not decisive for Eurasian cooling, they could modulate the cooling magnitude. Cooling rates differ significantly between P1 and P3 (−1.90 vs. −2.77 °C decade^−1) despite their comparable SCAND⁺ trends (1.36 vs. 1.35 decade^−1), coinciding with opposing NAO trends (0.55 vs. −0.82 decade^−1). With a strong SCAND⁺ signature and negligible or neutral NAO trend, the P2 episode presents intermediate cooling (–2.20 °C decade^–1 after excluding the extreme 1969 outlier of SAT_Eurasia < –6 °C). Correspondingly, the increase in Ural blocking frequency in P3 is more significant than that in P1 and P2 (Fig. 2d–f). This offers potential perspectives that NAO⁺ (NAO^−) trends may attenuate (amplify) the SCAND⁺-driven cooling, redefining NAO’s role as a modulator rather than an independent driver. Next, we will validate this hypothesis by synthesizing hundreds of SCAND⁺ cases (“Methods”) with different NAO configurations in LENS2.

Based on LENS2-simulated SCAND time series from 1850 to 2023, we identify 357 significant SCAND⁺ cases (“Methods”). On average, they exhibit a substantial positive SCAND trend of 1.25 decade^–1, which is comparable to the SCAND trends in P1, P2 and P3, and a moderate Eurasian cooling of −1.00 °C decade^–1 (Fig. 4a). Meanwhile, the cooling rates in individual cases exhibit a wide spread, strongly suggesting modulating factors in addition to SCAND⁺.

Fig. 4: Eurasian cooling intensity jointly determined by trends of SCAND and NAO indices in LENS2.

a Scatter plot of JF 15-year trends between SAT_Eurasia vs. SCAND index for all 357 SCAND⁺ cases (gray dots) and their mean (black dot). Reanalysis trends during P1, P2, and P3 are shown as gold, blue, and red crosses, respectively. b Same as (a), but for SAT_Eurasia vs. NAO index for all 357 SCAND⁺ cases, with gold dots standing for 59 SCAND⁺&NAO⁺ cases, blue dots for 233 SCAND⁺&NAO⁰ cases, and red dots for 65 SCAND⁺&NAO^− cases. c Probability distribution and corresponding fitting curve of 15-year SAT_Eurasia trends for SCAND⁺&NAO⁺ (gold), SCAND⁺&NAO⁰ (blue), and SCAND⁺&NAO^– (red) cases. d–f Composite HGT500 trend patterns (contours) of d SCAND⁺&NAO⁺, e SCAND⁺&NAO⁰, and f SCAND⁺&NAO^− group. g–i Same as (d–f), but for composite SAT trends (colors). Gray shadings in (d–f) and white hatches outlined by white contours in (g–i) indicate more than 80% of cases in corresponding groups agree with the composite signs of HGT500 and SAT trends, respectively. Composite trends of SCAND and NAO indices as well as SAT_Eurasia are given on top of each plot.

The trends of SCAND and NAO indices are independent in LENS2 simulations [r = −0.02 (p = 0.22) ± 0.18]. Thus, the NAO trends in LENS2-simulated SCAND⁺ cases are normally distributed relative to zero (Fig. 4b). To explore the additional influence of NAO on the cooling rate caused by SCAND⁺, we classify the 357 SCAND⁺ cases into three categories according to the simultaneous NAO trends: 59 SCAND⁺&NAO⁺ (NAO trends > 0.5 decade^–1; gold), 233 SCAND⁺&NAO⁰ (NAO trends between –0.5 and 0.5 decade^–1; blue), and 65 SCAND⁺&NAO^− (NAO trends <–0.5 decade^–1; red) cases. Their composite circulation patterns uniformly demonstrate a SCAND⁺ structure and differ markedly over the North Atlantic (Fig. 4d–f). Generally, three composite HGT500 trend patterns over SCAND⁺&NAO⁺, SCAND⁺&NAO⁰, and SCAND⁺&NAO^– cases highly correlate with the atmospheric circulation regimes in P1, P2, and P3 (Supplementary Fig. 4a–c; r = 0.75, 0.67, and 0.69, p < 0.01), respectively. Strong SCAND⁺ alone (SCAND⁺&NAO⁰ cases; Fig. 4h) can generate moderate Eurasian cooling of −1.00 °C decade^–1; whereas NAO⁺ weakens the cooling by 50% to −0.48 °C decade^–1 (SCAND⁺& NAO⁺ cases; Fig. 4g) and NAO^− amplifies it to −1.54 °C decade^–1 (SCAND⁺&NAO^− cases; Fig. 4i).

The role of NAO^– as a cooling amplifier is also underscored by distinct probabilities of strong Eurasian cooling (below −1.50 °C decade^–) in the three categories of SCAND⁺ cases (Fig. 4c). By further analyzing LENS2 daily HGT500, we find that the Ural blocking increasing trend is much more remarkable in SCAND⁺&NAO^– compared with SCAND⁺&NAO⁺ or SCAND⁺&NAO⁰ cases, possibly owing to the larger trends of positive HGT tendency induced by transient eddies related to NAO^–83 (Supplementary Fig. 6). With extra more Ural blocking plus reduced warm advection from the North Atlantic (Supplementary Fig. 7), over half of the Eurasian cooling in SCAND⁺&NAO^− cases reach the strong level that bears serious influences, double (quadruple) the strong cooling ratios of SCAND⁺&NAO⁰ (SCAND⁺&NAO⁺) configurations.

Projecting likelihoods of future occurrence

The above results have revealed the dependence of Eurasian cooling on SCAND⁺ and the modifying effect of NAO on cooling magnitude. Given the quasi-40-year periodicity of Eurasian cooling episodes both in reanalysis and simulations, a critical question emerges: will another significant cooling episode follow the 1998–2012 one in the coming decades? To address this question, we assess the likelihood of next IV-driven cooling occurring by 2050 using the framework of SCAND-NAO configurations, and then evaluate the masking effect of EX warming on the cooling under different emission scenarios.

We built a Bayesian multilevel logistic regression model, which quantifies the probability (contours in Fig. 5) of strong Eurasian cooling (i.e., <−1.50 °C decade^–1) based on specific SCAND and NAO trends. This model was trained by LENS2-simulated SCAND, NAO, and SAT_Eurasia trends from 1850 to 2023 (details in “Methods”), and it reflects the influence of SCAND-NAO configurations on the probability of strong Eurasian cooling. According to the model, a one-standard-deviation increase in SCAND trend raises the likelihood by 14.78% (95% credible interval: 13.14−16.51%), while a one-standard-deviation decrease in NAO trend increases it by 7.68% (6.58−8.87%).

Fig. 5: Recurrence probability estimation of strong Eurasian cooling episode by 2050.

Probability (contours) of strong Eurasian cooling (i.e., <−1.50 °C decade^–1) under various SCAND&NAO configurations based on the Bayesian multilevel logistic regression model (“Methods”). Circle (triangle) refers to the unconstrained (constrained) probability of at least one strong Eurasian cooling episode occurring before 2050, which is estimated by using the maximum 15-year SCAND trend during 2025–2050 and concurrent NAO trend within each LENS2 member. While the unconstrained probability indicates the 100-member ensemble mean, the constrained probability indicates the 6-member subsample mean, with the 6 members generally replicating both the recent negative SCAND and positive NAO trends during 2009–2023. Posterior means are given, with 95% credible intervals in brackets.

To maximize early warning of strong cooling risk in coming decades, we only consider the maximum SCAND trend out of all rolling 15-year SCAND trends during 2025–2050 in each LENS2 member, given the dependence of strong cooling on substantial SCAND⁺ trends. The 100-member ensemble mean of maximum SCAND trends during 2025–2050 is 0.80 decade^–1, yielding a very low probability of strong cooling before 2050 at 11.29% (circle, Fig. 5). However, the probabilities spread from near 0% to 90% in individual simulations. To obtain a more reliable estimation, we constrain the prediction by selecting out 6 members that exhibit SCAND decline and NAO increase during 2009–2023 same as in reanalysis (Fig. 3b, d). This constraint elevates the strong cooling probability to 32.57% (triangle, Fig. 5).

The above estimate is based on the IV component alone. Next, we used the Coupled Model Intercomparison Project Phase 6 (CMIP6) projections to obtain the time-dependent EX warming rates across various Shared Socioeconomic Pathways (SSPs; solid lines, Fig. 6). Note that the EX SAT change involves both anthropogenic [e.g., greenhouse gas (GHG) and anthropogenic aerosols] and natural (e.g., volcanic eruptions and solar incoming radiation) signals. The steady and long-term EX warming over Eurasia throughout the end of the current century could be, to a large extent, viewed as GHG-driven or anthropogenically dominated warming.

Fig. 6: Future projections of strong Eurasian cooling episodes.

EX components of JF SAT_Eurasia anomalies (right coordinate; relative to 2015–2025 mean; low-pass filtered; solid lines) in three CMIP6 scenarios of SSP245 (yellow), SSP370 (orange), SSP585 (red), as well as in SSP370-fixed LENS2 (gray) from 2015 to 2099. Light shadings in corresponding colors denote the inter-member spread of ± one standard deviation (29 members from 16 CMIP6 models, Supplementary Table 2; 100 members from LENS2). Reduction rate of strong Eurasian cooling (left coordinate; bars), defined as the fraction of IV-driven strong IV cooling episodes (i.e., <–1.5 °C decade^–1) that become weaker than –1.5 °C decade^–1 due to EX Eurasian warming in the near future (2025–2049), mid-term future (2050–2074), and far future (2075–2099) under corresponding scenarios.

The intrinsic characteristics (frequency and intensity) of IV cooling episodes are comparable between historical (Supplementary Fig. 2) and future climates (Supplementary Fig. 8), but the probability of strong IV cooling is significantly reduced when incorporating the EX component. In the far future (2075–2099), only 6% of IV cooling episodes persist under SSP585 with intense EX warming, compared to 81% under SSP245 with muted EX warming. This sharp diversity provokes a thought of tradeoff: if low-carbon heating alternatives remain limited in the affected region, more rigorous mitigation efforts might seemingly risk increasing fossil fuel demand, but this is beyond the scope of the current study. In contrast, for the near (2025–2049) and mid-term (2050–2074) future, the moderate Eurasian EX warming masks about 40–60% of IV cooling episodes, similar across scenarios. Thus, despite the scenario uncertainty, the strong Eurasian cooling driven by internal variability is likely reduced by approximately 50% owing to the anthropogenic warming.

Discussion

The 1998–2012 Eurasian winter cooling, sparking heated debate^14,57, is not an isolated episode; instead, here we show that episodes with analogous cooling intensity periodically occur every 40 years or so, both in reanalysis (1917–1931, 1958–1972, and 1998–2012) and LENS2 simulations. Such dramatic drop in seasonal mean temperature might be associated with more and stronger cold extremes or less and weaker warm extremes (Supplementary Fig. 9). Importantly, the ability of LENS2 to reasonably reproduce rare, high-impact cooling episodes like the strong 1998–2012 one (–2.77 °C decade^–1) provides a vital tool for their much-needed investigation. Our findings support the conclusion that Eurasian decadal cooling episodes are principally governed by internal variability^{12,46,48,49,50,51,55,60}, and we have disclosed two important atmospheric regimes (SCAND and NAO) that efficiently modulate Eurasian cooling. Specifically, potent positive SCAND trends favorably anchor the occurrence of Eurasian cooling, on top of which NAO trends jointly modify the eventual cooling rate. However, although most of the simulated Eurasian cooling episodes can be explained by our SCAND-NAO framework, the existing small portion of cases that occur with neither SCAND⁺ nor NAO^– also call attention to other atmospheric teleconnections, such as the negative phase of Arctic Oscillation⁶⁰ and East Atlantic pattern⁸⁵.

Interestingly, the 1917–1931 and 1998–2012 Eurasian cooling coincide with the early 20th-century Arctic warming and post-1990 rapid Arctic warming, respectively, the only two epochs with substantial Arctic warming since 1901^86,87. This indicates that, albeit with internal variability playing the dominant role, influences of Arctic warming and attendant sea-ice loss cannot be neglected, especially on the latest Eurasian cooling during 1998–2012 with a relatively elongated cooling area spreading into central Europe^54,55,88,89. Our results show that three historical cooling episodes are accompanied by divergent Arctic sea-ice conditions: the 1998–2012 period shows dramatic BKS sea-ice loss while the 1917–1931 and 1958–1972 periods exhibit no significant decrease (Supplementary Fig. 10), demonstrating sea-ice loss is not the prerequisite for Eurasian cooling. LENS2 simulations show evidence that more pronounced external forcings after 1970 potentially strengthen the coupling between BKS warming or sea-ice loss and Eurasian cooling, probably owing to distinct stratospheric and tropospheric pathways related to the changing external forcings⁹⁰. The relative contribution of Arctic sea-ice change in different cooling episodes should also be separately estimated. However, due to limitations in data accuracy, the comparison of sea ice changes (especially before the satellite era) presented here should be considered for reference only.

Based on the dynamic framework of SCAND and NAO, we have established a logistic regression predictive model, enabling earlier identification of extreme cooling episodes. While LENS2 projections show inherent atmospheric chaos (namely, SCAND/NAO trends varying substantially across members), seeking consensus from individual ensemble realizations could leave society unprepared for impending extremes⁹¹. To address this predictability gap, we applied observational constraints to SCAND and NAO indices, yielding an approximately 30% recurrence probability of strong Eurasian cooling during 2025–2050. One constrained member even projects the highest probability (90%) where strong SCAND⁺ (1.78 decade^–1) and NAO^– (–1.47 decade^–1) trends truly coincide. This warning-level risk—though partially attenuated by EX warming—demands proactive adaptation planning given the potential societal impacts of abrupt cooling events.

We have investigated whether the periodically occurring Eurasian cooling leaves discernible imprints in oceanic forcings like sea surface temperature (SST) anomalies that could enhance predictability. Previous studies have established connections between Eurasian temperatures and both the Atlantic Multidecadal Oscillation (AMO)^{13,49,62,77,92,93,94,95,96,97} and Pacific Decadal Oscillation (PDO)^63,98,99,100, yet these relationships remain inconsistent. While some reports associated positive AMO phases with Eurasian warming^92,93,96,97, others documented cooling effects^{49,62,63,94,95}, likely reflecting competing thermodynamic and dynamical influences of AMO. Our analysis reveals no coherent SST patterns across three cooling periods (Supplementary Figs. 11 and 12). However, we identify a robust linkage between SCAND⁺ trends and enhanced tropical eastern Atlantic convection [reanalysis: r = 0.40, p = 0.13; LENS2: r = 0.31 (p < 0.01) ± 0.16; Supplementary Fig. 13]. This connection may operate through diabatic heating-induced Rossby wave activity^101,102,103, which amplifies the SCAND⁺ pattern and consequently reinforces Eurasian cooling. In addition, cooling-favorable decadal changes in atmospheric regimes like SCAND⁺ may also be shaped by the variations of midlatitude jet stream (acting as a waveguide^78,104), which requires further study.

Methods

Reanalysis and observational datasets

We use two reanalysis products: ECMWF Reanalysis V5 (ERA5¹⁰⁵; regridded to 1° latitude × 1° longitude; from 1940 to present) and NOAA-CIRES-DOE 20th Century Reanalysis V3 (20CRV3¹⁰⁶; 1° latitude × 1° longitude; from 1901 to 2015). Equivalent monthly variables from both ERA5 and 20CRV3 are analyzed, including air temperature at 2 meters (SAT), geopotential height (HGT), horizontal wind (U and V), sea ice concentration (SIC), and precipitation. To reconstruct ERA5-20CRV3 combined records continuously from 1901 to 2023, we first align their climatologies by adjusting 20CRV3 data using ERA5’s long-term mean (over the overlapping period of 1940–2015) as the reference, then merge the ERA5 (1940–2023) and the adjusted 20CRV3 (1901–2015) reanalysis to extend the temporal coverage.

We also employ observational datasets, including HadISST¹⁰⁷ SST (1870-present) and multi-source SAT (Supplementary Fig. 1): the Climatic Research Unit at the University of East Anglia (CRU), the Met Office Hadley Centre (HadCRUT), NOAA (NOAAGlobalTemp), NASA Goddard Institute for Space Studies (GISTEMP), Berkeley Earth, and the University of Delaware (UDel).

CESM2 Large Ensemble and CMIP6 simulations

To validate reanalysis results and provide projections, we analyze 100-member CESM2 Large Ensemble (LENS2⁶⁶; 1850–2099) simulations under the CMIP6 historical (before 2015) and SSP370 (after 2015) scenarios. Additionally, future projections from 16 CMIP6 models (SSP245/SSP370/SSP585, 2015-2099; Supplementary Table 2) are also used. Note that the ensemble mean explicitly refers to the multi-member average across 100 members from LENS2 or across 29 members from 16 CMIP6 models.

Internal variability and externally forced components

To enable comprehensive analysis and robust quantification, SAT and SST anomalies (relative to the 1961–1990 mean) in reanalysis and LENS2 simulations are first decomposed into internal variability (IV) and externally forced (EX) components using distinct methodologies.

We first use the LENS2 multi-member ensemble mean (MMEM) of globally averaged SAT/SST anomalies to represent the externally forced signal ($X$). The SAT/SST anomaly in reanalysis at each grid point $i$ (${Y}_{i}$) is decomposed through simple linear regression^{87,98,108,109} over 1901–2023: ${Y}_{i}={b}_{i}X+{intercept}$, where ${b}_{i}$ is the regression coefficient. Then, the externally forced component (${{Y\_EX}}_{i}={b}_{i}X$) is subtracted from ${Y}_{i}$ to obtain the internal variability component of SAT/SST (${{Y\_IV}}_{i}={Y}_{i}-{{Y\_EX}}_{i}$, also known as the residual). Time series of IV SAT anomalies averaged over a small domain of Eurasia (gray bars in Fig. 1a) is similar to the detrended SAT_Eurasia time series through linear detrending (not shown). Given weaker trends in other variables in reanalysis (geopotential height, zonal and meridional wind, etc.), we directly apply linear detrending to those fields before analysis.

In LENS2 simulations, the separation of IV and EX components is more straightforward. The MMEM directly represents the EX response, while deviations of individual ensemble members from the MMEM capture the IV components. This decomposition is consistently applied to all LENS2 variables across both historical simulations and future projections.

Trend, correlation, and two-tailed Student’s t-test (effective degree of freedom)

Decadal trends over a 15-year window are estimated by regressing the unfiltered time series of interest onto time, and a two-tailed Student’s t-test (degree of freedom, N = 13) is used to evaluate the statistical significance. For significance tests of correlation coefficients calculated between two time series of rolling 15-year trends, we introduce effective degree of freedom (N_e; e.g., p values in Fig. 3e–g) to account for their respective high autocorrelation and the artificial increase in intercorrelation between them. N_e is given by the following equation¹¹⁰:

$${N}_{{\rm{e}}}=N{\left[1+2\mathop{\sum }\limits_{\tau =1}^{\frac{N}{2}}\frac{N-\tau }{N}{\rho }_{\tau }^{A}{\rho }_{\tau }^{B}\right]}^{-1}$$

(1)

where $N$ is the sample size of two variables ($A$ and $B$), while ${\rho }_{\tau }^{A}$ and ${\rho }_{\tau }^{B}$ are the autocorrelation coefficients at the lag of $\tau$ for $A$ and $B$, respectively.

15-year rolling trend analysis: identification of cooling episodes and SCAND⁺ cases

To identify independent decadal episodes from any time series, we first compute linear trends over successive 15-year windows. Take SAT_Eurasia as an example (Fig. 1b), cooling episodes are detected as local minima in this derived trend series where the trend value is (1) statistically significant at the 90% confidence level, and (2) lower than at least 20 surrounding data points (10 years before and after), ensuring each identified cooling episode represents an independent climatic event. The 15-year window is adopted to keep in line with previous literature^12,14,111, and the results have been tested insensitive to different lengths of windows (for example, 11, 13, 17, and 19 years; not shown). Note that an identified cooling episode refers to a significant downtrend of seasonal mean SAT, which is different from the definition of cold winters or cold extremes; while the former might be a favorable background for the latter.

Likewise, SCAND⁺ cases (Fig. 4) are identified using the same detection algorithm above, instead searching for local maxima in the derived SCAND trend series.

Detection of blocking

Daily 500-hPa geopotential height (HGT500) is used to detect instantaneous blocking days according to the one-dimensional blocking algorithm¹¹². First, two geopotential height meridional gradients, GHGS (Eq. 2) and GHGN (Eq. 3), both of fixed latitudinal span of 20° are calculated along each longitude as follows:

$$\mathrm{GHGS}=\frac{{{HGT}500}_{{\varphi }_{0}}-{{HGT}500}_{{\varphi }_{S}}}{{\varphi }_{0}-{\varphi }_{S}}$$

(2)

$$\mathrm{GHGN}=\frac{{{HGT}500}_{{\varphi }_{N}}-{{HGT}500}_{{\varphi }_{0}}}{{\varphi }_{N}-{\varphi }_{0}}$$

(3)

where ${\varphi }_{N}=80^\circ N+\Delta$, ${\varphi }_{0}=60^\circ N+\Delta$, ${\varphi }_{S}=40^\circ N+\Delta$, and $\Delta =-4^\circ ,0^\circ ,4^\circ$. An instantaneous blocking day at a given longitude is identified if at least one value of $\Delta$ exists satisfying both (1) $\mathrm{GHGS}$ > 0 and (2) $\mathrm{GHGN} < -10\,{\rm{m}}{\left({\rm{^\circ }}\mathrm{lat}\right)}^{-1}$. Blocking frequency in the winter season (in the unit of %) refers to the ratio of total instant blocking days to all winter days (59 days in January and February).

Wave activity flux (WAF)

To illustrate the wave energy propagation associated with atmospheric teleconnections, the horizontal wave activity flux (WAF) is calculated at the 200-hPa level as follows¹¹³:

$${\mathrm{WAF}}_{\mathrm{zonal}}=\frac{p\cos \varphi }{2\left|\vec{U}\right|}\left[U\left({{\psi }_{x}^{{\prime} }}^{2}-{\psi }^{{\prime} }{\psi }_{{xx}}^{{\prime} }\right)+V\left({\psi }_{x}^{{\prime} }{\psi }_{y}^{{\prime} }-{\psi }^{{\prime} }{\psi }_{{xy}}^{{\prime} }\right)\right]$$

(4)

$${\mathrm{WAF}}_{\mathrm{meridional}}=\frac{p\cos \varphi }{2\left|\vec{U}\right|}\left[U\left({\psi }_{x}^{{\prime} }{\psi }_{y}^{{\prime} }-{\psi }^{{\prime} }{\psi }_{{xy}}^{{\prime} }\right)+V\left({{\psi }_{y}^{{\prime} }}^{2}-{\psi }^{{\prime} }{\psi }_{{xy}}^{{\prime} }\right)\right]$$

(5)

where ${\psi }^{{\prime} }$ is the 15-year trend (or in most cases referring to the anomaly relative to the climatology) of 200-hPa streamfunction, with ${\psi }_{x}^{{\prime} }$/${\psi }_{y}^{{\prime} }$ and ${\psi }_{{xx}}^{{\prime} }$/${\psi }_{{xy}}^{{\prime} }$/${\psi }_{{yy}}^{{\prime} }$, respectively, referring to its first- and second-order partial derivatives of x or y. The vector $\vec{U}=(U,V)$ is the winter-mean horizontal wind, and $|\vec{U}|=\sqrt{{U}^{2}+{V}^{2}}$ stands for the velocity. The pressure $p$ is normalized by 1000 hPa, and $\varphi$ is the latitude in radians.

Rotated empirical orthogonal function (REOF) analysis: SCAND and NAO indices

Being more statistically stable and simpler to interpret compared with the conventional empirical orthogonal function (EOF) analysis¹¹⁴, the method of rotated EOF (REOF, also named rotated principal component analysis, RPCA) is employed in this study to extract distinct patterns of atmospheric circulation in the Northern Hemisphere following previous studies^67,103,115. First, the EOF analysis is performed on winter monthly detrended anomalies of 500-hPa geopotential height (HGT500) over the whole Northern Hemisphere north of 20°N, and 10 leading EOFs are yielded. Then, we rotate these 10 EOFs using Kaiser row normalization and the varimax criterion (employing the NCL function called eofunc_varimax, see details at https://www.ncl.ucar.edu/Document/Functions/Built-in/eofunc_varimax.shtml). The resultant 10 REOFs can one-to-one correspond to the individual teleconnection patterns provided by the Climate Prediction Centre (CPC; Supplementary Table 3; see details at https://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml), with the extracted REOF2 referring to NAO (explained variance of 10.77%) and REOF6 to SCAND (6.97%) of our concern (Supplementary Fig. 14).

For both reanalysis and LENS2 simulations, we derive the SCAND and NAO indices by projecting the respective HGT500 anomalies onto their corresponding REOF patterns: (1) SCAND index: projection onto REOF6 (Fig. 3a) over 30°–90°N, 60°W–150°E; (2) NAO index: projection onto REOF2 (Fig. 3c) over 30°–90°N, 120°W–90°E.

Bayesian multilevel logistic regression

We model the probability of strong Eurasian cooling episodes (i.e., <−1.5 °C decade^–1) across LENS2 members using hierarchical logistic regression, accounting for shared structure among ensemble realizations from the same model:

$$\Pr \left({{cooling}}_{{ij}}=1\right)={{logit}}^{-1}\left({\beta }_{0}+{\beta }_{1}* {tren}{d}_{{{SCAND}}_{{ij}}}+{\beta }_{2}* {{tren}{d}_{{NAO}}}_{{ij}}+{u}_{j}\right)$$

(6)

where ${{cooling}}_{{ij}}$ indicates whether ensemble member j exhibits strong cooling below −1.5 °C decade^–1 in any 15-year period i, with ${{cooling}}_{{ij}}=1$ for true and ${{cooling}}_{{ij}}=0$ for false; ${u}_{j} \sim {\rm{{\rm N}}}(0,{\sigma }_{u}^{2})$ represent deviations from global intercept ${\beta }_{0}$ for member j; ${\beta }_{1}$ and ${\beta }_{2}$ quantify effects of SCAND and NAO trends (log-odds scale).

Trained by 15-year trends from the 100-member LENS2 historical simulations (1850–2023; using the PyMC python package), our model reveals strong evidence for teleconnection influences on strong Eurasian cooling episodes. In terms of log-odds scale effects, posterior means (95% credible intervals) are 2.121 (1.994, 2.259) for trend_SCAND and −1.507 (−1.623, −1.382) for trend_NAO. Converting into probability-scale impacts yields that a 1-unit trend_SCAND increase raises the likelihood of strong Eurasian cooling episode by 14.78% (13.14−16.51%), while a 1-unit trend_NAO decrease increases it by 7.68% (6.58−8.87%).

The small but significant inter-member dispersion [${\sigma }_{u}=0.39$ (0.29, 0.49)] confirms that internal variability systematically alters cooling probabilities even under identical forcing scenarios and model physics. This variability implies that certain ensemble members exhibit higher cooling recurrence risks than the ensemble mean, with important implications for regional climate preparedness.