The Laurentian Great Lakes Basin (GLB) is prone to snowfall events developed from extratropical cyclones or lake-effect processes. Monitoring extreme snowfall trends in response to climate change is essential for sustainability and adaptation studies because climate change could significantly influence variability in precipitation during the 21st century. Many studies investigating snowfall within the GLB have focused on specific case study events with apparent under examinations of regional extreme snowfall trends. The current research explores the historical extremes in snowfall by assessing the intensity, frequency, and duration of snowfall within Ontario’s GLB. Spatiotemporal snowfall and precipitation trends are computed for the 1980 to 2015 period using Daymet (Version 3) monthly gridded interpolated datasets from the Oak Ridge National Laboratory. Results show that extreme snowfall intensity, frequency, and duration have significantly decreased, at the 90% confidence level, more so for the Canadian leeward shores of Lake Superior than that of Lake Huron, for the months of December and January. To help discern the spatiotemporal trends is snowfall extremes, several trend analyses for lake-induced predictor variables were analysed for two cities, Wawa and Wiarton, along the snowbelts of Lakes Superior and Huron, respectively. These variables include monthly maximum and minimum air temperature, maximum wind gust velocity, lake surface temperature, and maximum annual ice cover concentration. Resultant significant increase in December’s maximum and minimum air temperature for the city of Wawa may be a potential reason for the decreased extreme snowfall trends.

1. Introduction

During the winter season, heavy snowfall is a prominent meteorological phenomenon in the Great Lakes Basin (GLB) and is derived from either extratropical cyclones or lake-induced snowfall processes. Extratropical cyclones are generated by quasigeostrophic forcing from positive temperature or vorticity advection [1] and are associated with a low-pressure centre tracking zonally following the jet stream. Thus, frequent wintertime extratropical storms, such as the Alberta Clipper, Colorado Low, and Nor’easter, track from west to east, affecting surface-atmosphere conditions within the GLB region. Contrary to these large-scale synoptic systems are shallow meso-beta scale lake-effect snowfall (LES). LES ranges from approximately 5 km to a few hundred km in length and is triggered by turbulent energy and moisture fluxes off lakes [1].

The shear spatial extent and geographic location of the Laurentian Great Lakes make LES a prominent meteorological phenomenon during the late autumn and winter months [1]. LES forms when boundary layer convection is initiated as a result of a cold and dry continental air mass advecting over the relatively warm lake, generating turbulent moisture and heat fluxes that destabilize the lower part of the planetary boundary layer (PBL). The increase of moisture into the lower atmosphere enhances cloud formation and precipitation along the leeward shores of the Great Lakes [112]. The lake-effect vertical atmospheric profile often features a moist-neutral or unstable convective boundary layer that extends 1 to 4 km above lake surfaces, where a capping stable layer or inversion limits the vertical extent of convection [1216]. In addition to LES, lake-enhanced snowfall can occur when synoptic scale systems, driven by quasigeostrophic forcing, move over the GLB. The synoptic system can increase the altitude of the capping inversion, producing deeper cloud convection, and as a result, increase snowfall. Thus, the term lake-induced snowfall is used to describe both lake-effect and lake-enhanced snowfall processes, as defined by [17]. Lake-induced snowfall can organize into a spectrum of LES morphologies ranging from discrete and disorganized cells to organized mesoscale bands [16] and include widespread coverage, shoreline bands, midlake bands, and mesoscale vortices [1, 11, 1820]. Convective bands are typically long and narrow, between 50 and 300 km in length and 5 to 20 km in width [1, 10], and can produce highly localized snowfall in cities that are along the snowbelts of the GLB.

The highly populated GLB is home to over 33 million people including 90% of Ontario’s population. Over 1.5 million people reside along the Canadian shores of Lake Huron and 200,000 along Lake Superior’s [21]. The Laurentian Great Lakes support 40% of Canada’s economic activities, 25% of its agricultural capacity, and 45% of its industrial capacity. These lakes also support approximately 400 million dollars in cumulative recreational and commercial fishing industries and 180 billion dollars in Canada and US trade [22]. Thus, heavy lake-induced snowfall in the GLB can impact local residents and numerous industries [2327]. The GLB’s provision on the surrounding communities and its climatologically favourable snowfall location gives credence to assessing the historical spatiotemporal trends in snowfall extremes within the context of climate change.

According to the Intergovernmental Panel on Climate Change, Fifth Assessment Report, (IPCC AR5), climate change, driven by either natural or human forcings can lead to changes in the occurrence or strength of extreme weather and climate events such as extreme precipitation [28]. Therefore, assessing and monitoring spatiotemporal trends in snowfall extremes to a historically changing climate can provide knowledge as to future behavioural spatiotemporal trends in snowfall extremes within the densely populated Canadian GLB domain and can be useful for sustainability and adaptation studies.

However, only a few lake-induced snowfall studies focused on the Canadian GLB while many studies, including [25, 27, 2932], concentrated on the United States snowbelts. Furthermore, underexaminations of historically observed LES extremes compared to LES case studies are apparent. The majority of LES investigations have analysed specific LES events and associated features as opposed to LES climatology, for example, [9, 10, 3342].

