Introduction

Eurasia, a critical population center in the Northern Hemisphere, is highly vulnerable to climate change1. Recent decades have seen fluctuating warm and cold extremes in Eurasian surface air temperatures (SAT), with profound impacts on society and ecosystems2,3. In addition to anthropogenic warming, the Atlantic Multidecadal Oscillation (AMO) significantly contributes to SAT variations and extremes in Eurasia through both thermodynamic and dynamic processes4,5,6,7. The thermodynamic component involves SAT changes due to the release of heat fluxes from the North Atlantic8, directly warming the atmosphere and being transported downstream by westerly winds6,9,10. Meanwhile, the dynamic effect entails large-scale atmospheric circulation anomalies, including Rossby wave trains from the North Atlantic that propagate towards Eurasia, influencing SAT and SAT extremes there11,12.

Previous research has primarily emphasized the influence of the AMO on winter SAT in Eurasia through dynamic atmospheric circulation patterns11,12,13,14. Specifically, the Rossby waves generated by the AMO have been linked to contrasting winter SAT anomalies and extreme temperature events in Eurasia11,12,14,15. However, characteristics of the Rossby waves described by different studies vary, raising the question about the importance of the dynamic process affecting winter SAT in Eurasia. Furthermore, there has been relatively little exploration of the thermodynamic impact of the AMO on winter SAT co-variability across mid- to high- latitudes in Eurasia16,17,18. And it remains unclear whether dynamic circulation patterns dominate over thermodynamic heating in shaping this co-variability.

To unravel these complexities, we conducted comprehensive analyses using unforced pre-industrial experiments from the Coupled Model Intercomparison Project, Phase 6 (CMIP6), including the pre-industrial control (piControl) experiment and idealized pre-industrial experiments from the Decadal Climate Prediction Project (DCPP), which introduced positive/negative AMO anomalies (details of CMIP6 models provided in Materials and Methods and Supplementary Table 1). Our investigation aimed to discern the relationship between AMO and SAT over northern Eurasia without an influence of anthropogenic forcing. Strikingly, robust positive associations between the AMO and decadal SAT variations in northern Eurasia (40-70°N, 10°W-140°E) were observed across all piControl simulations (Supplementary Figure 1), with the AMO displaying decadal to multi-decadal periodicity (as shown in Supplementary Table 2). This suggests that the thermodynamic heating effect of the AMO, may exert a stronger influence on winter SAT over northern Eurasia than dynamic circulations on decadal and longer timescales. Thus, our study aims to elucidate the relative importance of thermodynamic and dynamic mechanisms of the AMO on SAT and SAT extremes in northern Eurasia. The findings of the study are expected to improve the understanding of the inherent impact of the AMO on climate variability and change in northern Eurasia.

Results

Weak dynamic wave train

We started by examining the influence of the AMO on geopotential height at 500 hPa (Z500) anomalies over northern Eurasia, aiming to understand the dynamic impact of the AMO. Normally, if the AMO triggers Rossby waves, Z500 anomalies should display a wave-like pattern across Eurasia. However, the analysis presented in Fig. 1a reveals that, in the 9-winter average of Z500 differences between positive AMO and negative AMO in the DCPP experiments, positive anomalies predominate in northern Eurasia. This pattern is particularly evident in the multi-member average using the EC-Earth model (Supplementary Fig. 2a). For other DCPP models like HadGEM3-GC31-MM (Supplementary Fig. 3a), and the IPSL-CM6A-LR (Supplementary Fig. 4a), positive Z500 anomalies are evident in most areas, despite the presence of wave train patterns over northern Eurasia.

Fig. 1: 9-winter mean of the Z500 differences between positive and negative phases of the Atlantic Multidecadal Oscillation (AMO) in the Decadal Climate Prediction Project (DCPP) experiments.
figure 1

a multi-member mean of the 9-winter mean Z500 difference (units: m) for all the DCPP models, with colored boxes representing subregions of northern Eurasia and region numbers displayed at the center. b cumulative sample means of the 9-winter mean Z500 differences averaged over subregions, accumulated from sample sizes ranging from 1 to 107. c box-and-whisker plots for sample means of the 9-winter mean Z500 differences averaged over subregions. In c, boxes and whiskers represent ranges from -1 to 1 standard deviation and from the 1st to the 99th percentile, respectively. These ranges are estimated through bootstrap random sampling (100,000 times). Red dashed boxes indicate the detection of robust signals, signifying a 99% confidence level in the results.

