Impact of Eurasian autumn snow on the winter North Atlantic 1 Oscillation in seasonal forecasts of the 20 th century.

Abstract


Introduction
As the leading climate variability pattern affecting winter climate over Europe, the North Atlantic Oscillation (NAO) has been extensively studied over the last decades (Wanner et al., 2001;Hurrell and Deser, 2010;Moore and Renfrew, 2012;Deser et al., 2017).The NAO state strongly impacts the hydroclimate as well as the ecological and socioeconomic conditions over major population clusters of Europe and North America.In its positive state, the NAO projects onto strong pressure gradients over the North Atlantic, strong westerly winds and mild but wet conditions for Central Europe.A negative winter NAO is connected to a southwardly displaced Atlantic jet stream, weaker westerlies and cold, dry conditions for Central Europe.The NAO also shows a distinct quadrupole signature in surface temperature straddling the Atlantic, with two opposite poles over northern Europe and Greenland /Labrador and an opposite pair further south over southern Europe/North Africa and the US East Coast.
Recent cases of extreme negative NAO states (Wang and Chen, 2010;Lü et al., 2020), including the winter 2020/2021, coincided with several extreme weather events across the Northern Hemisphere, including cold air outbreaks with record snowfall at locations over Southern and Northern Europe, as well as eastern parts of Canada and the United States (Blunden and Boyer, 2021).
The causal chain behind the snow impact is hypothesized as follows: due to the radiative and thermodynamical properties of snow (Cohen and Rind 1991;Vavrus 2007;Dutra et al., 2011;Thackeray et al., 2019), a thicker snowpack is associated with coherent surface cooling.Cohen et al.
(2007; see also Cohen et al., 2014;Henderson et al., 2018 for reviews) proposed a multi-step mechanism whereby this surface cooling leads to raised isentropic surfaces, triggering increased Rossby wave activity propagating upward and being absorbed in the stratosphere, warming it and subsequently weakening the polar vortex.The negative stratospheric Northern Annular Mode signal eventually propagates down into the troposphere and to the surface where it projects onto a negative NAO.
Investigating the robustness of this mechanism is challenged by several elements.Observational studies analyzing statistical links are restricted by the relatively short length (a few decades) of comprehensive and complete snow cover observations.Using long-term reanalyses, recent studies showed substantial Formatted: Font colour: Auto Formatted: Font: 11 pt, Font colour: Auto, English (US) Formatted: Font colour: Auto Deleted: " year and a length of 4 months has been used in several studies on the predictability of the NAO and other climate patterns (e.g., O'Reilly et al., 2017;Parker et al., 2019;Weisheimer et al., 2019;Weisheimer et al., 2020;O'Reilly et al., 2020).To investigate the influence of land surface conditions, in this case snow cover, on the evolution of the atmospheric state throughout the season, we use a novel, 21-member twin set of the ASF-20C forecasts with perturbed initial land conditions.This dataset was used as a pilot experiment in the context of the Land Surface, Snow and Soil moisture Model Intercomparison Program LS3MIP ( Van den Hurk et al., 2016), aimed at reproducing land surface potential predictability experiments as described by Dirmeyer et al. (2013).We aim to address the question of causality, pathway, stationarity and seasonal evolution of the proposed mechanism of the snow-stratosphere-troposphere linkage over decadal to centennial time scales.This paper is organized as follows.Sect. 2 describes the data and methods used.In Sect.3, we show winter evolution of climate anomalies for the different initialization runs and contrast them with observed anomalies.The results are discussed in Sect. 4 and finally summarized in Sect. 5.