Several studies have examined trends in snowfall patterns across the GLB over the 20th century and found a snowfall trend reversal, in which there was an increase in snowfall prior to the 1980s, followed by a decrease in snowfall over the past few decades. A study by [43] investigated snowfall trends over North America from 1948 to 2001 and observed a snowfall increase over a narrow band from Colorado to the lee of Lake Erie and Lake Ontario. Furthermore, [44] found that the area-averaged snowfall across Northern Canada (55°N to 88°N) significantly increased at a rate of 8.8 cm/decade in the late 20th century. However, in Southern Canada (below 55°N), a negative trend of 0.65 cm/decade during this period was discovered, with the most decrease occurring in the 1980s [45]. Study [17] suggests that while some studies have shown a general increase in snowfall over the GLB [8, 23, 25, 29, 46, 47], recent studies have seen a decline in LES through the later half of the 20th century and early 21st century. For example, [48] showed a decrease in snowfall along the leeward shores of Lake Michigan between 1980 and 2005 and [32] showed a decrease in Central New York between 1971 and 2012. While these studies outline trends in North American snowfall, there is still a lack in the examination of climatological snowfall extremes over the GLB. Little is still known about the physical processes influencing the past changes in daily snowfall extremes on the global and hemispheric scales. On a regional scale, observational studies show large interdecadal variations in snowfall extreme measures; however, long-term trends remain unclear [4951].

The objective of this study was to assess historical (1980–2015) spatiotemporal trends in snowfall extremes for the Canadian snowbelt zones of Lake Superior and Lake Huron-Georgian Bay (Figure 1). This will be conducted by examining snowfall intensity, frequency, and duration and provide potential explanations of the results in the context of lake-induced predictor variables, including monthly extreme maximum and minimum air temperature, maximum wind gust velocity, lake surface temperature (LST), and annual maximum ice cover concentration. The study is divided into two folds: the first will assess the statistical extremes of snowfall intensity, frequency, and duration over the 36-year period, while the second will examine the trends in these extremes over the given duration. The term “snowfall extremes” will be used when referring to all extremes in intensity, frequency, and duration, unless otherwise specified. In this paper, Section 1 presents the introduction that provides a background on the development of snowfall within the GLB region; Section 2 outlines the data and methodology; Section 3 presents the results; Section 4 provides discussions and analyses on the results; and Section 5 summarizes the study and suggests future potential research avenues.

2. Data and Methods

2.1. Datasets

The datasets used in this study include gridded interpolated precipitation and temperature data from Daymet. Predictor lake-enhanced meteorological variables, such as temperature and winds, are from weather observation stations for the cities of Wawa and Wiarton Ontario and are provided by Environment and Climate Change Canada (ECCC). In addition, the National Oceanic and Atmospheric Administration (NOAA) provides lake ice concentration and lake surface temperature (LST) datasets. Table 1 lists the datasets and sources used in this study.

2.1.1. Daymet

Daymet (Version 3) provides meteorological data used in this study and constitutes a collection of software and algorithms that produce gridded estimates of daily weather observations, including daily maximum and minimum temperature and precipitation. Daymet inputs contain land mask and the North American subset of the National Aeronautics and Space Administration (NASA) Shuttle Radar Topography Mission (SRTM) near-global 30-arc second digital elevation model, Version 2.1. It also includes spatially referenced ground-based observations obtained from NOAA's National Centers for Environmental Information’s Global Historical Climatology Network (GHCN), Version 3.22 [52, 53].

An interpolation method at each prediction point employs an iterative estimation of location density that uses spatial convolution of a truncated Gaussian filter to derive 1 km × 1 km grids of each meteorological parameter. In the algorithm, the search radius of stations is reduced in data-rich regions and increased in data-poor regions. The dataset was developed from the Environmental Sciences Division at the Oak Ridge National Laboratory [54, 55]. A detailed description of this dataset can be found at http://daymet.ornl.gov. Daily minimum and maximum temperature and precipitation over North America, for the years 1980 to 2015, were obtained online at http://daac.ornl.gov/cgi-bin/dsviewer.pl?ds_id=1328 and used to derive gridded daily snowfall.

It is also acknowledged that precipitation and temperature variables from Daymet V3 are the estimated values. Extensive cross-validation statistics for Daymet V3 data were computed by station-based daily observations and predictions over 2° × 2° tiles for North America [55]. The tiles provide information on the period-of-record mean absolute error (MAE) and bias statistics for input weather observations of maximum and minimum temperature and precipitation for each year, 1980–2015. A detailed validation summary can be found at https://daac.ornl.gov/DAYMET/guides/Daymet_V3_CrossVal.html#revisions. Daymet data are ideal for this study after considering other dataset options, such as the North American Regional Reanalysis (NARR). While NARR offers 3 hourly precipitation products, it comes at the cost of a coarser spatial resolution of 32 km. The 1 km daily Daymet products provide a high spatial and temporal resolution for producing and delineating highly localized snowsquall bands over the Ontario snowbelts.

2.1.2. Environment and Climate Change Canada Information Archive

Meteorological datasets for the cities of Wiarton and Wawa are provided by ECCC. The geographic coordinates and associated snowbelts are shown in Table 2. The observational data from ECCC’s historical archive data can be found at the National Climate Data and Information Archive, http://climate.weather.gc.ca. Historical ground-based weather observation stations from ECCC provide maximum and minimum air temperature and are defined as the highest and lowest observed air temperature at the specific location. The direction of maximum wind gust for each month is also attained and defined as the direction of maximum wind gust from which the wind blows, with the value of 360° indicating winds advecting from true north. Values are only reported for gust speeds that exceed 29 km/h. The speed of maximum wind gust is also analysed for both sites. The gust is defined as the peak instantaneous or single reading from the weather station anemometer. The elapsed time corresponding to the duration of the gust is typically between three to five seconds. Apart from the atmospheric meteorological observational variable, lake variables are also analysed.