Full size image

The behavior of Z500 in response to the AMO varies depending on the sample size, with different wave train patterns observed among various model members and years. Therefore, to determine whether the AMO stimulates or modifies stationary Rossby waves, it is crucial to have a sufficiently large sample size that can account for the phases of random synoptic waves and ensure an accurate selection of wave centers. To address this, the study initially conducted 9-winter average of the Z500 differences between positive AMO and negative AMO, and collected a total of 107 samples by considering combinations of different model members (i.e., member 1, EC-Earth model; …member 32, EC-Earth model; member 1, HadGEM3-GC31-MM model…member 25, HadGEM3-GC31-MM model; …). Furthermore, northern Eurasia was divided into 12 subregions, each spanning 30° of longitude and 30° of latitude. For instance, region 1 covered 0-30°E and 40-70°N, while region 2 extended 10° of longitude eastward, ranging from 10-40°E and 40-70°N (outlined as region 3-14 in Fig. 1a). The study also identified two wave train centers, located at 30-60°E, 55-70°N, and 60-90°E, 40-50°N, based on observation12 (outlined as region 1-2 in Fig. 1a) to further diagnose the existence of the wave center. Results indicated substantial fluctuations in the mean Z500 until a sample size of 30 was reached. However, beyond this point, all subregions consistently exhibited positive Z500 anomalies or insignificant Z500 anomalies, with values close to zero (e.g., regions 1, 6, 7, as shown in Fig. 1b). For individual DCPP model, samples were collected by considering combinations of different years and model members (i.e., member 1, year 1; …member 1, year 9; member 2, year 1…). Steady positive Z500 anomalies are evident in the EC-Earth model when reaching a sample size of 150 (Supplementary Figure 2b) and in the HadGEM3-GC31-MM model when reaching the sample size of 90 (Supplementary Figure 3b). For the IPSL-CM6A-LR model (Supplementary Figure 4b), the Z500 does not converge when accumulated to 450 samples. Nevertheless, positive Z500 anomalies can be observed across northern Eurasia, except for region 1, 6, 7.

To mitigate any bias introduced by resampling order, the study employed the Bootstrap resampling method19, which involved random sampling 100,000 times to assess the dependence on ensemble size. The sample size was systematically varied from 10 to 100. Regardless of the sample size, the mean Z500 values consistently showed positive anomalies across all subregions, except for region 1, where Z500 anomalies were approximately zero (as depicted in Fig. 1c and Supplementary Fig. 2-4 with sample size varying from 10 to 400). Considering the robustness of signal detection and its dependence on ensemble size, it was found that robust signals (defined as the 1st percentile value reaching above 0 or the 99th percentile falling below 0) could be detected within a sample size of less than 100 in regions 2-4, and 9-14 (as shown in Fig. 1c). This approach effectively helped rule out the influence of stochastic internal variability by accumulating a substantial sample size, confirming that the dynamic effect of wave train is relatively weak and is likely overshadowed by the thermodynamic effect, particularly on decadal and longer timescales.

Predominant thermodynamic effect

In Fig. 2a, the study examined the spatial patterns of SAT differences over northern Eurasia between warm and cold AMO phases in the pre-industrial control (piControl) simulations. During warm AMO phase, compared to cold phase, there is a notable overall increase in winter SAT across northern Eurasia, with warming exceeding 0.6 degree Celsius (°C) in the Europe region and relatively weak warming of more than 0.4°C in the Asian region.

Fig. 2: Winter surface air temperature (SAT) anomalies over northern Eurasia associated with the Atlantic Multidecadal Oscillation (AMO) in the piControl experiment.
figure 2

