Abstract

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.

1. Introduction

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.

Fig. 1.

Map of stations included in the Lake Michigan regional analysis, including symbols to indicate subregions.

Fig. 1.

Map of stations included in the Lake Michigan regional analysis, including symbols to indicate subregions.

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:

 
formula

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:

 
formula

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).

d. GIS

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).

Fig. 2.

(left) Time series and (right) decadal box-and-whisker plots of (a),(b) November–March and (c),(d) December–February Lake Michigan composite temperature from winter season 1950/51 through 2012/13. The decadal plots omit the three winters subsequent to 2009/10.

Fig. 2.

(left) Time series and (right) decadal box-and-whisker plots of (a),(b) November–March and (c),(d) December–February Lake Michigan composite temperature from winter season 1950/51 through 2012/13. The decadal plots omit the three winters subsequent to 2009/10.

Table 1.

Linear trends for selected variables, starting in the 1950/51 season, along with sample means for the first three decades and later period, respectively. For trend assessments, an asterisk indicates that we may reject the null hypothesis that the population trend is zero, based on a 95% confidence interval. For Welch’s t test and (for snowfall variables) bootstrapped early and late period comparisons, an asterisk indicates that we can reject the null hypothesis of equal early and late period population mean, using a 95% confidence interval. The p values are provided for all insignificant Welch’s t-test results, while the bootstrap values refer to the percentage of 10 000 realizations for which the early period sample mean exceeds the later period mean.

Linear trends for selected variables, starting in the 1950/51 season, along with sample means for the first three decades and later period, respectively. For trend assessments, an asterisk indicates that we may reject the null hypothesis that the population trend is zero, based on a 95% confidence interval. For Welch’s t test and (for snowfall variables) bootstrapped early and late period comparisons, an asterisk indicates that we can reject the null hypothesis of equal early and late period population mean, using a 95% confidence interval. The p values are provided for all insignificant Welch’s t-test results, while the bootstrap values refer to the percentage of 10 000 realizations for which the early period sample mean exceeds the later period mean.
Linear trends for selected variables, starting in the 1950/51 season, along with sample means for the first three decades and later period, respectively. For trend assessments, an asterisk indicates that we may reject the null hypothesis that the population trend is zero, based on a 95% confidence interval. For Welch’s t test and (for snowfall variables) bootstrapped early and late period comparisons, an asterisk indicates that we can reject the null hypothesis of equal early and late period population mean, using a 95% confidence interval. The p values are provided for all insignificant Welch’s t-test results, while the bootstrap values refer to the percentage of 10 000 realizations for which the early period sample mean exceeds the later period mean.

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).

Fig. 3.

(left) Time series and (right) decadal box-and-whisker plots of (a),(b) seasonal Lake Michigan composite snowfall and (c),(d) the fraction of December–February to seasonal snowfall from winter season 1950/51 through 2012/13. The decadal plots omit the three winters subsequent to 2009/10.

Fig. 3.

(left) Time series and (right) decadal box-and-whisker plots of (a),(b) seasonal Lake Michigan composite snowfall and (c),(d) the fraction of December–February to seasonal snowfall from winter season 1950/51 through 2012/13. The decadal plots omit the three winters subsequent to 2009/10.

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.

Fig. 4.

As in Fig. 3, but for the number of days with snowfall accumulation (a),(b) exceeding a trace and (c),(d) at least 15 cm.

Fig. 4.

As in Fig. 3, but for the number of days with snowfall accumulation (a),(b) exceeding a trace and (c),(d) at least 15 cm.

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).

Fig. 5.

As in Fig. 3, but for (a),(b) the number of November–March days with precipitation exceeding a trace and (c),(d) the fraction of the number of seasonal snow days exceeding a trace to the number of November–March days exceeding a trace of precipitation.

Fig. 5.

As in Fig. 3, but for (a),(b) the number of November–March days with precipitation exceeding a trace and (c),(d) the fraction of the number of seasonal snow days exceeding a trace to the number of November–March days exceeding a trace of precipitation.

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.

