The potential for anthropogenic climate change to impact patterns of seasonal snow cover has motivated numerous studies seeking trends in its extent and duration. Many have been based on the NOAA-Rutgers record of Northern Hemisphere snow cover. Several studies have found augmented early-season snow identified from this archive to be anomalous, and related it to the introduction of higher-resolution imagery and a more automated interpretation process in 1999. This study contributes to the discussion by describing in greater detail the spatial and temporal distributions of trends in the onset of seasonally snow-dominated conditions between 1972 and 2017, and relationships to their physiographic and climatological contexts. It also identifies changepoints between negative and positive onset-date anomalies, and relates these to corresponding meteorological patterns. Most trends identified indicated earlier onset, and were associated with midlatitudes, low to moderate elevations, and colder, drier climates. These were situated largely northeast of major topographic chains, southwest of increasingly ice-free Arctic waters, and to the east of areas associated with blocking systems. Onset-date anomalies switched from positive to negative in approximately 70% of the affected points before 1997. These changepoints generally occurred earlier at higher elevations to the south and west, and later at lower elevations to the north and east. Overall temporal trajectories correspond broadly to shifts in temperature and precipitation over the same areas. In contrast, positive (later) onset trends were found over much smaller areas, associated with warmer, wetter climates and higher elevations, particularly on west-facing slopes; temporal variations in anomalies of their onset dates and associated meteorological conditions were distinct from those having negative trends.
The annual transformation of terrestrial physical properties and conditions by seasonal snow cover (SSC) modulates a wide range of climatological, hydrological, ecological, and geomorphological processes, and controls both the provision of water resources and the risk of flood or drought for large human populations (Mankin et al. 2015). These influences are particularly important over the Northern Hemisphere (NH), which hosts approximately 98% of the total global area of SSC (Armstrong and Brodzik 2001). The extensive area of NH SSC even at relatively low latitudes means that snow contributes a greater negative top-of-atmosphere radiative forcing than sea ice in every month outside May–August (Flanner et al. 2011; based on analysis from 1979 to 2008). Investigations of the potential for disruption of hitherto reliable distributions of NH SSC by atmospheric warming and/or changing precipitation patterns are therefore highly relevant to research into trajectories of anthropogenic climate change.
The description and quantification of snow-related trends has been a key topic of interest among these studies, supported by a range of lengthening observational series. Remotely sensed imagery has proven particularly useful for monitoring the spatial extent of SSC, as it supports the frequent inference of conditions over wide areas, including remote and complex terrain, and long time periods. Among such resources, the U.S. National Oceanic and Atmospheric Administration’s climate data record of NH snow cover, compiled and curated by the Rutgers University Snow Laboratory [hereafter the NOAA-Rutgers Snow Archive (NRSA)], offers the longest period of record of all remotely sensed environmental datasets (Robinson and Estilow 2012; Estilow et al. 2015). It has therefore formed the foundation for a series of SSC-related investigations since the mid-1990s (Robinson et al. 1993; Groisman et al. 1994; Frei et al. 1999; Dye 2002; Rikiishi et al. 2004; Déry and Brown 2007; Choi et al. 2010; Brown and Robinson 2011; Cohen et al. 2012; Liu et al. 2012; Hernández-Henríquez et al. 2015; Allchin and Déry 2017).
However, the credibility of the NRSA has been debated in recent years. One point of discussion has centered on significant positive trends in snow-covered area identified from its records during the seasonal transition to boreal winter. These were first reported by Frei et al. (1999) and Brown (2000), based on annual series ending no later than 2000. Rikiishi et al. (2004), using the same time span, suggested that changing atmospheric circulations might have caused the widespread negative trends in SSC onset they described across the NH midlatitudes. Having extended the period of interest to 2008, Choi et al. (2010) also noted signs of earlier onset, and observed that warming need not preclude additional snowfall while temperatures remain below freezing. Déry and Brown (2007), Hernández-Henríquez et al. (2015), and Allchin and Déry (2017) all reported increases in NH snow-dominated area between October and December. However, several studies have found these early-season gains to be anomalous when compared to metrics derived from a variety of other modeled and observed datasets (Brown and Derksen 2013; Mudryk et al. 2017; Hori et al. 2017). The explanation most frequently suggested has been that they may be artifacts relating to the introduction of imagery at higher temporal and spatial resolutions, in conjunction with a more automated interpretation process, in 1999 [as described by Ramsay (1998)]. A further complication arises from recognition that calculating snow-covered area based on the conventional assumptions of 100% (0%) actual fractional snow cover within the NRSA’s spatial subunits when classified in the binary state of “snow-covered” (“snow-free”) is likely to exaggerate estimates of trend magnitudes, and that this may be particularly problematic during transitions to and from the snow-dominated season (Allchin and Déry 2017).
Nevertheless, others have incorporated these representations of more widespread early-season snow within conceptualizations of feedbacks between rapid rates of warming and more ice-free waters in the Arctic, leading to perturbations in atmospheric circulations (relating particularly to the boreal polar jet). These in turn have been identified as potential drivers of extreme weather in the midlatitudes, with some areas experiencing colder winter conditions as a consequence of meridional outflows from higher to lower latitudes, while others have been anomalously warm due to heat advection in the opposite direction (Overland et al. 2011; Cohen et al. 2014; Francis and Skific 2015; Wegmann et al. 2015; Francis et al. 2017; Kretschmer et al. 2018; Cohen et al. 2018).
This study was therefore conceived with the goal of generating more detailed information relating to the NH snow-onset season, without depending on estimates of SSC extent. Its principal aim was to describe the spatial distributions of areas which, according to the representation of snow cover provided by the NRSA, have experienced significant negative (earlier) or positive (later) trends in the onset date of SSC during recent decades, their associated physiographic and climatological characteristics, and the times of year in which these shifts occurred. A secondary goal was to determine the timings of any clear transitions during the period of record from one sign of snow onset date anomaly to the other and, if so, to compare these with equivalent trajectories in metrics of meteorological activity.
As well as enabling contributions to the discussion surrounding the veracity of early-season SSC trends identified from the NRSA, use of this archive also provided the benefit of a long period of record, and an opportunity for comparison of results with prior research relating to shifting patterns of SSC over the NH. The dataset’s spatial framework comprises a grid of 89 rows and columns in the NH stereographic polar projection, so that the areas associated with each dataset grid point (DGP) vary with latitude, from ~10 700 km2 at the lowest latitudes to ~41 800 km2 near the geographic North Pole (Estilow 2013). DGPs associated with primarily oceanic areas, or terrestrial areas of (hitherto) perennial cryospheric cover, together with an additional 33 previously associated with questionable snow cover records (Déry and Brown 2007), were excluded using masks supplied with the archive or developed in previous studies (Déry and Brown 2007; Hernández-Henríquez et al. 2015; Allchin and Déry 2017). The remaining spatial domain comprised 1980 DGPs, associated with a total area of 5.9 × 107 km2.
Each weekly NRSA granule comprises a series of binary values, one per DGP, signifying whether or not at least 50% snow cover existed over the associated terrestrial area on the last unobscured day of the prior seven (from Tuesday to Monday). While the archive’s first observations date from October 1966, gaps totaling 37 weeks between July 1968 and September 1971 have guided the selection of October 1971 as a start date for the period adopted in previous analyses (e.g., Déry and Brown 2007; Hernández-Henríquez et al. 2015; Allchin and Déry 2017). In this case, for reasons clarified below, week 32 (7–13 August) of 1972 was chosen as the start date. In the version of the archive kindly supplied by the Rutgers Snow Laboratory, records were available to week 31 (31 July–6 August) of 2018.
Through the long history of the NRSA, advances in hardware and software capabilities have driven a series of improvements to the sensors and processing algorithms used in its production. The greatest such shift occurred in 1999, when the original (mainly human-dependent) interpretation process was replaced by the more automated Interactive Multisensor Snow and Ice Mapping System (IMS). At the same time, the imagery’s spatial resolution improved from ~190 to 24 km, and its temporal resolution changed from weekly to daily (Ramsay 1998; Estilow et al. 2015). While considerable care has been taken to ensure overall consistency throughout the entire period of record (Romanov et al. 2000; Helfrich et al. 2007; Estilow 2013; Estilow et al. 2015), this upgrade has nevertheless been identified by several authors as a potential causative factor behind disagreements with snow-related metrics derived from other datasets (e.g., Brown and Derksen 2013; Mudryk et al. 2017; Hori et al. 2017), including the early-season trends of interest here.
The long-term climatological contexts of DGPs were derived from the 1961–90 mean temperature and precipitation datasets published by the University of East Anglia’s Climatic Research Unit (CRU) (version 2.0, released in 2002: http://crudata.uea.ac.uk/cru/data/hrg/tmc). The constituent values, supplied at a spatial resolution of 10 arc-minutes, were averaged within polygons representing the DGP surface footprints. The usual caveats relating to the uncertainties associated with interpolating ground observations over wide separations between stations (particularly in remote and/or complex terrain) are noted.
Meteorological activity during the period of interest was represented by monthly mean terrestrial temperature and precipitation datasets published (on a 0.5° grid) by Willmott and Matsuura of the University of Delaware (version 5.01, released August 2018; http://climate.geog.udel.edu/~climate), processed in the same manner.
a. Derivation of metrics of seasonal snow onset
To distinguish between the initial (often ephemeral) onset of snow cover during the transition from autumn to winter, and its greater persistence through the subsequent colder months, two metrics were tracked for every DGP. The first date in the latter part of each calendar year on which snow cover switched from <50% to ≥50% was named SnowOn1, and the beginning of the first period of uninterrupted snow-dominated conditions continuing for at least 4 weeks was termed SnowOn4. Mean SnowOn4 dates were found to correlate strongly with mean start dates of the longest continuous runs of snow-dominated conditions (R = 0.99, p < 0.0001, slope = 0.9), so it is reasonable to equate them with the beginning of the main winter season. SnowOn4 occurs within 1 week of SnowOn1 in ~50% of DGPs in which both dates were identifiable in at least 30 years: this fraction rises to 80% within 2 weeks, and 93% within 3 weeks.
For consistency with prior studies, this analysis adopted a snow-year beginning on 1 October. However, initial exploration of the dataset revealed that, of the 1980 DGPs considered, mean SnowOn1 dates from 1971 to 2017 were earlier than this for 185 of them (associated with a total area of 6.2 × 106 km2, comprising 10.5% of the spatial domain). Over the same period, mean SnowOn4 occurred before 1 October in 134 DGPs (total area 4.4 × 106 km2; 7.5% of the spatial domain). The earliest mean value of either metric was 30 August. To capture these transitions as comprehensively as possible, they were therefore sought in granules captured in or after week 32 of each year, as this begins between 2 and 9 August (depending on the year). Onset dates before 1 October in any year were therefore represented by negative values. This in turn dictated that the first year in which to seek onset dates would be 1972 (rather than 1971, as a data gap exists to September of that year), and the last 2017, a total of 46 years. (Note that all statistics quoted hereafter pertain to this 1972–2017 period of interest, unless explicitly stated otherwise.)
b. Identification and estimation of trends
The principal goal of this study was to identify DGPs associated with snow-onset time series manifesting significant monotonic trends. A common confounding factor in trend analysis is the existence of temporal autocorrelation within the series of interest (e.g., Weatherhead et al. 1998). Where this takes the form of a lag-1 autoregressive structure, it is possible to mitigate its effects by “pre-whitening” (Yue et al. 2002). Evidence of this phenomenon was therefore sought among the time series of snow-onset dates, using the full and partial autocorrelation functions described by Anderson (1977). However, no such patterns were identified.
Trends were identified from the SnowOn1 and SnowOn4 series using Mann–Kendall trend analysis (Mann 1945; Kendall 1975). With the goal of reporting conservative results, the series for each DGP was considered only if it comprised a full set of 46 annual values. The magnitudes of trends identified as significant at the 5% level were estimated using the Theil–Sen method (Theil 1950a,b,c; Sen 1968).
c. Illustration of spatial, temporal, and climatological distribution of trends
Having identified DGPs associated with significant snow-onset trends, descriptions were generated of their spatial and temporal distributions. The latter is a key point, as mean snow-onset dates vary with contextual attributes such as latitude, elevation, and proximity to an upwind coastline (Fig. 1), so the drivers of any trend identified for a given DGP must be active at (and/or potentially before) the corresponding time of year.
The principal physical characteristics of the spatial subdomains associated with each type and sign of trend were compared in cumulative spatial frequency curves, generated by ranking the corresponding subsets of DGPs by metrics representing their latitude, elevation, and mean day of snow-year (where 0 = 1 October) of SnowOn1 and/or SnowOn4.
To support the derivation of additional spatial statistics, the terrestrial area associated with each DGP was estimated by superimposing a layer of points drawn at 10-km intervals (identified previously as being located over land) over another containing polygonal features representing the surface footprints associated with the DGPs. Each of the raw DGP surface footprint areas (as supplied with the dataset) was then scaled by the corresponding fractional coverage of terrestrial points. This enabled illustrations integrating the spatial incidence and strength of trends to be generated, using the product of DGP terrestrial area (km2) and trend magnitude (days decade−1), termed “area-magnitude” plots (km2 days decade−1).
An additional aspect of interest was the range of climatological settings in which SnowOn1 and SnowOn4 trends were identified. To explore this, the CRU 1961–90 mean monthly precipitation and temperature within the spatial footprint associated with each DGP were retrieved for the month(s) containing that DGP’s 1972–2017 mean dates of SnowOn1 and SnowOn4. This permitted the preparation of cumulative spatial frequency curves describing the distribution of climatological characteristics among DGPs associated with each type and sign of trend.
d. Assessment of the temporal evolution of trends
The temporal evolution of shifting snow-onset dates giving rise to significant trends was explored using a simple form of changepoint analysis. While a range of relatively complex techniques has been developed for this purpose (e.g., Ducré-Robitaille et al. 2003; Reeves et al. 2007; Aminikhanghahi and Cook 2017), the large number of series to be considered here dictated that a more straightforward approach was required, and that described by Pettitt (1980) was selected. This method generates a cumulative summation (CUSUM) series Si from a time series of values Xi by adding successive anomalies relative to the series mean, thus: , where S0 = 0, and is the series arithmetic mean. In the resultant series, subsections with generally positive (negative) gradients correspond to periods during which the same sign of anomaly prevails, while approximately constant sections imply conditions close to the average: it follows that maxima and minima indicate changepoints from one sign of anomaly to the other, while inflections mark shifts in rates of change.
CUSUM series were generated for the time series of SnowOn1 and SnowOn4 dates associated with DGPs for which significant trends were identified, and for the series of monthly mean temperatures and precipitation depths (University of Delaware data) in the month containing the mean date of SnowOn1/SnowOn4. The month of mean onset was used because tracking values in the month of actual onset each year risks—for the example of temperature—jumping from a warmer (colder) month to a colder (warmer) one where a positive (negative) trend is identified. Series were also constructed from 3- and 5-yr running means of the CUSUM values. The overall form of these series, and the years in which their extreme values occurred, were used to identify pattern changes in snow-onset dates and the associated meteorological data.
a. Snow-onset trends
A SnowOn1 (SnowOn4) date was registered in each of the 46 years from 1972 to 2017 for 1228 (944) DGPs. Significant trends were identified for 364 (261) of the points in each of these subdomains, associated with approximately 29% (27%) of their respective areas (Table 1). The two signs of trend were found to occur within broadly distinct and contiguous regions: the overall spatial patterns were similar for both metrics, with relatively few exceptions and additions (Fig. 2). Variations of trend sign and magnitude with longitude, latitude, and elevation are depicted in Figs. 3–6. Their cumulative spatial frequency distributions in relation to latitude, elevation, mean onset date, mean monthly temperature, precipitation (in the month containing the mean onset date), and magnitude are illustrated in Fig. 7. Associations between magnitude and latitude, longitude, elevation, and mean onset date are shown in Fig. 8, with regression metrics in Table 2. Distributions of the product of area and trend magnitude with latitude, longitude, elevation, and mean onset date are provided in Fig. 9. The following sections summarize the distributions of positive and negative SnowOn1 and SnowOn4 trends in relation to these key contextual characteristics.
1) Significant SnowOn1 trends
SnowOn1 trends represent significant shifts in the date of the first transition from <50% to ≥50% snow cover, during the onset of cooler conditions in the latter stages of the calendar year. The summary spatial metrics shown in Table 1 indicate that significant negative trends (and thus progressively earlier snow onset) outnumber significant positive trends (later onset) by a ratio exceeding 5:1. When expressed in terms of their associated areas, this ratio rises to 6:1 (driven by more widespread occurrence of negative trends at higher latitudes, where the areas associated with individual DGPs are larger than at lower latitudes.)
The four main clusters of positive trends are located principally on the southwestern slopes of the Himalayan Massif northwestward to the Pamir, between the Tien Shan and Altai Mountains, in central to northern Norway, and along the Western Cordillera of North America, extending north and west to both sides of the Bering Strait (Fig. 2a). While they are distributed approximately evenly through NH latitudes (Fig. 7a), they are longitudinally much more limited, being associated primarily with the upper and western (and, at these latitudes, therefore windward) slopes of the Himalayas, central to northern parts of the Western Cordillera, and Norwegian mountains (Fig. 3). Note that 55% (33%) are found above elevations of 1000 m (1500 m) (Figs. 7b, 8e, 9c). Most of the lower-altitude DGPs experiencing positive trends are distributed on either side of the Bering Strait.
Negative trends are found over extensive swaths of northeastern and central Eurasia, central-eastern Europe, northeastern Canada, and western North America. The three largest spatially coherent clusters are found to the north and east of the Himalayas, European Alps, and Western Cordillera of North America, and occupy areas roughly proportional to the uplands’ lateral extents (Figs. 2a, 3, and 4). Smaller patches are seen in the eastern Himalayas, the Hindu Kush, and through Iran and Anatolia to the Balkans. There is a particularly clear association with the eastern slopes of major mountain ranges, except at relatively low latitudes in western North America, where they occur on the Pacific-facing upper western slopes of the coastal ranges (Fig. 3). They are more strongly associated with midlatitudes (40°–60°N) than those farther to the north or south (Figs. 7a and 9a). Notably, 79% (92%) of the associated area is found below 1000 m (1500 m) (Figs. 7b, 8e, 9c).
Snow onset in the areas affected by positive trends has historically occurred early in the seasonal transition, mainly within the first 2 weeks of October. In contrast, those experiencing negative trends have not generally seen their first major snowfall until sometime between mid-October and late November (Fig. 7e). In the month containing the mean date of initial snow onset, areas experiencing positive trends are warmer by ~5°C than those experiencing negative trends (Fig. 7c). Temperatures are below freezing over only ~53% of the area associated with positive trends, with 29% below −5°C and 5% below −10°C. In comparison, 86% of the area seeing negative trends is below freezing, with >60% below −5°C, and >24% below −10°C. Approximately 50% of the area affected by positive trends receives appreciably greater precipitation during the month of mean onset than those experiencing negative trends (Fig. 7d).
Positive trends have considerably larger absolute magnitudes than negative trends (Fig. 7f). Nearly 60% of the area associated with positive trends has been experiencing progressive delay of onset dates at rates greater than 5 days decade−1. Over 15% of the total affected area, the rate is in excess of 10 days decade−1, and over 3% it is greater than 20 days decade−1. Negative trends are more moderate, with 73% of the associated area seeing rates of change no greater than 5 days decade−1, and fewer than 1% above 10 days decade−1.
While the magnitudes of both signs of trend increase at higher elevations and lower latitudes, positive trends correlate more strongly with, and are more sensitive to, both apparent influences (Table 2). These associations persist within the meridional strips shown in Fig. 4, where sufficient DGPs with negative or positive trends exist (Table S1 in the online supplemental material). Much weaker relationships exist between climatologically cooler temperatures in the historical mean onset month and stronger negative trends (R = −0.15, p < 0.05, slope = −0.05 days decade−1 °C−1), and climatologically wetter conditions in the historical mean onset-month and stronger positive trends (R = −0.30, p < 0.05, slope = −0.03 days decade−1 mm−1; Fig. S1).
2) Significant SnowOn4 trends
SnowOn4 trends represent significant shifts in the date of the first transition from <50% to ≥50% snow cover, when this subsequently persists for at least 4 weeks. As established above, this approximates the onset of the longest continuous period of snow cover each winter. As evident from Table 1, the ratio of negative to positive trends is considerably higher than the equivalent for SnowOn1 trends, at 9:1 by the number of DGPs in each group, and 10:1 by area (again driven by more widespread occurrence of negative trends at higher latitudes, where the areas associated with individual DGPs are larger than at lower latitudes).
The overall spatial distribution of positive and negative SnowOn4 trends shadows that of the SnowOn1 trends, with several exceptions and one addition (Fig. 2b). There are considerably fewer instances of positive trends, with clusters involving more than a single DGP occurring only at higher elevations of the Himalayan Massif’s southwestern extremities, either side of the Bering Strait, and in a relatively small patch of central-western Russia, north of the Caspian Sea. Negative trends are also diminished in several regions, particularly in central-eastern Europe, northeastern Canada, and most of the more southerly parts of central-western North America affected by SnowOn1 trends. A possible anomaly exists in the Tien Shan (where mainly positive SnowOn1 trends were identified), in the form of a small group of negative SnowOn4 trends: together, these imply later initial snow onset, but earlier persistent onset, among DGPs in reasonable proximity to one another (although not coincident: no DGP yields trends of both types with opposite signs).
As positive SnowOn4 trends occur only within a few small patches, their latitudinal distribution appears as a series of discrete “steps” (Figs. 7a and 8b). The distribution for the negative trends is more continuous. A considerably larger fraction of the area associated with positive SnowOn4 trends occurs at lower elevations, compared to that for SnowOn1: this corresponds to the persisting incidence of positive trends either side of the Bering Strait. The distribution of elevations among negative SnowOn4 trends is identical to that for negative SnowOn1 trends. Thus, positive SnowOn4 trends are again more commonly associated with lower latitudes and higher elevations, while negative trends are found mostly at lower elevations in the midlatitudes. This is made additionally clear by the combination of area and magnitude in Figs. 9a and 9c.
It is noteworthy that the distributions of mean SnowOn4 onset date associated with the DGPs for which positive and negative trends were identified are essentially identical (Fig. 7e). This implies that—when considered over the entire domain—the influences driving both earlier and later arrival of more persistent snow cover are occurring at the same time of year in different areas.
The distributions of climatological mean precipitation expected during the mean onset-month are also identical for the DGPs associated with both signs of trend (Fig. 7d). Profiles of mean temperature at onset remain separated, with the areas seeing negative trends approximately 5°–7°C cooler than those with positive trends (Fig. 7c). Thus, among areas seeing positive trends, 86% have monthly mean temperatures in the mean month of the onset of persistent snow cover at or below freezing, 48% below −5°C, and 16% below −10°C, while for those with negative trends the corresponding values are 100% at or below freezing, 83% below −5°C, and 44% below −10°C. No significant correlations exist between either climatological metric and the magnitudes of either sign of SnowOn4 trend (Fig. S1).
The cumulative spatial frequency profiles of (absolute) trend magnitude are again similar, over the 80% of their respective associated areas with relatively weak values (~2–4 days decade−1) (Fig. 7f). However, the remaining 20% by area of positive trends are stronger than their negative equivalents, with ~13% greater than 7 days decade−1, and 10% above 9 days decade−1, the maximum being ~13 days decade−1. These extremes are all located at altitudes above 3500 m and latitudes south of 35°N.
3) Summary: Contextual comparison of SnowOn1 and SnowOn4 trends
Overall spatial patterns of the incidence of SnowOn1 and SnowOn4 trends are broadly similar, with positive trends associated with the upper western slopes of major mountain ranges, and negative trends found at more moderate elevations, almost always to the north and east of upland areas. Positive trends are more commonly associated with lower latitudes and higher elevations than negative trends, except in extreme eastern Siberia and western Alaska. Latitudinal distributions of both signs of SnowOn4 trend indicate slight northward shifts of a few degrees of latitude relative to those of SnowOn1 trends. During the month containing their respective mean onset dates, SnowOn4 trends of both signs are associated with generally cooler conditions (by ~5°C) than those seeing SnowOn1 trends, and somewhat (for negative trends) to considerably (for positive trends) drier conditions. The magnitudes of positive SnowOn1 trends are stronger than (the absolute values of) all other categories. Overall, magnitudes of both signs of trend strengthen at higher elevations, lower latitudes, and with later mean onset dates. These correlations are essentially identical for SnowOn1 and SnowOn4 trends of the same sign (Fig. 8, Table 2). The metrics plotted in Fig. 9, illustrating variation in both the rate of change in onset date and the area affected, further emphasize the distinctions between the two trend signs.
It is interesting to note that mean SnowOn1 dates among the DGPs experiencing progressively later onset of initial snow cover (positive SnowOn1 trends) are all earlier than those experiencing later onset of more persistent snow cover (positive SnowOn4 trends). However, the reverse is true among negative trends (i.e., mean onset dates of DGPs with negative SnowOn1 trends are later than those associated with DGPs having negative SnowOn4 trends). The patchy spatial distribution of positive trends severely limits the potential for inferences to be drawn about these from this observation. However, the more widespread incidence of negative trends is more useful, supported further by the strong metrics of correlation between latitude and onset date in the same subset of DGPs (Table 2). Thus, while there may be some possible sign of a link between progressively later onset (at high altitudes) each year in more southerly areas and later onset farther to the north, there is much stronger evidence that earlier onset at higher latitudes is followed by earlier onset at lower latitudes.
b. Are changepoints discernible in temporal trajectories of snow onset and meteorological conditions?
Where significant trends exist, it is useful to explore the timings of shifts from anomalies of one sign to the other within the associated series of snow-onset dates, and to compare these with patterns in meteorological data relating to the same areas. This was achieved using CUSUM series generated from the records of SnowOn1 and SnowOn4 onset dates, and also from series of mean temperature and precipitation (based on the datasets made available by Willmott and Matsuura of the University of Delaware) in each year during the month containing the mean onset date for the corresponding DGPs. In principle, any distinct extreme values within the CUSUM series mark changepoints in the overall trajectory of the metrics to which they relate, and thus provide information about the corresponding trends’ temporal development. [Note that meteorological records were not available for 10% (8%) of the DGPs for which negative SnowOn1 (SnowOn4) trends were identified, nor for 9% (12%) of those with the respective positive trends.] Given the number of points involved and the need for brevity, the comparison offered here is necessarily generalized, and it is not possible to consider the nuances of relationships between subsets of patterns among the various series. Nevertheless, it provides a foundation for subsequent more detailed analysis.
1) Snow-onset changepoints
Among the CUSUM profiles generated from snow-onset dates associated with significant negative SnowOn1 and SnowOn4 trends, the extreme values are all clear maxima, indicating a distinct temporal changepoint from positive to negative anomalies. These suggest that snow generally arrived later than the series mean until the mid-1990s, and earlier after that point (Figs. 10a,b, Figs. 11a,b). The modal years for the annual counts of peak 1-yr, 3-yr mean, and 5-yr mean CUSUM values for both SnowOn1 and SnowOn4 lie predominantly between 1994 and 1997: maxima occur before 1997 over approximately 70% (66%) of the total area in which negative SnowOn1 (SnowOn4) trends were identified (Fig. S4). Mapping the peak years, and plotting them against latitude, longitude, and elevation (Figs. S7–S10) reveals that the earlier SnowOn1 trends located inland from the northern Pacific coast of North America, in central-eastern Europe, and in isolated pockets of the Middle East and eastern Eurasia correspond to areas where the switch from positive to negative anomalies occurred in the 1980s. The comparatively small clusters of negative SnowOn1 trends situated on higher, west-facing slopes are members of this subset. Transitions from later to earlier onset-date anomalies throughout most of eastern Eurasia occurred during the 1990s, whereas in central and eastern North America they generally took place somewhat later, until the mid-2000s.
Overall, highly significant (p < 0.001) but weak (absolute values of R ~0.25–0.30) positive (negative) correlations exist between the year of the extreme 3-yr mean value of each CUSUM series and latitude (elevation) among DGPs with significant negative SnowOn1 trends, and with elevation only for those with negative SnowOn4 trends. However, stronger positive associations exist on each continental landmass between the year of peak 3-yr mean SnowOn1 CUSUM and longitude, with R = 0.52 in the Western Hemisphere and R = 0.36 in the Eastern Hemisphere (p < 0.0001 for both). These findings imply that (in highly simplistic terms) transitions to earlier onset developed initially at higher elevations and lower latitudes, and “migrated” downward and poleward in a generally northeasterly direction.
The extreme CUSUM values associated with DGPs manifesting significant positive trends are all (with the exception of a single SnowOn1 instance) minima, indicating clear shifts from earlier- to later-than-average onset. However, the timings of these transitions are generally less distinct and coherent than the equivalents for negative trends. The annual mean and 95th percentile curves show clear minima in 1998 for both SnowOn1 and SnowOn4 (Figs. 10b and 11b), and the modal frequencies of extreme CUSUM values in the 1-yr, 3-yr mean, and 5-yr mean series coincide with these (Fig. S4). However, between 1985 and the early 2000s the curve representing the median of annual SnowOn1 CUSUM values is more or less flat (indicating that onset dates during this interval were near the series mean), suggesting that outliers during this period may be exaggerating the means. Given that nearly half of the 1-yr, 3-yr mean, and 5-yr mean CUSUM changepoints fall between 1997 and 2001 (Fig. S4), and in light of the temporal proximity of these years to the introduction of major technological improvements to the NRSA with implementation of the IMS in 1999, this in particular demanded investigation of associated patterns of temperature and precipitation.
2) Temperature changepoints
Trajectories of temperature during the period of interest in the areas associated with DGPs yielding significant SnowOn1/SnowOn4 onset-date trends were explored using CUSUM series based on mean temperatures in the months containing the corresponding mean onset dates. (Note again that this approach was preferred over using the month of actual onset in each year, as there are considerable differences in monthly means, so using these would have exaggerated the results.)
The extreme values of these CUSUM series were found to be minima for ~93% of those associated with both SnowOn1 and SnowOn4 negative (i.e., earlier) snow-onset trends, and ~85% of those having positive (later) SnowOn1/SnowOn4 trends. This suggests that the great majority of areas seeing all four categories of trend have experienced distinct shifts from negative to positive temperature anomalies, and thus overall warming, during the period of interest (Figs. 10c,d and 11c,d). While it is perhaps somewhat counterintuitive that this applies to points for which negative, as well as positive, snow-onset trends were identified, there are distinct contrasts between their detailed trajectories.
Among those experiencing negative trends, extreme minima occurred principally in three main temporal clusters, in the latter 1980s (more dominant for points having SnowOn4 than SnowOn1 trends), latter 1990s, and early 2000s (Fig. S5). Over approximately 40% (55%) of the total area experiencing negative SnowOn1 (SnowOn4) trends, temperature anomalies switched from negative to positive before 1990. Warming progressed at a relatively slow rate from 1987 to 2003, although it was faster (apparently with a brief hiatus in the late 1990s) for areas seeing SnowOn4 than SnowOn1 trends. Since the early 2000s, warming has accelerated considerably in these areas. When compared with the corresponding snow-onset CUSUM series (Figs. 10a and 11a), the peak of the latter may be seen to occur approximately 5–10 years after the first major shift from negative to slightly positive temperature anomalies.
Among DGPs for which positive trends were identified, a tighter grouping of changepoints from negative to positive temperature anomalies is evident through the latter 1990s (slightly earlier for SnowOn4) to early 2000s, with a sharp minimum in 2001. This follows the clear minimum in the snow-onset CUSUM series (Figs. 10b and 11b) by 3–5 years. Mean overall cumulative cooling until ~2000 in areas associated with positive SnowOn4 trends was approximately twice that for DGPs having negative SnowOn4 trends (Figs. 11c,d). These profiles demonstrate that areas seeing later onset, although generally warmer than those with negative trends, experienced a longer initial period of negative temperature anomalies, followed by a sharper transition to the rapid rates of warming observed since ~2000.
3) Precipitation changepoints
The CUSUM series generated from records of mean monthly precipitation in the mean onset month within each category of trend display considerably greater spreads of annual values between 1972 and 2017, and comprise more subdecadal variability in the implied anomalies, than those seen in the equivalent series for snow onset and temperature (Figs. 10e,f and Figs. 11e,f). This translates into less temporal consensus in the incidence of changepoints, which often involve a mixture of minima and maxima of appreciable magnitudes in individual series (Fig. S6). (As is so often the case, information is unfortunately not available to describe the precipitation phase.)
However, among DGPs experiencing negative trends and thus progressively earlier SnowOn1 (SnowOn4) dates, the extreme values of their precipitation CUSUM series are minima for 62% (68%) of the total number of each type of trend, suggesting an overall shift from mainly negative to positive anomalies, and thus a general increase in precipitation during the period of interest. Where positive trends were identified, the extreme values for 52% (64%) of the series associated with SnowOn1 (SnowOn4) trends are maxima, implying the opposite (Fig. S6).
Generally, areas associated with DGPs for which negative SnowOn1 and SnowOn4 trends were identified exhibit mainly negative precipitation anomalies until the early 1980s, followed by 15–20 years of comparatively low positive anomalies (greater, but more variable, for SnowOn4 than SnowOn1), before more rapid increases during the 2000s. Three principal clusters of minima are evident, in the late 1980s, mid-1990s, and late 2000s (Fig. S6). These correspond reasonably closely (but not precisely) with the years in which minima occur in the temperature CUSUM series associated with DGPs having negative snow-onset trends (Fig. S5). The first two clusters occur during the period of positive snow-onset anomalies, and the third aligns with a small cluster of maxima in the latter 2000s (Fig. S4).
The relationship between trajectories of precipitation and positive snow-onset trends is the least distinct of all combinations considered. Nevertheless, the excess of maximum extreme values over minima suggests an overall switch from positive to negative anomalies, and thus a slight overall drying tendency. The plots of mean and median CUSUM values (Figs. 10f and 11f; see also Figs. S2f and S3f) imply a transition from moderate negative anomalies in the early 1970s to increases until the late 1990s, followed by a steep drop to 2013. The local maximum around 1998 corresponds closely with the clear minimum in the snow-onset CUSUM series (and is also near the switch from negative to positive temperature anomalies) associated with DGPs having positive trends. Although only relatively few DGPs are involved in these subsets, it would be useful, given the complexity of the underlying fine detail, to explore these patterns further.
5. Discussion and conclusions
This analysis explored monotonic trends in the onset date of seasonal snow cover across the NH between 1972 and 2017. The first major snowfall of the season was termed SnowOn1, and the first persisting subsequently for at least four weeks was called SnowOn4. Widespread statistically significant negative trends, indicating shifts toward progressively earlier onset dates, were identified, while positive trends (shifts to later onset) were fewer in number and more localized. Negative SnowOn1 (SnowOn4) trends were found to occur over areas larger than those associated with positive trends by a factor of 5 (9). Trends of each sign were associated with clearly distinct contexts, in terms of their physiography, latitude, longitude, and climatology.
Perhaps the most striking spatial pattern is the clear longitudinal alternation between areas affected by negative and positive trends (Figs. 3, 4, 8c,d, and 9b). Those experiencing a delay in snow onset are almost all located to the west (often southwest) of appreciable uplands, in northwestern North America, Scandinavia, to the immediate southwest of the Urals, and in southwestern parts of the Himalayan Massif. Meanwhile, earlier onset is associated largely with areas of low-to-moderate elevation to the north and east of major uplands between 40°N and 60°N. The three principal examples of the latter (in decreasing order of extent) extend from the Tien Shan, through the Altai, Sayan, Yablonoi, and Stanovoy chain to northeastern Siberia; from the Rocky Mountains through central-western to central-northern North America; and from the Alps through central-eastern Europe. The northeastward continuation of each of these alignments leads to a circumpolar marine area in which autumn sea ice extent has been receding particularly rapidly, namely the Chukchi Sea (Serreze and Stroeve 2015), Hudson Bay (Andrews et al. 2018), and the Barents/Kara Sea (Wegmann et al. 2015; Onarheim and Årthun 2017), respectively. Indeed, the largest increases in ocean-to-atmosphere surface turbulent heat fluxes in the entire Arctic have occurred in the first and last of these (Cohen et al. 2018).
Together with latitude, elevation, aspect, and proximity to an Arctic coastline, the other key modulating influence on trend development is the spatiotemporal evolution of atmospheric circulations, and their interactions with surface relief and other characteristics. While detailed consideration of these factors lay outside the scope of this study, it is noteworthy that the longitudinal positions of positive trends identified here coincide closely with those in which persistent blocking systems have been found to occur relatively frequently from September to November, while negative trends align with areas immediately to the east of these meridional bands (e.g., Barriopedro et al. 2006; Tyrlis and Hoskins 2008; Cheung et al. 2013). The pattern described here over eastern Eurasia conforms closely with current consensus on the most robust mechanism by which amplification of climate change in the Arctic is affecting extreme meteorological events at midlatitudes—namely, that ice loss in the Barents and Kara Seas is contributing to the intensification of the Siberian high and consequent troughing to its east, resulting in additional humidity released from increasingly ice-free areas of the Arctic Ocean falling as snow farther to the south (Cohen et al. 2018). Given the juxtaposition of climatologically frequent anticyclonic conditions and more open Arctic waters “upstream” of the two other largest regions affected by earlier snow onset (to the north of Canada and central Europe), it is reasonable to speculate that comparable processes may also explain these trends. Together with the generally cool climatological settings in which they occur, and the potential for some degree of orographic enhancement of precipitation by uplands farther to the south and/or southwest in each instance, these factors help to explain how snow could have been arriving earlier in the year over these areas. Meanwhile, the positive trends could relate to a reduction in the availability of precipitable water, as systems are prevented from traveling eastward by the blocking systems themselves, and/or the northward advection of warmth from areas farther to the south. It would thus be useful to explore in more detail how these circulations have evolved within the areas affected by both signs of trend over the same period, particularly in relation to any changes in the spatial patterns of dominant pressure systems.
A key finding pertinent to the wider discussion relating to the credibility of the NRSA is that approximately 70% (66%) of DGPs associated with significant negative SnowOn1 (SnowOn4) trends had shifted from positive to negative snow-onset anomalies by 1997, 2 years before introduction of the IMS and imagery of higher spatial and temporal resolutions. These technical improvements would have resulted in the more distinct resolution of snow-covered from snow-free land than the prior system, particularly over areas with greater spatial and/or temporal frequency of variation in snow cover, such as in complex mountainous terrain or at the fringes of snow-prone areas. It might therefore be expected that snow-onset changepoints around 1999 would cluster preferentially over such settings. However, those DGPs for which changepoints between 1997 and 2000 were identified are distributed predominantly over the interior shield regions of North America and Eurasia (Fig. 12). Of the relatively few exceptions, most—particularly in the southwest Himalayan ranges and central Rocky Mountains—are positive trends. While these are less contentious in the context of widespread warming, it would be informative to explore their associations with shifting patterns of precipitation, temperature, and meteorological systems in more detail.
While the temporal trajectories of snow-onset anomalies in each subset of trends have been shown to correspond generally to shifts in temperature and precipitation anomalies, this analysis does not claim to offer a conclusive explanation for each individual DGP’s negative or positive onset trend. The association between negative onset trends and switches from mainly negative to positive temperature anomalies is also counterintuitive, particularly given current understanding of the sensitivity of snow cover to temperature (e.g., Mudryk et al. 2017). However, these results offer the basis for a possible explanation, by demonstrating that the areas affected are found in relatively cool climatological settings, immediately to the east of regions in which blocking systems develop frequently (and which are therefore prone to troughing), and are set between patches of the southern Arctic Ocean in which freeze-up has been delayed to the northeast, and uplands of appreciable height and extent to the south and/or southwest. Together, these influences appear to have resulted in a moderate increase in precipitation falling as snow earlier in the transition from late summer to winter. It seems probable, though, given the accelerated rate of warming observed since the early 2000s, that this augmented early snowfall will change to rain in the near to medium term. In contrast, some combination of a later but sharper shift to warming, together (in some areas) with a downward swing in precipitation, led to snowfall arriving later in the year over DGPs associated with positive onset trends, largely within regions prone to frequent and/or persistent anticyclonic conditions at the same time of year.
It is again emphasized, though, that the extreme interannual variability of the precipitation patterns severely limits the reliability of the conclusions drawn here. A useful approach for more detailed analysis might be to identify clusters of DGPs with broadly similar CUSUM profiles generated from onset dates, temperature, and precipitation on the basis of curve shape, using an appropriate multidimensional classification method. This would in turn support examination of relationships between the various changepoints within subsets having comparable trajectories (and thus, notwithstanding equifinality, exogenous influences), and potentially also in more localized spatial contexts. It would also be useful to identify and quantify temperature and precipitation trends within the same spatial framework as that on which the onset trends were based, and to compare temporal patterns of Arctic sea ice (particularly in the areas of key interest) in the months corresponding and prior to mean snow onset with those of the onset dates.
While these results seek to augment understanding of trends in the onset of NH snow cover derived from the NRSA by identifying their temporal trajectories and relating them to meteorological activity, other explanations must yet be considered. For example, widespread clearances over large areas of the Russian and North American boreal forests during recent decades—caused by fires, major pest outbreaks, and the resource-extraction industries [such as those described since 2000 by Hansen et al. (2013)]—will have resulted in a proliferation of more open areas. It seems possible that snow may previously have existed at a given time of year in such areas, but was previously hidden by the dense needleleaf forest (once it had been shed from the canopy). The new clearings may have provided more visual cues for the IMS, potentially resulting in the identification of such pre-existing snow cover as “new.” Therefore, if—within the footprint associated with a given DGP for which a negative snow-onset trend was identified—the year in which substantial clearance occurred is found to correspond closely with a shift in its onset-date anomalies, this would warrant further investigation. On the other hand, it is at this stage by no means clear that the spatial dimensions of these clearances would be sufficiently large for any newly revealed snow to be detectable within the imagery from which the NRSA is derived.
This study thus yields a range of new information relating to snow-onset trends across the NH, and relates these details to the physiographic and climatological characteristics of the affected areas. Although no claims are made that these findings represent conclusive explanations for either sign or type of the onset-date trends derived from the NRSA, it is hoped that these new perspectives on their spatial and temporal incidence and evolution will provide a useful contribution to the continuing exploration and discussion of this topic.
This study was inspired by and extends prior research supported by the Canadian Sea-Ice and Snow Evolution Network (principal investigator Dr. Paul Kushner), funded by the Climate Change and Atmospheric Research initiative of the Natural Sciences and Engineering Research Council of Canada. The authors are particularly grateful to Dr. David Robinson and Mr. Thomas Estilow of Rutgers University for supplying the NH snow cover Climate Data Record, and for their helpful technical assistance. They also thank the University of East Anglia’s Climatic Research Unit and the University of Delaware’s Department of Geography for making their climate datasets publicly available. Thanks are also due to Dr. Richard Kelly and Dr. Chris Fletcher of the University of Waterloo for their helpful comments relating to this work, and to three anonymous reviewers for their insightful and constructive contributions. M.I.A. also gratefully acknowledges the financial support provided by the University of Northern British Columbia in the form of a Doctoral Dissertation Completion Award.
Supplemental information related to this paper is available at the Journals Online website: https://doi.org/10.1175/JCLI-D-18-0686.s1.