Multi-model means of the differences in (a) full SAT anomalies, (b) dynamic SAT anomalies, and (c) thermodynamic SAT anomalies (units: °C) between warm and cold phases of the AMO (when the AMO anomalies are greater than 1 or less than -1 standard deviation). Hatched areas in (a-c) indicate significance at the 95% confidence level for SAT in more than 5 models. d correlation coefficients between the AMO and decadal variations in the decomposed thermodynamic SAT anomalies (represented as dots) and dynamic SAT anomalies (represented as triangles). Significant correlations are denoted by solid dots/triangles, while non-significant correlations are shown as hollow dots/triangles. Significance is determined at the 95% confidence level using a Student’s t-test, taking into account the effective degree of freedom (details in Materials and Methods).

Full size image

To gain a deeper understanding of the mechanisms responsible for the AMO’s effect, the study decomposed the full SAT anomalies into dynamic (Fig. 2b) and thermodynamic (Fig. 2c) components using a constructed circulation analog approach (with details provided in Materials and Methods)9,10. The warming observed over northern Eurasia in Fig. 2a is primarily attributed to the thermodynamic contribution, as indicated by the spatial patterns and the magnitude of warming in Fig. 2c. Conversely, the dynamic circulations make less conspicuous positive contributions (less than 0.2 °C) and, in some cases, even negative contributions, as shown in Fig. 2b. Despite variations in the spatial patterns of winter SAT over northern Eurasia among different models, a consistent and dominant thermodynamic effect of the AMO emerges (as observed in Supplementary Figs. 5−20a-c). Specifically, the warm phase of the AMO consistently leads to warmer winters over northern Eurasia, particularly in the Europe region. This analysis underscores the dominant role of thermodynamic processes in driving SAT changes during the warm phase of the AMO, while the influence of dynamic circulation is relatively minor.

The study extended its analysis by investigating temporal correlations between the AMO and decadal variations in winter SAT anomalies, which were calculated using a 10-winter moving average and adjusted for an effective degree of freedom (as detailed in Materials and Methods). This analysis is depicted in Fig. 2d and Supplementary Figures 5-20d, e. The results reveal strong positive correlations between the AMO and thermodynamic SAT anomalies in all models. For example, in the EC-Earth3 (r1i1p1f1) model, the AMO can explain more than 80% of the variance in the thermodynamic SAT with the correlation coefficient exceptionally high at r = 0.94 (as illustrated in Supplementary Figure 10e). In contrast, the correlations between the AMO and dynamic SAT anomalies are much weaker and instability. All the correlation coefficients are less than 0.5, and more than one-third of them are statistically insignificant (7 out of 16 models). Moreover, the positive correlations between the AMO and thermodynamic SAT components (as shown in Fig. 2d) tend to be stronger than the correlations with the full SAT (as depicted in Supplementary Figure 1). This finding confirms the prominent role of the thermodynamic effect of the AMO in driving long-term temporal variations in winter SAT anomalies over northern Eurasia, while the connection between the AMO and dynamic SAT anomalies is notably weaker and less consistent across the models.

The examination of upward surface sensible and latent heat fluxes (referred to as USLH) from the North Atlantic Ocean sheds light on the thermodynamic impact of the AMO. The multi-model average of the differences in USLH between warm AMO and cold AMO phases in the piControl simulations reveals substantial surface heat fluxes being released from the North Atlantic into the atmosphere, with magnitudes exceeding 20 W·m-2 (as illustrated in Fig. 3a and Supplementary Figure 21). In certain models from the CMCC and EC, which exhibit greater AMO variations, more heat releases are observed, surpassing 40 W·m-2 (as shown in Supplementary Figure 21c-d, g). Similarly, in the DCPP experiments, significant USLH releases are evident over the North Atlantic Ocean when comparing positive and negative AMO phases (Fig. 3b and Supplementary Figure 22a, b, c). The advection of this warmer air from the North Atlantic, carried by the climatological westerlies9 (as indicated by positive climatological zonal wind at 850 hPa over the North Atlantic in Fig. 3a), results in pronounced increases in winter SAT across the northeastern hemisphere. Notably, the most prominent anomalies are observed over northern Eurasia (Fig. 3c and Supplementary Figure 22d, e, f). Additionally, we sought to understand the reasons behind the relatively poor performance of two specific models in representing the relationship between the AMO and winter SAT over northern Eurasia during certain periods (Supplementary Figure 23). The result confirms that the variations in USLH, influenced by the AMO, appear to drive changes in SAT, highlighting the critical role of ocean-atmosphere interactions and heat transport in understanding the thermodynamic impact of the AMO on winter SAT over northern Eurasia.