2.1.3. National Oceanic and Atmospheric Administration Ice Atlas and Coast Watch

NOAA provides lake-wide annual maximum ice cover concentration and monthly average LST for Lakes Superior and Huron. LST datasets can be found at https://coastwatch.glerl.noaa.gov. Ice chart data are a blend of observations from ships, shore stations, aircraft, and satellites to estimate ice cover data for the entire Great Lakes. The ice charts were digitized and made available at https://www.glerl.noaa.gov/data/ice/atlas/index.html.

2.2. Snowfall Derivation

Snowfall is derived from daily precipitation, P1 (mm), and daily midpoint 2 m air temperature, Tmid (°C), based on the method used by [56]. Because Daymet does not provide daily mean 2 m air temperature, the midpoint value is used and calculated by taking the sum of the maximum and minimum daily temperature and dividing by 2, as shown in the following equation:

The Tmid is used to estimate the ratio (R), which denotes the daily precipitation falling as snowfall to that of the total daily precipitation, where , equation (2) [57]. The following function was derived using a logistic curve and fitted to monthly data with a reported absolute deviation of 0.06 mm per month [58]:

Multiplying P1 and R yields the water equivalent of the new daily snowfall, Pn (mm), as shown in the following equation:

A density value is then computed as a function of Tmid (4). It is noted that when Tmid ≤−15°C, the density of snowfall is assigned to 0.05 gcm−3:

Finally, the estimated height of new snowfall, Hn (cm) is given by

Additional uncertainty analysis was conducted by [56] to quantify uncertainties in derived snowfall. Considering that Hn is a function of R, P1, , and , first, annual uncertainty values were calculated for each of these variables. To calculate the initial uncertainty of precipitation and , two 2° × 2° tiles that contained Wawa, a city within the snowbelt of Lake Superior, and Wiarton, a city within Lake Huron-Georgian Bay’s snowbelt, were selected from Daymet’s cross-validation data resources. Both tiles provided the MAE of precipitation, which were converted to the error of standard deviation and used as the precipitation uncertainty. Furthermore, annual uncertainty measures of were calculated by computing the annual standard deviation between both Daymet and an independent dataset, NARR, for the cities of Wawa and Wiarton. The uncertainty estimate of density was given as 0.005, 10% of its initial value. The calculated annual uncertainties of R, P1, , and were then inputted into the derivative of the above equations to determine the annual sensitivity of Hn to , Hn to P1, R to , to , Hn to R, Hn to (if ≥−15°C), and Hn to (if <−15°C).

Uncertainty analysis shows that the overall 36-year annual uncertainties in Hn relative to are 1.7 cm for Wiarton and 6.2 cm for Wawa. The sensitivity of Hn to P1 is 0.03 cm for Wiarton and 0.12 cm for Wawa. The sensitivity of snowfall to R averages 1.4 cm for both Wiarton and Wawa. It is expected that the tile containing Wawa would have higher snowfall uncertainties than that of Wiarton because it is farther north and contains less numbers of weather observation stations. Furthermore, when predicting snowfall accumulation, meteorologists usually provide a 5 cm predictive range in snowfall, for example, 1–5 cm or 10–15 cm of snow. Thus, an uncertainty value of 1.7 cm for Wiarton and 6.2 cm for Wawa are within a reasonable uncertainty. These uncertainties are fairly low and provide confidence in the derived snowfall method used in this study. The derived gridded daily snowfall is computed over the GLB and is bounded by the coordinates 95° W, 75° W, 40° N, and 45° N. Over this region, the extremes in snowfall (intensity, frequency, and duration) for each month (November, December, January, February, and March) over the 36-year duration (1980 to 2015) are computed.

2.3. Defining Snowfall Extremes

Snowfall extreme analyses are evaluated in two folds: first by exploring the monthly spatiotemporal extremes in snowfall intensity, frequency, and duration over the 36-year period (from here on referred to as the 36-year extreme) and secondly by determining the trends in these extremes. Firstly, the 36-year extremes are explained. Snowfall intensity is defined as the rate of snowfall over time. The 99th percentile of daily snowfall accumulation is computed over all time steps, that is, all days between 1980 and 2015 for each individual month. The 99th percentile over the 36-year period is then spatially mapped to determine the extreme values of snowfall intensity and its location for each month of the cold season (November, December, January, February, and March). Based on each month’s extreme intensity snowfall, a snowfall midpoint value is chosen that best represent the extreme snowfall threshold intensity (S_threshold) for the overall cold months.

The S_threshold value is then used to evaluate the 36-year extreme in snowfall frequency and duration. Frequency is defined as the number of snowfall days over a particular time period. In this study, the extreme snowfall frequency is determined by calculating the number of snowfall days that the given S_threshold is met or exceeded for each month over all time steps of the 36-year period. Monthly maps are generated to determine regions within the GLB of high and low extreme snowfall frequencies during the 36-year period. Snowfall duration is described as the length of time it snows at a particular intensity. The 36-year duration extreme is presented by mapping the maximum number of consecutive snowfall days, for which events equaled or exceeded the acquired S_threshold value over this time period. This procedure is repeated for each month of the cold season to determine spatial patterns in extreme snowfall duration throughout the season.

