This study explores the wintertime extratropical atmospheric response to Siberian snow anomalies in fall, using observations and two distinct atmospheric general circulation models. The role of the quasi-biennial oscillation (QBO) in modulating this response is discussed by differentiating easterly and westerly QBO years. The remote influence of Siberian snow anomalies is found to be weak in the models, especially in the stratosphere where the “Holton–Tan” effect of the QBO dominates the simulated snow influence on the polar vortex. At the surface, discrepancies between composite analyses from observations and model results question the causal relationship between snow and the atmospheric circulation, suggesting that the atmosphere might have driven snow anomalies rather than the other way around. When both forcings are combined, the simulations suggest destructive interference between the response to positive snow anomalies and easterly QBO (and vice versa), at odds with the hypothesis that the snow–North Atlantic Oscillation/Arctic Oscillation [(N)AO] teleconnection in recent decades has been promoted by the QBO. Although model limitations in capturing the relationship exist, altogether these results suggest that the snow–(N)AO teleconnection may be a stochastic artifact rather than a genuine atmospheric response to snow-cover variability. This study adds to a growing body of evidence suggesting that climate models do not capture a robust and stationary snow–(N)AO relationship. It also highlights the need for extending observations and/or improving models to progress on this matter.
Teleconnections (i.e., long-distance relationships between two or more climate phenomena) have been the subject of extensive research since the pioneering works of Blanford (1884) and Walker (1910). One of the most intriguing proposed teleconnections is the relationship between the extent of snow cover (SNC) in fall over Siberia and the subsequent winter atmospheric circulation in the Northern Hemisphere (NH), especially in the North Atlantic sector. Cohen and Entekhabi (1999) first identified a statistically significant correlation between the fall Siberian snow extent derived from satellite data and the first mode of atmospheric variability of the NH atmosphere in winter, that is, the Arctic Oscillation (AO), also referred to as the northern annular mode (NAM) (e.g., Barnston and Livezey 1987; Thompson and Wallace 1998; Kushner 2010). Several observational studies followed, focusing on the correlation with the North Atlantic Oscillation (NAO) (e.g., Bojariu and Gimeno 2003), exploring possible physical mechanisms for the teleconnection (e.g., Saito et al. 2001) or suggesting a variety of snow-based statistical models to predict the North Atlantic Oscillation/Arctic Oscillation [(N)AO] (Fletcher and Saunders 2006; Cohen and Fletcher 2007; Cohen and Jones 2011).
Numerous modeling studies have also investigated the influence of Siberian snow on the wintertime atmospheric circulation (Gong et al. 2003; Fletcher et al. 2009; Orsolini and Kvamstø 2009; Allen and Zender 2011; Peings et al. 2012; Henderson et al. 2013). Despite a lack of consensus in their results, these studies helped identify a robust physical mechanism consistent with the observed teleconnection, which is described in detail in Cohen et al. (2007) and referred to as the C07 mechanism hereafter. According to this mechanism, increased fall Siberian snow influences NH atmospheric variability in winter through a stratospheric pathway, whereby stratospheric circulation anomalies provide the memory to communicate the effect of snow anomalies several weeks later. Siberia is a source of atmospheric planetary waves due to the presence of orography and land–sea thermal contrast (Held 1983). The largest-scale planetary waves may propagate in the stratosphere and break, inducing eddy heat/momentum flux anomalies that influence the strength of the stratospheric polar vortex (e.g., Polvani and Waugh 2004; Limpasuvan et al. 2004). The aforementioned observational and modeling studies have shown that a greater Siberian snow-cover extent is associated with an increase in upward planetary waves in the Eurasian sector, which propagate into the stratosphere where they may decelerate the westerlies of the polar night jet. The weaker stratospheric polar vortex is then associated with tropospheric anomalies in the following weeks due to stratosphere–troposphere coupling mechanisms (e.g., Baldwin and Dunkerton 2001; Matthewman and Esler 2011; Kidston et al. 2015; Hitchcock and Haynes 2016). This tropospheric response typically projects onto the negative phase of the (N)AO, which has a well-known signature in terms of temperature and precipitation over Eurasia and North America (Hurrell and van Loon 1997).
However, most of the aforementioned studies prescribed unrealistic snow forcing to the AGCMs in order to isolate the snow-driven signal from noise due to internal atmospheric variability [an exception is, e.g., Gong et al. (2003), who found a large response to observed snow anomalies in their model]. In contrast with these idealized sensitivity experiments, the snow–(N)AO teleconnection is not found in coupled ocean–atmosphere simulations from the Coupled Model Intercomparison Project (CMIP3 and CMIP5) (Hardiman et al. 2008; Smith et al. 2011; Furtado et al. 2015). Unlike snow-oriented numerical experiments, the CMIP runs are fully coupled GCM simulations that include other sources of variability and a freely evolving Siberian snow cover. This discrepancy between results from CMIP models and sensitivity experiments has been attributed to GCM biases, including a lack of variance in Siberian snow extent in fall and a misrepresentation of some stratosphere–troposphere coupling processes (Hardiman et al. 2008; Furtado et al. 2015). But it also challenges the robustness of the snow–(N)AO relationship in observations. For instance, Furtado et al. (2015) noticed that significant snow–(N)AO correlations episodically emerge in the CMIP5 simulations over periods of time corresponding to the length of the observational record (~40 years). Such nonstationarity in the snow–(N)AO teleconnection is also found in observations. Peings et al. (2013) and Douville et al. (2017) have made use of recently available extended atmospheric reanalyses to explore the snow–(N)AO teleconnection over a longer period of time. Peings et al. (2013) noticed that Siberian snow cover in fall is remarkably realistic in the NOAA Twentieth Century Reanalysis (20CR; Compo et al. 2011) when compared with satellite and in situ data over the overlap period. After extending the correlation analysis between Siberian snow indices and the (N)AO over the entire twentieth century, Peings et al. (2013) found large multidecadal variability in the snow–(N)AO teleconnection. Only recent decades show significant correlations, a finding confirmed by Douville et al. (2017) using a different twentieth-century reanalysis from the European Centre for Medium-Range Weather Forecasts (ECMWF) (ERA-20C; Poli et al. 2016) to document the atmospheric variability. Douville et al. (2017) also found that over the course of the twentieth century, the influence of the fall Siberian snow on the wintertime circulation switches from significant correlations with the Pacific–North American (PNA; Wallace and Gutzler 1981) pattern to negative correlations with the (N)AO in recent decades. This lack of robustness and the apparent nonstationarity in the snow–(N)AO teleconnection, in nature as well as in the GCMs, is not well understood and calls into question the physical mechanisms of the teleconnection, as well as the relevance of snow-based statistical models for empirical forecasts of the NAO (Cohen and Jones 2011).
Various hypotheses may explain the nonstationarity of the snow–(N)AO teleconnection and help reconcile past results with recent studies: (i) running correlations between pairs of stochastic time series are typically characterized by a low-frequency evolution that may just arise from noise (e.g., Gershunov et al. 2001; Douville et al. 2017); (ii) one or both time series are contaminated by another source of interannual variability, which can alter their relationship (e.g., Peings et al. 2013; Douville et al. 2017); and (iii) multidecadal climate variability alters the atmospheric mean state and thereby the planetary wave propagation that sustains the snow–(N)AO relationship (e.g., Peings et al. 2012). For example, low-frequency sea surface temperature patterns such as the Atlantic multidecadal variability (AMV; Kerr 2000) or the Pacific decadal oscillation (PDO; Mantua et al. 1997) influence the atmospheric state at decadal/multidecadal time scales (e.g., Knight et al. 2006; Peings and Magnusdottir 2015; Kren et al. 2016).
The focus here is on the second hypothesis. At interannual time scales, a potential source of modulation of the snow–(N)AO teleconnection is the quasi-biennial oscillation (QBO). The QBO is a zonally symmetric, descending oscillation of equatorial wind in the stratosphere, from easterly to westerly winds, with a period of about 28 months. The QBO originates from the propagation of mesoscale to planetary-scale atmospheric waves (gravity, inertia–gravity, Kelvin, and Rossby–gravity waves) into the stratosphere, where they transport easterly and westerly zonal momentum. The amount of momentum that is deposed by each type of wave depends on the zonal wind profile and thus differs with the altitude, resulting in zonal wind anomalies that originate in the middle stratosphere and descend with time toward the tropopause (e.g., Lindzen 1987; Haynes 1998). For a comprehensive review of the QBO and associated physical mechanisms, see Baldwin et al. (2001). The QBO in turn influences the propagation of extratropical large-scale wave activity by modifying the location of the zero wind line. Through this mechanism, known as the Holton–Tan mechanism (Holton and Tan 1980), the QBO modulates the strength of the polar vortex, the easterly phase of the QBO being associated with a weaker polar vortex and increased frequency of sudden stratospheric warmings (SSWs), and vice versa (e.g., Anstey and Shepherd 2014). In their studies using twentieth-century reanalyses, Peings et al. (2013) and Douville et al. (2017), point out that the sign of fall Siberian snow anomalies (positive or negative) in recent decades has concurred with the phase of the QBO (easterly or westerly) that promotes similar stratospheric and NAM anomalies. Since the 1970s, high snow-cover years have concurred with easterly QBO (QBO-E) years and low snow-cover years have concurred with westerly QBO (QBO-W) years (this is discussed later in the manuscript and illustrated in Fig. 13). We refer to years with such a concurrence as synchronous years. Before the 1970s, nonsynchronous years predominated; that is, high snow was mostly associated with QBO-W and conversely low snow with QBO-E. Stochastic multidecadal variability of synchronous snow/QBO anomalies could therefore have contributed to the nonstationarity of the snow–(N)AO teleconnection. This is consistent with results from Garfinkel et al. (2010) that reported a drop in correlation between October Eurasian snow and the midwinter polar vortex when the variance explained by other parameters, including the QBO, is removed. By its influence in the extratropics, the QBO potentially preconditions the polar vortex and makes it more or less sensitive to snow-driven planetary wave anomalies. As illustrated in Peings et al. (2012), the climatological state of the polar stratosphere can therefore influence the nature of the atmospheric response to Siberian snow. The potential modulation of the snow–(N)AO teleconnection by the QBO is the motivation of the present paper. Our objective is to determine if such a modulation is found when forcing two atmospheric general circulation models (AGCMs), one of which has a spontaneous QBO and the other a prescribed QBO, with a realistic snow anomaly derived from a recent land surface reanalysis.
Section 2 describes the data, the two AGCMs, and the experimental design. Results are presented in section 3, and the main conclusions are drawn in section 4. They emphasize the difficulty in reconciling the model sensitivity to realistic snow anomalies with the observed snow–(N)AO relationship and, thereby, do not rule out our first hypothesis whereby the recent strengthening of this relationship might be a pure stochastic effect. They also show that the QBO signal dominates the snow signal in the polar stratosphere, in line with our second hypothesis that at least one additional source of climate variability (here the QBO) contaminates the snow–(N)AO teleconnection in observations.
a. Data and models
Two high-top AGCMs are used in this study: the Whole Atmosphere Community Climate Model, version 4 (WACCM4; Marsh et al. 2013), from the National Center for Atmospheric Research (NCAR, United States) and ARPEGE-Climat, version 6, from the Centre National de Recherches Météorologiques (CNRM, France). The two AGCMs are forced with the 1979–2008 climatological annual cycle of sea surface temperature (SST) and sea ice concentration (SIC) derived from the HadISST dataset (Rayner et al. 2003). This setup has been chosen to improve the signal-to-noise ratio by removing all oceanic sources of variability that can otherwise interfere with the snow signal and/or with the QBO variability in the equatorial stratosphere. In line with the C07 mechanism, it assumes that the lagged tropospheric response to the snow forcing is mediated by the stratosphere and not by the ocean, a hypothesis that will be further discussed in this manuscript.
WACCM is a high-top chemistry–climate model with 66 vertical levels from the surface to 5.1 × 10−6 hPa (approximately 140 km). Version 4 includes parameterizations from the low-top CAM4 model (Neale et al. 2013), as well as a parameterization of nonorographic gravity waves that improves the frequency of Northern Hemisphere (NH) SSWs (Richter et al. 2010). It is used with a 1.9° latitude by 2.5° longitude resolution in this study. WACCM4 does not spontaneously simulate the QBO, so the QBO is prescribed by relaxing equatorial zonal winds between 86 and 4 hPa to observed radiosonde data (~28-month period).
ARPEGE-Climat is a high-top climate general circulation model. It is used in this study with an equivalent 1.4° latitude by 1.4° longitude horizontal resolution and with 91 vertical levels from the surface to 0.01 hPa (approximately 80 km). Version 6 used here includes a set of new physical parameterizations: the vertical diffusion scheme is a prognostic turbulent kinetic energy scheme (Cuxart et al. 2000), where the microphysics is the detailed prognostic scheme of Lopez (2002), used for both the large-scale and convective precipitation. The shallow and deep convection are those of the prognostic condensates microphysics transport (PCMT) scheme described in Guérémy (2011). The model also includes a parameterization of nonorographic gravity waves based on the stochastic parameterization described in Lott et al. (2012), which allows the simulation of a spontaneous QBO (see Fig. 2d).
Model results are systematically compared with observations and discussed in light of observational results by Douville et al. (2017). Concerning atmospheric fields, sea level pressure from the ERA-20C (1900–2010) and 20CR (1851–2014) reanalyses are used since they are constrained with surface pressure observations (plus surface winds in the case of ERA-20C). For upper-level fields, where the twentieth-century reanalyses are poorly constrained, the NCEP–NCAR reanalysis is used (1948–2015; Kalnay et al. 1996). Observed daily snow mass anomalies are derived from the recent ERA-Interim land surface reanalysis, available over 1979–2014 (Balsamo et al. 2015). This reanalysis is produced by a land surface model forced with the ERA-Interim atmospheric reanalyses. Two versions of the reanalysis are available, with and without precipitation correction toward monthly precipitation values of the Global Precipitation Climatology Project (GPCP; Adler et al. 2003). The uncorrected precipitation version is used in this study since it is in better agreement with observations concerning Siberian snow variability (E. Dutra 2016, personal communication). The snow water equivalent (SWE; kg m−2) of the ERA-Interim land surface reanalysis is used to prescribe October–November (ON) snow mass anomalies in the perturbation experiments (see below).
b. Observational indices and diagnostic tools
A snow-cover index (SCI) is defined as the averaged snow cover in ON over Siberia [35°–60°N, 40°–180°E], using the 20CR snow cover over 1901–2014. An evaluation of the fall snow cover in 20CR and a comparison with in situ/satellite products is made in Peings et al. (2013). This index is referred to as the “snow index” or SCI in the following. For composite analyses, high and low snow-cover years are defined as years with greater than ±0.5 standard deviation, respectively. The influence of snow is then estimated by averaging years in each category and then subtracting them, and the statistical significance of the composites is assessed using a two-tailed Student’s t test of the null hypothesis.
The QBO is defined as the averaged zonal wind in the 5°S–5°N equatorial band. We use a QBO index that represents the phase of the QBO at 30 hPa in October–November. The ON QBO is used in order to be consistent with the snow index and to assess the predictive potential of the QBO on the winter climate. The QBO reconstruction from S. Brönnimann (Brönnimann et al. 2007) provides an “observed” monthly mean QBO index over 1901–2010. It is extended to 2015 using the NCEP data. Years are split into three categories based on the 30-hPa QBO (westerly, easterly, neutral) using the upper and lower terciles as thresholds. Similar to snow, composite analyses based on the QBO are constructed by subtracting the averaged easterly years from the averaged westerly years.
SSWs are defined as a reversal of the 10-hPa daily zonal wind at 60°N, following the definition of Charlton and Polvani (2007). Events are detected over November–April and must be separated by at least 20 days to be counted as two SSWs, with the final warming occurring at the end of winter being discarded. Following this methodology, we obtain a frequency of 6.1 SSWs decade−1 in NCEP over 1958–2014, in agreement with Butler et al. (2015). Note that the 1948–57 data are removed from this analysis owing to the lack of stratospheric observations assimilated in NCEP prior to 1958. We detect the SSWs events over December–March (DJFM) in the present study, which explains slightly lower frequencies shown in Table 2 compared to using the whole November–April period.
The wave activity in the extratropics is evaluated using the Eliassen–Palm (EP) flux (Edmon et al. 1980) and its three-dimensional generalization, the Plumb flux (Plumb 1985). The horizontal component of the EP flux is proportional to the meridional eddy momentum flux, and the vertical component is proportional to the meridional eddy heat flux. The divergence of the EP flux is a measure of the total nondiabatic forcing on the zonal mean zonal flow with the assumptions of quasigeostrophic theory and linear perturbations [see section 2a in Edmon et al. (1980)]. The zonal mean zonal flow is accelerated where there is divergence of the EP flux and decelerated where there is convergence of the EP flux. The Plumb flux is computed at 850 and 150 hPa in order to geographically localize sources and sinks of wave activity near the surface and near the tropopause.
c. Set of experiments
The set of simulations that are used in this study is described in Table 1. For each AGCM, a 200-yr control simulation is run, and then 200 branched seasonal runs perturbed with a Siberian snow anomaly are started on each 1 October of the control run. The 200 perturbed simulations run from 1 October to 31 March. From 1 October to 30 November, evolving daily snow mass anomalies are imposed over Siberia (40°–80°N 35°–180°E domain) either by replacing (in WACCM) or by nudging (in ARPEGE) the prognostic snow mass. In each grid cell, the prescribed positive snow mass anomalies correspond to two daily standard deviations of the SWE from the ERA-Interim land surface reanalysis (1979–2014 period). After 30 November, the snow constraint is removed and the snowpack evolves freely. Figure 1 shows the averaged snow forcing that is imposed in both models over October–November (Figs. 1a and 1b), as well as its daily evolution (Fig. 1c). Our snow mass perturbation methodology ensures a progressive growth of the anomaly and thus a smooth evolution of the snowpack. The SWE anomaly does not melt during the course of the run; it remains present over Siberia until the end of March (Fig. 1c). Unlike some previous studies that have used fixed and unrealistic snow forcings (e.g., Fletcher et al. 2009; Peings et al. 2012), here the amplitude, spatial pattern, and temporal evolution of the snow anomalies are more realistic, although still of relatively large amplitude.
For each model, the simulated response to the snow anomaly is isolated by subtracting the 200-yr ensemble mean of the perturbation experiment from the 200-yr mean of the control. As for observations, the QBO index is defined as the averaged ON zonal mean zonal wind in the 5°S–5°N tropical latitudinal band at 30 hPa. The 200 years of simulations are then split into westerly and easterly QBO based on the upper/lower terciles of the QBO index, such that 67 years are selected for each phase. Differences in ensemble mean of westerly versus easterly subsets isolate the Holton–Tan effect as simulated by the AGCMs. Like for the observations, the statistical significance of the results is assessed using a two-tailed Student’s t test of the null hypothesis.
a. Evaluation of the Holton–Tan mechanism in the models
A prerequisite for capturing the Holton–Tan mechanism in a model is a good representation of the QBO and of the polar vortex, as well as of planetary waves. Note first that model evaluation is potentially hampered by the lack of interannual SST variability in our AGCM experiments. Yet this is not a major issue since the focus here is on the wintertime extratropical variability, which is mainly driven by internal atmospheric processes. For instance, the interannual variability of the Siberian snow cover, the stratospheric polar vortex, and the midtroposphere midlatitude circulation is realistic in our control experiments, and not much different from the variability obtained in AMIP-type simulations driven by observed rather than climatological monthly mean SST (not shown).
The QBO from NCEP and the models is evaluated in Fig. 2 by comparing it to direct observations from radiosonde measurements in the tropics (http://www.geo.fu-berlin.de/en/met/ag/strat/produkte/qbo/). A subset of 20 years is shown for each dataset (1981–2000 for observations/reanalysis, first 20 years of the control run for the models). The QBO in NCEP has a good timing, but the amplitude of the zonal wind anomalies is lower than in radiosonde observations (Fig. 2b vs Fig. 2a; note different vertical axis in Fig. 2a vs Figs. 2b–d). By design, the QBO in CTLQ is close to radiosonde data since these observations are used to prescribe the QBO in WACCM. ARPEGE-Climat simulates the QBO spontaneously but with a shorter period (~25 months), a smaller amplitude, less downward propagation in the lower stratosphere, and a weaker easterly regime than in observations.
The polar vortex in the two models is compared to NCEP in Fig. 3. The daily zonal mean zonal wind at 10 hPa (U10) is averaged at 60°N in Fig. 3a to compare the seasonal cycle of the polar vortex in NCEP and the experiments. Figures 3b–d, respectively, show a time–latitude Hovmöller plot of U10 in NCEP and the corresponding biases in CTLQ and STCTRL. WACCM lacks seasonal variability in the polar vortex, with westerlies that are too weak in early winter and too strong in fall and late winter (Figs. 3a and 3c). In comparison, ARPEGE-Climat underestimates the strength of the polar vortex during the entire season (Figs. 3a and 3d). Consistent with this mean bias, the daily variability of the polar vortex is larger in ARPEGE-Climat than in NCEP (envelopes of Fig. 3a). This leads to a higher frequency in DJFM SSWs in ARPEGE-Climat (Table 2; 8.8 SSWs decade−1 in STCTRL compared to 5.7 SSWs decade−1 in NCEP). Conversely, WACCM underestimates the frequency of SSWs (4.2 SSWs decade−1), in line with too-strong westerlies in late winter. The polar vortex variability is driven by the dissipation of planetary waves in the stratosphere through wave–mean flow interactions. A good representation of large-scale waves in the atmosphere is therefore important for representing a realistic polar vortex. Figure 4 shows the October–November climatology of the vertical component of the stationary wave activity (Plumb flux WAFz) at 850 and 150 hPa. The main sources of stationary waves originate from orography and land–sea thermal contrast in the North Atlantic, Siberia, and the northwest Pacific. Siberia exhibits the largest amplitude in WAFz, especially at 150 hPa (Figs. 4d–f), making it a key region for generating anomalous planetary waves, in particular due to snow-cover anomalies. The models give a good representation of the climatological WAFz pattern, but WACCM underestimates its amplitude while ARPEGE-Climat overestimates it, especially in the Siberia–North Pacific region at 150 hPa. These biases are persistent throughout winter (not shown), and they are consistent with the biases in the polar vortex since at the interannual time scale increased WAFz at 150 hPa is associated with a weaker polar vortex in both NCEP and the models. This is illustrated in Fig. S1 in the supplement with time–longitude Hovmöller plots of differences in daily 150 hPa WAFz for years of strong minus weak polar vortex in DJF (respectively defined as years with greater than ±0.5 standard deviation of the average 50-hPa geopotential height over the polar cap, north of 65°N). In NCEP as in the models, a weak polar vortex in winter is preceded by a pulse of WAFz at 150 hPa in December over the Eurasian sector. Note that in NCEP this signal originates in November over the western Siberian sector, a key region discussed in section 3b when we evaluate the response to Siberian snow anomalies.
We now examine whether the Holton–Tan mechanism (i.e., easterly QBO promoting a weaker polar vortex, and vice versa) is represented by the models. Figure 5 shows time–latitude Hovmöller plots of differences in daily 30-hPa zonal wind (U30) for QBO-E minus QBO-W years in NCEP (23 yr per QBO cycle), CTLQ, and STCTRL (67 yrs per QBO cycle). By construction, the easterly QBO pattern is apparent in both NCEP and the models, with large negative U30 anomalies in the equatorial band that persist throughout winter. Compared to NCEP (Fig. 5a), the latitudinal extent of the QBO is slightly wider in WACCM (Fig. 5b) and narrower in ARPEGE-Climat (Fig. 5c). In NCEP, the Holton–Tan effect appears in November with negative zonal wind anomalies in the northern high latitudes. This weakening of the polar night jet westerlies persists until the end of winter. The Holton–Tan effect is found in the models, but it appears later on in December and has a smaller amplitude, especially in ARPEGE-Climat in which the weakening of the polar vortex is restricted to midwinter and is not statistically significant.
Figure 6 illustrates the QBO effect on the DJF zonal mean zonal wind and EP flux in a latitude–pressure cross section (vectors for EP flux, red contours for its divergence). Composites from NCEP (Fig. 6a) suggest that there are two mechanisms at play in the Holton–Tan mechanism (illustrated as the difference between QBO-E and QBO-W): a convergence of anomalous upward planetary wave activity north of 50°N and a horizontal divergence of EP flux between high and low latitudes in the stratosphere. The convergence of the EP flux in the polar stratosphere represents a deposit of momentum that decelerates the westerly winds (black contours). The easterly zonal wind anomalies are statistically significant in the stratosphere but not in the troposphere, suggesting that QBO forcing alone is not sufficient to drive negative seasonal (N)AO anomalies at the surface (as discussed at the end of this section). WACCM simulates a weaker polar vortex but does not show the increase in upward planetary wave activity from the surface into the stratosphere (Fig. 6b). It only captures one part of the Holton–Tan mechanism (i.e., the horizontal divergence of EP flux between the subtropics and the pole in the stratosphere). This suggests that the Holton–Tan effect in WACCM is not driven by the meridional eddy heat flux (vertical component of EP flux) but rather by the eddy momentum flux (horizontal component of the EP flux) between the tropics and the polar region in the stratosphere. Not surprisingly given the smaller amplitude of the QBO in this model, ARPEGE-Climat exhibits a smaller weakening of the polar vortex, although statistically significant (shading). The horizontal divergence signal in the midlatitude stratosphere is missing in this model; it only exhibits a weak convergence of upward wave activity in the high latitudes.
The Holton–Tan effect is expected to modulate the frequency of SSWs, with an increase in the likelihood of such events during the easterly versus westerly QBO due to reduced strength of the polar vortex (e.g., Richter et al. 2011). This is found in NCEP, with a doubling in SSW frequency during QBO-E years compared to QBO-W (Table 2). Figure 7 complements results from Table 2 by showing the daily evolution of geopotential-height, polar-cap (Zcap) anomalies along the vertical, for QBO-E versus QBO-W years. This diagnostic is a good proxy of the strength of the polar vortex and is useful for assessing the downward extent of the NAM anomalies that originate in the stratosphere (Baldwin and Thompson 2009). Positive Zcap anomalies denote a weak polar vortex and negative NAM anomalies, and vice versa. In NCEP (Fig. 7a), the weakening of the polar vortex seen during QBO-E is maximum in early December and mid-January. In line with Fig. 6a, this response is only found in the stratosphere and never extends into the troposphere and to the surface. WACCM is consistent with NCEP in simulating positive winter Zcap anomalies (Fig. 7b) and an increase in SSW with QBO-E (~50% increase from QBO-E to QBO-W; Table 2). ARPEGE-Climat shows less dependence on the QBO, with significant Zcap anomalies in midwinter only (Fig. 7c) and consequently little change in SSW frequency depending on the state of the QBO (Table 2).
In summary, the Holton–Tan mechanism is found in both NCEP and the models but with some differences. Not surprisingly given the weak QBO amplitude in this model, the Holton–Tan effect is too small in ARPEGE-Climat to significantly impact the state of the polar vortex climatology/variability throughout the entire winter. The Holton–Tan effect is larger in WACCM, highlighting that a realistic QBO is critical for simulating its effect in the polar stratosphere. However, the physical mechanism seems to differ from NCEP with a different pattern of EP-flux anomalies and the absence of vertical convergence from the troposphere into the stratosphere in high latitudes. These biases in large-scale wave mean-flow interactions have to be kept in mind since they are a limitation for capturing the influence of surface forcings in our models.
The Holton–Tan effect has been shown to project at the surface with significant difference in the (N)AO depending on the phase of the QBO (Watson and Gray 2014). Figure S2 shows seasonal SLP anomalies for QBO-E versus QBO-W, in reanalyses and the models. Note that here we use ERA-20C to create QBO-E versus QBO-W composites (for the longer time span 1901–2010), but results are quite similar when using NCEP (1948–2015). Results from observations and the two models are different and exhibit little statistical significance. The most notable signals are a negative NAO pattern in ERA-20C in December–January (DJ) (Fig. S2b) and a deepening of the Aleutian low in WACCM in DJ (Fig. S2e). The anomalies differ in ON and February–March (FM), which does not support a robust influence of the QBO on the surface atmospheric circulation. However, recall that our QBO index is defined over October–November. The SLP composites are more consistent throughout winter when the corresponding synchronous QBO index is used for each 2-month period (DJ QBO for DJ composites, FM QBO for FM composites; not shown).
b. Response to the Siberian snow anomaly
This section focuses on the influence of Siberian snow on the atmospheric circulation, regardless of the phase of the QBO. For observations, we present composites based on the snow-cover index of 20CR (over 1901–2010 for ERA-20C fields, 1948–2014 for NCEP fields). For the models, the influence of snow is shown by taking the difference of the ensemble mean of the snow perturbation experiment and the control run (SNOWQ-CTLQ for WACCM, STQBO-STCTRL for ARPEGE-Climat). To present a complete overview of the observed composites and of the simulated responses, we show bimonthly (ON/DJ/FM) anomalies of various fields, for high versus low ON Siberian snow cover.
In fall, snow anomalies mostly modify the surface energy budget through changes in land surface albedo (radiative effect of snow); therefore, in this study anomalies in SNC are more relevant than SWE or snow depth anomalies. Figure 8 presents high versus low snow SNC anomalies in 20CR, WACCM, and ARPEGE-Climat. In 20CR, SNC anomalies are present over southern Siberia in ON, with a maximum over central Siberia (Fig. 8a). As snow cover reaches 100% throughout winter, only the southern part of the anomalies persist (Figs. 8b and 8c). The SNC anomalies are much larger in the models, especially in ON, due to the fact that the daily snow mass anomalies that are imposed represent two standard deviations at each grid point (while the observed snow-cover anomalies are constructed based on the averaged SNC anomaly in the Siberian domain). Note, however, that the radiative effect of snow depends not only on its extent but also on the amount of incoming solar radiation. Thus, north of 60°N, snow-cover anomalies that are present in the models but not in observations have little impact on the energy budget due to low incoming solar radiation in fall/winter. High snow anomalies induce a surface cooling over Siberia, as shown in Fig. S3 with 2-m temperature (T2M) anomalies. The cooling is particularly pronounced in ON and persists during winter. Note that ERA-20C 2-m temperatures (as well as HadISST sea surface temperatures) suggest a downstream propagation of the Siberian cold anomaly over the North Pacific, which was also found in snow sensitivity experiments using a slab-ocean model (Henderson et al. 2013). Not surprisingly this downstream cooling is not captured by the models given the use of prescribed climatological SST. However, according to the C07 mechanism, this oceanic response is secondary in comparison to the land surface cooling that is sufficient to alter the atmospheric circulation through a reinforcement of the Siberian high. The latter signal is clearly found in the models, especially in fall when SNC and T2M anomalies are the most pronounced (Figs. 9d and 9g). However, the observed pattern of ON SLP anomalies differs in a remarkable way (Fig. 9a). Instead of local high pressure anomalies, our composite analysis identifies an upstream high to the northwest of the SNC anomalies, as well as low pressure anomalies in the Bering Strait region. A similar dipole pattern, consisting of an upstream high over eastern Europe and a downstream low in the North Pacific, has been identified as a precursor of weak polar vortex due to constructive interferences with the climatological stationary waves 1 and 2 (Garfinkel et al. 2010). The western expansion of the Siberian high is also found in November SLP composites, following October snow anomalies, suggesting a lagged influence of snow on the atmosphere (Cohen et al. 2014). However, whether the pattern in Fig. 9a is truly a response to SNC anomalies is dubious. First, the model responses differ greatly from observations (Figs. 9d and 9g), even though the models present no obvious reason to misrepresent the local atmospheric response to the snow forcing. Second, it extends well beyond the latitude of the SNC anomaly over the Arctic. This anticyclonic pattern advects polar air and leads to increased snowfall over Siberia, such that it is possibly causing the Siberian snow anomaly, rather than being forced by it [supporting similar results by Kryjov (2015)]. A similar high pressure anomaly has been found in response to decreased sea ice in the Barents and Kara (BK) Seas (Honda et al. 2009), which is also associated with cold Eurasian temperature and polar vortex weakening through anomalous upward wave activity (e.g., Peings and Magnusdottir 2014; Kim et al. 2014; Kug et al. 2015). However, when compositing ON SLP over a September Barents–Kara sea ice index, we do not find signals that resemble Fig. 9a, suggesting that this atmospheric pattern is not driven by sea ice anomalies (Fig. S4). How the SLP/T2M anomalies of Fig. 9 and Fig. S3 connect with Siberian snow cover and Arctic sea ice anomalies is beyond the scope of the present study and will be the focus of future works. But clearly, the discrepancy between the observed and simulated fall SLP anomalies is a warning flag that the observed anomalous atmospheric patterns should not be interpreted as a response to snow only. It may also reveal limitations of using AGCM experiments as oceanic and sea ice feedbacks that have the potential to reinforce the response are neglected.
In winter, as previously shown in Douville et al. (2017), the SLP pattern resembles the negative NAM pattern, although it is stronger over the North Pacific (Fig. 9b). This signal masks some large multidecadal variability. The atmospheric pattern associated with Siberian snow has shifted from a PNA pattern prior to the 1970s to a (N)AO pattern only in recent decades (Douville et al. 2017). The signal is predominant in midwinter [DJ; Fig. 9b] but not statistically significant in late winter (FM; Fig. 9c). In comparison, the two models exhibit very little SLP response, except for maintaining local high SLP anomalies over the forcing region. The absence of a larger hemispheric atmospheric response in the models can be explained by the absence of phasing between the snow-driven anomalous stationary waves and the climatological stationary waves. As stated before, snow is expected to exert a perturbation of the stratosphere and an associated remote response when it forces stationary wave anomalies that constructively interfere with the climatological stationary waves (Smith et al. 2011). Figure 10 compares the response to forcing to the climatology in stationary wavenumber 1, estimated by Fourier transforming the geopotential height at 60°N in a pressure–longitude plot. A similar figure is shown in SI for all wavenumbers (Fig. S5, only removing the zonal mean of geopotential height without Fourier transforming). In NCEP (Fig. 10b), large constructive interferences occur in midwinter between the anomalous (contour) and climatological (shading) wave pattern, as indicated by a spatial pattern correlation of 0.81. Such constructive interference is absent in the models, which rather exhibit negative correlations (i.e., destructive interferences that prevent the surface anomaly from communicating upward and into the stratosphere). Note that the good correspondence between observations and models in ON (Figs. 10a, 10d, and 10g) is only found for wavenumber 1. It is absent when all wavenumbers are considered (Fig. S5). The absence of upward propagation of the signal in the models is striking when comparing WAFz anomalies at 850 (Fig. S6) and at 150 hPa (Fig. S7). Anomalous wave activity is present in fall near the surface in the models (Figs. S6d,g); although similar to the SLP, the pattern of anomalous WAFz in the models differs from observations (Fig. S6a). Once again, the fact that WAFz anomalies are located above the snow forcing in the models, but not in observations, calls into question the attribution of the ON observed composites to snow only. The pulse in WAFz near the surface does not propagate at 150 hPa in the models (Figs. S7d,g), while this is the case in observations (Fig. S7a). Reasons for the absence of constructive interferences in the models are unclear. Either the signal is not only driven by snow in observations such that the imposed snow forcing is not sufficient to drive significant upward wave activity propagation in the models (in this case, the models are not wrong). Or the models lack sensitivity to the forcing and thus underestimate the response. In the latter case, model biases in terms of zonal wind profile and in the mean state of the polar vortex (discussed in section 3a) may be involved, since the zonal circulation acts as a waveguide for planetary waves (Edmon et al. 1980). In DJ, observations exhibit positive WAFz anomalies over the SNC anomalies while the signal vanishes in the models (Figs. S6b,e,h and S7b,e,h). Unlike the fall signal discussed before that is ambiguous, this midwinter pulse in planetary wave activity is more likely induced by snow anomalies that might have a delayed influence on the polar vortex in midwinter. Note, however, that a large anomaly of WAFz is also found over the northeast Pacific sector that may not be driven by snow and possibly contaminates the observational composites (we will discuss this signal further at the end of the section).
The EP-flux anomalies of Fig. 11 further illustrate the lack of upward anomalous wave propagation in the models compared to NCEP (vectors). Consistently, the models show small zonal wind anomalies that have little statistical significance in the stratosphere (black contours; note that the contour interval is 5 times smaller for the models). Unlike the models, the NCEP composites exhibit a large increase in upward planetary wave activity that propagates and converges into the stratosphere, where it decelerates the westerlies through wave–mean flow interactions (Figs. 11a,d,g). This perturbation of the stratosphere is reflected in Zcap anomalies (Fig. 12), where a weaker polar vortex is identified in NCEP. The stratospheric response is already present in October–November, before it gets reinforced during early winter with the formation of the polar vortex. In line with previous results, Zcap anomalies are small in the models (Figs. 12b,c), especially compared to the Holton–Tan effect (Fig. 7). It is interesting to also measure the stratospheric response in terms of SSWs. Changes in frequency of SSW for high Siberian snow years are given in Table 2. The change in SSW is also given for low snow years in NCEP but not for the models since our simulations only examine the response to an excess of snow. The findings are consistent with the Zcap anomalies. In NCEP, the response in SSW is quite linear between high and low snow years, with an increase of ~35% during high snow years and a decrease of ~40% during low snow years. The change is comparable to the change associated with QBO-E versus QBO-W, suggesting that the QBO is as likely to perturb the frequency in SSW as is Siberian snow. WACCM simulates a small increase in SSWs that is also on the same order of amplitude as the QBO signal, while in ARPEGE-Climat the frequency of SSW decreases, consistent with the late-winter reinforcement of the polar vortex that is found in this model (Figs. 11i and 12c).
To sum up this section, our numerical experiments do not capture any significant snow–(N)AO teleconnection. They only reproduce the first link of the C07 mechanism, which is an increase in upward planetary wave activity in the vicinity of the enhanced snow cover and of the associated surface temperature and sea level pressure anomalies. The anomalous wave activity does not propagate to upper levels, resulting in a smaller effect (or even slightly opposite in the case of ARPEGE-Climat) on the polar vortex and SSW than in NCEP. Yet the observational results should be interpreted with caution, since our composite analyses reveal atmospheric precursors that can hardly be attributed to snow alone. From our findings we cannot rule out that the observed snow–(N)AO relationship is purely stochastic or is the result of some other atmospheric forcing on both snow and (N)AO. For instance, a positive PDO signal is found in composites of SST based on the ON SCI, both in the pre-NCEP and NCEP era (1901–47 and 1948–2014; Figs. S8a–f). Composites of SLP and 150-hPa WAFz based on an ON PDO index (defined as the first EOF of SST north of 20°N in the Pacific) exhibit some resemblance with similar composites constructed using the snow index (Figs. S8g and S8h vs Figs. 9a and S6b, respectively), especially over the North Pacific sector. Although Siberian snow anomalies in fall potentially impact the downstream atmospheric circulation and SST in the North Pacific (e.g., Henderson et al. 2013), the PDO has been shown to be related to various other climate phenomena, including El Niño–Southern Oscillation (ENSO). Therefore, the PDO-like signals in the snow composites are unlikely to be driven by snow alone and are probably contaminated by other climate processes. Furthermore, in the polar stratosphere the effect of the QBO appears to be at least as large as the effect of snow (comparable amplitude in NCEP, larger amplitude in the models). Therefore, the potential role of the QBO in canceling or reshaping the response to snow is also substantial. This is investigated in the following section.
c. Response to Siberian snow depending on the phase of the QBO
Peings et al. (2013) and Douville et al. (2017) have shown that over recent decades there has been a remarkable concurrence of QBO and snow anomalies that promote the same sign of perturbation on the polar vortex (i.e., SCI+/QBO-E and SCI−/QBO-W have dominated, in years that we refer to as “synchronous years”). Conversely, before the 1970s, “nonsynchronous years” have predominated; that is, SCI+ years rather concurred with QBO-W, while SCI– concurred with QBO-E. This is illustrated in Fig. 13, which shows the 31-yr moving window correlation between the ON SCI and the DJF AO (defined as the first EOF of DJF SLP in the NH) over 1901–2010. Synchronous snow–QBO years (green dots) prevail after the 1970s when the snow–(N)AO correlation is statistically significant. Before the 1970s, nonsynchronous snow–QBO years (red dots) predominated and no significant correlation is found. The change in synchronous versus nonsynchronous snow–QBO years is a possible stochastic source of multidecadal variability in the snow–(N)AO teleconnection. It is not necessarily associated with any physical mechanism but can simply arise from random fluctuations in the concurrence between snow and QBO anomalies.
The snow–QBO coexistence is apparent in the zonal wind cross sections of Figs. 11a,d,g. The presence of the easterly QBO signal contaminates the observational composites (as does the PDO pattern mentioned at the end of section 3b) and makes the attribution of the signal to snow alone questionable. We make the two following hypotheses: either 1) the QBO promotes a greater response to snow by preconditioning the polar vortex, such that both forcings interfere constructively to induce significant stratospheric and (N)AO perturbations (which was hypothesized in Peings et al. 2013), or 2) the QBO dominates the response to snow, such that the apparent influence of snow in observations is artificially inflated by its concurrence with the QBO. While it cannot be ruled out from observational analyses such as Douville et al. (2017), the assumption that snow dominates and induces an artificial effect of the QBO is not considered here since it is not supported by our simulations. To test these hypotheses, we have plotted in Fig. 14 daily Zcap anomalies for different cases that allow us to summarize the main findings of this study:
Figures 14c and 14h show the combined response to the QBO and to Siberian snow in the two models. It is computed as the difference between high snow/QBO-E years and low snow/QBO-W years (i.e. QBO-E years from the perturbation snow experiment minus QBO-W years from the control). As stated before, in WACCM snow anomalies reinforce the response in early and late winter, but they are of secondary importance compared to the Holton–Tan effect. In ARPEGE-Climat, snow has a negligible effect except for a strengthening of the polar vortex in late winter that reinforces the similar QBO-driven signal.
Figures 14d and 14i show the response to Siberian snow under the easterly phase of the QBO. It is computed as the difference between QBO-E years from the perturbation run and the control run (67 years per phase).
Similar results from NCEP are shown in Fig. S9 for the sake of completeness, but they are subject to high uncertainties given the small sample sizes after compositing on both snow and QBO (less than 10 years).
A larger weakening of the polar vortex is found in WACCM during the QBO-W (Fig. 14e) than during the QBO-E (Fig. 14d). Therefore in WACCM, snow mitigates the strengthening of the polar vortex forced by QBO-W, rather than amplifying its weakening during QBO-E. This is confirmed when looking at changes in SSW frequency (Table 2). The increase in SSW during high snow years is more pronounced under QBO-W (3.1 to 4.6 SSWs decade−1) than under QBO-E for which it is already relatively high (4.5 to 4.9 SSWs decade−1). Destructive interferences between snow and QBO thus dominate in WACCM. ARPEGE-Climat also exhibits a larger response to snow under QBO-W conditions (Fig. 14j vs Fig. 14i), although it does not lead to an increase in SSW frequency (Table 2) in average over winter. These results are at odds with the hypothesis that the QBO promotes the effect of snow through constructive interferences, since snow mitigates, rather than amplifies, the Holton–Tan effect in our models. They do not support the Peings et al. (2013) hypothesis that synchronization in snow and QBO anomalies has enhanced the influence of snow on the stratosphere in recent decades. Rather, they suggest that the apparent influence of Siberian snow on the polar vortex is artificially inflated by its concurrence with the QBO. However, recall that this applies to the stratospheric signal only since the effect of the QBO is not significant at the surface (Fig. S2). Therefore, the QBO is not sufficient to explain multidecadal fluctuations of the snow–(N)AO teleconnection, which could potentially originate from atmospheric stochastic variability only or other climate phenomena such as the PDO (Douville et al. 2017).
This study has investigated the influence of fall Siberian snow anomalies on the wintertime extratropical atmospheric circulation. Realistic (albeit rather strong in magnitude and spatial extent) October–November snow mass anomalies have been prescribed in two state-of-the-art AGCMs with prescribed or spontaneous QBO, WACCM and ARPEGE-Climat, respectively. The effect of snow anomalies has been decomposed depending on the phase of the QBO, in order to test whether the QBO can modulate the atmospheric response to fall Siberian snow anomalies, as hypothesized in Peings et al. (2013). Our findings and their implications are summarized as follows:
Despite some differences, the models qualitatively capture the Holton–Tan mechanism that is identified in the NCEP reanalyses (i.e., QBO-E is associated with a weaker polar vortex, and vice versa). However, the amplitude is underestimated, especially in ARPEGE-Climat, which also lacks in amplitude of the QBO. Analyses of wave activity anomalies reveal that the models only partially reproduce the mechanisms behind the Holton–Tan effect. In WACCM, it is only associated with eddy momentum rather than eddy heat flux anomalies in the stratosphere, while ARPEGE-Climat only simulates a weak vertical convergence of planetary waves in the polar stratosphere.
The effect of snow anomalies is weak in the models, in both the stratosphere and the troposphere. WACCM simulates a weaker polar vortex with higher Siberian snow, but the signal is small compared to the Holton–Tan effect. ARPEGE-Climat shows very little sensitivity to the snow forcing, except for a local cooling and reinforced Siberian high at the surface. In contrast, in NCEP the effect of Siberian snow and of the QBO on the polar stratosphere are of comparable amplitude. However, we point out that the composites based on the Siberian snow index are ambiguous owing to their contamination by the QBO and an SST pattern that resembles the PDO. Added to the fact that the local surface atmospheric response in our simulations greatly differs from the observational composites, we question the causal relationship between snow and the atmosphere since in nature the atmospheric circulation might have driven snow variability, rather than the other way around. In the case that snow is actually responsible for the observed fall SLP anomalies, reasons for the failure of our models in capturing this direct effect of snow would have to be identified. One possibility is that the indirect influence of snow through its downstream influence on North Pacific SST (Henderson et al. 2013) is missing in our AGCM-only study. SST variability might amplify the snow influence in the real world, although this view is challenged by the PDO-like signal identified in the SST composites of Fig. S8, a basinwide mode of variability hardly attributable to snow variability alone.
Our simulations show a larger influence of snow in the stratosphere under QBO-W than under QBO-E. Yet in recent decades, when the snow–(N)AO teleconnection is significant in observations, high snow has predominantly concurred with QBO-E and conversely low snow with QBO-W. This contradiction between observations and simulations further challenges our understanding of the snow–(N)AO teleconnection, especially concerning the stratospheric pathway.
Returning to the hypotheses that were listed in the introduction to explain the nonstationarity of the snow–(N)AO teleconnection, our modeling experiments rather support the first hypothesis; that is, the snow–(N)AO teleconnection potentially arises from stochastic variability alone. The absence of a significant negative (N)AO pattern in response to snow forcing only does not support the existence of a robust snow–(N)AO teleconnection. This statement is not necessarily in conflict with observational results. As discussed by Douville et al. (2017), multidecadal fluctuations similar to the observed snow–(N)AO correlation over the whole twentieth century can be found in random white noise time series. Such multidecadal variations in the correlations are also found in long preindustrial control runs of CMIP5 models (Furtado et al. 2015) as well as in the ARPEGE-Climat control run despite the use of climatological SSTs (Fig. S10). Regarding the stratospheric component of the snow–(N)AO teleconnection, the second hypothesis (i.e., the time series are influenced by another mode of interannual variability, which alter their relationship) is also supported since random concurrence of snow and QBO anomalies in observations artificially inflates the influence of snow on the polar vortex variability. Although not considered in the present study, our last hypothesis [i.e., multidecadal climate variability alters the atmospheric mean state and thereby the planetary wave propagation that sustain the snow–(N)AO relationship] cannot be ruled out. In particular, a PDO signal is found in the snow composites that, similarly to the QBO, contaminates the observed signals. Whether the PDO amplifies or mitigates the influence of Siberian snow on the atmosphere at long time scales will be the object of future studies.
These conclusions apply to our AGCMs only and one should not generalize the results. As shown in section 3a, the models have biases, both in terms of polar vortex variability and in terms of their representation of the QBO (prescribed in WACCM, spontaneous but too weak and symmetric in ARPEGE-Climat). As discussed before, the use of climatological SST, although useful to isolate the response to snow, is a simplification that might lead to underestimating the response to snow. We also used an idealized snow forcing (even though it is more realistic than most previous studies), and it is possible that the atmospheric response depends on the pattern of the snow-cover anomalies rather than on its amplitude, especially for interfering constructively with the climatological background. To what extent these model biases and setup affect the result of the study is hard to address, but a similar set of experiments with different models and/or setups would be useful to assess the possible limitations of our study.
Recent studies have reported encouraging skill in predicting the NAO up to more than a year ahead, using both dynamical and statistical forecast schemes (Scaife et al. 2014; Dunstone et al. 2016; Wang et al. 2017). However, the predictive skill of the NAO fluctuates at multidecadal time scale, for reasons that are still to be understood (Weisheimer et al. 2017). The present study is an illustration of the challenges that we face for understanding nonstationarities in large-scale teleconnections and their associated predictability. This does not apply to the snow–(N)AO teleconnection only; for instance, Koenigk and Brodeau (2017) report similar nonstationarity in the Arctic sea ice–(N)AO teleconnection in quasi-equilibrium GCM experiments. The relative influence of Siberian snow, of the QBO, and of other potential drivers of the interannual variability of the (N)AO is hard to extract from the limited observational record owing to large internal variability of the atmosphere. Further, the model biases in representing the polar vortex or the QBO still limit their relevance for examining the processes, making it hard to draw robust conclusions. The availability of extended observations and reanalyses in the future, as well as new model developments, will hopefully help to reconcile the observational with the model results. In the meantime, our results question whether the recent strength of the snow–(N)AO relationship may be a stochastic artifact rather than a genuine atmospheric response to snow-cover variability.
Thanks are due to Emanuel Dutra for providing us with the snow field from the ERA-Interim Land surface reanalysis. We thank Paul Kushner, Judah Cohen and an anonymous reviewer for their comments that helped to improve this manuscript. GM and YP are supported by NSF Grants AGS-1407360 and AGS-1624038. We acknowledge high-performance computing support from Yellowstone (ark:/85065/d7wd3xhc) provided by NCAR’s CISL, sponsored by the NSF.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-17-0041.s1.