Fig. 3: Thermodynamic heating source from the North Atlantic Ocean and its influence on winter surface air temperature (SAT) over the northern Hemisphere.
figure 3

a multi-model mean of differences in upward surface sensible and latent heat fluxes (USLH, shadings, units: W·m^-2) between warm and cold phases of the Atlantic Multidecadal Oscillation (AMO) in the piControl simulations. It is accompanied by multi-model average for the climatological (years 100−149) mean state of the zonal wind at 850 hPa (contours, units: m·s^-1, with positive in an eastward direction shown in solid lines and negative values in dash lines). Multi-model composition of differences in winter USLH (b, units: W·m^-2) and SAT (c, units: °C) between positive and negative AMO phases in the DCPP models. Hatched areas in (b) and (c) indicate significance at the 95% confidence level for winter USLH and SAT, respectively. The black box in (c) represents the northern Eurasia region selected in this study.

Full size image

Equivalent but opposite impact on warm and cold extremes

We analyzed frequency histograms of winter SAT over northern Eurasia in multi-model piControl simulations, employing the generalized extreme value (GEV) distribution fitted using the Fischer-Tippett Theorem20. These analyses were conducted for both warm and cold phases of the AMO, as shown in Fig. 4a. Remarkably, the SAT distributions during warm and cold AMO phases exhibited similar GEV shapes, with shape parameters falling below 0, indicating Weibull distributions. The steepness parameters were 1.57 for warm AMO and 1.68 for cold AMO phases. Notably, the mean SAT shifted from approximately -0.66 °C during cold AMO to 0 °C during warm AMO, leading to a distinct contrast in the location parameters (-1.11 for warm AMO and -0.49 for cold AMO). Furthermore, the bilateral GEV tails, which represent SAT extremes21, displayed a balanced increase or decrease as the AMO transitioned. This indicates that the shift in mean SAT has equivalent effects on both warm and cold extremes during AMO phases.

Fig. 4: Decadal variations in the frequencies of winter surface air temperature (SAT) extremes over northern Eurasia in response to mean SAT shifts induced by the Atlantic Multidecadal Oscillation (AMO) in the piControl simulations.
figure 4

a multi-model composition of frequency histograms (bars) and their generalized extreme value fittings (curves) for winter SAT anomalies (units: °C) over northern Eurasia during warm and cold phases of the AMO. b decadal variations in the frequencies for the sum of winter warm and cold extremes over northern Eurasia in different models. Values in the brackets of the legend represent the difference in the total SAT extremes between warm and cold AMO phases. c difference in decadal variations in the frequencies of winter warm (dots) and cold (asterisks) extremes between warm and cold AMO phases. Correlation coefficients between the AMO and decadal variations in warm/cold extremes are labeled near the dots/asterisks for different models (shown as the x-axis, in the same color as the dots/asterisks).

Full size image

The sums of frequencies for warm and cold extremes, computed as detailed in Materials and Methods, exhibited minimal variations of less than 2 days per season over a period spanning more than 180 years across all models, as depicted in Fig. 4b. The differences in total SAT extremes between warm and cold phases of the AMO in all piControl simulations, as indicated in the legend of Fig. 4b, were negligible, with values falling below 0.15. This underscores the equivalent but reverse impact of the AMO on decadal variations in warm and cold extremes over northern Eurasia, as also evident in Fig. 4c. Statistical analysis confirmed significant positive correlations between the AMO and decadal variations in the frequencies of warm extremes, along with negative coefficients for cold extremes (as depicted in Supplementary Figure 24). Moreover, the frequencies of SAT extremes exhibited dependencies on the intensities of the AMO, as shown in Fig. 4c and Supplementary Figure 24. For instance, greater AMO variations were accompanied by more SAT extremes in models such as CMCC-ESM2 and the EC-Earth3 models (r1i1pif1 and r2i1p1f1 members). In summary, the primary thermodynamic effect of warm (cold) AMO has relatively equivalent but reverse impacts on increasing (decreasing) warm extremes and decreasing (increasing) cold extremes over northern Eurasia, all mediated through the shift in mean SAT.