Secondly, the trends in 36-year extremes are calculated. Trends in extreme intensity are computed by applying a filter to extract daily snowfall events greater than or equal to the S_threshold value for each month and for each year. The filtered values are aggregated and then divided by the total number of snowfall days in the month, yielding the monthly average snowfall intensity. Note a snowfall day is defined as a full 24-hour day, for which snowfall exceeds 0 cm of snow. The procedure is repeated for each of the 36 years. The Mann–Kendall (MK) test is then applied to this time series to compute the trend in monthly extreme snowfall intensity. The trend is computed for each month of the cold season, and monthly spatiotemporal extreme intensity trends are produced. Trends in snowfall frequency are computed by counting the number of monthly extreme events (filtered values) for each year. The MK test is then applied over this time series to evaluate the spatiotemporal trends in monthly snowfall extremes. The trend in extreme snowfall duration is calculated by computing the MK test over each monthly time series. The time series comprises the monthly maximum number of consecutive snowfall days, for which events equaled or exceeded the S_threshold.

The MK test is applied to the extreme snowfall time series to determine whether there are monotonic upward or downward trends in intensity, frequency, and duration for each grid cell over the 36-year period. If grid cells do not have a complete 36-year time series of extreme values, then no trend analysis is computed for that grid cell, and it is assigned “not a number” and shaded gray. The MK test evaluates the slope of an estimated linear regression line nonparametrically and does not require the residual from a fitted regression line to be normally distributed, like that of a parametric linear regression [59, 60]. In this study, slope is calculated over the 36-year period for each grid cell in order to generate spatiotemporal trend maps. For the purpose of this study, the significance of the extreme snowfall trends is calculated at the 90% confidence level. This is because a higher confidence level may not be able to capture the highly episodic and disorganized spatial patterns of lake-induced trends seen in the results. The predictor variable trends are also estimated using the MK approach for the cities of Wawa and Wiarton Ontario.

3. Results

3.1. Snowfall Intensity Extreme Values

A predefined extreme intensity snowfall value is required in order to conduct extreme snowfall analyses. For this study, the S_threshold value is computed by calculating the 99th percentile of daily snowfall for each month over the 36-year period (Figures 2(a)2(e)). Within the Canadian GLB domain, the leeward shores of Lake Superior and Lake Huron-Georgian Bay show the greatest value in snowfall intensity for all months. Along these shores, the lower snowfall intensity values range between 5 cm/day and 15 cm/day for the months of November and March (Figures 2(a) and 2(e)) and the higher values range between 15 cm/day and 30 cm/day for December and January (Figures 2(b) and 2(c)).

The S_threshold value of 15 cm/day of snowfall is adopted to evaluate extreme snowfall over the 36-year study period. It is noted that preliminary work also investigated smaller threshold values (5 cm/day and 10 cm/day) to determine whether there were any month-to-month spatiotemporal variations. However, results indicated similar spatiotemporal trends to that of 15 cm/day, but with smaller magnitudes of the extreme values. Due to these findings, only the 15 cm/day S_threshold values are plotted to reduce repetitiveness in the results.

Furthermore, the 15 cm/day S_threshold is ideal because it is the midpoint value among the cold months being analysed. This S_threshold is also in agreement with Environment and Climate Change Canada’s warning criteria for Ontario. Environment Canada issues a snowfall warning for Ontario when the alerting snowfall parameter reaches or exceeds 15 cm of snowfall within 12 hours. The alerting parameter for Ontario’s snowsquall warnings is similar but also factors in a reduced visibility criterion of 400 m with or without the presence of persistent blowing snow [61]. Since this study is limited to a daily temporal resolution, 15 cm/day is used as the S_threshold value. Furthermore, 15 cm/day is a reasonable S_threshold value. For example, when separately considering historical 1981 to 2015 monthly snowfall average for the city of Barrie, which lies within Lake Huron’s snowbelt, Barrie experiences its highest snowfall in January with an average of 66 cm. This suggests an average of approximately only 2 cm/day with assumed continuous daily snowfall [62]. Therefore, a snowfall rate of 15 cm/day would be considered an extreme snowfall event of high intensity. This S_threshold is now applied to evaluate extremes in snowfall frequency and duration.

3.2. Snowfall Frequency Extreme Values

Results show that extremes in snowfall intensity are predominant along the leeward shores of Lake Superior and Lake Huron-Georgian Bay during the cold season. Next, the spatiotemporal frequency in extreme snowfall events is explored by determining the number of days that snowfall events equaled or exceeded S_threshold of 15 cm/day over the 36-year period. At the start of the LES season, in November, the leeward shores of Lake Superior begin to experience higher frequencies in extreme snowfall events compared to other regions of the GLB (Figure 3(a)). The greatest spatial extent with highest frequencies upward of 50 days is seen on the leeward shores of Lake Superior and Lake Huron for the months of December (Figure 3(b)) and January (Figure 3(c)). The spatial and temporal onsets of high-frequency snowfall behave similar to that of lake-induced snowfall spatiotemporal patterns. In February, (Figure 3(d)) the frequency is still high along Lake Superior’s Canadian snowbelt, but becomes highly localized, while both the spatial extent and higher frequencies along the snowbelt of Lake Huron-Georgian Bay have decreased. By March, both the high-frequency values and spatial extent have decreased for both snowbelts with frequencies less than 10 days (Figure 3(e)). Similar spatiotemporal results are exhibited when examining extreme snowfall duration.