Fig. 6.

Scatterplots of November–March composite temperature with (a) composite seasonal snowfall, (b) number of snow days, (c) number of 15-cm snow days, and (d) fraction of seasonal snow days to November–March precipitation days. Correlation coefficients r are provided, with an asterisk indicating a correlation that is significant at the 95% threshold.

Fig. 6.

Scatterplots of November–March composite temperature with (a) composite seasonal snowfall, (b) number of snow days, (c) number of 15-cm snow days, and (d) fraction of seasonal snow days to November–March precipitation days. Correlation coefficients r are provided, with an asterisk indicating a correlation that is significant at the 95% threshold.

Fig. 7.

(a) Time series and (b) decadal plots of the observed minus predicted seasonal snowfall using a regression of composite seasonal snowfall onto composite November–March temperature.

Fig. 7.

(a) Time series and (b) decadal plots of the observed minus predicted seasonal snowfall using a regression of composite seasonal snowfall onto composite November–March temperature.

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).

Fig. 8.

As in Fig. 3, but for the number of days with snow cover of (a),(b) at least 2.5 and (c),(d) at least 15 cm.

Fig. 8.

As in Fig. 3, but for the number of days with snow cover of (a),(b) at least 2.5 and (c),(d) at least 15 cm.

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).

Fig. 9.

GIS maps for the 1950–2013 (a) mean and (b) standard deviation of seasonal snowfall, (c) correlation of seasonal snowfall and November–March composite temperature, and (d) linear trend of seasonal snow days per decade. The dots indicate the locations of the observing stations.

Fig. 9.

GIS maps for the 1950–2013 (a) mean and (b) standard deviation of seasonal snowfall, (c) correlation of seasonal snowfall and November–March composite temperature, and (d) linear trend of seasonal snow days per decade. The dots indicate the locations of the observing stations.

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).

Table 2.

The PC loadings for the first two PCs of wintertime variability in the Lake Michigan region, using regionwide composite data. The first two PCs explain 66.8% and 14.7% of the variance, respectively.

The PC loadings for the first two PCs of wintertime variability in the Lake Michigan region, using regionwide composite data. The first two PCs explain 66.8% and 14.7% of the variance, respectively.
The PC loadings for the first two PCs of wintertime variability in the Lake Michigan region, using regionwide composite data. The first two PCs explain 66.8% and 14.7% of the variance, respectively.
Fig. 10.

Time series of the first two principal components as based on regionwide wintertime variables. The PC time series are each scaled to the variance, such that the units are standard deviations from the mean.

Fig. 10.

Time series of the first two principal components as based on regionwide wintertime variables. The PC time series are each scaled to the variance, such that the units are standard deviations from the mean.

Fig. 11.

November–March composite differences in Northern Hemisphere Reanalysis 1 (a),(c) 500-hPa geopotential height and (b),(d) 300-hPa zonal wind between the five highest and lowest (top) PC1 and (bottom) PC2 values, as based on PCA of the Lake Michigan composite variables. NCEI reanalysis data provided by the NOAA Earth System Research Laboratory/Physical Sciences Division from their website (http://www.esrl.noaa.gov/psd/).

Fig. 11.

November–March composite differences in Northern Hemisphere Reanalysis 1 (a),(c) 500-hPa geopotential height and (b),(d) 300-hPa zonal wind between the five highest and lowest (top) PC1 and (bottom) PC2 values, as based on PCA of the Lake Michigan composite variables. NCEI reanalysis data provided by the NOAA Earth System Research Laboratory/Physical Sciences Division from their website (http://www.esrl.noaa.gov/psd/).

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).

Fig. 12.

(a) Scatterplot of November–March composite temperature with December–February NAO. (b) As in (a), but for the temperature and PNA. (c) Box-and-whisker plot of composite November–March temperature with the December–February ENSO category. (d) As in (a), but for the temperature and PDO. Correlation coefficients r are provided, with p values in parentheses for insignificant correlations; MEI values were used for the ENSO-related correlation. Welch’s t tests were performed for ENSO group pairings, allowing us to reject the null hypothesis of equal population mean temperature for the El Niño and neutral comparison.

Fig. 12.

(a) Scatterplot of November–March composite temperature with December–February NAO. (b) As in (a), but for the temperature and PNA. (c) Box-and-whisker plot of composite November–March temperature with the December–February ENSO category. (d) As in (a), but for the temperature and PDO. Correlation coefficients r are provided, with p values in parentheses for insignificant correlations; MEI values were used for the ENSO-related correlation. Welch’s t tests were performed for ENSO group pairings, allowing us to reject the null hypothesis of equal population mean temperature for the El Niño and neutral comparison.

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).

