A 21st century shift in the mechanisms of the early-winter United States snowfall variability

Snowfall is a critical element of natural disasters to the United States (US) with strong climatic and socioeconomic influences. Meanwhile, snowfall acts as a driving force to the US water supplies for agriculture, drinking water and hydropower. However, so far, what factors influence the US snowfall variations and how these factors change under global warming remain unclear. Here, we found that large-scale influences of the early-winter US snowfall experienced a shift from the Pacific to the Atlantic side around 2000, through observational analysis and climate model simulations. The Pacific/North American pattern was identified as a dominant driver of the early-winter US snowfall before 2000, but its impact became much weaker in the 21st century as its associated western North American cell shifted northward away from the US. Instead, the tropical and subpolar North Atlantic surface temperature has been influencing the early-winter US snowfall variations via teleconnections after 2000. This changed influence of US snowfall around 2000 is demonstrated to be related to the observed global warming pattern since the 1950s. Our study provides new perspectives in understanding large-scale snowfall pattern and variability and its connection to the global warming pattern.


Introduction
Snowfall invariably influences both human society and natural ecosystems with positive or negative impacts (Schmidlin and Kosarik 1999, Smith and Matthews 2015, Stoy et al 2022, Zhu et al 2022, Do et al 2023).In December 2022, a sequence of winter storms hit the majority of the United States (US), which killed at least 87 people and caused billions of dollars in damage (National Centers for Environmental Information 2023).Meanwhile, more than one sixth of the global population relies on the snowmelt water supply (Barnett et al 2005), and the western US has emerged as a hot spot of snow drought in the past four decades (Huning and AghaKouchak 2020).It is therefore of great importance to understand why the US snowfall fluctuates from year to year and how it will change in a warming climate.
The formation of snowfall involves complex physical processes and requires favorable climatic conditions, including moisture convergence and low temperature (Kapnick et al 2014, Yang et al 2019).
Various climate factors, depending on the region of interest, have been proposed to potentially influence the North American snowfall, including the Pacific/North American pattern (PNA) (Serreze et al 1998, Coleman and Rogers 2003, Notaro et al 2006), the Pacific sea surface temperature (SST) variability associated with the El Niño-Southern Oscillation (ENSO) or the related Pacific decadal oscillation (PDO) (Kunkel and Angel 1999, Neal et al 2002, Newman 2007, Sobolowski and Frei 2007, Seager et al 2010),the North Atlantic oscillation (Notaro et al 2006, Ghatak et al 2010, Seager et al 2010), and the Arctic sea ice variability (Liu et al 2012).However, it remains ambiguous whether there is a dominant driver of the US snowfall variations on the continental scale and if so, how it may respond to global warming.These questions will be addressed in this study.We will focus primarily on the early-winter season (i.e.October-December or OND) when the US snowfall typically starts to accumulate and could cause severe societal impacts although the total amount of snowfall may be less than that for winter (Schmidlin and Kosarik 1999, Stoy et al 2022, Yao et al 2023).As we will show later, this is also the season when a prominent teleconnection pattern shift relevant to US snowfall variability is identified.The interannual variations of OND-averaged snowfall are directly linked to the first snowfall date, across the entire US (figure S1), which holds significant implications for individuals and areas affected by early snowfall.