3.3. Snowfall Duration Extreme Values

Extreme snowfall duration plots show the maximum number of consecutive days when snowfall intensity equaled or exceeded 15 cm/day for the 36-year period. As early as November, highly localized areas along the leeward shores of Lake Superior experienced high duration of extreme snowfall events upwards of 3 days, while Lake Huron-Georgian Bay’s snowbelt shows no maximum snowfall durations (Figure 4(a)). In December and January, most of Northeastern Ontario experienced extreme events greater than 2 days. However, the duration of extreme snowfall events greater than 5 days was spatially confined to the snowbelt regions. By February, the duration of extreme snowfall events has spatially decreased along the Lake Huron-Georgian Bay snowbelt, but continue to stay dominant along the leeward shores of Lake Superior. By March, both snowbelts show a spatially less extreme in snowfall duration (Figure 4(e)).

3.4. Trends in Snowfall Extremes

This section evaluates the 1980 to 2015 spatiotemporal trends in extreme snowfall intensity, frequency, and duration for each month of the cold season, with S_threshold of 15 cm/day. The study [62, 63] showed clear, consistent, and contiguous spatiotemporal trends in monthly snowfall totals for Lake Superior’s and Lake Huron-Georgian Bay’s snowbelts. However, spatiotemporal trends in monthly mean snowfall extremes are less coherent and less obvious to delineate.

3.4.1. November

Extreme snowfall intensity, frequency, and duration experience no significant changes over the Great Lakes Basin for the month of November. Although not significant, the snowfall intensity trends can provide indication of extreme snowfall evolution, shown in Figure 5. For instance, the northern leeward shores of Lake Superior show a decrease in extreme snowfall intensity, while its southern leeward shores present an increase, indicating the inconsistency of coherent spatial distribution along this snowbelt. Furthermore, there is a strong negative trend within the region of Northeastern Ontario, which is farther inland from the leeward shores of Lake Superior. This geographic location implies that the strong decrease in extreme snowfall trend may not be due to LES but perhaps to synoptic and extratropical scale snowstorms. The spatiotemporal trends for the month of November suggest that an increase in snowfall extremes in Northern Ontario is perhaps due to lake-induced processes, but a decrease in extratropical and synoptic scale-driven snowstorms.

In Southern Ontario, most of this region does not experience daily extreme snowfall greater than 15 cm/day in November. Recall from Figure 2, the monthly mean extreme snowfall for November does not exceed 12 cm/day for most grid cells in Southern Ontario and is indicated by gray pixelated regions. Furthermore, extreme snowfall frequency and duration show no changes for the month of November and are not included in the results. December results provide more substantive evidence of trends in snowfall extremes for both snowbelt regions.

3.4.2. December

December exhibits the most coherent significant spatiotemporal trends in snowfall extremes for all months of the cold season. In December, Lake Superior’s snowbelt displays two distinct trends. The southern leeward shore reveals a significant decrease in extreme snowfall intensity, while the northern leeward shore shows a highly localized significant increase (Figure 6(a)). Both the northern and southern leeward shores exhibit a significant decrease in extreme snowfall frequency (Figure 6(b)). However, the northern shore displays a weaker negative trend. The duration also significantly decreases for the southern shore, with no change in extreme snowfall duration for the northern leeward shore of Lake Superior’s snowbelt (Figure 6(c)). During this month, the snowbelt of Lake Huron-Georgian Bay shows a highly localized decrease in snowfall intensity but no changes in extreme snowfall frequency and duration. It is evident that not only are spatiotemporal trends in snowfall extremes different for both snowbelts but also, within the same snowbelt, trends in snowfall extremes are not spatially coherent.

3.4.3. January

The month of January shows less spatial coherency than December for both snowbelt regions. The northern leeward shore of Superior shows a significant decrease in snowfall intensity (Figure 7(a)), frequency (Figure 7(b)), and duration (Figure 7(c)). Southern Ontario shows spatial incoherent negative trends in extreme snowfall intensity, while there are no changes in frequency and duration. Compared to December, the spatiotemporal trends in extreme snowfall are highly localized and overall less spatially coherent.

3.4.4. February and March

The month of February shows incoherent spatial positive and negative trends in extreme snowfall intensity for both snowbelts (Figure 8(a)). Extreme snowfall frequency results show a highly localized decrease for the Superior snowbelt (Figure 8(b)), while there are no changes in extreme snowfall duration, results are not shown. There are also no changes in frequency and duration for Lake Huron-Georgian Bay’s snowbelt.

March exhibits a significant decrease in extreme snowfall intensity along the leeward shores of Lake Superior and farther inland, outside of the snowbelt zone (Figure 9). The geographic location and timing of this negative intensity trend suggest that the decrease in extreme snowfall intensity may not solely derive from LES but, perhaps as well from synoptic scale snowstorms. Frequency and duration in snowfall extremes exhibit no changes for this snowbelt, results are not shown. Lake Huron-Georgian Bay’s snowbelt also shows negligible changes in snowfall extremes.