a. Climate reanalysis and reconstruction
We use the European Centre for Medium-Range Weather Forecasts (ECMWF) product ERA-20C (ERA20C;Poli et al., 2016) to investigate pre-conditions and the initialization of the seasonal predictions, to compute the DJF NAO index as well as to create a Eurasian snow dipole index.ERA-20C only assimilates surface pressure and marine wind observations, with sea surface temperature (SSTs) boundary conditions taken from the HadISST2.1.0.0 datasets (Rayner et al., 2003).ERA-20C was found to represent interannual snow variations over Eurasia remarkably well.For an in-depth discussion of its performance and the technical details concerning snow computation, see Wegmann et al., (2017b).Due to the aforementioned statistical impact for the winter NAO evolution, we focus on the November Eurasian snow dipole index as a predictor for the following NAO state (Gastineau et al., 2017;Han and Sun 2018;Santolaria-Otín et al., 2021).Following Han and Sun (2018) who explicitly selected western and eastern domains because of the high co-variance with DJF NAO, we calculate the index over the period 1901-2010 by averaging snow depths over the western domain (30°-60°E, 48°N-58°N) and the eastern domain (80°-130°E, 40°-56°N), eventually subtracting the western domain from the eastern domain to derive the west-east snow depth gradient.Hence, a positive snow index indicates higher snow depths in the eastern domain and a positive longitudinal snow depth gradient.The index is normalized and linearly detrended with respect to the overall time period.To comply with the initialization date of 1 st of November for the seasonal hindcasts, we calculate the index for 1 st of November instead of November mean snow (Figure 1a).Even though Han and Sun (2018) calculated the dipole index using snow cover, we used snow depth since ERA-20C provides snow depth as the actual prognostic variable.We hence refrained from using empirical rules to convert snow depth to snow cover.We found the index based on snow depth to be virtually identical (also see Supplementary Figure S1) to the index using snow cover (see also Wegmann et al., 2020 for more insights).
To compute the winter NAO index, we normalize the first Empirical Orthogonal Function of ERA-20C DJF sea level pressure (SLP) for the region (90°-50°E, 20°-80°N).We use the same approach to calculate the NAO DJF index based on the seasonal hindcasts and compare those with the reconstructed, independent DJF NAO index by Jones et al., (1997) from the Climate Research Unit (CRU).

b. Seasonal prediction experiments
Additionally, we use atmospheric seasonal retrospective hindcasts covering the 110-year period 1901-2010 of ERA-20C with 51 ensemble members of the ASF-20C hindcasts (hereafter ASF-20C CTL) (Weisheimer et al., 2017).The atmospheric model used for the 4-month hindcasts is the ECMWF Integrated Forecast System Version CY41R1 and is initialized at four start dates per year (1 st of Feb, May, Aug and Nov) with ERA20C land and atmospheric conditions.It uses the same lower boundary conditions for SST and sea ice as ERA-20C.Here, we only use hindcasts initialized on the 1 st of November.The horizontal spectral resolution of the model of T255 is similar to ECMWF's previous Deleted: -operational system System 4 (Molteni et al., 2011) and corresponds to a grid length of approximately 80 km.The model has 91 vertical levels and a top at 0.01 hPa.The ensemble has been created by perturbing each member through the stochastic physics schemes to represent model uncertainties in a similar way as the aforementioned System 4.
To investigate the impact of Eurasian snow depth we use an additional set of perturbed hindcasts, based on a 21-member subset of the ASF-20C CTL experiment (hereafter, the "Experiment" or ASF-20C EXP).Each member run is initialized with different land surface conditions, sampled from the neighbouring 20 years.For example, the range of land surface conditions for the 21-member ensemble forecast initialised on 1 st of November 1950 spans the land surface conditions of the years 1940-1960: member 01 is initialized with the land surface conditions of 1940, member 02 of 1941, member 03 of 1942 and so forth.For the beginning and ending ten years of the hindcast dataset, the land surface conditions are sampled from the closest 21 neighbouring years within the dataset.Here, land surface conditions include the entire land state, including soil moisture, snow depth and soil temperatures.We argue that for investigating northern hemisphere climate anomalies of the 1 st of November initialization, snow depth has by far the largest impact on atmospheric dynamics compared to soil moisture and soil temperatures, thus allowing us to attribute the differences to snow changes.The main bulk of the experiment data has a monthly resolution, daily data is only available for selected variables and three tropospheric levels.After initialization, oceanic components like SSTs and sea ice are nondynamical and based on observations among all members in ASF-20C CTL and ASF-20C EXP.
Taking advantage of the shuffled initial land conditions of ensemble members in ASF-20C EXP, we subsample members with positive or negative initial Eurasian snow dipole (Figure 2).This conditional sampling approach has been used when testing the sensitivity of extended range forecasts to soil moisture (Koster et al., 2011;van den Hurk et al., 2012) or to snow initial conditions (Li et al., 2019;Garfinkel et al., 2020).For each start date, we can identify those members with positive or negative initial snow dipole indices, corresponding to different years of the shuffled land initialisation.We further proceed with compositing these two selected sets.Due to the decadal variability in the November snow cover, the amount of "high snow dipole members" (positive dipole index) and "low snow dipole members" (negative dipole index) varies throughout the 110 years.There might be periods when a majority of the neighbouring 20 years shows a positive snow dipole index and other periods when a minority does.To avoid this variation of the composited ensemble size across the years, we only use the five ensemble members with the most positive and most negative initial Eurasian snow dipole, creating two ensemble means (each of size N=5), namely a high snow dipole ensemble-mean and low snow dipole ensemble mean, for each winter through the 110-year period.
It should be noted that the absolute magnitude of the ensemble-mean snow differences is still changing from year to year.For example, the most positive snow dipole for the period 1910-1930 might be lower than in the time window 1980-2000, and the same applies for negative dipole indices.Due to the definition of the ASF-20C EXP, this setup is unavoidable, but it also allows for realistic magnitudes of snow forcings and for incorporating a realistic natural variability into the experiment.The (5-member) ensemble-mean difference (Figure 3a) displays a snow depth increase of 1-2 cm over Central and Eastern Siberia, together with a 0.2-1 cm snow depth decrease over Western Russia, as expected from the snow dipole definition.Concomitant negative anomalies (1-2 cm snow depth) nevertheless extend outside of the dipole definition domains to more northern latitudes, e.g., over Western Russia and the Russian Far East, or over the coastal mountain ranges of the North American Pacific Northwest.Note that the two domains forming the dipole are in snow transition zones, where the snow cover is rare on November 1 st and shows some variability (Figure 3b-c).The dipole positive phase corresponds to anomalously high snow depths over Eastern Eurasia, where the ERA20C snow depth climatology indicates a few centimeters of snow.It also corresponds to anomalously low snow over the west of Russia, in regions with no to rare snow cover in the ERA20C November 1 st climatology.The eastern domain partly covers the Mongolian Plateau region which was shown to exert a strong impact of the wintertime wave fluxes in the stratosphere (White et al., 2017).
If not stated otherwise we compute differences between the 5-member ensemble means of the "high snow dipole" and "low snow dipole" in ASF-20C-EXP as well as differences of each ensemble mean relative to the ensemble mean of ASF-20C CTL.We compute significance using a two-sided Student´s t-test.