Discussion

The study’s findings suggest that the internal variability of the AMO primarily influences winter SAT and its extremes over northern Eurasia through thermodynamic heating rather than dynamic wave trains. Previous studies have investigated the physical processes underlying the influence of the AMO-related USLH on winter Eurasian temperatures8. The Atlantic warm sea surface temperature (SST) anomalies increase local upward latent heat flux, thereby leading to anomalous net moist static energy (MSE) flux input into the atmosphere. Additionally, the advection of this warmer air from the North Atlantic, can be transported by the climatological westerlies9 and through water vapor feedback or Planck feedback22,23, resulting in pronounced increases in winter SAT over northern Eurasia.

It is noteworthy that apparent wave train patterns are observed in the difference in geopotential height at 500 hPa (Z500) between positive and negative AMO phases in the DCPP experiments of the EC-Earth model (hereinafter referred to as AMOp - AMOn DCPP, Supplementary Figure 25a-i) and in certain warm AMO years in the CMCC-CM2-SR5 model of the piControl experiment (Supplementary Figure 26a, b). Interestingly, when the zonal mean of Z500 across the Atlantic and Eurasia (90°W-180°E) is removed, potentially induced by the thermodynamic effect of the AMO, similar Rossby wave patterns are observed in the Z500 difference of the 9-winter in the AMOp - AMOn DCPP (Supplementary Figure 25j) and 10-winter mean of the warmer AMO phase in the piControl simulation (Supplementary Figure 26d). This suggests that the AMO could indeed modulate dynamic stationary Rossby waves, but the dynamic effect of these waves is relatively weaker than the thermodynamic effect in influencing winter SAT over northern Eurasia.

Recent observations have indicated inter-decadal variations in SAT extremes over northern Eurasia16,18, which cannot be explained solely by the equivalent but opposite impact on warm and cold extremes induced by the thermodynamic effect of the AMO or anthropogenic forcings24. However, along with anthropogenic warming25,26 and the northward heat transport facilitated by the AMO27, the rapid Arctic warming, known as “Arctic amplification”, affects the meridional thermal gradient. This, in turn, leads to abnormal atmospheric circulations, such as intensified and prolonged Ural blockings, which increase the likelihood of cold extremes over northern Eurasia24,28,29. These processes may partially explain the observed increase in the cumulative occurrences of both of warm and cold extremes in recent decades. Additionally, the “cold Eurasia” phenomenon associated with the “warm Arctic cold Eurasia” pattern is more likely induced by cold extremes rather than changes in the mean SAT. Even when external forcing signals and SAT extremes are removed from the observational SAT variation (as demonstrated in Supplementary Figure 27 using a linear regression method30, with further details provided in Materials and Methods), the AMO still exhibits a consistent warming effect over northern Eurasia. This highlights the ongoing prominence of the thermodynamic heating in influencing the region’s climate.

When comparing Eurasia SAT during positive and negative AMO years in the piControl experiment, other internal variabilities of the climate system, such as the Interdecadal Pacific Oscillation (IPO), may be included when it is in the same phase as the AMO in Fig. 2. However, the IPO signal can be eliminated when conducting the multi-model mean. In contrast, the DCPP experiment can be regarded as solely capturing the AMO effect when focusing on the difference between positive and negative AMO in the DCPP experiments. Additionally, the AMO’s substantial impact on winter SAT and frequencies of SAT extremes over northern Eurasia still holds even in historical simulations that encompass all external forcings (Supplementary Figure 28, Supplementary Figure 29). And if we decompose the winter SAT anomalies in northern Eurasia utilizing CMlP6 AMO pacemaker experiment of the hist-resAMO models (details provided in Materials and Methods) which also encompass all external forcings, our conclusion regarding the dominant role of the AMO’s thermodynamic processes of in driving SAT changes remains valid (Supplementary Figure 30), notwithstanding the extensive debates about the nature of the AMO31,32,33,34,35,36,37, Consequently, the AMO cycle holds substantial importance for predicting decadal variations in future SAT and SAT extremes. Looking ahead, the North Atlantic subpolar gyre is currently cold and is highly likely to transition into a cold-phase AMO in the 2020s38. This transition presents two potential scenarios for the future SAT and SAT extremes in northern Eurasia. In the first scenario, if the cold AMO effect outweighs anthropogenic warming, northern Eurasia may experience a cooler mean SAT, fewer warm extremes, and more cold extremes. In the second scenario, assuming that anthropogenic warming prevails over the cold AMO, SAT may become warmer, and there may be an increase in both warm and cold extremes, influenced by Arctic amplification.