3.4.5. Predictor Variables

Lake-induced predictor variables are also analysed for each month of the cold season for the cities of Wawa (Table 3) and Wiarton (Table 4). Both locations and all months show a rise in maximum and minimum air temperature, with the exception of February (Wawa) and November (Wiarton). There is significant warming in maximum monthly air temperature for December (Wawa) and January (Wiarton). There is a significant warming in minimum monthly air temperature in both December and February (Wawa) and for the month of March (Wiarton).

Both cities do not show a significant change in direction of maximum wind gust. For Wawa, the trends in maximum wind gusts are positive for all months, except for March, when there is no change. For Wiarton, all trends are positive except for the month of December, which shows a negative trend. Furthermore, the speed of maximum wind gust shows no changes for the Wawa, except for a decrease occurring in November. In Wiarton, the speed of maximum wind gusts has decreased for each month of the cold season, with a significant decrease occurring in January.

LST has also warmed for Lakes Superior and Huron for all months of the cold season, with a significant increase occurring for Lake Superior in January. Refer to Tables 3 and 4 for Sen slope values and significance. The final variable analysed is maximum annual ice cover concentration for both lakes. These lakes both show a decrease in maximum annual ice cover concentration, though not significant (Table 5).

4. Discussion

4.1. Snowfall Intensity, Frequency, and Duration Extreme Values

The 36-year snowfall intensity, frequency, and duration, as determined by the daily 15 cm/day S_threshold, are assessed over the Canadian domain of the GLB. The resultant behavioural patterns in the spatial and temporal snowfall extremes suggest that they are predominantly attributed to lake-induced snowfall rather than extratropical storms. This is because the snowfall extremes exhibit similar temporal and spatial patterns to that of LES. Spatially, these extremes are predominant along the Canadian leeward shores of Lakes Superior and Huron, where lake-induced snowfall occurs. Temporally, these resultant extremes in snowfall are dominant during the lake-effect months. According to [4], these snowbelt regions experience high annual snowfall totals, exceeding 250 cm during the autumn and winter months due to LES. These spatial patterns in extreme snowfall are also seen for all cold months in this study, with the highest extreme snowfall intensities occurring during December and January. This is in agreement with the study [64] that suggests that LES within the GLB is most intense during December and January when there is maximum open water and cold Arctic air masses that facilitate the exchange of moisture and heat between the lake and atmosphere.

Results show that snowfall extremes in February and March are lower than the other cold months. This is in agreement with the study [65] that suggests that only approximately 14% of LES occurs downwind of Lakes Superior and Huron during February and March. This is attributed to the fact that ice cover is most extensive over the Great Lakes in February and March [66, 67]. As ice cover formation becomes more prominent, during the winter season, the dampening of the exchange in energy and moisture fluxes occurs, diminishing the production of snowfall [9]. Thus, lake ice is a regulator of LES [68], and the change in ice cover fraction can also influence the spatial distribution and extent of LES.

From a spatial perspective, when compared to the rest of the Canadian GLB, the leeward shores of Lakes Superior and Huron-Georgian Bay exhibit the highest snowfall extremes, with a decrease in the snowfall extremes farther inland. The spatial patterns of the extreme snowfall resemble that of lake-induced snowfall near the leeward shores of Lakes Superior and Huron. LES does not fall very far inland and is dependent on wind shear at upper levels to determine the intensity and organization of snow bands. LES can produce the heaviest amounts of snow within 50 to 150 km of the leeward shore. Most moisture supplied from the lake to the upper atmosphere is removed through precipitation within this distance [1, 69, 70].

Furthermore, in February, the spatial extent of extreme snowfall has substantially decreased for Lake Huron-Georgian Bay’s snowbelt more than Lake Superior’s. In November, the onset of extreme snowfall has started to develop for Lake Superior’s snowbelt but not for Lake Huron’s. These differences in extreme snowfall behaviours are attributed to the influence of lake-induced predictor variables, which are influenced by bathymetry and topography, giving additional credence that the resultant extremes in snowfall are mostly driven by lake-induced snowfall.

The Great Lakes have different geographic locations and lake bathymetry such as size, depth, and latitudinal extent that influence the onset of the unstable season and trigger LES along each snowbelt region. The unstable season refers to the period when the lake is warmer than the ambient air temperature. Lake Superior, for example, is larger, deeper, and farther north than the other Great Lakes, and as a result, the warming season for Lake Superior begins relatively later in the spring and rarely ever attains the high LST values that are found for the lower Great Lakes. The prolonged unstable season over Lake Superior results in a longer LES season, lasting from mid November to early April [1], and is indicative of what is being observed in the current study. Lake depth is also a controlling factor when considering ice freeze-up and ice break-up dates at high latitudes [71, 72]. In shallow lakes, LST rises faster than deeper lakes, allowing for strong evaporation to occur earlier in the thaw season [73]. In the fall and winter, shallower lakes permit rapid and early cooling of LST [4]. These lakes also promote faster ice growth due to shorter thermal turnover rates, in the order of a week, which stores less heat. Thus, the spatiotemporal behaviours in snowfall extremes may vary among snowbelts because of lake bathymetry and other important factors.