a. DJF NAO comparison
Figure 1b shows the time series of the normalized reconstructed (i.e., based on station data), reanalysed and predicted winter NAO state for the period 1901-2010.Unsurprisingly, the ensemble means of the ASF-20C CTL and ASF-20C EXP hindcasts show reduced temporal variance compared to the observation-based NAO datasets.However, single realizations and member spread of the CTL and EXP runs cover the whole range of variability displayed by the observation-based product (see also Supplementary Figure S2).
The correlation between the ERA-20C and CRU NAO index is 0.83, indicating that the EOF approach is a good approximation of the station-based index.It should be noted that the DJF average has a higher correlation between hindcasts and reanalyses than the individual months within the season (see Supplementary Table S1).
The ASF-20C CTL ensemble mean DJF hindcasts achieves an overall correlation of 0.33 with the CRU NAO reconstruction for the complete time period, with ASF-20C EXP having a nearly identical correlation (0.34).This near-identical correlation is expected given that the land state perturbations Contrasting the (initial) high-dipole and low-dipole composites constructed from the ASF-20C EXP ensemble, we see decadal variability in the difference of winter-mean NAO (Figure 1c&d).The first two decades of the 20 th -century are characterized by rather strong negative NAO responses to a strong positive snow dipole.This is followed by two decades spanning the early twentieth century Arctic warming (Polyakov et al., 2003), which shows the opposite response: A strongly positive west-east snow depth gradient, as depicted in Figure 3, leads to more positive NAO-like states compared to a weak west-east snow depth gradient.After several decades with changing responses to the snow dipole between the two ensembles, eventually the 21 st -century starts with a weak negative NAO response to a strong positive snow dipole.Averaged over the whole period, the high snow dipole ensemble shows a slightly stronger negative NAO response: 51 cases of positive NAO response versus 59 cases of negative NAO response.For more extreme NAO states (1 SD exceedance), the difference is more pronounced: 18 versus 29, and for 2 SD exceedance, 2 versus 9. Possible reasons for the decadal response to the snow forcing will be considered in the discussion section.

