An NAO-dominated mode of atmospheric circulation drives large decadal changes in wintertime surface climate and snow mass over Eurasia

The leading mode of wintertime atmospheric variability over the North Atlantic-North Eurasia sector is dominated by the North Atlantic Oscillation (NAO) and accounts for more than one third of the total variability. This study explores the influences of the leading mode on decadal climate variability of Northern Eurasia. We focus on the little-explored decadal covariations of surface air temperature (SAT), snowfall, snow water equivalent (SWE) and snow cover over the region, using extensive model output from the Coupled Model Intercomparison Project sixth phase. Recent decadal trends (−0.92σ per decade) in the leading mode identified, are found to be largely a manifestation of internal climate variability (at least two thirds from the most conservative estimate). These internally-generated decadal trends strongly contributed to recent trends in SAT, snowfall, SWE and snow cover over Eurasia. External forcings should have played a minor role over Eurasia as they usually suggest opposite decadal trends to those observed. An exception is found for snowfall and SWE in east Eurasia, for which external forcings may have driven a large part of the recent upward trends, equally as important as the NAO-dominated mode. This points to a complex interplay between internally-generated and externally-forced climate variability over Northern Eurasia. Model discrepancies are identified in reproducing the linkages between the leading mode and the Eurasian surface climate variability. The internally-generated variability of this leading mode thus represents a large source of uncertainty in future decadal climate projections over Eurasia and, due to the memory effects of snow, also in modelling springtime climate variability.