Observational data
For snowfall, we used National Operational Hydrologic Remote Sensing Center (NOHRSC) snowfall dataset provided by the National Oceanic and Atmospheric Administration, and European Centre for Medium-Range Weather Forecasts (ERA5) snowfall (Hersbach et al 2020).NOHRSC uses all the available snowfall observations to produce a best estimate of snowfall characteristics over the US with 1 km spatial resolution and hourly temporal resolution (Carroll et al 2006).While this data might not perform optimally in high-resolution and complex terrain (Gillan et al 2010), our study concentrates on examining snowfall variability across a broad scale throughout the US.The study also incorporates the widely used ERA5 snowfall data (Nouri and Homaee 2021) due to the limited time span of NOHRSC, which commences only from 2008.For the PNA index, we mainly used the PNA index constructed by the modified pointwise method from Climate Prediction Center (CPC).Besides, we also used the PNA index calculated based on the rotated principal component analysis from CPC to verify the result.For SST, we used the Hadley Centre Sea Ice and SST data set (Rayner et al 2003).For surface air temperature, precipitation, water vapor convergence, we used ERA5 dataset.For 500 hPa geopotential height (Z500), we mainly used ERA5 dataset and also used NCEP-NCAR Reanalysis 1 dataset (Kalnay et al 1996) for validation.All data have a monthly resolution, except for ERA5 snowfall with a daily resolution.ERA5 data, NCEP-NCAR Reanalysis 1 data, SST and PNA indices have a temporal coverage for 1959-2022, and NOHRSC covers from October 2008 to December 2022.

Coupled model intercomparison project phase 6 (CMIP6) outputs
Monthly data from 30 CMIP6 models (Eyring et al 2016) (table S1) are used in this research (we use all models with available monthly data of snowfall, 500 hPa geopotential height and SST).In particular, we use each model's 'r1i1p1f1' member (to weigh all models equally), from two experiments: historical (between 1959 and 2014), SSP5-8.5 (between 2015 and2022).

Atmosphere-only general circulation model (AGCM) simulations
In this study, we conduct sensitivity model experiments using the general circulation model (GCM) developed by the National Center for Atmospheric Research (NCAR), the Community Atmosphere Model version 5 (Tilmes et al 2015) to investigate the SST-US snowfall link.We perform atmosphere-only GCM experiments forced by various SST boundary conditions with the radiative forcing fixed at the values in 2000.The horizontal resolution of the atmospheric model is about 1.9 • latitude × 2.5 • longitude (i.e.f19).This spatial resolution may not be able to accurately resolve the snowfall behavior in mountainous regions on a fine scale but should be sufficient to serve our purpose as the main focus of this study is on the large-scale climatic control of snowfall variability across the US.We conduct two pairs (control and forcing scenarios) of large-ensemble AGCM simulations in the Pacific SST forcing experiments and Atlantic SST forcing experiments, respectively.The details of these experiments can be found in section 3. We start the model simulations on September 1st to minimize the influence of initial conditions on the OND snowfall change.For each set of simulations, we conduct 50 ensemble members that start from slightly perturbed atmospheric initial conditions (on the order of 10 −14 K) to isolate the influence of SST forcing on the snowfall change.The total OND response to the SST forcing is the ensemble mean difference in each experiment, and statistical significance is evaluated by the Student's t test of the difference between ensemble means (Wilks 2006).

Definitions of the PNA index in CMIP6
Consistent with the observation, the monthly PNA index is constructed by following the modified pointwise method (Wallace and Gutzler 1981): where Z500 * denotes normalized monthly Z500 anomaly from each CMIP6 model, which has been interpolated to the same horizontal resolution (1 • latitude × 1 • longitude).

Interannual variations of early-winter US snowfall
We analyze the NOHRSC snowfall dataset covering 2008-2022 and find that the recent devastating winter storms mentioned above were in fact part of the widespread, strong snowfall in early winter 2022 across the US, especially over the Sierra Nevada, the Cascade Range, and the Rocky Mountains (figure 1(a)).Accompanied with that, the northeastern US and the central US to the southeast of the Rocky Mountains experienced weak, negative snowfall anomalies.Although the high-resolution and broad-coverage NOHRSC dataset is useful, its limited length of duration precludes long-term statistical analysis.Therefore, we investigated the atmospheric reanalysis, ERA5, and found that it well captured this observed anomalous snowfall pattern (figure 1(b)).ERA5 will be the main dataset for our following analysis, given its longer temporal coverage .
Based on empirical orthogonal function (EOF) analysis, the dominant mode (EOF1; 40% variance explained) of the US OND snowfall in 1959-2022 (figure 1(c)) largely resembles the pattern of 2022 OND snowfall anomaly with a spatial correlation of 0.78 (figure 1(b)).It implies that the early-winter storms in 2022 were not just coincidental or regional extreme events, but rather reflected spatially-coherent climate fluctuations potentially with large-scale influences.Furthermore, we find that the principal component associated with this dominant mode (PC1) is tightly correlated with the averaged OND snowfall in the entire US (r = 0.80; figure 1(d)).Therefore, understanding this dominant pattern of OND US snowfall and its associated temporal variations is relevant for both the snowfall onset in the western US but also the total US snowfall amount on a continental scale.