b. Regression analysis
Previous studies showed that regressing observed boreal winter zonal-mean temperature and zonal wind anomalies onto an observed Eurasian autumn snow index reveals a significant stratospheric warming and slow-down of the polar vortex starting in November, migrating down towards the tropopause until February (Wegmann et al., 2020).A similar relation between Eurasian snow and the polar stratosphere can be found in the dataset used here.
Figure 4 shows a strongly reduced polar vortex for the ERA20C autumn to winter climate anomalies regressed on the November snow dipole index.The zonal wind anomalies in the troposphere highlight a weakened polar jet and an increased subtropical jet, especially in January and February.The concurrent polar stratospheric warming signal moves towards the upper troposphere throughout the winter months, with peak warming at around 100 hPa in February.
Spatially, pressure anomalies regressed onto the November snow dipole index reveal that the geographical center of the stratospheric warming is located over the Canadian Arctic (Figure 5).

c. Spatial anomalies in the experiment
In the following paragraphs we investigate the spatial differences in the atmospheric response associated with the high and low snow dipole ensemble means of ASF-20C EXP, focusing on the initial response in December as well as the average DJF response, as the November response is not yet significant for almost all climate variables.
Figure 6a&b shows stratospheric geopotential heights anomalies at 10 hPa.In December, a significant negative anomaly formed above Eurasia, corresponding to a polar vortex displacement toward the Eurasian sector and a high over Alaska (albeit not significant), as commonly found during stratospheric warming events.Over the course of the winter, this pattern develops into increased geopotential heights over the Arctic with significantly reduced geopotential heights over the extratropics, albeit only significant over Southern Europe and the Caucasus.
To better understand the wave activity flux into the stratosphere, we investigated the meridional eddy heat flux at 100 hPa, which is proportional to the vertical component of the wave activity flux (Figure 6c): it highlights a wave train of circumpolar anomalies in December (hence, following the surface signal forcing in November) with significant positive anomalies over the Ural mountains, eastern North Pacific and the European part of the North Atlantic and negative anomalies over Central and Northern Europe and along the North American Pacific coast.The average DJF response highlights a circumpolar wave-train but shows significant anomalies only for the increased northward heat flux over the northern North Atlantic.
Tropospheric circulation anomalies are depicted for geopotential heights at 500hPa in Figure 6e&f.In December, a strong positive anomaly is located over the Barents-Kara Sea sector, with significantly negative anomalies up-and downstream.A second region of positive anomalies emerges at the Canadian Atlantic coast.Both regions match the significant positive anomalies in the 100 hPa heat flux well.The averaged DJF anomalies highlight a negative mid-tropospheric NAO signal with significantly increased geopotential heights above Greenland and Iceland.
Sea level pressure anomalies largely mirror the 500 hPa geopotential height anomalies.The averaged DJF pattern only shows significant positive anomalies over the northern North Atlantic, but still projects onto a meridional pressure gradient characteristic of a negative NAO phase (Fig. 6h).It is important to note, that the absolute difference is rather small compared to interannual SLP variability.Anomalies between the two ensemble-means are less than 1 hPa.Even though this number can be assumed to be Deleted: sampling random smaller than in observational datasets due to the ensemble averaging process, it only constitutes about 15% of the average 1901-2010 DJF SLP standard deviation over the Euro-Atlantic sector.
Due to its large variability, composites of the near-surface temperature are largely non-significant (Figure 6i&j).Yet, in December a clear cooling signal emerges over Central and Eastern Eurasia, as expected from the location of the positive snow anomalies at the time of forecast initialization.At the same time, eastern North America and south-eastern Europe show significant positive temperature anomalies, a result of northward heat advection at the eastern flanks of low-pressure anomalies (Figure 6g).Averaged DJF 2m temperatures are significant only for Greenland and Eastern Eurasia, with the cooling over the latter a direct result of the persistence of the anomalously high initial snowpack.