Methods

Warm and cold extremes

SAT extremes in this study are defined according to the guidelines provided by the Expert Team on Climate Change Detection and Indices (ETCCDI)39. Specifically, a cold extreme is identified if SAT anomaly falls below the 10th percentile of the SAT anomaly distribution, while a warm extreme is identified if the SAT anomaly lies above the 90th percentile of the SAT anomaly distribution. The SAT distribution is determined using all 5-day intervals centered on each calendar day during years 100-149 of the entire period in the piControl and historical experiments. Winter in this study encompasses the period from November in the current year to February in the next year, with the exclusion of February 29 due to its occurrence in leap years.

CMIP6 model simulations

This study utilizes data from four types of simulations within the CMIP6 framework: piControl, DCPP, historical, and hist-resAMO simulations. The piControl experiment involves 18 ensemble members from 16 models and serves as a control simulation with non-evolving pre-industrial conditions, providing insights into the unforced variability of climate system. The conditions in the piControl simulations are chosen to be representative of the period prior to the onset of large-scale industrialization, with a reference year of 185040. It’s worth noting that the piControl experiments have varying time lengths across different models, with some models having durations of less than 500 years due to data limitations. To ensure consistency and for analytical purposes, this study redefines the years starting from year 1. The AMO index is defined as the area average of the 10-year low-pass filtered North Atlantic SST anomalies within the region spanning 0°-65°N and 80°W-0°. Different piControl models exhibit diverse AMO periods (Supplementary Table 2).

The DCPP experiments consist of a total of 107 members from three different models (32 members from the EC-Earth3 model, 25 members from the HadGEM3-GC31-MM model, and 50 members from the IPSL-CM6A-LR model). These experiments are designed to be idealized, involving the imposition of positive/negative AMO anomalies. The primary objectives are to investigate the climate impacts of these AMO anomalies and to understand how models respond to slowly evolving SST anomalies in the Atlantic. In these experiments, SSTs are restored to the model’s climatology plus a positive AMO anomaly over the North Atlantic region (10°-65°N). Outside of this restored region, the model is allowed to evolve freely, enabling a full climate system response41. The duration of the DCPP experiment spans 10 years, from which data for 9 years of winters can be derived and analyzed.

The historical experiment comprises a total of 341 members from 44 different models. These experiments are designed as all-forcing simulations of the recent past and involve different realizations (r), initialization (i) schemes, physics (p), and forcing (f) indices. They are used to evaluate model performance in simulating both present climate conditions and observed climate change40.

The hist-resAMO experiment, containing 5 members of 3 models, is initialized from the historical run year 1870 and integrated up to the year 2014 with historical forcings and is a pacemaker-coupled historical climate simulation that includes all forcings but with SST restored to the model climatology plus observational historical anomaly in the AMO domain (0-70°N, 70°W-0°) using the same model resolutions as the CMIP6 historical simulation42.

Daily outputs from these experiments include data on near-surface air temperature (SAT), sea level pressure (psl), upward sensible heat (hfss), upward latent heat (hfls), geopotential height at 500 hPa (Z500), and zonal wind at 850 hPa (U850). To ensure consistency and comparability, all model outputs from the piControl and historical simulations were interpolated to a common resolution of 2.5°×2.5°. A detailed list of the models used in these experiments, along with their specific configurations, can be found in Supplementary Table 1 for reference.

Significance test

A 10-year running mean was conducted to represent the decadal variation in SAT and SAT extremes. This may reduce the effective degree of freedom (EDF) when calculating statistical significance of correlation coefficients. Thus, EDF is estimated based on the following equation43.

$$\frac{1}{{N}^{\ast }}\,\approx \frac{1}{N}+\frac{2}{N}\mathop{\sum }\limits_{j=1}^{N-2}\frac{N-j}{N}\beta xx(j)\beta yy(j)$$
(1)

