Regionalization of snowfall frequency and trends over the contiguous United States

INTERNATIONAL JOURNAL OF CLIMATOLOGY Int. J. Climatol. 35: 4348–4358 (2015) Published online 19 February 2015 in Wiley Online Library (wileyonlinelibrary.com) DOI: 10.1002/joc.4292 Regionalization of snowfall frequency and trends over the contiguous United States Daria Kluvera* and Daniel Leathersb a Department of Earth and Atmospheric Sciences, Central Michigan University, Mount Pleasant, MI, USA b Department of Geography, University of Delaware, Newark, DE, USA ABSTRACT: This study examines the regional variations in the frequency of snowfall across the conterminous United States from 1930 to 2007. Principal components analysis and cluster analysis are used to group stations together based on the main modes of variation in snowfall frequency. Results indicate the existence of seven unique snowfall regions, which correspond to predominant storm tracks across the United States. These are the southeast, the south central Plains and southwest, the Ohio River Valley and mid-Atlantic, the Pacific Northwest, and three sub-regions in the Upper-Midwest. Quantile regression reveals that the distribution functions of each region’s snowfall frequency are different and in some regions, changing over time. The northern part of the Upper-Midwest is experiencing increasing trends in all percentiles of snowfall frequencies, the Pacific northwest is experiencing declines in greater than median snowfall frequencies, and the southeast is seeing a decline in extreme frequency years. Correlation analysis between large-scale teleconnection patterns and regionally averaged snowfall frequencies corroborate previous findings and indicate specific forcing mechanisms for snowfall frequency in each region. KEY WORDS snowfall; snow frequency; regionalization; principal components analysis; cluster analysis; quantile regression; trends Received 16 September 2014; Revised 12 January 2015; Accepted 19 January 2015 1. Introduction Snowfall, as a variable to study, has practical importance because of its immediate impact on human activities, such as transportation (Changnon, 1979; Hanbali and Kuemmel, 1993; Schmidlin, 1993; Norrman et al., 2000; Strasser, 2008). It is also responsible for the formation of a snow pack, which alters the energy balance (Gray and Male, 1981) and water budget (Hartmann, 1994; Dingman, 2008) of the climate system. Fundamentally, total snowfall varies through time due to changes in amount per event (based on temperature and moisture availability) and frequency of events. Examination of either of these components of total snowfall will shed light on the causes of change over time. However, snowfall frequency is more directly influenced by large-scale atmospheric forcing mechanisms via changes to storm track positions. To better understand the impacts of these large-scale forcing mechanisms on snowfall and therefore better predict its occurrence, it is important to understand the spatial distribution of snowfall events within the climate system. In this study a snowfall frequency regionalization is presented to contribute towards that goal. Several earlier studies have examined the spatial variability in United States snowfall, but are either not current or not spatially extensive enough for a robust snowfall * Correspondence to: D. Kluver, Department of Earth and Atmospheric Sciences, Central Michigan University, 314 Brooks Hall, Mount Pleasant, MI 48859, USA. E-mail: kluve1db@cmich.edu © 2015 Royal Meteorological Society regionalization (Harrington et al., 1987; Leathers et al., 1993; Groisman and Easterling, 1994; Hughes and Robinson, 1996; Hartley and Keables, 1998; Serreze et al., 1998; Smith and O’Brien, 2001; Bradbury et al., 2003; Patten et al., 2003; Morin et al., 2008; Ghatak et al., 2010). Some of these studies have described snowfall regions based on harmonic analysis (Harrington et al., 1987) or on principal component loading patterns (Leathers et al., 1993; Hughes and Robinson, 1996). Notable among these, Leathers et al. (1993) used principal components analysis to identify regions in which snowfall varied in a coherent manner from 1945 to 1985 at stations east of the Cascades. The resulting regions were the Great Lakes/Upper Midwest region, the Central Plains and Southern Rockies, the Eastern Mid-Atlantic region through Southern New England, the Southern Mid-Atlantic to Central Plains, the Northern Mid-Atlantic region, and New England. The changes in each region were attributed to the movement of storm tracks and frequency of synoptic patterns. The importance of cyclone tracks on snowfall in the United States has also been documented in earlier studies and highlights a mechanism by which larger-scale climate variability in the form of teleconnection patterns can influence regional snowfall (Kunkel and Angel, 1999; Bradbury et al., 2002; Bradbury et al., 2003; Changnon et al., 2008). Therefore, to delineate areas based on the large-scale circulation pattern responsible for the snowfall regime, snowfall frequency, rather than total snowfall amounts, is used in this study to characterize the homogeneous snowfall regions. REGIONALIZATION OF SNOWFALL FREQUENCY OVER THE CONTIGUOUS UNITED STATES 4349 Figure 1. Location of 440 USHCN stations from Kunkel et al. (2009) and their average seasonal frequency of events greater than or equal to 2 inch. The regionalization in this study is accomplished by using Principal Components Analysis (PCA) and Cluster analysis. PCA reduces the data into their main modes of variation; resulting clusters isolate spatially coherent regions in which snowfall frequency varies in a similar manner; and snowfall frequency in each region is examined through Quantile Regression, Correlation analysis, and composite analysis. The snowfall data and methods used in this study are described in Section 2. In Section 3 the results are presented and discussed. A summary and conclusions are presented in Section 4. 2. Data and methods The snowfall data used in this study are a subset of 440 stations from the United States Historical Climatology Network. These data are described in Kunkel et al. (2009) as high quality and were deemed suitable for trend analysis after undergoing inspection by a panel of snow data experts. With the exception of the area of North Dakota, South Dakota, Wyoming, and Montana, there is even spatial coverage of snowfall stations across the United States (Figure 1). Daily data are used over the snow seasons 1930/31 to 2006/2007, where the snow season is defined as 1 July to 30 June. Both Leathers et al. (1993) and Bradbury et al. (2003) indicate in their results that regions highlighted by the main principal components of total snowfall amounts are attributable to the location of common storm tracks. To more specifically define the snow regions by storm tracks, the current study uses the frequency of snowfall events rather than snowfall amounts. A threshold is set for the minimum size of snowfall events included in the analysis, with the purpose of considering only events that are large enough to cause disruption to normal human activities and to elicit road maintenance actions in the affected regions. As the levels of service vary among states and between countries and are often not explicitly defined, © 2015 Royal Meteorological Society events that are greater than or equal to 2 inch are used in this study, as these disrupt driving and require action by transportation authorities (Bremner, 1977; Gray and Male, 1981; Katko, 1993). These data are standardized (converted into standardized anomalies or z-scores) by z= x−x sx where x is the annual snowfall frequency, x is the annual snowfall frequency mean from 1930 to 2007, and sx is the standard deviation. This is done to remove the influence of each station’s mean and standard deviation by reducing the data to anomalies from a mean of zero, with a standard deviation of one (Wilks, 2006). PCA (or Empirical Orthogonal Function Analysis) is used to identify the main modes of variation among the 440 snowfall stations used in this study. PCA is a data reduction method that uses eigenvectors of the covariance matrix to define the patterns of simultaneous variation within the dataset. This method is considered a standard multivariate analysis technique in climate science and has also been used specifically in snow studies (Leathers et al., 1993; Cayan, 1996; Hughes and Robinson, 1996; Serreze et al., 1998; Frei and Robinson, 1999; von Storch and Zwiers, 1999; Derksen et al., 2000; McCabe and Dettinger, 2002; Bradbury et al., 2003; Jin et al., 2006; Wilks, 2006; Sobolowski and Frei, 2007; Morin et al., 2008; Ge et al., 2009). As the snow frequency data are standardized, the PCA is applied to the covariance matrix of the values, as opposed to using the correlation matrix with non-standardized data (Navarra and Simoncini, 2010). The resulting orthogonal component vectors are not rotated, to preserve as much of the variance explained as possible (Jolliffe, 1989; von Storch and Zwiers, 1999). The Principal Component score time series at each station are used in a cluster analysis to identify unique clusters of stations that vary in the same manner (Hughes and Int. J. Climatol. 35: 4348–4358 (2015) 4350 D. KLUVER AND D. LEATHERS Robinson, 1996). Within-groups clustering is used because it optimizes homogeneity within the clusters by minimizing within-cluster variance. Several cluster solutions are tested. When more than seven clusters are allowed, no new snowfall regions are identified. Regional time series of the resulting snowfall frequency clusters are examined using Quantile Regression. Quantile Regression is a method originally developed in the field of econometrics (Koenker and Bassett, 1978; Koenker and Hallock, 2001) for estimating conditional quantile functions and to examine how an entire distribution changes, rather than just the mean as is done in Least-Squares linear regression. It has been employed a number of times in climate research (Baur et al., 2004; Barbosa, 2008, Elsner et al., 2008; Timofeev and Sterin, 2010, Lee et al., 2013; Tareghian and Rasmussen, 2013). Similar to Least-Squares linear regression, QR solves for equations of the form p p Y (p|x) = 𝛽0 + 𝛽1 x + 𝜀 where a vector, 𝛽 is found with coefficients for each quantile (p). 𝜀 is the error term, with the expectation of zero. Unlike the Least-Squares method, which minimizes the squared errors, n ∑ ( )2 yi − 𝜇 , min i=1 where (yi − 𝜇) are the residuals, QR requires the minimization of the sum of absolute residuals, but with a weighting function 𝜌𝜏 , n ∑ ( )| | min 𝜌𝜏 |yi − ̂ yp xi |, | | i=1 where ŷp (xi ) is the model for the pth quantile of Y. The weighting function, 𝜌𝜏 , is (1 − p) if the observation is less than the quantile regression line and is weighted by p if above the line. This ensures that observations above the quantile regression line have positive residuals, while those below the line have negative residuals. There are no assumptions on the error distribution. Not only does QR provide information on the changes in the entire distribution of the data, but it is also non-parametric and less vulnerable to influence by outliers compared to Least-Squares linear regression. Regional snowfall frequency is correlated with several December, January, and February teleconnection indices to verify relationships among these patterns and winter storm tracks bringing snowfall to each region. The teleconnection indices used are selected due to their impact on either the location of precipitation or temperature anomalies that influence snowfall in the United States. They are also selected because they were identified by Kluver and Leathers (2014) as having significant correlations with individual snowfall stations. Included in the analysis are: the Arctic Oscillation (AO) (Thompson and Wallace, 2000), the North Atlantic Oscillation (NAO) (van Loon and Rogers, 1978; Wallace and Gutzler, 1981; Hartley and Keables, 1998; Bradbury et al., 2002; Morin et al., 2008; Climate Prediction Center, 2010; Ghatak et al., 2010), Nino 3.4 region sea surface temperatures (Rasmusson and Wallace, 1983, Ropelewski and Halpert, 1986, Kunkel and Angel, 1999; Smith and O’Brien, 2001; Patten et al., 2003; Climate Prediction Center, 2010), the Pacific North American index (PNA) (Wallace and Gutzler, 1981, Leathers et al., 1991; Serreze et al., 1998; Notaro et al., 2006; Coleman and Rogers, 2007; Morin et al., 2008; Climate Prediction Center, 2010; Abatzoglou, 2011), and the Pacific Decadal Oscillation (PDO) [Mantua and Hare, 2002; Joint Institute for the Study of the Atmosphere and Ocean (JISAO), 2010]. Northern Hemisphere annual temperatures are also considered (Jones et al., 2013) to identify areas where a change in temperature during snowfall events is changing precipitation from the solid to liquid phase (Hartley, 1996; Serreze et al., 1998). Spearman Rank Correlation is used to estimate the strength of the relationship as it is a non-parametric measure (Wilks, 2006). Composite analysis is conducted using seasons (November through April) with snowfall frequencies greater than positive one standard deviation (+1 s.d.) and smaller than negative one standard deviation (−1 s.d.) from the mean for each region. Variables examined are 500 hPa Figure 2. Results of cluster analysis on principal components score time series of snowfall frequency. © 2015 Royal Meteorological Society Int. J. Climatol. 35: 4348–4358 (2015) REGIONALIZATION OF SNOWFALL FREQUENCY OVER THE CONTIGUOUS UNITED STATES (a) (b) (c) (d) (e) (f) 4351 (g) Figure 3. Time series of regional average standardized frequency (events 2 inch or greater) for each of the seven regions given in Figure 2. Quantile Regression lines for the 10th, 25th, 50th, 75th, and 90th percentile are overlaid on the graphs. Bold dashed lines indicate statistically significant (p ≤ 0.05) trends. © 2015 Royal Meteorological Society Int. J. Climatol. 35: 4348–4358 (2015) 4352 D. KLUVER AND D. LEATHERS geopotential height (m) anomalies, Sea Level Pressure (SLP) (hPa) anomalies, and Sea Surface Temperature (∘ C) anomalies available from the NCEP/NCAR Reanalysis (Kalnay et al., 1996). As colinearlity is a common problem in climatological analyses, an evaluation of conditions during extreme snowfall frequency years can aid with interpretation. 3. 3.1. Results 3.2. Regional trends in snowfall frequency Regionalization A cluster analysis is performed on the first 8 principal components, which explain 46.2% of the variation in snowfall frequency. They are selected via North’s rule of thumb (Von Storch and Zwiers, 1999). The resulting snowfall frequency clusters are shown in Figure 2. The seven snowfall regions are: the southeast (Region 1, shown in red plus signs), the south central plains and southwest (Region 2, shown in orange x’s), the Ohio River Valley and mid-Atlantic states (Region 3, shown in yellow boxes), the pacific Northwest (Region 4, shown in green triangles), and the Upper Midwest which is separated into three clusters. The Upper Mid-west region is sensitive to the position of the ‘Alberta Clipper’ type storm tracks from the northwest (Thomas and Martin, 2007) and ‘Colorado Lows’ coming from the southwest (Changnon et al., 2008). Stations impacted by the more northward tracks create a cluster in the Dakotas, Minnesota, and Wisconsin (Region 7, shown in purple filled circles). Another cluster is made up of stations that receive more frequent snowfall associated with southerly storm tracks such as ‘Colorado Lows’. This region extends from Iowa, Indiana, southern Wisconsin and Michigan to northern New England (Region 6, shown in dark blue open circles). A third cluster is formed between these two in Nebraska, Iowa, southern Minnesota and Wisconsin (Region 5, shown in light blue stars). The other clusters also clearly show spatial footprints tied to the typical tracks of low pressure systems, such as the Ohio River Valley and central Great Plains snowfall associated with ‘Colorado Lows’, and the southeast with ‘Nor’ Easters’ and ‘Gulf Lows’ (Whittaker and Horn, 1981). (a) Even though the current study is based on a larger data set (both spatially and temporally), the seven regions are fairly consistent with previous findings that use other snow variables. In particular, Region 6, which extends from Iowa to New England matches the Leathers et al. (1993) PC 6 based on mean snowfall amounts. Also, splitting the Upper-Midwest into multiple regions is similar to the results from Hughes and Robinson (1996) for snow cover duration in the Great Plains. (b) Each of the seven regions is examined by plotting the time series associated with the average regional snowfall frequency (Figure 3). Quantile regression analysis is conducted on the time series, and the 10th, 25th, 50th, 75th and 90th percentile regression lines are overlaid on the graphs. Trends in the 50th percentile (or the median) are also known as the least absolute-deviations regression, and show that only the Pacific Northwest (Region 4, decreasing) and the northern Upper Mid-west (Region 7, increasing) have statistically significant (p ≤ .01) trends in the frequency of snowfall events over the period of record. Both these regions correspond to large areas of homogeneous station trends of total seasonal snowfall found in Kunkel et al. (2009). However, if the trends in the rest of the data distribution are considered, there is a more complete picture of how snow frequency is changing which would be missed by only examining the mean or median. Regions that display statistically significant trends in any of the quantiles are also plotted as process-diagrams in Figure 4. This type of figure shows the quantile regression results for taus = 0.1 to 0.9 at 0.1 increments on the x-axis, and the trends in standard deviations per year on the y-axis. The grey shading around the trend values is the band of 90% confidence computed via bootstrap estimates of standard error. Figure 5 displays box and whisker plots of the estimated data distribution based on the QR results, plotted for 1930 and 2007 to help visualize the change in the distribution over the period of record. Region 1 has a statistically significant (p < 0.01) decreasing trend in the 90th percentile, indicating a decline (c) Figure 4. Process-diagram of quantile regression trend coefficients. 90% confidence band generated using bootstrap method. © 2015 Royal Meteorological Society Int. J. Climatol. 35: 4348–4358 (2015) 4353 REGIONALIZATION OF SNOWFALL FREQUENCY OVER THE CONTIGUOUS UNITED STATES (a) (b) (c) Figure 5. Estimated box-and-whisker plots of 1930 compared to 2007 based on the quantile coefficients from the quantile regression. Whiskers denote the estimated 10th and 90th percentile values. in the number of seasons exceeding that quantile. The process-diagram of the regression trend coefficients is given in Figure 4. For Region 1, the last black dot on the diagram represents the slope of the tau = 0.9 quantile regression line. The trend of this line is −0.015 standard deviations per year, or −1.5 standard deviations per 100 years. The Pacific Northwest (Region 4) not only has significant declines in the median snowfall frequency over time, but the 75th and 90th quantiles also show statistically significant (p < .05) decreasing trends. These negative trends at several quantiles can be seen in Figures 3 and 4 where quantiles greater than or equal to the 50th percentile are statistically significant. This results in a narrowing of the distribution, which can be seen in the theoretical box-and-whisker plots in Figure 5. The 1930 inter-quartile range is 1.21851 standard deviations and decreases dramatically to 0.47623 standard deviations in 2007. This results in a snowfall frequency that would have been ‘average’ at the beginning of the period of record ranking as an ‘extreme’ event at the end of the period of record. Region 7, defined by the Dakotas, northern Minnesota, and northern Wisconsin, has experienced an almost complete shift in the snowfall frequency distribution. Figure 3 shows that all percentiles have statistically significant (p < 0.01) increasing trends. The process-diagram of quantile coefficients in Figure 4 shows that trends are positive at all quantiles. The result of these trends on the distribution can be seen in the box and whisker plots for 1930 and 2007 in Figure 5. As the trends are larger at higher quantiles, the range has increased with time. The inter-quartile range doubles from 0.48 in 1930 to 0.92 in 2007. Using these estimated values, it can be seen that what is a 90th percentile frequency snowfall year in 1930 falls below the 10th percentile of the distribution in 2007. 3.3. Atmospheric forcing mechanisms In order to determine which atmospheric forcing mechanisms are associated with annual snowfall frequency, Table 1. Spearman correlation coefficients for regional average snowfall frequency and monthly teleconnection indicesa. December NAO January NAO February NAO December PDO January PDO February PDO December PNA January PNA February PNA December Nino 3.4 January Nino 3.4 February Nino 3.4 December AO January AO February AO Annual NH temperature a Region 1: Southeast Region 2: south central Plains and Southwest Region 3: Ohio River Valley and mid-Atlantic Region 4: Pacific Northwest Region 6: IA, IN, WI, MI and New England −0.182 (0.111) −0.331 (0.003) −0.375 (0.001) 0.093 (0.419) 0.085 (0.461) 0.098 (0.397) −0.236 (0.074) −0.142 (0.289) 0.034 (0.802) 0.160 (0.161) 0.153 (0.185) 0.152 (0.186) −0.243 (0.066) −0.423 (0.001) −0.374 (0.004) −0.328 (0.003) 0.047 (0.680) −0.048 (0.681) −0.223 (0.051) 0.218 (0.055) 0.100 (0.389) 0.057 (0.625) −0.135 (0.314) 0.077 (0.567) 0.127 (0.341) 0.231 (0.042) 0.264 (0.020) 0.305 (0.007) 0.084 (0.531) 0.064 (0.634) −0.235 (0.076) −0.068 (0.554) −0.210 (0.065) −0.168 (0.143) −0.132 (0.253) 0.121 (0.291) 0.177 (0.124) 0.128 (0.266) −0.046 (0.730) 0.002 (0.990) 0.076 (0.573) −0.158 (0.168) −0.170 (0.139) −0.157 (0.172) −0.369 (0.004) −0.294 (0.025) −0.235 (0.076) −0.084 (0.464) −0.115 (0.316) 0.033 (0.773) −0.124 (0.284) −0.364 (0.001) −0.518 (0.000) −0.584 (0.000) −0.540 (0.000) −0.717 (0.000) −0.487 (0.000) −0.381 (0.001) −0.374 (0.001) −0.363 (0.001) 0.003 (0.982) 0.067 (0.618) −0.109 (0.417) −0.433 (0.000) −0.195 (0.088) −0.177 (0.123) −0.063 (0.588) −0.182 (0.112) −0.156 (0.175) −0.255 (0.025) −0.296 (0.024) −0.388 (0.003) −0.387 (0.003) −0.308 (0.006) −0.279 (0.014) −0.267 (0.019) −0.224 (0.091) −0.105 (0.433) −0.080 (0.552) −0.188 (0.099) Two tailed p values given in parentheses and bold values indicate statistically significant. © 2015 Royal Meteorological Society Int. J. Climatol. 35: 4348–4358 (2015) 4354 D. KLUVER AND D. LEATHERS (a) (b) Figure 6. Composite anomaly of Sea Level Pressure for Region 1 for November through April snowfall frequencies (a) smaller than 1 s.d and (b) larger than 1 s.d. From NCEP/NCAR reanalysis (Kalnay et al., 1996). Spearman correlation coefficients are computed between the time series of the regional mean frequencies and several concurrent teleconnection patterns. If statistically significant relationships are indicated composite anomaly maps are constructed using seasons (November through April) with snowfall frequencies greater than positive one standard deviation (+1 s.d.) and smaller than negative one standard deviation (−1 s.d.) from the mean. The coefficients are displayed in Table 1 for all regions with the exception of Regions 5 and 7, which had no significant correlations. In Region 1, the southeast, there are significant (p ≤ 0.05) correlations with the NAO (January and February), AO (January and February), and Northern Hemisphere annual temperature. The correlation coefficients are all negative, indicating that a positive phase of the NAO is correlated with less frequent snowfall events in this region. © 2015 Royal Meteorological Society Figure 6 shows the composite SLP anomalies for seasons with frequencies greater or less than one standard deviation from the mean. During the low frequency years the pressure pattern resembles a relatively weak positive/warm AO and NAO with low pressures over the Arctic and higher pressure in the mid-latitudes (van Loon and Rogers, 1978; Thompson and Wallace, 2000). The high frequency years are associated with a strong negative/cold AO and NAO, with high pressure over the Arctic and a strong pressure dipole in the North Atlantic which results in more frequent and stronger storms moving along the US east coast (Hartley and Keables, 1998). There is also a negative correlation with Northern Hemisphere annual temperature, indicating that increased temperatures correspond to fewer snowfall events in the southeast. This can be seen in the quantile regression plot for this region, where the 90th percentile regression line shows decreases in large snow frequency years over time. Region 2, the south central Plains and the southwest, has a clear connection with Pacific Ocean teleconnection patterns. It is positively correlated with Nino 3.4 region sea surface temperatures in December, January, and February (p < 0.05). Figure 7 displays the composite of SST anomalies for seasons with frequencies greater than +1 s.d. minus those years with frequencies less than −1 s.d. During high frequency seasons, central Pacific SSTs are warmer than during low frequency seasons. A warm equatorial Pacific indicates warm ENSO conditions during high snow frequency seasons when the jet stream is shifted farther south (Kunkel and Angel, 1999). In the Ohio River Valley and mid-Atlantic (Region 3) there are negative correlations with December and January AO (p < 0.05). In Figure 8 the composite SLP anomalies show a weak positive NAO like dipole pattern in the eastern Atlantic during low snow frequency years and a strong negative AO/NAO signal during high frequency years. During a positive AO year the higher-than-normal pressure located over the mid-latitudes shifts storm tracks more northward. This prohibits many Canadian arctic air masses from entering the central United States, leading to warm temperatures and less frequent snowfall (Thompson and Wallace, 2000). The average regional frequency in the Pacific Northwest (Region 4) has significant correlations with the PDO (December, January, and February), PNA (December, January, and February), and Nino 3.4 (December, January, and February). The highest correlation coefficients are in January with the PNA (−.717, p = 0.000), which shows up most clearly in the composite analysis in Figure 9. During a positive PNA there is enhanced meridional flow resulting in a 500 hPa ridge over the western United States (Wallace and Gutzler, 1981; Leathers et al., 1991), which has been linked to decreases in snow water equivalent (Jin et al., 2006; Sobolowski and Frei, 2007), snow depth (Ge and Gong, 2009; Ge et al., 2009; Ghatak et al., 2010) and in this study low snowfall frequency seasons (Figure 9(a)). During high frequency seasons, a 500 hPa trough is found over the region, bringing colder temperatures to the Pacific Northwest (Figure 9(b)). A Int. J. Climatol. 35: 4348–4358 (2015) REGIONALIZATION OF SNOWFALL FREQUENCY OVER THE CONTIGUOUS UNITED STATES 4355 Figure 7. Region 2 composite of tropical SST anomalies, snowfall frequency years greater than 1 s.d. minus snowfall frequency years less than 1 s.d. From NCEP/NCAR reanalysis (Kalnay et al., 1996). negative relationship also exists with the PDO, where during a positive PDO the enhanced Aleutian low results in decreased winter precipitation in the Pacific Northwest (Mantua and Hare, 2002). During a winter with positive Nino 3.4 SST anomalies there is a southward shift in the jet stream, resulting in a decrease in snowfall frequency in the Pacific Northwest (McCabe and Dettinger, 2002; Lute and Abatzoglou, 2014). This region also has a negative correlation with Northern Hemispheric annual temperature indicating that precipitation is likely falling more frequently in liquid rather than solid form during warm years. Region 6, including portions of Iowa, Indiana, southern Wisconsin, and Michigan, has significant correlations with February PDO, the PNA (December, January, and February) and Nino 3.4 SSTs (December, January, and February). The PNA affects this region’s snowfall frequency by shifting the location of the polar front jet, resulting in a southward shift during a positive PNA and more dry air advection (Serreze et al., 1998; Notaro et al., 2006). The El Nino Southern Oscillation has been documented as impacting stations in the same area as Region 6 through a weakened polar jet stream during the warm phase and reduced moisture advection from the Gulf of Mexico (Smith and O’Brien, 2001; Patten et al., 2003), as well as a reduction in the frequency of ‘Alberta Clipper’ type storms (Kunkel and Angel, 1999). © 2015 Royal Meteorological Society Composite analysis for this region showed no strong anomaly pattern. 4. Summary and conclusions This study develops a regionalization over the conterminous United States based upon snowfall frequency to better understand the distribution of this important societal and hydrological variable. This snowfall regionalization updates previous snowfall and snow cover regionalization work (Leathers et al., 1993; Hughes and Robinson, 1996) and adds a unique perspective as it is based upon snowfall frequency. Several snowfall frequency clusters are in similar regions and are corroborated by the previous snowfall amount and snow cover studies. Quantile Regression is used to determine how the distribution of snowfall frequency is changing over the period of record in each region. The influence of hemispheric and global-scale forcing mechanisms is also investigated through correlation analysis and composite analyses. A combined PCA/Cluster analysis technique results in seven unique regions in which snowfall frequency varies in a similar manner: • Region 1 – the southeastern United States, which resembles the spatial footprint of of ‘Nor’Easter’ and ‘Gulf Low’ cyclones. The time series for this region Int. J. Climatol. 35: 4348–4358 (2015) 4356 D. KLUVER AND D. LEATHERS (a) (a) (b) (b) Figure 8. Composite anomaly of Sea Level Pressure for Region 3 for November through April snowfall frequencies (a) smaller than 1 s.d. and (b) larger than 1 s.d. From NCEP/NCAR reanalysis (Kalnay et al., 1996). Figure 9. Composite anomaly of 500 hPa Geopotential Height for Region 4 for November through April snowfall frequencies (a) smaller than 1 s.d. and (b) larger than 1 s.d. From NCEP/NCAR reanalysis (Kalnay et al., 1996). is correlated with the NAO, AO, and Northern Hemispheric temperature. Quantile regression shows that large snowfall frequency years (90th percentile) have a significantly decreasing slope with time. • Region 2 – the south central Plains and southwestern United States. This region corresponds to the area of maximum cyclogenesis in the United States (Whittaker and Horn, 1981) and is most highly correlated with the Nino 3.4 sea-surface temperatures. • Region 3 – the Ohio River Valley and mid-Atlantic states, is most strongly associated with the AO, which controls the availability of arctic air masses in those areas. • Region 4 – the Pacific Northwest. This region is associated with the North Pacific storm track, most active in the winter (Paciorek et al., 2002). It is most highly correlated with the PDO, PNA, Nino 3.4, and Northern Hemispheric annual temperature and has statistically significant decreasing trends in snowfall frequency. The quantile regression shows that years with median and higher snow frequencies are occurring less often later in the period of record. • Regions 5, 6, and 7 –these regions split up the Upper-Midwest and extend to stations in New England. This area is influenced by ‘Alberta Clippers’ and is sensitive to the exact location of the jet stream as these cyclones move through the area. The northern region (North and South Dakota, Minnesota and Wisconsin) has experienced statistically significant increasing trends in all quantiles of snowfall frequency during the period of record. This indicates an upward shift in the entire snow frequency distribution for this area. The southern region has significant correlations with the PDO, PNA, and Nino 3.4. © 2015 Royal Meteorological Society Int. J. Climatol. 35: 4348–4358 (2015) REGIONALIZATION OF SNOWFALL FREQUENCY OVER THE CONTIGUOUS UNITED STATES The snowfall frequency regions defined in this study clearly reflect the prominent storm tracks across the United States. Correlations between regional frequency time series and teleconnection patterns as well as composite analyses indicate that several concurrent large-scale forcing mechanisms may be used for regional forecasts at short time scales. This is helpful for any user attempting to better prepare for the impact of snowfall events on human activities, such as snow removal (Cohen, 1982). Properly utilizing and removing snowfall is an important aspect of winter maintenance and water resource management in many regions of the United States. As the change in distributions in this study show, relying on average snowfall frequency or old or anecdotal snow frequency distributions for planning purposes could leave areas unprepared, or with wasted resources. These results also identify regions where future research is needed to better document the cause for the changing distributions. For example, the Upper-Midwest is experiencing increases in all quantiles of snowfall frequency, which could be due to a change in moisture transport (such as a shift of the mean location of ‘Alberta Clipper’ storm tracks) or temperature (such as increased temperature during snowfall events leading to a higher atmospheric water vapour capacity). Future work to answer these questions will help those interested in snowfall to anticipate future changes in the frequency distributions. References Abatzoglou JT. 2011. Influence of the PNA on declining mountain snowpack in the Western United States. Int. J. Climatol. 31: 1135–1142, doi: 10.1002/joc.2137. Barbosa SM. 2008. Quantile trends in Baltic Sea Level. Geophys. Res. Lett. 35: L22704. Baur D, Saisana M, Schulze N. 2004. Modelling the effects of meteorological variables on ozone concentration – a quantile regression approach. Atmos. Environ. 38: 4689–4699. Bradbury JA, Keim BD, Wake CP. 2002. U.S. east coast trough indices at 500 hPa and New England winter climate variability. J. Clim. 15: 3509–3517. Bradbury JA, Keim BD, Wake CP. 2003. The influence of regional storm tracking and teleconnections on winter precipitation in the Northeastern United States. Ann. Assoc. Am. Geogr. 93: 544–556. Bremner RM. 1977. Report of City of Toronto Winter Services. Department of Public Works: Toronto, Canada. Cayan DR. 1996. Interannual climate variability and snowpack in the western United States. J. Clim. 9: 928–948. Changnon SA Jr. 1979. How a severe winter impacts on individuals. Bull. Am. Meteorol. Soc. 60: 110–114. Changnon D, Merinsky C, Lawson M. 2008. Climatology of surface cyclone tracks associated with large central and eastern U.S. snowstorms, 1950–2000. Mon. Weather Rev. 136: 3193–3202. Climate Prediction Center. 2010. Northern Hemisphere Teleconnection Patterns. http://www.cpc.ncep.noaa.gov/data/teledoc/telecontents. shtml (accessed 1 February 2010). Cohen SJ. 1982. User oriented climatic information for planning a snow removal budget. J. Appl. Meteorol. 20: 1420–1427. Coleman JSM, Rogers JC. 2007. A synoptic climatology of the central United States and associations with Pacific teleconnection pattern frequency. J. Clim. 20: 3485–3497. Derksen C, LeDrew E, Walker A, Goodison B. 2000. Winter season variability in North American Prairie SWE distribution and atmospheric circulation. Hydrol. Processes 14: 3273–3290. Dingman SL. 2008. Physical Hydrology, 2nd edn. Waveland Press Inc: Long Grove, IL, 656 pp. © 2015 Royal Meteorological Society 4357 Elsner JB, Kossin JP, Jagger TH. 2008. The increasing intensity of the strongest tropical cyclones. Nature 455: 92–95, doi: 10.1038/nature07234. Frei A, Robinson DA. 1999. Northern hemisphere snow extent: regional variability 1972–1994. Int. J. Climatol. 19: 1535–1560. Ge Y, Gong G. 2009. North American snow depth and climate teleconnection patterns. J. Clim. 22: 217–233. Ge Y, Gong G, Frei A. 2009. Physical mechanisms linking the winter Pacific-North American teleconnection pattern to spring North American snow depth. J. Clim. 22: 5135–5148. Ghatak D, Gong G, Frei A. 2010. North American temperature, snowfall, and snow-depth response to winter climate modes. J. Clim. 23: 2320–2332. Gray DM, Male DH. 1981. Handbook of Snow: Principles, Processes, Management and Use. The Blackburn Press: Caldwell, NJ, 776 pp. Groisman PY, Easterling DR. 1994. Variability and trends of total precipitation and snowfall over the United States and Canada. J. Clim. 7: 184–205. Hanbali RM, Kuemmel DA. 1993. Traffic volume reductions due to winter storm conditions. Transportation Research Record No. 1287, TRB, National Research Council, Washington, DC, 159–164. Harrington JA, Cerveny RS, Dewey KF. 1987. A climatology of mean monthly snowfall for the conterminous United States: temporal and spatial patterns. J. Clim. Appl. Meteorol. 26: 897–912. Hartley S. 1996. Atlantic sea surface temperatures and New England snowfall. Hydrol. Processes 10: 1553–1563. Hartley S, Keables MJ. 1998. Synoptic associations of winter climate and snowfall variability in New England, USA, 1950–1992. Int. J. Climatol. 18: 281–298. Hartmann DL. 1994. Global Physical Climatology. Academic Press: San Diego, CA, 409 pp. Hughes MG, Robinson DA. 1996. Historical snow cover variability in the Great Plains region of the USA: 1910 through to 1993. Int. J. Climatol. 16: 1005–1018. Jin J, Miller NL, Sorooshian S, Gao X. 2006. Relationship between atmospheric circulation and snowpack in the western USA. Hydrol. Processes 20: 753–767. JISAO (Joint Institute for the Study of the Atmosphere and Ocean). 2010. The Pacific Decadal Oscillation (PDO). http://jisao.washington.edu/pdo/ (accessed 1 February 2010). Jolliffe IT. 1989. Rotation of ill-defined principal components. Appl. Stat. 38: 139–147. Jones PD, Parker DE, Osborn TJ, Briffa KR. 2013. Global and hemispheric temperature anomalies – land and marine instrumental records. In Trends: A Compendium of Data on Global Change. Carbon Dioxide Information Analysis Center, Oak Ridge National Laboratory, U.S. Department of Energy: Oak Ridge, TN, doi: 10.3334/CDIAC/cli.002. Kalnay E, Kanamitsu M, Kistler R, Collins W, Deaven D, Gandin L, Iredell M, Saha S, White G, Woollen J, Zhu Y, Leetmaa A, Reynolds R. 1996. The NCEP/NCAR 40-year reanalysis project. Bull. Am. Meteorol. Soc. 77: 437–470. Katko K. 1993. Goals and methods of winter maintenance in Finland. Transportation Research Record: Journal of the Transportation Research Board, No. 1387, TRB, National Research Council, Washington, DC, 8–11. Kluver DB, Leathers DJ. 2014. Winter snowfall prediction in the United States using multiple discriminant analysis. Int. J. Climatol., doi: 10. 1002/joc.4103. Koenker R, Bassett G. 1978. Regression quantiles. Econometrica 46: 33–50. Koenker R, Hallock KF. 2001. Quantile regression. J. Econ. Perspect. 15: 143–156. Kunkel KE, Angel JR. 1999. Relationship of ENSO to snowfall and related cyclone activity in the contiguous United States. J. Geophys. Res. 104: 19425–19434. Kunkel KE, Palecki M, Ensor L, Hubbard KG, Robinson D, Redmond K, Easterling D. 2009. Trends in twentieth-century U.S. snowfall using a quality-controlled dataset. J. Atmos. Oceanic Technol. 26: 33–44. Leathers DJ, Yarnal B, Palecki MA. 1991. The Pacific/North American teleconnection pattern and United States climate. Part I: Regional temperature and precipitation associations. J. Clim. 4: 517–528. Leathers DJ, Mote TL, Kuivinen KC, McFeeters S, Kluck DR. 1993. Temporal characteristics of USA snowfall 1945–1946 through to 1984–1985. Int. J. Climatol. 13: 65–76. Lee K, Baek HJ, Cho CH. 2013. Analysis of changes in extreme temperatures using quantile regression. Asia-Pac. J. Atmos. Sci. 49: 313–323. Int. J. Climatol. 35: 4348–4358 (2015) 4358 D. KLUVER AND D. LEATHERS van Loon H, Rogers JC. 1978. The seesaw in winter temperatures between Greenland and Northern Europe. Part I: general description. Mon. Weather Rev. 106: 296–310. Lute AC, Abatzoglou JT. 2014. Role of extreme snowfall events in interannual variability of snowfall accumulation in the western United States. Water Resour. Res. 50: 2874–2888, doi: 10.1002/ 2013WR014465. Mantua NJ, Hare SR. 2002. The Pacific decadal oscillation. J. Oceanogr. 58: 35–44. McCabe GJ, Dettinger MD. 2002. Primary modes and predictability of year-to-year snowpack variations in the western United States from teleconnections with Pacific Ocean climate. J. Hydrometeorol. 3: 13–25. Morin J, Block P, Rajagopalan B, Clark M. 2008. Identification of large scale climate patterns affecting snow variability in the eastern United States. Int. J. Climatol. 28: 315–328. Navarra A, Simoncini V. 2010. A Guide to Empirical Orthogonal Functions for Climate Data Analysis. Springer: New York, NY, 200 pp. Norrman J, Eriksson M, Lindqvist S. 2000. Relationships between road slipperiness, traffic accident risk, and winter road maintenance activity. Clim. Res. 15: 185–193. Notaro M, Wang WC, Gong W. 2006. Model and observational analysis of the Northeast U.S. regional climate and its relationship to the PNA and NAO patterns during early winter. Mon. Weather Rev. 134: 3479–3505. Paciorek C, Risbey J, Ventura V, Rosen RD. 2002. Multiple indices of Northern Hemisphere cyclone activity: winters 1949–99. J. Clim. 15: 1573–1590. Patten JM, Smith SR, O’Brien JJ. 2003. Impacts of ENSO on snowfall frequencies in the United States. Weather Forecast. 18: 965–980. Rasmusson EM, Wallace JM. 1983. Meteorological aspects of the El Nino/Southern oscillation. Science 222: 1195–1202. Ropelewski CF, Halpert MS. 1986. North American precipitation and temperature patterns associated with the El Nino/Southern Oscillation (ENSO). Mon. Weather Rev. 114: 2352–2362. © 2015 Royal Meteorological Society Schmidlin TW. 1993. Impacts of severe winter weather during December 1989 in the Lake Erie snowbelt. J. Clim. 6: 759–767. Serreze MC, Clark MP, McGinnis DL. 1998. Characteristics of snowfall over the eastern half of the United States and relationships with principal models of low-frequency atmospheric variability. J. Clim. 11: 234–250. Smith SR, O’Brien JJ. 2001. Regional snowfall distributions associated with ENSO: implications for seasonal forecasting. Bull. Am. Meteorol. Soc. 82: 1179–1191. Sobolowski S, Frei A. 2007. Lagged relationships between North American snow mass and atmospheric teleconnection indices. Int. J. Climatol. 27: 221–231. von Storch H, Zwiers FW. 1999. Statistical Analysis in Climate Research. Cambridge University Press: Cambridge, UK, 496 pp. Strasser U. 2008. Snow loads in a changing climate: new risks? Nat. Hazards Earth Syst. Sci. 8: 1–8. Tareghian R, Rasmussen P. 2013. Analysis of Arctic and Antarctic sea ice extent using quantile regression. Int. J. Climatol. 33: 1079–1086, doi: 10.1002/joc.3491. Thomas BC, Martin JE. 2007. A synoptic climatology and composite analysis of the Alberta Clipper. Weather Forecast. 22: 315–333. Thompson DWJ, Wallace JM. 2000. Annular modes in the extratropical circulation. Part I: month-to-month variability. J. Clim. 13: 1000–1016. Timofeev AA, Sterin AM. 2010. Using the quantile regression method to analyze changes in climate characteristics. Russian Meteorol. Hyrdol. 35: 310–319, doi: 10.3103/S106837391005002X. Wallace JM, Gutzler DS. 1981. Teleconnections in the geopotential height field during the Northern Hemisphere winter. Mon. Weather Rev. 109: 784–812. Whittaker LM, Horn LH. 1981. Geographical and seasonal distribution of North American cyclogenesis, 1958–1977. Mon. Weather Rev. 109: 2312–2322. Wilks DS. 2006. Statistical Methods in the Atmospheric Sciences, 2nd edn. Academic Press: Burlington, MA, 648 pp. Int. J. Climatol. 35: 4348–4358 (2015)