d. Vertical anomalies in the experiment
To get a better understanding on how the different land initial conditions impact the vertical distribution of temperature and zonal wind, Figure 7 shows meridional cross-section of the zonal-mean anomalies of zonal wind and temperatures from November to February.
While November anomalies (Figure 7) are overall insignificant, a strong snow dipole is associated with an increased polar vortex and cooler stratosphere.In December, zonal wind anomalies are indicative of the tropospheric subtropical jet shifted northward concurrent with a weak Arctic surface warming.
Changes are substantial in January, when the stratospheric polar vortex is significantly weakened, with a slight increase in westerlies in the mid-troposphere.The corresponding temperature anomalies show a widespread stratospheric warming and negative anomalies in the lower Arctic troposphere.Eventually in February, the slow-down of westerlies is predicted to reach all the way down from the stratosphere into the troposphere.On the southern flank of these negative zonal wind anomalies, westerly winds are increasing, especially so in the stratosphere.The stratospheric warming signal migrates downwards to the lower stratosphere and tropopause layer.As the warming has migrated down, a stratospheric cooling is occurring aloft.
As a further confirmation, polar cap heights (Supplementary Figure S3) reveal a development of positive anomalies from the surface in December up to the stratosphere in January, migrating back to the troposphere in February.Note that the development of these anomalies is delayed compared to the one shown in the ERA20-C reanalyses (compare Figures 4 and 7) since initial atmospheric conditions are identical in the perturbed ensemble members.

e. Daily evolution of anomalies in the experiment
To investigate the temporal evolution and importance of tropospheric anomalies, Figure 8  November going into December and is preceding the development of the North Atlantic ridge, which is the main component of the negative NAO-like feature in our results.It should also be noted that the absence of meaningful anomalies during the first ten days of the composite difference again reflects the identical tropospheric anomalies arising from the pre-conditions.The anomalies generated by the end of November do indeed arise from the impact of snow cover differences and snow-atmosphere feedbacks.

f. Non-linearities in the snow forcing impact
Two distinct non-linearities need to be considered.First, a non-linearity in the physical snow feedback: adding a few centimetres of snow in a snow-covered region will not change the radiative and thermodynamic properties of the already snow-covered land surface substantially (due to a saturation effect) but, by contrast, removing a few centimetres of snow might remove the snow layer altogether, changing drastically the albedo and thermodynamics of the surface-atmosphere boundary.This nonlinearity may be important for the Rossby wave generation as air flows over the uplifted isentropes above the snow-covered area.The non-linear effect of snow cover saturation and the impact of the relative magnitude of regional surface cooling in our experiments is addressed by Figure 9.In years when the high minus low dipole anomalies preceeded a negative NAO anomaly (see Figure 1d for indication of years), the December cooling anomaly over Eastern Eurasia is much stronger than for the opposite case when it preceeded a positive NAO anomaly.Concurrently, the formation of a Ural ridge anomaly is much more pronounced, flanked by troughs up and downstream, with positive eddy heat fluxes into the stratosphere over the Barents-Kara Sea and polar stratospheric warming.This supports the notion that adding an absolute amount of snow (in either of the two longitudinal domains) is not sufficient for the causal chain to be triggered.Instead, it is a large (in magnitude and extent) relative surface temperature impact of the additional snow that triggers the initial anomalous Rossby wave generation part of the hypothesized causal chain.
A second non-linearity is the asymmetrical role of the eastern and western domains of the snow dipole.
Our subsampling of the ASF-20C EXP simulation allows to estimate the respective roles of these two domains.Interestingly, the difference between the low snow dipole ensemble mean and the CTL ensemble mean for DJF sea level pressure (Figure 10) reveals a much stronger response to a negative snow dipole (i.e., with high snow depths over Western Russia and low snow depths over Eastern Eurasia) than to the positive snow dipole (i.e., with high snow depths over Eastern Eurasia and low snow depths over Western Russia).The reason behind this is a combination of study design, the nonlinearity of snow cover and the snow climatology of Eurasia.In Figure 10a, we compare the effect of a very non-climatological snow depth gradient to the impact of a climatological snow depth gradient and as such get a very strong response in SLP anomalies.This comparison equalizes or reduces the snow depth gradient and as a result very zonal flow occurs over high-latitudes.In Figure 10b, we compare the effect of a slightly increased (to the climatology) snow depth gradient to the impact of a