Dominant role of PNA before 2000
What influences the interannual variability of earlywinter US snowfall?PNA is the dominant mode of winter atmospheric variability over the North America (Wallace andGutzler 1981, Barnston andLivezey 1987) with strong impacts on surface temperature and precipitation variations (Leathers et al 1991).We find that the PC1 associated with the dominant mode of US OND snowfall variability is negatively correlated with the PNA index during 1959-2022 with a correlation of −0.55 (figure 2(a)).Their connection was particularly strong before 2000 (r = −0.68),but surprisingly, their connection disappeared at the beginning of the 21st century (r = −0.03).Our 17 year running correlation analysis further confirms that their negative correlation remained statistically significant since the 1960s until it sharply declined in the early 2000s (figure S2(a)).The reduction of correlation can potentially be due to the weakened variabilities of both the PNA and the US snowfall (i.e.reduced signal-to-noise ratio).To test that, we isolated the years during 1959-2000 that had the normalized values of PNA and PC1 both between −1 and 1, similar to those after 2000.In total, 19 years was selected during 1959-2000, and the correlation coefficient of the PNA index and the PC1 of US snowfall was −0.52, still significantly more negative than that for 2001-2022.It implies that the reduction of PNA-PC1 correlation after 2000 is not solely attributed to the reduced variability.
Before 2000, a typical negative phase of the PNA is characterized by a wave-train structure extending from the North Pacific to the North America, resulting in a cyclonic circulation over the northwestern North America and an anticyclonic circulation over the southeastern North America (figure 2(c)).
The resultant westerly wind anomalies over the western US and the southeasterly wind anomalies in the eastern US carry additional water vapor towards the continent, leading to the convergence of water vapor flux and thus increased precipitation almost across the entire US (figure S3).In the meantime, the cold advection on the west and the warm advection on the east lead to a cold-warm dipole structure across the US (figure S3).Taken together, the combined effects of cold temperature and enhanced precipitation anomalies contribute to a prominent increase in snowfall in the western and north-central US (figure 2(c)) as is also seen in the EOF1 pattern (figure 1 The second pair, named Pacific-Post2000, is similar to the first pair, except that we use the global climatological SST over 2001-2022.Each set contains 50 ensemble members, and the ensemble-mean difference of the US snowfall response between the two sets highlights the net impact of a La Niña-like SST forcing (figure 2(e)).Our modeling results confirm that a La Niña-like SST pattern can indeed induce a negative PNA-like circulation response over the North Pacific-North American sector, a cold west-warm east dipole, and an enhanced precipitation across the US (figure S4), which together lead to an increased snowfall in the western and north-central US (figure 2(e)).This simulated snowfall response pattern looks similar to the observed pattern (figure 1(c), and its magnitude also agrees reasonably well with the observations despite a spread across the ensemble members (red line in figure 2(g)).
So why did the US snowfall-PNA connection break after 2000?Interestingly, the negative phase of PNA during the post-2000 period exhibits a different pattern with the pre-2000 cyclonic circulation originally centered over the US-Canada border now shifting northward to the northern Canada (figure 2(d)).This PNA pattern change around 2000 is a robust feature across different reanalysis datasets and different PNA indices (figure S5).Consequently, the PNA-associated cooling is found only in the far western US and is much weaker than the pre-2000 period (figure S3).Meanwhile, the nationwide moisture convergence found for a negative PNA before 2000 now becomes much weaker after 2000 (figure S3).As a result, the impact of PNA on the US OND snowfall is rather weak during the 21st century (figure 2(d)).Consistent with that, our sensitivity simulations demonstrate that the same Pacific SST anomaly pattern imposed onto the 2001-2022 climatological SST barely affects the US snowfall in the ensemble mean (figures 2(f) and (g)).

Effects of global warming pattern
The simulated distinct impacts of La Niña-like SST on the US snowfall before and after 2000  (figures 2(e) and (f To further confirm that the historical warming pattern may modulate the PNA-US snowfall relation, we next investigate the spread among 30 CMIP6 models that can well simulate both the dominant mode of early-winter US snowfall variability and the PNA pattern (figures S6 and S7).We find that most of the models meeting those two criteria can reasonably capture the negative correlation between the PNA and the dominant mode of snowfall variability during 1959-2022 with a correlation coefficient ranging from −0.24 to −0.79; on average, the correlation is about −0.56, close to the observations (figure 3(c)).For each individual model, the correlation coefficient can vary for the two periods before and after 2000 but mostly within a smaller range than the observations.In terms of multi-model average, the PNA-snowfall correlation keeps almost unchanged for those two periods.These results imply that climate models may have underestimated the influence of internal variability or have misrepresented the global warming pattern, among other possibilities.After regressing the OND mean SST warming magnitude (i.e. 2001-2022 minus 1959-2000) of each model against their PNAsnowfall correlation during 2001-2022, we identify a pattern that somewhat resembles the observed warming pattern with a warmer Northern Hemisphere and a warmer Indo-western Pacific (figure 3(d)).These results, together with our model experiments (figures 2(e)-(g)), indicate that the observed drop of the PNA-snowfall correlation from −0.68 during 1959-2000 to −0.03 during 2001-2022 is likely linked to the observed warming pattern.

Increasingly important role of Atlantic SST after 2000
When the PNA lost its influence on the US snowfall after 2000, what might have taken the turn to drive the snowfall variability instead?As illustrated Our linear regression analysis indicates that after 2000 the tropical and subpolar North Atlantic warming is associated with a strong cyclonic circulation over the western North America (figure 4(d)).The westerly wind anomalies at the southern flank of this cyclonic circulation transport additional water vapor inland and cause the moisture convergence and precipitation increase over the western and the northern US (figure S8).Dominated by the precipitation effect, the snowfall in the region significantly increases as seen in the dominant pattern of snowfall variability (figure 1(c)), despite only a limited area of cooling near the western coast and the northern US (figure S8).Before 2000, however, the cyclonic circulation is displaced to the west coast of the North America and also much weaker in magnitude, limiting its effect on the US snowfall increase (figure 4(c)).
To verify the Atlantic impact on the US snowfall after 2000, we performed another two sets of 50-member AGCM simulations forced by the 2001-2022 climatological SST and the Atlantic SST anomalies (shown in figure 4(a)) added to the climatology, named Atlantic-Pre2000 and Atlantic-Post2000, respectively.We find that the imposed tropical and subpolar North Atlantic warming induces a cyclonic circulation anomaly in the western North America (figure S9) and a snowfall increase in the western US (figure 4(f)), similar to the observations.The detailed mechanisms of the snowfall increase induced by the Atlantic warming differ between the observations and the model simulations due to the exact location of the anomalous cyclonic circulation.More specifically, the snowfall increase is dominated by the precipitation increase as the anomalous southwesterlies transfer more moisture inland in the observations (figure S8) but is dominated by the temperature decrease associated with the anomalous advection of cold air from the north in the model simulations (figure S9).Despite this discrepancy, both the observations and the model simulations exhibit an anomalous cyclonic circulation in the Northeastern Pacific-North America that leads to a snowfall increase over the western US accompanied.When the background climatological SST is taken from 1959-2000, the same tropical and subpolar North Atlantic warming has a much weaker impact on the US snowfall (figure 4(e)).The probability distribution function of western US snowfall is also clearly different for the two sets of simulations with the post-2000 one being closer to observations (figure 4(g)).These results imply that the historical global warming pattern, while weakening the PNA impact, strengthens the Atlantic impact on the US snowfall.

Summary and discussions
Using atmospheric reanalysis, sensitivity model experiments, and CMIP6 historical simulations, we report a dominant mode of early-winter US snowfall and its changing influence from the PNA to Atlantic SST around 2000.The change of influencing factor for the US snowfall is shown to be related to the northward shift of the western North American cell of the OND PNA pattern after 2000.It is noteworthy that this PNA pattern shift is a robust feature that can be seen in different reanalysis, and with different definitions of PNA (figure S5).Although our paper is focused on the US snowfall, this PNA pattern shift presumably may also affect other aspects of the climate system.Why state-of-the-art coupled climate models fail to reproduce this observed PNA pattern shift and whether the PNA pattern shift will continue into the future also need further investigations.In the meantime, if the tropical and subpolar North Atlantic SST variability is intensified under global warming (Yang et al 2021), whether it will keep serving as the dominant influence on the US snowfall in the coming decades remain to be explored.
We propose a potential role of the historical warming pattern in modulating the influencing factor of early-winter US snowfall.We speculate that, after 2000, the weakened PNA impact is related to the Indo-western Pacific warming and the central-eastern Pacific cooling, while the strengthened Atlantic impact is related to the enhanced warming in the North Atlantic (figure 3(a)).Although the factor that changes the US snowfall influence from the Pacific to the Atlantic may be thought to be associated with the transition of the PDO or the Atlantic multidecadal oscillation (AMO) phase, the correlation coefficient between these natural variabilities with US snowfall is not significant.This implies that the 21st century shift in the influencing mechanism of the earlywinter US snowfall cannot be solely attributed to the natural variability of the PDO or AMO and that anthropogenic warming can be critical.Follow-up work on a more definitive attribution is currently underway.It remains under debate to what extent this observed warming pattern is externally forced by radiative agents or internally generated by natural climate variability (Seager et al 2022, Rezaei et al 2023, Tian et al 2023), directly relevant to the future projection of US snowfall variability (O'Gorman 2014).Our study demonstrates that the spatial structure of global ocean warming not only influences, e.g.climate sensitivity (Armour et al 2013), tropical rainfall pattern (Xie et al 2010), inter-ocean teleconnections (Jia et al 2016), and tropical cyclone distribution (Vecchi and Soden 2007), but also significantly affects regional and continental climate variability and the cryosphere, pointing to an urgent need of accurately projecting global warming pattern.

Figure 1 .
Figure 1.(a), (b) 2022 October to December (OND) US snowfall anomaly from (a) National Operational Hydrologic Remote Sensing Center (NOHRSC) and (b) ERA5 over the period of 2008-2022.(c), (d) The first empirical orthogonal function (EOF) mode of the OND US snowfall (c) and the associated normalized principal component (PC1) time series (d, red line), and the normalized OND US average snowfall index (d, black line) over the period of 1959-2022.The correlation coefficient (r) between the PC1 and US average snowfall is 0.80, statistically significant at the 95% confidence level.

Figure 2 .
Figure 2. (a) The normalized PC1 of the OND US snowfall (red line), and the OND PNA index (black line).The correlation coefficient (r) between the PC1 and PNA is −0.68 during 1959-2000, statistically significant at the 95% confidence level, and only −0.03 during 2001-2022.(b) The regression of Pacific SST against the normalized PC1 during 1959-2000, and the dotted areas indicate the correlation coefficients between PC1 and SST at the 95% confidence level.(c), (d) OND Z500 (contours at 4 m intervals) and snowfall (shaded) regressed against the negative PNA index over 1959-2000 (c) and 2001-2022 (d), respectively.The zero-contour line is omitted for clarity.(e), (f) The OND snowfall responses over the US in the Pacific-Pre2000 experiment (e) and Pacific-Post2000 experiment (f), respectively.The dotted areas indicate the responses at the 95% confidence level.(g) Probability density function (PDF) for the average of snowfall responses over the western US (orange box in e or f) in the Pacific-Pre2000 experiment (red curve) and Pacific-Post2000 experiment (blue curve).The black line corresponds to the observed magnitude of snowfall anomaly in (c).
(c)).Since the PNA can often be triggered by theENSO (Lau 1997, Trenberth et al  1998), we attempt to further link the US snowfall to ENSO by regressing the global SST anomalies onto the PC1 of US OND snowfall.Our results show that the dominant pattern of US snowfall variability (figure 1(c)) is associated with a La Niña-like SST structure (figure 2(b)), consistent with previous studies(Simmons et al 1983).To confirm this La Niña-PNA-US snowfall link before 2000, we conducted two sets of large-ensemble AGCM simulations forced by various SST fields, named as the Pacific SST forcing experiments.In the first pair, named Pacific-Pre2000, one set uses the global climatologicalSST over 1959SST over  -2000   as the boundary condition, and the other uses the regressed pattern of Pacific SST onto the snowfall PC1 over 1959-2000 (figure 2(b)) plus the global climatological SST over 1959-2000 as the boundary condition.
)) imply a potential role of varying background SSTs in modulating the US snowfall influencing factors.Indeed, the early-winter background SST change from 1959-2000 to 2001-2022 exhibit notable patterns on top of global warming (figure 3(a)).The Indian Ocean, North Atlantic, western Pacific, and extratropical North Pacific all experience enhanced warming as compared to the global average, while the equatorial central-eastern Pacific, southwest subtropical Pacific, South Atlantic, and the Southern Ocean undergo a suppressed warming or even a cooling.Historical simulations from the CMIP6 can reproduce some of the features (e.g.hemispheric asymmetry), suggesting an external origin, but the model-observation mismatches are also apparent (figure 3(b)).For example, the observed SST change shows an Indo-western Pacific warming and a La Niña-like pattern, which is absent in the multi-model average, implying a model bias or a role of internal variability (Lehner et al 2020, Tokarska et al 2020) or more likely a combination of both.

Figure 3 .
Figure 3. (a), (b) Difference between the average of1959-2000 and 2001-2022 (2001-2022 minus 1959-2000)  in OND SST for observation (a) and multi-model mean of the CMIP6 historical simulations (b).The dark green contour corresponds to the global average.(c) The correlation coefficient between the PC1 and the corresponding PNA index of each CMIP6 model at different stages; the black asterisks correspond to 1959-2022, the red circles correspond to 1959-2000, and the blue boxes correspond to 2001-2022.The two rightmost columns correspond to the average correlation coefficient of multi-model and the observed correlation coefficient, respectively.The dotted lines of different colors represent the 95% confidence level of each stage.(d) The difference of OND mean SST (i.e.2001-2022 minus 1959-2000) regressed against the correlation coefficient of PNA and PC1 during 2001-2022 across 30 CMIP6 models (table S1 in supporting information S1).

Figure
Figure (a) The regression of the Atlantic SST against the normalized first principal component (PC1) during 2001-2022, and the dotted areas indicate the correlation coefficients between PC1 and SST at the 95% confidence level.(b) The normalized PC1 of the OND US snowfall (red line), and the OND Atlantic SST index (ATL, black line) based on the average of the subpolar North Atlantic (top orange box in a) and the tropical Atlantic (bottom orange box in a).The correlation coefficient (r) between the PC1 and ATL is only −0.11 during 1959-2000, but reaches 0.74 during 2001-2022, statistically significant at the 95% confidence level.(c), (d) OND Z500 (contours at 4 m intervals) and snowfall (shaded) regressed against the ATL over 1959-2000 (c) and 2001-2022 (d), respectively.The zero-contour line is omitted for clarity.(e), (f) The OND snowfall responses over the US in the Atlantic-Pre2000 experiment (e) and Atlantic-Post2000 experiment (f), respectively.The dotted areas indicate the responses at the 95% confidence level.(g) Probability density function (PDF) for the average of snowfall responses over the western US (orange box in e or f) in the Atlantic-Pre2000 experiment (red curve) and Atlantic-Post2000 experiment (blue curve).The black line corresponds to the observed magnitude of snowfall anomaly in (d).