where N is the sample size, \(\beta xx(j)\) and \(\beta yy(j)\) are the autocorrelation coefficients of the two time series at lag j.

Dynamical adjustment using the constructed circulation analog approach

In this study, a method based on the constructed circulation analog approach9,10 was employed to estimate the dynamical and thermodynamical components of the SAT anomalies. Here’s how this method was implemented:

  1. 1.

    For each winter month in each piControl simulation, we selected a set of \({N}_{a}\)=50 closest sea level pressure (SLP) analogs. These analogs were chosen based on their similarity to the target SLP field, with the selection process involving ranking them by their Euclidean distance from the target SLP field. To ensure that neighboring years were not included, analogs occurring within 5 years of the target month were excluded9,10.

  2. 2.

    From the \({N}_{a}\) analogs, we randomly subsampled \({N}_{s}\)=30 analogs. We then computed the optimal linear combination of these \({N}_{s}\) analogs that best matched the target SLP field using multiple linear regression.

  3. 3.

    The dynamically induced SAT anomaly field was defined as the optimal linear combination of the SAT anomalies associated with the Ns selected SLP analogs. This was achieved by multiplying the SAT anomaly for the corresponding month by the regression coefficients obtained in the previous step.

  4. 4.

    We repeated the random sampling procedure described in step 2 a total of 40 times.

  5. 5.

    Finally, we averaged the Nr optimal sets of SLP analogs to obtain a “best estimate” of the circulation-induced component of the SAT anomalies.

  6. 6.

    To calculate the residual thermodynamical anomaly, we subtracted the dynamical anomaly (obtained in step 5) from the full anomaly.

  7. 7.

    Additionally, we estimated the uncertainty associated with the 40 random subsamples used to construct analogs. For each month in a season, one of the subsamples was randomly selected, and this process was repeated 1000 times. The uncertainty was measured based on the resulting distribution, and the shading in Supplementary Figures 5-20d, e indicates the 10%–90% range of this distribution for each year.

  8. 8.

    When comparing the SAT difference between the warm and cold AMO phases, 10-winter moving averages were conducted for the undecomposed SAT, thermodynamical SAT and dynamical SAT.

A linear regression method to obtain the internal SAT variations

The response of the global-mean SAT to historical external forcing changes can be calculated by the ensemble mean of the historical simulations. This estimated variations associated with forced SAT time series were then removed from the observed SAT series at each grid point via a linear regression method30. The method is briefly described below.

Generally, each \(\Delta T\) series are considered as consisting of an externally-forced component (denoted by subscript F, due to solar cycles, GHGs, aerosols and other external forcings) and an unforced component (denoted by subscript I, due to internal climate variabilities). For the global-mean \(\Delta T\) from observations,

$$\overline{T}(n)=\overline{{T}_{F}}(n)+\overline{{T}_{I}}(n)$$
(2)

where \(T(n,i)\) is the winter SAT anomaly (\(\Delta T\)) at grid point i for year n from observations, where n = 1, …, n = 55 for 1960-2014. The over bar denotes area-weighted average over the globe, and \(\overline{T}(n)\) is the global-mean of \(T(n,i)\) from observations, and \({\overline{T}}_{m}(n)\) is the global-mean SAT anomaly for year n from the multi-model ensemble mean of CMIP6 historical simulations. Since the internal variations among individual realizations are usually uncorrelated and thus are averaged out over a large number of ensemble members, the global mean \(\overline{{T}_{m}}(n)\) from the large historical ensemble contains mostly externally-forced response with little internal variation. Then, \(\overline{{T}_{F}}(n)\) is estimated using

$$\overline{{T}_{F}}(n)=\overline{{b}_{F}}\overline{{T}_{m}}(n)$$
(3)

where \(\overline{{b}_{F}}\) is the regression slope in \(\overline{T}(n)=\overline{{b}_{F}}\overline{{T}_{m}}(n)\) + ε (residual). Thus, the global-mean internal signal that can be estimated as

$$\overline{{T}_{I}}(n)=\overline{T}(n)-\overline{{b}_{F}}\overline{{T}_{m}}(n)$$
(4)