559
generation by much and we do not find positive the 560 geopotential height anomalies in its vicinity.¶ 14 member" approach resulting in different magnitudes and shorter response times, the seasonal evolution of the ASF-20C EXP high-low snow dipole anomalies similarly indicates a weakened polar vortex.It also supports the notion of a surface cooling over the Eastern domain anchoring a Ural ridge anomaly on its western flank in December (Figure 6e).This Ural ridge triggers an increased northward heat flux in the lower stratosphere, thereby reducing the polar vortex strength and increasing polar stratospheric temperatures.In January and February, the signal moves downwards into the troposphere where it evolves into a negative NAO anomaly.In general, these results agree with the framework proposed by Cohen et al., (2007) and the experiments with the ECMWF seasonal prediction model by Orsolini et al. (2016).The physical causal chain in our experiment is also in line with recent model studies investigating the impact of Eurasian snow on stratospheric warmings and possible surface climate anomalies (Cohen et al., 2021).However, it should be highlighted that the absolute ensemble-mean, time-average SLP signal, diagnosed as the conditional composite difference in ASF-20C EXP, is very small, less than 1 hPa.As mentioned before, this represents only a small fraction of the interannual SLP variability in the Euro-Atlantic region.Nevertheless, for single realisations of winter forecasts, this impact can be much higher.By design, we excluded the impact of sea ice on the NAO evolution, since SSTs and sea ice stay the same through time in all EXP members.We checked for significant tropical precipitation patterns in the high minus low anomalies for November and December as potential covariates in driving an DJF NAO signal but found no coherent significance across tropical latitudes.As such, we exclude tropical rainfall as first order driver behind the EXP NAO response.
The role of the Ural ridge in the snow cover à NAO causal chain has been discussed and analysed in several recent studies (Peings 2019;Santolaria-Otín et al., 2020).Here we find that the Ural ridge is a pre-condition of predicted negative NAO winters in ASF-20 CTL (Supplementary Figure S6), together with a cold 2m temperature anomaly in Eastern Russia and a cold stratospheric polar vortex displaced over Eurasia, downstream of the Ural ridge.However, these initial conditions are subtracted out in the ASF-20C EXP high--low snow dipole composite difference, and we find that the composite difference indicates a reinforced Ural ridge (Figure 6e).We find the mid-troposphere Ural ridge is reinforced only at the end of November going into December, which precedes the formation of a North Atlantic ridge that prevails until February (Figure 8).This result indicates that the snowpack does indeed play a feedback role (see also Orsolini et al., 2016).Thus, we propose that the relation between the Ural ridge and Eurasian snow cover consists of a mutual interaction: the circulation anomaly associated to a preexisting Ural ridge shovels cold polar air southwards along its eastern flank, allowing for an deeper snowpack to form over Eastern Eurasia (eastern domain of the dipole).In addition to this process (Figure 9c), our analysis reveals that the snow cover anomaly reinforces the Ural ridge, allowing for increased wave flux into the stratosphere.This location of a tropospheric ridge interferes constructively with climatological stationary wave-1 and wave-2 patterns (Garfinkel et al., 2010) and seems to be key for a skilled forecast of the polar winter stratosphere (Portal et al., 2021).