Introduction
The mid-to-high latitudes have experienced dramatic climate changes since about 1979, the beginning of the satellite era. These changes are dominated by Arctic amplification (AA), Arctic sea-ice loss, and counterintuitive Eurasian winter cooling trends (Cohen et al 2014(Cohen et al , 2020.
Many studies have attempted to link Eurasian winter cooling trends since the early 1990s with AA and Arctic sea-ice loss, but with differing conclusions. Some studies have argued for a driving role of the Arctic sea-ice reduction (e.g. Honda et al 2009, Screen and Simmonds 2010, Mori et al 2014 in the cooling trends, while others have contested this (e.g. McCusker et al 2016, Sun et al 2016, Ogawa et al 2018, Blackport et al 2019, Ye and Messori 2020. Others yet have suggested an indirect role of Barents-Kara Seas (BKSs) sea-ice (e.g. Luo et al 2016, Matsumura and, mediated by atmospheric blocking or tropical ocean-driven teleconnections. However, Woollings et al (2014) have argued for a marginal impact of temperature variability in the BKS region on the occurrence of cold European winters. Therefore, there is still no consensus on the driving factors of the recent Eurasian winter cooling trends.
Other variables such as snowfall and snow accumulation during the winter season can modulate climate variability in Eurasia too. These variables are closely linked and reflect complicated dynamic and thermodynamic processes. Snow exerts significant influences on ecosystems (Nobrega andGrogan 2007, Sullivan et al 2008), surface energy balance and atmospheric circulation (Barnett et al 1989, Ye et al 2015, Ye and Lau 2018. Snow accumulation, represented by snow water equivalent (SWE), is an important integrator of winter seasonal atmospheric circulation variability and climate variability, triggering a memory effect (Ye 2018). This affects snowmelt in spring. The driving factors of SWE variability during cold months in Eurasia during recent decades have been studied in Sun et al (2019) and Ye (2018), with North Atlantic warming and the North Atlantic Oscillation (NAO)/Arctic Oscillation (AO) being important factors on decadal and interannual timescales, respectively. Bailey et al (2021) further noted that Arctic sea-ice loss may contribute to extreme European snowfall. Climate projections suggest that daily snowfall frequency across much of the Northern hemisphere will decrease, but frequency of daily heavy snowfall will increase in some regions (Danco et al 2016. However, concurrent variations in temperature, SWE, snow cover and snowfall in recent decades and the relative importance of internally-generated versus externally-forced variability have not been well studied. Here, we assess the role of the leading mode of wintertime atmospheric variability over the North Atlantic-North Eurasia sector (see Ye and Messori 2020) in driving the recent decadal trends in surface climate in Eurasia, including snowfall, SWE and temperature. We use data from the Coupled Model Intercomparison Project sixth phase (CMIP6; Eyring et al 2016) to better understand the role of internally-generated versus externally-forced variability. We quantify their relative contributions by making use of large-ensemble simulations (ensemble size ⩾25) from six CMIP6 models. Results from the analysis significantly improve our understanding of decadal climate variability in Eurasia by providing an assessment of internally-generated climate variability in recent decades in the region. While model simulations are subjected to uncertainties, the use of a multimodel, large-ensemble approach and physical reasoning-based analysis in this study are our attempts to mitigate these uncertainties. We further raise the issue of uncertainties by discussing the agreements and disagreements between models and suggesting possible limitations in model physics.

Data
We use monthly-mean geopotential height at 500 hPa (Z500) and surface air temperature (SAT) from both the fifth generation European Centre for Medium-Range Weather Forecasts atmospheric reanalysis (ERA5; Hersbach et al 2020) and the ERA-Interim reanalysis (Dee et al 2011). These two reanalysis data sets (0.5 • × 0.5 • latitude-longitude grid for the former and 0.75 • × 0.75 • latitude-longitude grid for the latter) provide near-identical results in our study and only results from the ERA5 reanalysis are shown. The monthly indices of NAO/AO are provided by the Climate Prediction Centre of the National Oceanic and Atmospheric Administration. Monthly-mean SAT and precipitation on a 0.5 • × 0.5 • latitudelongitude grid for the period 1901-2018 are provided by the Climatic Research Unit (CRU version 4.03) at the University of East Anglia (Harris et al 2014).
To obtain snowfall estimates (snowfall hereafter), a monthly snow accumulation and melt model (McCabe and Wolock 2010) was used with CRU monthly-mean SAT and precipitation as inputs. Monthly and daily GlobSnow-2 SWE data on a 25 km × 25 km equal-area scalable Earth grid over non-alpine areas are provided by the European Space Agency, covering the period 1979-2016. The SWE data was re-gridded to a 0.5 • horizontal resolution. The weekly Snow Cover Extent Dataset from the Global Snow Lab-Rutgers University Climate Lab is used to depict snow cover extent (Robinson and Estilow 2012); it is referred to as snow cover fraction (SCF) as in model output described below.
Model output from both the Historical and the Atmospheric Model Intercomparison Project (AMIP) experiments is analysed. For most models, the Historical experiment covers 1850-2014 and the AMIP experiment covers 1979-2014. A brief summary of the models is provided in table S1 (available online at stacks.iop.org/ERL/17/044025/mmedia). The original model output was re-gridded to a 2.5 • horizontal resolution. Only the first ensemble members (labelled as 'r1i1p1f1') are considered for multimodel mean results or the inter-model spread. Large ensemble historical simulations from CanESM5_p1, CanESM5_p2, CNRM-CM6-1, IPSL-CM6A-LR, MIROC6,and NorCPM1,with 25,40,30,32,50,and 30 realizations respectively, are used to study forced versus internal climate variability. For all data, the months of December, January, February, and March are referred to as 'winter' . In this study, internallygenerated and externally-forced variability refers respectively to climate variability generated owing to processes within the climate system without external forcings and variability generated owing to external forcings, such as increasing greenhouse gases.

Methods
Europe, Russia, and East Eurasia are defined as the regions bounded by 10-40 • E, 40-60 • N, by 50-130 • E, 60-75 • N, and by 90-140 • E, 40-60 • N, respectively (figure 1(d)). Note that these domains do not cover exactly the geographic areas that are suggested by these names. The Theil-Sen trend estimate was used to estimate trends, and their significance was assessed by the Mann-Kendall non-parametric test for monotonic trends (Mondal et al 2012).
Empirical orthogonal function (EOF) analysis was applied separately to both the reanalysis and CMIP6 model output (weighted by the square root of the cosine of latitude) for the available winters, to obtain the leading mode of geopotential height at 500 hPa (hereafter Z500-PC1) over the region 90 • W −160 • E, 20-90 • N. The results are not sensitive to reasonable domain changes (Ye and Messori 2020). The first leading modes in all the ensemble members in the CMIP6 data are qualitatively similar to those of the reanalysis, although the detailed structures differ. These leading modes in both reanalysis and model output are well separated from other modes according to the North et al (1982) criterion.
The procedures for estimating the contribution of Z500-PC1 to decadal trends in a particular variable (e.g. SWE) are described separately for reanalysis/ observational data and model output as follows. The time series of Z500-PC1 and the particular variable are denoted as TsA and TsB, respectively. For reanalysis/observations, a five year running average of the time series was subtracted from the raw time series to retain interannual variability for both TsA and TsB, and then the linear least-square trends were removed before obtaining regression coefficients. In computing the five year running average, reflective (symmetric) conditions were applied to the end points of the time series to avoid data points loss. Any trend may be introduced due to this was then removed by the linear least-square method. The trends in TsA were then multiplied by these regression coefficients to estimate the contribution of TsA to TsB. For model ensemble data, regression coefficients of decadal trends in TsB against decadal trends in TsA across ensemble members were calculated, and then multiplied by the decadal trends in TsA across ensemble members to obtain the decadal trends in TsB as contributed to by TsA for each ensemble member.

Observed impacts on surface climate and SWE
The wintertime leading mode of atmospheric variability over the North Atlantic-Northern Eurasia region (figure 1(a)) is dominated by the NAO, and represents the complex interactions between the NAO, downstream circulation, and other teleconnection patterns (Ye and Messori 2020). As suggested by Dommenget and Latif (2002), EOF analysis may lead to artificial patterns without physical soundness. To prevent this from happening, regression patterns of the Z500 field against the NAO indices (figures S1(a) and (c)) and area-mean Z500 over the eastern Asia region (figures S1(b) and (d)) are compared with figure 1(a). The comparison suggests that this NAO-dominated circulation pattern reasonably reflects physical interactions between the NAO and the downstream circulations. The NAO and its downstream developments have been proposed in a number of studies as important modulators of Eurasian climate (e.g. Hurrell 1996, Branstator 2002, Watanabe 2004, Sung et al 2011, Bollasina and Messori 2018, Ye and Messori 2020. The collaborative impacts of the NAO and downstream circulations on extreme cold weather have been highlighted recently in several studies (e.g. Yao et al 2016, Li et al 2020.
There are two main reasons for the use of the NAO-dominated EOF mode of the atmospheric circulation rather than a conventional NAO index. Firstly, this ensures consistency and comparability with the results of Ye and Messori (2020). Secondly, this allows a more robust representation of the downstream impacts of the NAO than focusing on a North-Atlantic based index. As has been demonstrated by previous studies, an important mechanism for such downstream impacts is the downstream propagation of Rossby waves from the North Atlantic to east Eurasia and their interaction with the surface climate in the region (see Sung et al 2011, Bollasina and Messori 2018, Ye and Messori 2020. A pronounced downward decadal trend is observed for the period from the early 1990s to the early 2010s (black lines in figure 1(b)), consistent with similar decadal downward trends observed in the NAO/AO (blue lines in figure 1(b)). The negative decadal trend in the wintertime NAO for a similar period (1989-2013) is a robust feature (Iles and Hegerl 2017). The Z500-PC1 trend during 1990-2013 (a 24 year window) is −0.95σ. This is the dominant feature during the satellite era, and has been discussed in previous studies (Cohen et al 2014;1990in Ye and Messori 2020. Such a downward trend has a strong link to winter cooling trends, and the number of extremely cold days in Eurasia (Ye and Messori 2020). Here, we further link it to snowfall and SWE and try to address covariations of these variables.
The wintertime SWE shows a pronounced tripole-like trend in Eurasia during the period of a strong decrease in Z500-PC1 (figure 1(c)), with upward trends over Europe and East Eurasia and downward trends over Russia. This trend pattern of wintertime SWE shows notable differences to that of wintertime SWE during a longer period 1980-2016 (not shown), and also to that of March SWE (annual maximum snow mass) computed for a longer period 1980-2018 based on GlobSnow v3.0 data (Pulliainen et al 2020). Thus, the trend pattern of winter SWE during 1990-2013 evidently seems to be a temporal fluctuation. The detrended regression patterns of SWE against Z500-PC1 closely resemble the SWE trends (cf figures 1(d) and (c)), particularly for Europe and Russia. This suggests that the trends of SWE and Z500-PC1 may be related. Their link can be further investigated by studying the impacts of Z500-PC1 on SAT, snowfall, and SCF (figure 2), because these are important factors for SWE variability (Ye 2018). During the negative phase of Z500-PC1, cooler SATs are found across Eurasia ( figure 2(a)). Snowfall shows a roughly south-north contrasting pattern, with positive anomalies over Europe and East Eurasia, and negative anomalies across Northern Eurasia ( figure 2(c)). The impact on snowfall is particularly large in Europe, where climatologically relatively high SATs are observed (orange contours in figure 2(a)). Note that the relatively weak link between Z500-PC1 and snowfall over the Urals region seems to explain the relatively weak link between Z500-PC1 and SWE in the same region ( figure 1(d)). Negative SAT anomalies and positive snowfall anomalies sustain larger SCF in East Eurasia and particularly Europe (figure 2(e)), which in turn favours snow accumulation. This explains the significantly more frequent high-SWE days during negative phases of Z500-PC1 (figures 2(b) and (f)). The reduced snowfall over Russia accounts for fewer high-SWE days ( figure 2(d)). Therefore, the link between the trends in SWE over Eurasia and in Z500-PC1, as suggested in figures 1 and 2, is supported by relevant atmospheric and surface processes that contribute to the SWE trends. The question then is: have the recent trends in Z500-PC1 been exceptional, and what are the impacts of internally-generated versus externally-forced Z500-PC1 trends on the surface climate in Eurasia? To address this question, the CMIP6 model output is used to study the 24 year-window decadal trends in Z500-PC1.

Externally-forced versus internally-generated decadal trends in Z500-PC1 in ERA5 data and CMIP6 models
The standardized principal components (PCs) derived from the CMIP6 Historical and AMIP simulations exhibit considerable inter-model variability during both the satellite era and the period of interest spanning 1990-2013 ( figure 3(a)). The mean trends of the AMIP and Historical simulations from 1990 to 2013 are relatively small compared to those from ERA5 (0.01σ and −0.18σ versus −0.95σ, respectively; figure 3(a)). We interpret these as an externally-forced response, assuming cancellation of internal variability effects. The strongest individual negative trend is −0.7σ (historical simulations in CanESM5_p2), which is still smaller than that in ERA5. The externally-forced response is stronger and the ensemble range is larger in historical simulations than in AMIP simulations. One may thus infer that atmosphere-ocean interactions have enhanced both these quantities. This discrepancy between both Historical and AMIP simulations and ERA5 suggests that the remarkable recent trend in Z500-PC1 in ERA5 is likely a manifestation of large internal climate variations, that is internally generated. Different model formulations (e.g. dynamical core and parameterizations) can also contribute to the CMIP6 inter-model variability.
A similar EOF analysis was repeated for ensembles of the six selected CMIP6 models ( figure 3(b)). All PCs are standardized. The large inter-model spread across the CMIP6 models, and inter-ensemble spread within the individual CMIP6 models (small red dots in the figure mark one standard deviation), suggest an important role of internal climate variations during the period 1990-2013. The ensemble-mean responses of the individual CMIP6 models suggest a weak to moderate external forcing of Z500-PC1 towards negative trends, especially for CanESM5_p1/p2. The CanESM5 family of models is known to respond strongly to CO 2 forcing (Swart et al 2019), and the two models produced fairly large global warming trends (figure S1(a)). Nevertheless, these ensemblemean responses are still much smaller than the trend revealed by ERA5. The trend in ERA5 is located at the lower end of the ensemble ranges for CanESM5_p2, MIROC6, and NorCPM1, and is close to the lower range of that in CanESM5_p1. Therefore, the recent notable downward trend in Z500-PC1 is considered primarily a result of large internal climate variations, which are nonetheless within the modelled range. The externally-forced responses are mostly of the same sign as-but much weaker than-the trend in ERA5.
A further issue that needs to be addressed, is whether the skewness of the trend distribution for the period 1990-2013 is caused by the rapid and consistent decadal warming since about 1980 (figures S2(a) and (b)). The statistics for the 24 year movingwindow trends in Z500-PC1, both for individual CMIP6 models and ERA5, are shown in figures 3(c) and (d), for the periods starting from 1980 and 1851, respectively. Most of the trends are negative, with a mean of about −0.5σ for the ERA5 during the most recent decades (figure 3(c)). The mean trends in most of the CMIP6 models are also negative, but are much weaker than the majority of the trends in ERA5, and the ensemble spread is large ( figure 3(c)).
For the entire historical period, externallyforced responses in the models are negligible, and the ensemble range is larger than that for the period after 1980 ( figure 3(d)). The recent rapid decadal warming thus seems on average to have forced mostly negative decadal trends in Z500-PC1, but has not systematically led to unprecedented trends. For CanESM5_p1/p2, IPSL-CM6A-LR, and MIROC6, some of the externally-forced ensemblemean decadal trends since 1980 are unprecedented, but not for NorCPM1 nor for CNRM-CM6-1 ( figure S2(c)). On the other hand, one standard deviation of the ensemble spread is not exceptional after 1980, except for NorCPM1 ( figure S2(d)), and this is also the case for internally-generated maximum/minimum trends across ensemble members (figure S3). In short, most models agree on the role of rapid decadal warming in forcing a moderate negative decadal trend in Z500-PC1, which leads to the skewness of the trend distribution (figures 3(b) and (c)), and on its marginal importance for modulating internally-generated decadal variability in Z500-PC1.

Impacts of internally-generated decadal trends in Z500-PC1 on surface climate and SWE in observations/reanalysis and CMIP6 models
As discussed above, the recent decadal trends in Z500-PC1 are largely internally-generated. In this section, the impacts of these internally-generated trends in Z500-PC1 on the Eurasian surface climate, including snow variability, are studied.
The modelled, externally-forced decadal trends in SAT exhibit consistent warming trends, with generally increasing magnitudes over the three regions (figures 4(a), (c) and (e)), largely consistent with the global and mid-latitude warming trends (figures S2(a) and (b)). There is a corresponding, externally-forced decadal snowfall increase in both Russia (figure 4(c)) and East Eurasia (figure 4(e)), and decreases in Europe ( figure 4(a)). In terms of snow cover, negative decadal SCF trends are generally modelled in the three regions, consistent with decadal warming trends. The CanESM5_p1/p2 models (red/ green dots) have particularly strong externally-forced responses in both SAT and snowfall over the three regions, and in SCF across both Europe and Russia. In response to the combined influences of SAT, snowfall, and SCF, SWE decadal trends are consistent with snowfall trends in Europe and East Eurasia and are opposite to snowfall trends in Russia. The decreasing trends in SCF may contribute to shorter snow accumulation periods, which overcompensate for the increasing trends in snowfall in Russia.
The externally-forced decadal trends in the three regions are generally either opposite to, or much weaker (SWE trends in Russia) than those seen in reanalysis/observations (figures 4(a), (c) and (e); blue triangles). In fact, the temporary slowdown of the decadal warming in the globe (figure S2(a); green line) and particularly midlatitudes (figure S2(b); green line) during recent decades contrasts sharply with the externally-forced decadal trends. The snowfall trends in East Eurasia are likely determined largely by external forcings, because they lie approximately in the middle of the model ensemble means (figure 4(e)). Most of the other observed decadal trends are likely not externally forced. Can these then be attributed to internally-generated variability in Z500-PC1? A comparison of the reconstructed decadal trends using Z500-PC1 (figures 4(a), (c) and (e); red triangles; see section 2.2 for details) and those from reanalysis/observations (figures 4(a), (c) and (e); green triangles), suggests that Z500-PC1 can largely explain these internally-generated decadal trends. The exception is the snowfall trends in East Eurasia, whose PC1 reconstructions are much weaker than most of the model ensemble means, as well as the reanalysis/observations for about half of the 24 year window periods. However, snowfall trends based on PC1 reconstructions seem to be important in some of the later 24 year window periods. The reconstructed SWE trends in the same region are comparable to those in CanESM5_p1 and MIROC6, and are much weaker than the reanalysis/observations. In both cases, the internally-generated variability in Z500-PC1 thus does not appear to play a dominant role over the external forcings. The decadal trends in both snowfall and SWE in East Eurasia are driven by a combination of externally and internally generated trends, without a dominant factor between the two.
The role of Z500-PC1 in driving internallygenerated decadal trends is further supported by evaluating the link between ensemble spread in the decadal trends and in Z500-PC1 within the large ensemble simulations (figures 4(b), (d) and (f); see section 2.2 for estimating the contribution of Z500-PC1 to the decadal trends in various variables). For decadal trends in Europe, five of the models (except CanESM5_p1) agree that internallygenerated decadal trends in Z500-PC1 exert significant impacts on the decadal trends, in at least some of the variables ( figure 4(b)). The MIROC6 model ( figure 4(b); green stars) exhibits the most consistent impacts of internally-generated decadal trends in Z500-PC1 on the four variables. A similar pattern is seen for decadal trends in Russia, where significant impacts on snowfall are observed in all the models ( figure 4(d)). The impacts on SCF trends in Russia are weak, because of the extremely weak SCF trends in most of the models. For decadal trends in East Eurasia, the impacts of internally-generated decadal trends in Z500-PC1 on SWE trends are virtually absent in all the models (figure 4(f)). Most of the models also do not exhibit consistent significant impacts on snowfall trends. These findings year-moving windows starting from 1980 until 1991 for regional-mean variables over (a) Europe, (c) Russia and (e) East Eurasia. Trends in SAT, snowfall, SCF and SWE are computed from reanalysis/observations (green triangles), model ensemble-means (dots and stars) and reconstruction by Z500-PC1 (red triangles). See text for details of the reconstruction of trends. (Right) Similar to left panels but for standard deviations of trends across the ensemble members over 24 year-moving windows, obtained through regression against Z500-PC1, for (b) Europe, (d) Russia and (f) East Eurasia. See text for details of obtaining trends of ensemble members through regression. Only trends computed from regression coefficients significant at the 5% level are shown. No SWE data are available for the NorCPM1 model. suggest that internally-generated decadal trends in Z500-PC1 have not dominated the decadal trends in either snowfall or SWE in recent decades over East Eurasia. This is consistent with the analysis using the reanalysis/observations. The dynamic and thermodynamic processes driven by external forcings may have strongly modulated those forced linked to the internally-generated Z500-PC1 trends in terms of driving snowfall, leading to the impacts of internallygenerated Z500-PC1 trends on snowfall and SWE in East Eurasia being relatively weak. On the other hand, significant and consistent impacts of internallygenerated decadal trends in Z500-PC1 on SAT and SCF are seen in most of the models, particularly for SAT trends. A key finding is that large-ensemble simulations do suggest that internally-generated decadal trends in Z500-PC1 have exerted strong impacts on winter climate variability during recent decades in the three regions studied, consistent with what the reanalysis/observations reveal.

Discussion
The possible impacts of AA and Arctic sea-ice loss on the midlatitude climate variability and atmospheric circulation have been hotly debated in the literature (see references in section 1). The 24 year moving decadal trends in Z500-PC1 correlate mostly weakly with those in global-mean SAT, Arctic seaice, and Arctic SAT across the model ensembles ( figure S4). The only exceptions are for Arctic sea-ice (SAT) for CanESM5-p1/p2 (IPSL-CM6A-LR). However, without concurrent significant correlations for both Arctic sea-ice and Arctic SAT appearing in the same model, it is hard to argue for an impact of Arctic sea-ice on Z500-PC1. Therefore, we find no evidence that internally-generated variability in the recent decadal trends in Z500-PC1 is significantly driven by the global-mean SAT variability and Arctic sea-ice/SAT. Atmosphere-ocean interactions may instead play an important role in driving the recent decadal trends in Z500-PC1. This also suggests that Z500-PC1 may have contributed to the recent decadal climate variability in Eurasia including the cooling trends largely independently from the Arctic sea-ice and Arctic warming. This is consistent with the results based on reanalysis/observations in Ye and Messori (2020). This also highlights that the recent decadal climate variability over Eurasia is not necessarily directly tied to AA/Arctic sea-ice loss. Our results are aligned with other studies that found no evidence of Arctic sea-ice driving of the recent Eurasian winter cooling trends (e.g. McCusker et al 2016, Sun et al 2016, Ogawa et al 2018, Blackport et al 2019. The dynamical leading mode in this study as mentioned in the introduction reflects the downstream impacts of the NAO on the atmospheric circulation in Northern Eurasia. Although this leading mode is largely independent of the Arctic sea-ice loss and warming, its interactions with the Arctic sea-ice and AA, and whether these are competing mechanisms, are interesting topics for future studies. For example, The Z500-PC1 has been shown to drive a large part of the warming trends in Northern Canada and Greenland for the period 1990-2010 (Ye and Messori 2020).

Conclusions
The leading mode of wintertime atmospheric circulation variability over the North Atlantic-Northern Eurasia sector (whose timeseries we designate Z500-PC1) exhibits robust influences on climate variability in Northern Eurasia. Analysis of CMIP6 model output suggests the recent decadal trends in Z500-PC1 are most likely due to internallygenerated climate variability. Internally-generated decadal trends in Z500-PC1 are generally found to have driven a significant portion of the recent decadal trends in variables such as SAT, snowfall, SWE, and SCF in Eurasia. External forcing is, on the other hand, not a major contributing factor in most cases, except for snowfall and SWE in East Eurasia wherein it is arguably equally important as internally-generated climate variability. Therefore, in a decadal perspective, internally-generated variability in Z500-PC1 and the associated impacts on Eurasian climate variability may oppose or amplify the externally-forced responses under strong background warming. This is a significant source of uncertainty and a challenge for future decadal climate projections in Eurasia. A better understanding of this leading mode of variability will be beneficial to both scientific studies and climate service for the region. Another important source of uncertainty in future decadal climate projections for Eurasia is associated with the model discrepancies in the linkages of Z500-PC1 to decadal trends in surface climate. Discrepancies in the sensitivity to external forcing, detailed atmospheric pathways for downstream impacts, rain/snow partitioning, land-atmosphere interactions, and land surface schemes are all possible causes. For example, the snow model has only one layer in CanESM5_p1/p1 (Swart et al 2019), while it has 12 layers in CNRM-CM6-1 (Voldoire et al 2019). Reducing these discrepancies would mitigate model discrepancies in capturing the linkages of Z500-PC1 to the surface climate in Eurasia.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).