Orographic lift is another factor that influences LES and can provide additional instability to an air parcel, affecting timing, and spatial distribution of precipitation [12]. Orographic lift can add approximately 13 to 20 cm of mean annual snowfall for every 30 m increase in elevation. It can also double the accumulation spawned from LES [27, 74]. The snowbelts of Lake Superior and the Tug Hill Plateau produce the highest annual snowfall precipitation although the Tug Hill Plateau is located farther south and experiences warmer air temperatures aloft than some of the other snowbelts [4]. The plateau rises 500 m above the leeward shores of Lake Ontario, and even with its modest altitude, it can influence the distribution of lake-effect precipitation off the long axis of the lake. Mean annual snowfall in this area can exceed 700 cm in Redfield, New York, on the western slope, which is more than twice the accumulation observed on the lowlands [12].

Similar behaviours in extreme snowfall found in this study, for both Lake Superior’s snowbelt and the Tug Hill Plateau, further affirm that the extreme snowfall is produced by lake-induced processes and not solely by extratropical synoptic storms. If the extremes were mainly a product of large-scale frontal systems, then spatial distributions in snowfall extremes would follow climatologically persistent frontal boundaries and squall lines as opposed to snowbelt zones. The spatiotemporal snowfall extremes have been assessed for the 36-year duration; and the next step is to discuss whether there were any trends in these snowfall extremes.

4.2. Trends in Snowfall Extremes

Previous studies have shown an increase in the frequency of winter circulation patterns (between 1951 and 1982) such as 850 mb westerly winds and superadiabatic air temperature, which are both required for the development of LES and have given rise to an increase in LES events [3, 23, 25, 29, 64]. Comparatively, the last two decades of the 20th century and the first part of the 21st century show a significant decrease in monthly snowfall totals over the Canadian region of the Great Lakes Basin [62, 63]. Recent trends in extreme snowfall over the Canadian domain of the Great Lakes are not as apparent and show less coherent spatiotemporal patterns. While, overall, the results show a decrease in extreme snowfall intensity, frequency, and duration, they are spatially inconsistent for each month and fluctuate between the snowbelts and nonlake-effect zones.

Therefore, to further identify potential changes in extreme snowfall, lake-enhanced predictor variable trends are analysed for Wawa and Wiarton. Overall, lake-wide annual maximum ice cover has been decreasing for both lakes. These results are in agreement with [62, 63], who investigated trends in monthly averaged ice cover concentration from 1980 to 2015 and showed significant decrease for Lake Superior for the months of January and February, while Lake Huron showed a significant decrease in January. Trends in monthly ice cover do not exactly align with the timing of trends in extreme snowfall. While significant reduction in ice cover occurs in January for Lake Superior, the most significant reduction in extreme snowfall occurs in December (Figure 6). Decrease in annual maximum lake ice can be attributed to the lake-wide monthly LST warming trends.

Increase in LST and subsequent decrease in maximum ice cover can exhibit an increase in energy fluxes into the lower planetary boundary layer, which can be conducive to the development of lake-enhanced snowfall. However, extensive work in [62, 63] suggests that as air temperatures warm, air parcels can hold more moisture, extending the resonance time of water vapour in the air and delaying the saturation and precipitation of an air parcel along the immediate leeward shores of the Great Lakes.

Results from the current study display a similar idea. Snowfall extremes have been decreasing within the past 36 years along the leeward shores of Lake Superior and Lake Huron, despite the increase in LST and decrease in maximum annual ice cover. A potential reason for the decrease in extreme snowfall could be attributed to the increase in both maximum and minimum monthly air temperature. As the observed weather stations’ maximum and minimum temperatures rise, the air temperature will become too warm for precipitation to fall as snow. Furthermore, an increase in both monthly maximum and minimum air temperatures suggest that the atmosphere is warming and, thereby, has the ability to retain moisture, as previously discussed.

The paper also presents trends in the direction and speed of maximum monthly wind gusts. Although there are no significant changes in direction of maximum wind gust, the cities of Wawa and Wiarton show an increase in direction of maximum wind gust for all months, with an exception to December for Wiarton. The predominant increase in these values suggests a clockwise shift in the origin location of the wind gusts. This clockwise shift is known as veering and can change the onset, location, frequency, intensity, and duration of lake-induced snowfall. For example, the Canadian leeward shores of Lake Superior and Huron are predominant snowbelt zones because of the prevalent cold and dry northwesterly winds that encounter maximum fetch across these lakes. However, if the origins of these wind directions begin to veer towards a northerly and easterly flow, the leeward snowbelt zones of Lake Superior will be predominant along the shores in Northern Michigan and Eastern Minnesota. The predominant snowbelt of Lake Huron would also shift to the shores along Eastern Michigan.

The speeds of the maximum wind gusts show no changes for all months, except for a decrease in November for the city of Wawa. All months of the cold season indicate a decrease in maximum wind gusts for the city of Wiarton, with a significant decrease in January. The speed of an air parcel will influence the fetch time over lakes. The fetch time of an air parcel will modify the air temperature over lakes and influence the development of lake-induced snowfall. A study by [75] indicates that it only takes approximately 10 minutes for an air parcel, at heights between one and 15 m above lakes, to undergo an increase in temperature modification. The longer the resonance time of an air parcel over lakes, the greater the increase in air temperature. If the air temperature warms too much (exceeds values above 0°C), this can decrease the production of lake-induced snowfall, as observed in this study.