Deleted: T
Furthermore, the high minus low composite highlights a subpolar North Pacific surface and midtropospheric low-pressure anomaly that appears first in December and remains throughout all of DJF (Figures 6,f,g,i and 10).The generation of this circulation feature was pointed out by previous studies (Orsolini and Kvamstø, 2009;Garfinkel et al., 2010;Garfinkel et al., 2020), and has been attributed to an enhanced vertical propagation of Rossby waves into the stratosphere and horizontal downstream of the cooled Eurasian land mass.
Subsampling of the experimental multi-decadal historical hindcasts (ASF-20C EXP) highlighted an interdecadal variability and non-stationarity of the snow dipole impact, despite the cancelling out of common boundary forcings such as SSTs in the composite difference.The configuration of our experiment does not allow to explain this behaviour completely; however, we can address some possible reasons.A potential influence on the decadal variability of the snow cover impact might be the precursory climate system state, promoting or counteracting the tendency for the (perturbed) snow forcing towards a given NAO state.
Surprisingly, the positive snow dipole forcing tends to favour a negative NAO signal when the climate system is "tuned" for a positive winter NAO in ERA20C, for example when high ERA20C Barents-Kara sea ice extent and La Niña SST conditions prevail (Supplementary Figure S7).This supports the idea of a clear and strong snow impact when the relative cooling anomaly in Eastern Eurasia is relatively strong and the climate state is preconditioned to a rather positive NAO-like condition.This might explain the strong positive NAO anomaly during the early twentieth century Arctic warming in Figure 1d: the period 1920-1940 was characterized by a strong positive mid-tropospheric high anomaly from Northern Europe to East Siberia (Wegmann et al., 2017a).We find that the 500 hPa anomalies between high and low snow composites show only a weak to non-existing Ural ridge for the period 1921-1940, when compared to e.g.1991-2010 (Supplementary Figure S8).On the contrary, increased snow in an already snow covered Eastern Eurasia will not provide the same response as the pre-existing anomalies favoured by other background conditions.Rather, strong non-linearities seem to occur, which is reasonable given the non-linear thermodynamic and radiative impacts of a deeper snowpack.
On that note, we find that the relative magnitude of regional cooling compared to the existing climate state in our experiments is of crucial importance.In years when the high-low snow dipole anomalies preceeded a negative NAO, the December cooling anomaly over Eastern Eurasia is much stronger than for the opposite case (Fig. 9).Moreover, we found that in our model experiment a negative snow dipole forcing leading to a positive NAO signal has a much larger relative impact compared to the positive snow dipole resulting in a negative NAO signal, which is due to the much stronger changes in the Earth System we impose with the negative snow dipole ensemble.Moreover, due to the Eurasian snow climatology, a similar level of snow depth variability in the western domain will have higher impacts on the tropospheric variability.In our experimental setup, a weak snow depth gradient from west to east allows for a rather zonal circulation in the following months, with no subsequent stratospheric warming Deleted: experiment dataset (created by Bart van den Hurk) can be made available on the ECMWF MARS system upon request.(a,c,e,g,j), and DJF (b,d,f,i,k): a&b) 10 hPa geopotential heights, c&d) 100 hPa meridional eddy heat flux, e&f) 500 hPa geopotential heights, g&h) sea level pressure and i&j) 2m temperature.Stippled areas represent 90% significance.For monthly anomalies see Supplementary Font colour: Auto, Not Highlight Deleted: The location of the sub-sampled snow forcing 266 however adds… 267 Deleted: towards across the 21 members are two-sided.Differences between the predicted NAO index of ASF-20C CTL and EXP ensemble means are generally small, with the NAO indices having the same sign during most winters.The correlation between CTL and EXP is 0.8 for the 110-year period.The slightly stronger variability of ASF-20C EXP can partly be attributed to the reduced ensemble size.
Tropospheric pressure differences highlight a strong ridging over Western Russia and the Ural Mountains in December, which subsequently over the course of winter is shifted more towards Greenland and the Northern North Atlantic region, reflecting a negative NAO-like atmospheric state.This state is further supported by negative DJF SLP anomalies over Southern Europe and the Mediterranean region.Downstream of the Eurasian snow signal, a negative SLP anomaly is found over extreme NAO states with 51( the Northern North Pacific.The question remains, if these patterns derived by the regression analysis are a result of co-variability, common climate drivers or causal physical processes. shows daily mean meridional mean 500 hPa GPH anomalies (high minus low snow dipole ensembles) averaged over 60-70°N.The Hovmøller diagram illustrates the Ural ridge developing only at the end of List Paragraph, Numbered + Level: 2 + Numbering Style: a, b, c, … + Start at: 1 + Alignment: Left + Aligned at: 1,9 cm + Indent at: 2,54 cm Formatted: Font: 11 pt, Bold Deleted: (?) 20C EXP simulation allows to estimate the 525 respective roles of these two domains.Interestingly, the 526 difference between the low-snow dipole ensemble mean and 527 the CTL ensemble mean for DJF sea level pressure (Figure 528 11) reveals a much stronger response to a negative snow 529 dipole (i.e., with high snow depths over Western Russia and 530 low snow depths over Eastern Eurasia) than to the positive 531 snow dipole (i.e., with high snow depths over Eastern Eurasia 532 and low snow depths over Western Russia).In other words, 533 an anomalously negative NAO signal in boreal winter in the 534 high-snow minus low-snow dipole anomalies is mostly a 535 result of an anomalously positive NAO signal in the low-536 snow ensemble.The negative dipole corresponds to lower 537 snow depths over the eastern domain (Mongolian Plateau and 538 surroundings areas), consistent with lessened wave fluxes into 539 the stratosphere over this region which is the important 540 orographic driver of climatological upward wave fluxes in 541 winter (White et al., 2017).On this note, additional snow 542 right around the Ural Mountains (a negative dipole index) 543 does not enhance the pre-existing role of the Ural mountains 544 in Rossby wave generation by much in our setup.¶ 545 A possible reason for this non-linear behaviour might be 546 found in the importance of the Rossby wave generation via 547 the Eastern Russia cold air dome over the snow-covered area.548 Based on the snow depth climatology of Eurasia (Figure 3), 549 the negative dipole forcing represents a much more severe 550 disruption to the climatological snow distribution than the 551 positive index forcing is.Removing The former with low 552 snowdepths over Eastern Eurasia seems to favour zonal flow 553 more than adding the latter with high snow depths is 554 favouring meridional flow.Again, the relative impact of the 555 snow forcing is key in this context.Interestingly, additional 556 snow upstream or just right around the Ural mountains 557 (during a negative dipole index) does not seem to enhance the 558 pre-existing role of the Ural mountains in Rossby wave

Figure 1 FormattedFigure 2 :
Figure 1: a) Normalized 1st of November Eurasian snow dipole index for the period 1900-2010 as derived from ERA20C.b) Normalized DJF NAO index in the CRU station-based reconstruction, ERA20C EOF-based index, ASF-20C CTL and ASF-20C EXP EOF-based index.Hollow points

