This study has investigated the spatiotemporal structure and changes in Lake Michigan snowfall for the period 1950–2013. With data quality caveats acknowledged, a larger envelope of stations was included than in previous studies to explore the data using time series analysis, principal component analysis, and geographic information systems. Results indicate warming in recent decades, a near-dearth of serial correlation, midwinter dependence on teleconnection patterns, strong sensitivity of snowfall to temperature, peak snowfall variability and dependence on temperature within the lake-effect belt, an increasing fraction of seasonal snowfall occurring from December to February, and temporal behavior consistent with the previously reported trend reversal in snowfall.
The Laurentian Great Lakes of North America have numerous impacts on the regional climatology, particularly during the winter. Augmented snowfall within lake-effect belts impacts socioeconomic conditions in numerous ways, ranging from recreational benefits to deleterious road conditions and other societal costs (Schmidlin 1993; Schmidlin and Kosarik 1999; Monmonier 2012). These impacts have motivated a substantial body of research focused on climatological studies (e.g., Braham and Dungey 1995), field experiments (e.g., Kristovich et al. 2000), forecasting (e.g., Niziol 1987), numerical simulations (e.g., Lavoie 1972), and snowband morphology (e.g., Laird et al. 2003). In addition, there remains a need to evaluate the interannual variability of snowfall characteristics and the sensitivity of regional snowfall to increasing wintertime temperature (Kunkel et al. 2002; Burnett et al. 2003; Bard and Kristovich 2012). The former informs seasonal prediction, while the latter helps frame climate projections for the region.
In discerning the impacts of climate change on snowfall, the difficulty in quantifying the lake-effect contribution to snowfall (Braham and Dungey 1995) adds ambiguity, as does the scale separation between mesoscale lake-effect events and larger scales important for climate variability. Nevertheless, there is a clear relationship between seasonal snowfall and Great Lakes composite temperature, with greater snowfall typically occurring during colder winters (Braham and Dungey 1984; Norton and Bolsenga 1993; Burnett et al. 2003; Bard and Kristovich 2012). Consistent with this sensitivity, Kunkel et al. (2002) examined favorable conditions for heavy Lake Erie events within global climate simulations and found that warming decreased the inferred event frequency by the late twenty-first century. Recently, Notaro et al. (2015) employed dynamic downscaling from two global climate simulations and found a projected decrease in lake-effect snowfall for most of the Great Lakes during the twenty-first century, as a concurrent increase in lake-effect precipitation was more than offset by an increase in rainfall at the expense of snowfall.
However, as summarized by Bard and Kristovich (2012), Great Lakes snowfall trend assessments have largely indicated an increase in snowfall during the twentieth century (e.g., Ellis and Johnson 2004; Burnett et al. 2003). A possible mechanism is through the warming lakes and resulting reduction in ice cover (Assel et al. 2003; Burnett et al. 2003; Wang et al. 2012), which can impact lake-effect snowfall through the effect of ice concentrations on turbulent fluxes of heat and moisture (Gerbush et al. 2008). The relationship is complex, but studies have demonstrated the sensitivity of individual events to lake ice and temperature (Cordeira and Laird 2008; Wright et al. 2013).
The apparent contradiction between increasing twentieth-century snowfall and the projected future decreases may represent a nonlinear response to warming, while also being affected by the time period considered. Later trend assessments include the decades of warming since the 1970s, with visual and statistical results depending on the region, data utilized, and encompassed time period (Burnett et al. 2003; Ellis and Johnson 2004; Kunkel et al. 2009a). For example, Burnett et al. (2003) focused on the eastern Great Lakes and found an upward trend throughout the century, while Ellis and Johnson (2004) found similar long-term results but minimal change since the 1970s. Studies have typically focused on seasonal snowfall, though Grover and Sousounis (2002) utilized a synoptic classification of fall precipitation events and did not find an upward trend in lake-effect frequency.
To foster trend assessments, Kunkel et al. (2009b) examined a myriad of issues with snowfall data, ranging from observer changes and station relocation to missing data. The results generated a temporally homogeneous, quality-controlled subset of snowfall locations across the United States. Kunkel et al. (2009a) evaluated data from the 19 stations within the Great Lakes basin and assessed trends within snow belts. They confirmed the twentieth-century increase for the aggregated Lake Michigan–Huron region, yet with a significantly smaller trend than previously estimated. Of interest to the present study, visual inspection of the Kunkel et al. (2009a) time series indicates stationary Lake Michigan–Huron north snowfall since the 1970s, while the Lake Michigan–Huron south composite has a downward trajectory during this period.
Important context for the present study comes from Bard and Kristovich (2012), who employed eight Kunkel et al. (2009b) locations near Lake Michigan and data from Valparaiso, Indiana, to evaluate transects. For a series of 5-yr periods, they assessed the total lake-contribution snowfall (LCS), which encompasses both lake-effect and lake-enhanced snowfall. Results indicated a trend reversal, with LCS increasing until the 1950–80 period and subsequently decreasing. Greater LCS was coincident with colder periods and latitudinal dependence was apparent, with an earlier and less amplified reversal for the southern transect. Future explorations were suggested regarding event-scale analysis, latitudinal dependence, and the sensitivity to teleconnection patterns.
Lake-effect snowfall variability has not previously been explored extensively in the context of Northern Hemisphere teleconnection patterns, but large-scale patterns have been shown to significantly affect United States snowfall (Serreze et al. 1998). The effects of El Niño–Southern Oscillation (ENSO) result from physical modes in the climate system and impact snowfall in the United States (Patten et al. 2003). Other modes of variance such as the North Atlantic Oscillation are (arguably) better contextualized as statistical reflections of Northern Hemisphere storm-track behavior (Gerber and Vallis 2005; Vallis et al. 2004) than as pure physical modes (Dommenget and Latif 2002). In any case, a combination of these teleconnections covaries with storm-track characteristics across North America (Grise et al. 2013), wintertime climate variability in the United States (Higgins et al. 2002), and lake ice concentration within the Great Lakes (Wang et al. 2012).
The current study focuses on the spatiotemporal variability of snowfall characteristics in the Lake Michigan basin and the sensitivity to temperature, precipitation frequency, and Northern Hemisphere teleconnection patterns. This approach provides additional context for the Bard and Kristovich (2012) trend reversal and the impact of warming since the 1970s. The envelope of stations was expanded from the Bard and Kristovich (2012) Lake Michigan transects to 47 locations to foster a broad-based analysis, while the study period encompasses winters 1950/51–2012/13to reduce missing data and station relocation issues. For long-term trend assessment, we recommend a best-practice use of the quality-controlled Kunkel et al. (2009b) data. However, the expanded number of stations herein fosters this complementary and detailed spatiotemporal analysis.
2. Data and methods
a. Daily data and monthly composites
Daily snowfall, precipitation, and maximum and minimum temperature data from the National Weather Service Cooperative Observer Program (COOP) stations were obtained from the National Centers for Environmental Information (NCEI); these stations are also part of the Global Historical Climatology Network. The goal was to include numerous stations on all sides of the lake, as well as extend the domain far enough downwind of Lake Michigan to ensure that the lake-effect belt was well defined for the spatial analysis. The study period was from 1950/51 through 2012/13 and stations were limited to those with less than 10% of missing snow days within the overlapped available record and our study period. The fraction of missing daily values is more directly relevant than the fraction of missing November–March values, since stations with missing values are redacted from the respective monthly composite calculation (discussed in more detail below), yet can still be utilized in other months within the season. The spatial distribution of stations is shown in Fig. 1, while a list of station information, subregion, period of availability, and fraction of missing days is provided in Table S1 of the online supplemental information. The 34 stations with less than 10% of daily November–March temperature data missing were also used for the temperature composite.
Checks for missing daily snowfall and precipitation data were undertaken using the archived monthly Record of Climatological Observations Form reports, available through NCEI. Through this process, missing dates were identified and a small fraction of missing values were recovered. Only the most straightforward cases were adjusted; these were nearly all cases in which snowfall values were left blank on days without snowfall. The snowfall was subsequently coded as missing, but can be set to zero given the apparent monthly report quality, lack of precipitation, or mild temperature. Much rarer replacements of missing data with nonzero values were made for instances of late monthly reports for which daily data were previously demarked as missing. For the more numerous cases in which the monthly report is nearly illegible or contains numerous missing values or other ambiguity, the missing values were not further examined for adjustment. This process did not meaningfully affect the results, and was largely undertaken for future work exploring event-level analysis; for the present study, most of the affected stations remain excluded for the given month because of at least one remaining missing value.
For snowfall (using summed daily snowfall), temperature (determined using maximum and minimum daily temperature), and the frequency of daily snowfall, precipitation, and temperature exceedances, composites of monthly values were calculated for each of the subregions (listed in Fig. 1), using standardized anomalies (Wilks 2006). First, the standardized anomaly z of each station for the ith monthly observation was expressed in terms of the distance in standard deviations σ from the mean μ for that month:
These standardized monthly anomalies were then averaged among stations without missing data within each subregion, yielding μz for observation i, then converted to the subregional monthly composite Csr using the subregional monthly mean μsr and standard deviation σsr:
The monthly subregion composite values were subsequently averaged to generate the Lake Michigan composite, with months summed to calculate November–March and (for the case of snowfall variables) seasonal results. The seasonal composites for snow variables include all months with snowfall, while the composites for temperature and the number of precipitation days are limited to the November–March period. While not eliminating the effect of missing data, this composite procedure reduces the influence of missing values for instances in which a typically high or low value location is missing. A simple average within each subregion was also generated for comparison, with results very similar to the composite results.
b. Teleconnection indices
The North Atlantic Oscillation (NAO), Pacific–North American index (PNA), multivariate ENSO index (MEI), and Pacific decadal oscillation (PDO) teleconnection indices were acquired online through the National Oceanic and Atmospheric Administration (NOAA) Physical Sciences Division and Climate Prediction Center. For determining the impact of ENSO in the Lake Michigan region, December–February MEI values were used to provide cutoff values for categorical treatment, employing threshold values of 1 and −1 for El Niño and La Niña, respectively. The Arctic Oscillation was considered as an alternative to the NAO, but it does not correlate as well with the Lake Michigan temperature and snowfall composites evaluated herein.
c. Statistical analysis
Temporal changes were assessed using a linear trend through the study period, as well as Welch’s t test (Crawley 2005) and (only for snow variables) bootstrap-estimated comparisons of inferred population means from early and late periods, employing a 95% confidence level to determine significance. For the case of bootstrap estimations, the sample means from 10 000 realizations of 20-member (i.e., year) random samples from the first three decades and latter period were compared to assess whether we can reject the null hypothesis of equal population means. The bootstrap procedure used sampling with replacement and the results are very similar to Welch’s t-test results, yet provide additional confidence given the non-Gaussian distribution of snowfall variables. To augment test assessments, the serial structure of the data was examined using autocorrelation analysis; this is a useful tool, since red noise processes are capable of producing spurious behavior within short series (Wunsch 1999).
The sensitivities of composite November–March temperature, seasonal snowfall, and seasonal snowfall frequency to teleconnection indices were assessed using data visualization, as well as Pearson correlation coefficients (with reported significance based on a 95% confidence interval). Linear approaches have limitations, yet this approach benefits from simplicity (Crawley 2005) and offers a baseline assessment of the impact of teleconnection patterns.
Principal component analysis (PCA) was utilized to examine the temporal covariability of Lake Michigan composite temperature and snowfall-related variables and develop a winter severity index for the region (the leading PC). Rotation was not utilized, and the correlation matrix was used since the temperature and snowfall variables are on widely different scales. PCA methods and terminology vary in the literature (Wilks 2006), with a rich history of application within meteorology (e.g., Lorenz 1956; Horel 1984; Thacker and Lewandowicz 1996).
To visualize regional spatial patterns, ArcGIS 10.2 software was used to interpolate data. The station points were projected to universal transverse Mercator zone 16N, with linear units measured in meters. The kriging method was applied for interpolation, using the spherical semivariogram option for spatial dependence. Search radius was set to “variable” and the search criterion was set to 47 nearest input sample points. Since the goal was to derive a smooth visual pattern, using the maximum possible measured points to calculate the value for each unmeasured point was deemed appropriate in this case. Kriging is an effective interpolation technique that minimizes mean-squared error and predicts values at unmeasured points using surrounding measured points by weighting locations by distance, as well as the overall spatial arrangement of the measured points using spatial autocorrelation statistics (Cressie 1993; Oliver and Webster 1990). The modest number of observations in this study limits the quantitative certainty in the spatial predictions but produces effective maps for visual examination of regional-scale spatial patterns within the data. Mapping with GIS has become an important tool for exploration of meteorological and climatological datasets (Wilhelmi and Brunskill 2003).
3. Results and discussion
a. Temporal and spatial variability
The time series and decadal plots of November–March and December–February Lake Michigan composite temperature reveal an upward trend, superimposed on substantial noise (Fig. 2). Despite the variance and limited series length, the assessed linear trends and comparison of early and late period means allow us to reject the null hypothesis of a stationary temperature, based on both measures for the case of November–March temperature and the latter method for December–February temperature (Table 1).
The composite seasonal snowfall time series exhibits substantial noise, along with peak decadal snowfall in the 1970s (Figs. 3a,b), consistent with the trend reversal found in Bard and Kristovich (2012). We cannot reject the null hypothesis of stationary seasonal snowfall during the study period, whether measured by a linear trend or Welch’s t test and bootstrap comparisons of early and later periods (Table 1). However, there has been a significant increase in the fraction of the snowfall that occurs from December to February, as measured by linear trend or Welch’s t test and bootstrapped comparison of early and late periods (Figs. 3c,d; Table 1).
Given the potential sensitivity of snowfall reports to changes in observer, it is also valuable to assess the frequency of days with measurable snowfall. The series is noisy and contains a downward trend since the 1970s (Figs. 4a,b), largely because of the most southern and eastern subregions. There has also been a significant temporal decrease in the enhancement of snowfall frequency within eastern subregions, relative to western subregions (not shown); this is found for both east–west and east–west measures of the snow frequency increase east of the lake.
The number of days with at least 15 cm of snowfall has instead increased over time (Figs. 4c,d). The trends are not large and the assessment of significance depends on the approach employed (Table 1), but the difference between snow days exceeding a trace and 15 cm, respectively, is consistent with the combination of stationary seasonal snowfall and the decreasing number of snow days. This warrants additional investigation given the conceivable roles of more intense precipitation systems and observer issues but is beyond the scope of the present study.
We also examined the frequency of precipitation days, since the interannual variability in precipitation frequency impacts the concurrent variability in snowfall. The number of days with precipitation exceeding a trace has been quite variable during the study period and peaked in the 1970s, but has not trended significantly (Figs. 5a,b). In contrast, the fraction of seasonal snow days to November–March precipitation days has decreased significantly (as indicated by all three approaches) in recent decades (Figs. 5c,d; Table 1).
The differing temporal trends among snowfall variables may be explained in part by a varying sensitivity to temperature. The seasonal snowfall and (especially) the number of snow days are quite sensitive to temperature (Figs. 6a,b), though the dependence on temperature is greatly reduced for 15-cm snow days (Fig. 6c). Furthermore, since the number of snow days is only weakly correlated with temperature (r = −0.31), the fraction of seasonal snow days to November–March precipitation days is strongly dependent on temperature (Fig. 6d). Additional analysis regressing seasonal snowfall onto temperature reveals that the observed snowfall in recent decades has increased significantly relative to the snowfall expected from the temperature-based prediction (Fig. 7; Table 1). In essence, warming since the 1970s has diminished the snowfall frequency and increased the fraction of snowfall that occurs between December and February, while minimally affecting seasonal snowfall.
Last, we consider the number of days with snow cover of at least 2.5- and 15-cm depth, respectively (Fig. 8). There has been a decrease in regional snow cover since the 1970s (Figs. 8b,d), though the change is only significant at the 95% level for the Welch’s t test on the 2.5-cm or greater snow cover (Table 1). Although a decrease is consistent with the decreasing number of snow days, the substantial amount of missing snow-cover data for several of the stations and the difficulty in interpreting blanks on monthly reports reduces our confidence in this analysis. The result is nonetheless consistent with decreasing springtime snow cover in recent decades (Brown and Robinson 2011) and the changing fraction of precipitation that falls as snow in a warming climate (Mankin and Diffenbaugh 2014), although regional snow-cover responses are greatest in coastal regions and sensitive to altitude and seasonal mean temperature (Brown and Mote 2009).
In addition to visual inspection and trend analysis, the temporal structure was inspected. There is minimal serial correlation in the temperature and snowfall data, as indicated by autocorrelation function analysis (provided in Fig. S1 of the online supplemental material). The correlation is positive for the initial lags but insignificant and insufficient for an autoregressive model fit, while the cumulative periodogram is consistent with independent, noisy data (Chatfield 2004; Venerables and Ripley 2002). This noisy structure is similar to the wintertime analysis found for many U.S. stations in Madden (1977) and is beneficial for statistical tests. This does not mean that climatic factors with low-frequency variability, such as multidecadal oceanic modes, do not impact the climatology; their effect is simply small enough that noise dominates given the length of the series. Detailed spectral analysis of a longer series could foster an assessment of the lower-frequency behavior, as well as cyclic behavior, but this is beyond the scope of the present study.
Visual inspection of the GIS maps (Fig. 9a) reveals the impact of lake-effect snowfall, with peak seasonal snowfall within the lake-effect belt. The sensitivity to latitude is also evident, with seasonal snowfall exceeding 300 cm in the northeast subregion and less than 70 cm in the southwest. The interannual variability of snowfall has a similar pattern, with peak standard deviation in the northeastern subregion (Fig. 9b). The sensitivity of snowfall to composite temperature peaks within the eastern subregion of the lake (Fig. 9c), while the downward trend in the number of snow days is maximized adjacent to the eastern shoreline (Fig. 9d).
To further explore the temporal covariability of temperature and snowfall metrics, principal component analysis was performed on the Lake Michigan composite data. The leading principal component explains 66.8% of the variance and essentially functions as a Lake Michigan winter severity index. The combination of loadings indicates that low values of the index correspond to warmer winters with fewer cold and precipitation days as well as a reduction in seasonal snowfall and snow frequency (Table 2). The time series is quite noisy (Fig. 10a) and indicates peak winter severity in the late 1970s and a minimum severity during the winter of 2011/12. Visual inspection of the difference between composite Northern Hemisphere 500-hPa geopotential height and 300-hPa zonal wind for the most and least severe winters (from PC1) demonstrates the role of large-scale patterns, including an upper-level trough and equatorward shift in 300-hPa zonal wind maximum within the eastern United States for the most severe winters (Figs. 11a,b).
In contrast to PC1, the second principal component has trended upward and produces higher index values for warmer temperatures, greater snowfall, and more frequent precipitation, yet only explains only 14.7% of the variance (Fig. 10b). Consistent with a large-scale pattern that fosters more frequent precipitation for high PC2 winters, visualization of the difference between composite Northern Hemisphere 500-hPa geopotential height and 300-hPa zonal wind for the highest and lowest PC2 winters reveals warmer conditions in the Great Lakes during the highest (relative to lowest) PC2 cases, along with an upper-level trough in the western United States and an associated displacement in 300-hPa zonal wind maximum (Figs. 11c,d).
b. Sensitivity to teleconnection patterns
The composite temperature is sensitive to Northern Hemisphere teleconnection patterns, with the (relatively) strongest linear sensitivity to the NAO (Fig. 12a). The sensitivity to the PNA is weaker and just outside the threshold for 95% significance (Fig. 12b), while there is no obvious sensitivity to the PDO (Fig. 12d). The sensitivity to ENSO was evaluated through the correlation with December–February MEI, as well as a categorical treatment of ENSO. While the temperature is poorly correlated with the December–February MEI values, categorical treatment of ENSO reveals substantial warming during El Niño winters, whereas La Niña winters are largely indistinguishable from neutral winters (Fig. 12c).
Snowfall is very sensitive to ENSO, with substantially lower composite snowfall associated with El Niño (Fig. 13c) winters. The sensitivities to the NAO (Fig. 13a) and PNA (Fig. 13b) are weaker, although the PNA has a significant negative correlation with composite snowfall. The sensitivity to the NAO (PNA) is furthermore stronger (weaker) to the east of Lake Michigan (maps are provided in Fig. S2 of the online supplemental material); this west–east change may simply reflect the role of lake-effect snow east of the lake. It may relate to the west–east pattern characteristic of the canonical PNA and associated storm-track position as well, but this was not evident near Lake Michigan in the snow-cover and teleconnection analysis of Ge and Gong (2009). As with the temperature, the composite snowfall is insensitive to the PDO (Fig. 13d).
The time series of seasonal snowfall and temperature in the Lake Michigan region are noisy, with minimal autocorrelation. Both temperature and snowfall are dependent on large-scale teleconnection patterns during the core winter months, as encapsulated by the combination of ENSO phase, NAO, PNA, and PDO. The ENSO impact is primarily within El Niño winters, which are significantly milder and less snowy. Milder conditions are also associated with the high values of the NAO, while snowfall is more sensitive to the PNA. Based on principal component analysis, the leading mode of wintertime variability is sensitive to teleconnection patterns and indicates peak winter severity in the late 1970s.
The seasonal snowfall and frequency of snowfall are quite sensitive to the composite temperature, as is the fraction of snow days to precipitation days. However, this sensitivity is greatly reduced for higher snowfall exceedances. The recent decades of warming have not produced a significant trend in seasonal snowfall, yet the frequency of snowfall and (especially) the fraction of snow days to precipitation days has decreased since the 1970s. Furthermore, the proportion of seasonal snowfall that occurs from December to February has increased.
Spatially, the peak snowfall, frequency, and variability of seasonal snowfall each occur within the northeastern subregion, with a secondary peak in the southeast subregion. The correlation of snowfall to temperature peaks within the eastern subregion, while the downward trend in the frequency of snow days is maximized all along the eastern shore and produced from months outside of the December–February period.
While seasonal snowfall has thus far been impervious to warming, this may not be true over a longer time horizon, once the core winter months are substantially warmer. However, the roles of lake temperature and ice concentration in a warming climate make this a complex question. The sensitivity to temperature and spatiotemporal variability of fall and spring transition season snowfall in the Lake Michigan region are also of interest, and are the subject of an ongoing complementary study.
The authors thank Valparaiso University for funding this work. We also thank student and faculty colleagues within the Department of Geography and Meteorology at Valparaiso University, especially Kevin Goebbert and Ryan Connelly, for their ongoing help with this project. We also thank Kristin Smedley for participating in early stages of the project, the National Oceanic and Atmospheric Administration (NOAA) and National Centers for Environmental Information (NCEI) for the online access to daily station data and scanned monthly COOP reports, the NOAA Physical Sciences Division and Climate Prediction Center for teleconnection indices, and the open-source statistical package R (available at http://cran.r-project.org/).
Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JAMC-D-15-0285.s1.