This study has two apparent findings. Firstly, both Canadian snowbelt zones exhibit different spatiotemporal trends in snowfall extremes. For example, Lake Superior’s Canadian snowbelt exhibits significant decreases in extreme snowfall intensity, frequency, and duration, predominantly, for the months of December and January, while there are no significant changes for Lake Huron’s Canadian snowbelt. As discussed in the previous section, this is attributed to the geographic locations of the lakes, topography, lake bathymetry, and lake orientation that will influence fetch, ice cover, vertical temperature, and instability differently, thereby resulting in various spatial and temporal behaviours. It is therefore important to assess LES within each snowbelt separately, as the Laurentian Great Lakes and their locations have their own unique properties that influence LES and its ingredients such as lift, instability, and moisture.

Secondly, within each snowbelt, the spatiotemporal trends are not contiguously spatially coherent. This is evident from the resultant increase and decrease in extreme snowfall intensity seen along the northern and southern Canadian leeward shores of Lake Superior, respectively. The spatial variability within each snowbelt can result from local effects such as small-scale surface-atmosphere meteorological factors, which include localized ice cover fraction, shifts in vertical and horizontal wind velocity in the lower PBL, and storm tracks from synoptic scale systems, leading to lake-enhanced snowfall. Furthermore, changes in the frequency and intensity of cold air outbreaks can affect the vertical temperature gradient and oscillation patterns, both of which influence the production of LES. Assessing trends in Arctic air outbreak will be worth examining in future lake-effect and climate change studies.

5. Conclusion

The current study focused on examining the monthly cold season spatiotemporal extremes in snowfall intensity, frequency, and duration over the historical period of 1980 to 2015 for the Canadian GLB. The study has two folds: the first examines the 36-year period snowfall extremes, while the second explores the trends in snowfall extremes over the given period. In the first fold, extremes in snowfall intensity are derived by evaluating the 99th percentile of all daily snowfall for each cold month of the 36-year period over the GLB domain. Based on the spatiotemporal results, a snowfall intensity of 15 cm/day was found as the ideal threshold to evaluate extremes in snowfall. Counting the number of snowfall days that equaled or exceeded this S_threshold, the frequency of extreme snowfall days was determined. The maximum number of consecutive snowfall days that equaled or exceeded 15 cm/day represented the duration of extreme snowfall events. Spatiotemporal results show that the extremes in snowfall intensity, frequency, and duration are spatially coherent along the leeward shores of Lake Superior and Lake Huron-Georgian Bay. Results also indicate that snowfall extremes are most predominant in the months of December and January, suggesting that extremes in snowfall are attributed to lake-enhanced processes rather than driven solely by extratropical storms.

In the second fold, there are two apparent findings in the spatiotemporal trends of snowfall extremes. Firstly, individual snowbelts behave differently when assessing spatiotemporal trends. For instance, while portions of the Canadian leeward shores of Lake Superior show significant decreases in extreme snowfall, no changes are observed for Lake Huron’s snowbelt. Secondly, within each snowbelt, there are spatial variability and inconsistent snowfall patterns. This is evident as Lake Superior’s Canadian snowbelt shows a decrease in extreme snowfall intensity along the southern leeward shores and an increase along its northern leeward shores. The lakes’ topographical, geographical, and bathymetric features influence these results, in addition to, mesoscale and synoptic scale shifts in lake-induced predictor variables.

Trends in lake-induced predictor variables were also analysed to provide potential explanations in the resultant negative trends in extreme snowfall intensity, frequency, and duration. The predominant variable influencing the decrease in extreme snowfall is suggested to be monthly maximum and minimum air temperature. This is because, in comparison with the other predictor variables, the maximum and minimum air temperature produces the most number of significant trends. The significant increases in these temperatures for both Wawa and Wiarton suggest that these variables can influence the decrease in snowfall, as previously discussed. The study [62, 63] noted that trends in monthly snowfall totals shifted from an increase in the first half of the 20th century to a decrease during the last part of the 20th and early 21st century. However, the current study suggests that not only are the totals in monthly snowfall decreasing along the Canadian leeward Shores of the Laurentian Great Lakes for the 1980 to 2015 period but also the extremes in snowfall intensity, frequency, and duration are also decreasing, albeit the spatial coherency is difficult to delineate.

Data Availability

The data used to support the findings of this study are available from the corresponding author upon request.


This research was presented at the 2017 European Geosciences General Assembly, in Vienna, entitled Climatological assessment of spatiotemporal trends in observational monthly snowfall totals and extremes over the Canadian Great Lakes Basin.

Conflicts of Interest

The authors declare that they have no conflicts of interest.


The authors would like to thank Michele Thornton, from the Oak Ridge National Laboratory Distributed Active Archive Center (ORNL DAAC), for support with the Daymet product, as well as Nan Jiang, from Beijing Normal University, and Dr. Sofia Antonova from Alfred Wegener Institute, Potsdam, Germany, for their computational assistance. In addition, gratitude is expressed to Dr. Andrea Scott for her assistance with the uncertainty analyses and Mike Lackner for his technical mapping support, both from the University of Waterloo in Ontario, Canada. This research was conducted as part of the Canadian Network for Regional Climate and Weather Processes (CNRCWP) and funded by the Natural Sciences and Engineering Research Council (NSERC) of Canada.