Figure 3 FormattedFigure 4 :
Figure 3: a) Average (1900-2010) 1 st of November snow depth difference between the high-dipole and low-dipole ensemble.b) Average (1900-2010) 1 st of November snow depth.c) Average (1900-2010) 1 st of November snow depth standard deviation.Hatched in green (olive) is the western (eastern) domain of the snow index.All 3 plots are based on ERA20C.

Figure 4 :
Figure 4: Time series of DJF NAO for the period 1901-2010.

.Figure 5 :
Figure 5: ERA20C anomalies of a&b)10 hPa geopotential heights, c&d) 500 hPa geopotential heights and e&f) Sea Level Pressure regressed onto the snow dipole index in November from ERA20C covering 1901-2010 for December and DJF mean.Shading indicates 95% significance level.For monthly anomalies see Supplementary Figure S9.

Figure 7 :
Figure 7: Zonal-mean cross-section of left) zonal wind anomalies and right) temperature anomalies for the period 1901-2010 between high-dipole and low-dipole ASF-20C EXP ensemble means.Shading indicates 90% significance level.

FormattedFigure 8 :
Figure 8: Hovmøller diagram of daily mean predicted 500 hPa geopotential height anomalies for the period 1901-2010 averaged for the latitude band 60°-70°N difference between high-dipole and lowdipole ASF-20C EXP ensemble means.Stippled areas represent 90% significance.Days from NOV 1 st are indicated on y-axis.

FormattedFigure 9 :
Figure 9: Climate anomaly composites of predicted December fields after which a positive snow dipole 1081

Figure 10 :
Figure 10: Mean sea level pressure [Pa] DJF anomalies for the period 1901-2010 between a) lowdipole ASF-20C EXP ensemble mean and ASF-20C CTL ensemble mean (subsampled from 21 CTL members) and b) high-dipole ASF-20C EXP ensemble mean and ASF-20C CTL ensemble mean (subsampled from 21 CTL members).Stippled areas represent 90% significance.c) Represents differences between high-dipole ASF-20C EXP ensemble mean and ASF-20C CTL ensemble mean DJF seasons after positive snow dipole November (see Figure 1a).d) as c) but for DJF seasons after negative snow dipole November.e) as c) but for the low-dipole ASF-20C EXP ensemble mean and f) as d) but for the low-dipole ASF-20C EXP ensemble mean.