Fig. 13.

As in Fig. 12, but for teleconnection indices and Lake Michigan composite seasonal snowfall. Welch’s t tests were also performed for ENSO group pairings, allowing us to reject the null hypothesis of equal population mean temperature for the El Niño–neutral and El Niño–La Niña comparisons.

Fig. 13.

As in Fig. 12, but for teleconnection indices and Lake Michigan composite seasonal snowfall. Welch’s t tests were also performed for ENSO group pairings, allowing us to reject the null hypothesis of equal population mean temperature for the El Niño–neutral and El Niño–La Niña comparisons.

4. Conclusions

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.

Acknowledgments

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/).

REFERENCES

REFERENCES
Assel
,
R. A.
,
K.
Cronk
, and
D.
Norton
,
2003
:
Recent trends in Laurentian Great Lakes ice cover
.
Climatic Change
,
57
,
185
204
, doi:.
Bard
,
L.
, and
D. A. R.
Kristovich
,
2012
:
Trend reversal in Lake Michigan contribution to snowfall
.
J. Appl. Meteor. Climatol.
,
51
,
2038
2046
, doi:.
Braham
,
R. R.
, Jr.
, and
M. J.
Dungey
,
1984
:
Quantitative estimates of the effect of Lake Michigan on snowfall
.
J. Climate Appl. Meteor.
,
23
,
940
949
, doi:.
Braham
,
R. R.
, Jr.
, and
M. J.
Dungey
,
1995
:
Lake-effect snowfall over Lake Michigan
.
J. Appl. Meteor.
,
34
,
1009
1019
, doi:.
Brown
,
R. D.
, and
P. W.
Mote
,
2009
:
The response of Northern Hemisphere snow cover to a changing climate
.
J. Climate
,
22
,
2124
2145
, doi:.
Brown
,
R. D.
, and
D. A.
Robinson
,
2011
:
Northern Hemisphere spring snow cover variability and change over 1922–2010 including an assessment of uncertainty
.
Cryosphere
,
5
,
219
229
, doi:.
Burnett
,
A. W.
,
M. E.
Kirby
,
H. T.
Mullins
, and
W. P.
Patterson
,
2003
:
Increasing Great Lake–effect snowfall during the twentieth century: A regional response to global warming?
J. Climate
,
16
,
3535
3542
, doi:.
Chatfield
,
C.
,
2004
: The Analysis of Time Series: An Introduction. 6th ed. Chapman and Hall, 333 pp.
Cordeira
,
J. M.
, and
N. F.
Laird
,
2008
:
The influence of ice cover on two lake-effect events over Lake Erie
.
Mon. Wea. Rev.
,
136
,
2747
2763
, doi:.
Crawley
,
M. J.
,
2005
: Statistics: An Introduction Using R. John Wiley and Sons, 327 pp.
Cressie
,
N. A. C.
,
1993
: Statistics for Spatial Data. Rev. ed. John Wiley and Sons, 900 pp.
Dommenget
,
D.
, and
M.
Latif
,
2002
:
A cautionary note on the interpretation of EOFs
.
J. Climate
,
15
,
216
225
, doi:.
Ellis
,
A. W.
, and
J. J.
Johnson
,
2004
:
Hydroclimatic analysis of snowfall trends associated with the North American Great Lakes
.
J. Hydrometeor.
,
5
,
471
486
, doi:.
Ge
,
Y.
, and
G.
Gong
,
2009
:
North American snow depth and climate teleconnection patterns
.
J. Climate
,
22
,
217
233
, doi:.
Gerber
,
E. P.
, and
G. K.
Vallis
,
2005
:
A stochastic model for the spatial structure of annular patterns of variability and the North Atlantic Oscillation
.
J. Climate
,
18
,
2102
2118
, doi:.
Gerbush
,
M. R.
,
D. A. R.
Kristovich
, and
N. F.
Laird
,
2008
:
Mesoscale boundary layer and heat flux variations over pack ice–covered Lake Erie
.
J. Appl. Meteor. Climatol.
,
47
,
668
682
, doi:.
Grise
,
K. M.
,
S. W.
Son
, and
J. R.
Gyakum
,
2013
:
Intraseasonal and interannual variability in North American storm tracks and its relationship to equatorial Pacific variability
.
Mon. Wea. Rev.
,
141
,
3610
3625
, doi:.
Grover
,
E. K.
, and
P.
Sousounis
,
2002
:
The influence of large-scale flow on fall precipitation systems in the Great Lakes basin
.
J. Climate
,
15
,
1943
1956
, doi:.
Higgins
,
R. W.
,
A.
Leetmaa
, and
V. E.
Kousky
,
2002
:
Relationships between climate variability and winter temperature extremes in the United States
.
J. Climate
,
15
,
1555
1572
, doi:.
Horel
,
J. D.
,
1984
:
Complex principal component analysis: Theory and examples
.
J. Climate Appl. Meteor.
,
23
,
1660
1673
, doi:.
Kristovich
,
D. A. R.
, and Coauthors
,
2000
:
The Lake-Induced Convection Experiment and the Snowband Dynamics Project
.
Bull. Amer. Meteor. Soc.
,
81
,
519
542
, doi:.
Kunkel
,
K. E.
,
N. E.
Westcott
, and
D. A. R.
Kristovich
,
2002
:
Effects of climate change on heavy lake-effect snowstorms near Lake Erie
.
J. Great Lakes Res.
,
28
,
521
536
, doi:.
Kunkel
,
K. E.
,
L.
Ensor
,
M.
Palecki
,
D.
Easterling
,
D.
Robinson
,
K. G.
Hubbard
, and
K.
Redmund
,
2009a
:
A new look at lake-effect snowfall trends in the Laurentian Great Lakes using a temporally homogeneous data set
.
J. Great Lakes Res.
,
35
,
23
29
, doi:.
Kunkel
,
K. E.
,
L.
Ensor
,
M.
Palecki
,
D.
Easterling
,
D.
Robinson
,
K. G.
Hubbard
, and
K.
Redmund
,
2009b
:
Trends in twentieth-century U.S. snowfall using a quality-controlled data set
.
J. Atmos. Oceanic Technol.
,
26
,
33
44
, doi:.
Laird
,
N. F.
,
D. A. R.
Kristovich
, and
J. E.
Walsh
,
2003
:
Idealized model simulations examining the meoscale structure of winter lake-effect circulations
.
Mon. Wea. Rev.
,
131
,
206
221
, doi:.
Lavoie
,
R. L.
,
1972
:
A mesoscale numerical model of lake-effect storms
.
J. Atmos. Sci.
,
29
,
1025
1040
, doi:.
Lorenz
,
E. N.
,
1956
: Empirical orthogonal functions and statistical weather prediction. Massachusetts Institute of Technology Department of Meteorology Statistical Forecasting Project Scientific Rep. 1, 49 pp.
Madden
,
R. A.
,
1977
:
Estimates of the autocorrelation and spectra of seasonal mean temperatures over North America
.
Mon. Wea. Rev.
,
105
,
9
18
, doi:.
Mankin
,
J. S.
, and
N. S.
Diffenbaugh
,
2014
:
Influence of temperature and precipitation variability on near-term snow trends
.
Climate Dyn.
,
45
,
1099
1116
, doi:.
Monmonier
,
M.
,
2012
: Lake Effect: Tales of Large Lakes, Arctic Winds, and Recurrent Snows. Syracuse University Press, 272 pp.
Niziol
,
T. A.
,
1987
:
Operational forecasting of lake-effect snowfall in western and central New York
.
Wea. Forecasting
,
2
,
310
321
, doi:.
Norton
,
D.
, and
S.
Bolsenga
,
1993
:
Spatiotemporal trends in lake effect and continental snowfall in the Laurentian Great Lakes, 1951–1980
.
J. Climate
,
6
,
1943
1956
, doi:.
Notaro
,
M.
,
V.
Bennington
, and
S.
Varvus
,
2015
:
Dynamically downscaled projections of lake-effect snow in the Great Lakes basin
.
J. Climate
,
28
,
1661
1684
, doi:.
Oliver
,
M. A.
, and
R.
Webster
,
1990
:
Kriging: A method of interpolation for geographical information systems
.
Int. J. Geogr. Inf. Syst.
,
4
,
313
332
, doi:.
Patten
,
J. M.
,
S. R.
Smith
, and
J. J.
O’Brien
,
2003
:
Impacts of ENSO on snowfall frequencies in the United States
.
Wea. Forecasting
,
18
,
965
980
, doi:.
Schmidlin
,
T. W.
,
1993
:
Impacts of heavy snowfall during 1989 in the Lake Erie snowbelt
.
J. Climate
,
6
,
759
767
, doi:.
Schmidlin
,
T. W.
, and
J.
Kosarik
,
1999
:
A record Ohio snowfall during 9–14 November 1996
.
Bull. Amer. Meteor. Soc.
,
80
,
1107
1116
, doi:.
Serreze
,
M. C.
,
M. P.
Clark
,
D.
McGinnis
, and
D. A.
Robinson
,
1998
:
Characteristics of snowfall over the eastern half of the United States and relationships with principal modes of low-frequency atmospheric variability
.
J. Climate
,
11
,
234
250
, doi:.
Thacker
,
W. C.
, and
R.
Lewandowicz
,
1996
:
Climatic indices, principal components, and the Gauss-Markov theorem
.
J. Climate
,
9
,
1942
1958
, doi:.
Vallis
,
G. K.
,
E. P.
Gerber
,
P. J.
Kusher
, and
B. A.
Cash
,
2004
:
A mechanism and simple dynamical model for the North Atlantic Oscillation and annular modes
.
J. Atmos. Sci.
,
61
,
264
280
, doi:.
Venerables
,
W. N.
, and
B. D.
Ripley
,
2002
: Modern Applied Statistics with S. 4th ed. Springer-Verlag, 495 pp.
Wang
,
J.
,
X.
Bai
,
H.
Hu
,
A.
Clites
,
M.
Colton
, and
B.
Lofgren
,
2012
:
Temporal and spatial variability of Great Lakes ice cover, 1973–2010
.
J. Climate
,
25
,
1318
1329
, doi:.
Wilhelmi
,
O. V.
, and
J. C.
Brunskill
,
2003
:
Geographic information systems in weather, climate, and impacts
.
Bull. Amer. Meteor. Soc.
,
84
,
1409
1414
, doi:.
Wilks
,
D. S.
,
2006
: Statistical Methods in the Atmospheric Sciences. 2nd ed. Elsevier, 627 pp.
Wright
,
D. M.
,
D. J.
Posselt
, and
A. L.
Steiner
,
2013
:
Sensitivity of lake-effect snowfall to lake ice cover and temperature in the Great Lakes region
.
Mon. Wea. Rev.
,
141
,
670
689
, doi:.
Wunsch
,
C.
,
1999
:
The interpretation of short climate records, with comments on the North Atlantic and Southern Oscillations
.
Bull. Amer. Meteor. Soc.
,
80
,
245
255
, doi:.

Footnotes

Supplemental information related to this paper is available at the Journals Online website: http://dx.doi.org/10.1175/JAMC-D-15-0285.s1.

Supplemental Material