Abstract

Lake ice was simulated for the province of Quebec, Canada, for both contemporary and future climate conditions using a one-dimensional thermodynamic ice model. The model was forced with NARR data (32 km) and both the daily IMS product (4 km) and the MODIS snow product (500 m) were assessed for their utility at determining lake ice phenology at the subgrid cell level (based on the 32 km NARR grid). Both products were useful for detecting ice-off; however, the MODIS product was advantageous for detecting ice-on, mainly due to the finer resolution and resulting spatial detail. The subgrid cell variability in ice-on/off dates was typically less than 2% of the mean, although it ranged up to 10% for some grid cells. The simulations were found to be within the satellite-detected subgrid cell variability: 62% of the time for ice-off and 80% of the time for ice-on. Forcing the model with future climate scenarios from the Canadian Regional Climate Model predicts the regional ice cover durations will decrease by up to 50 days from the current 1981–2010 means to the 2041–2070 means and decrease from 15 to nearly 100 days shorter from the current means to the 2071–2100 means.

1. Introduction

Lakes are an important part of the cryosphere. Lake ice cover both plays a role in and responds to climate variability and the presence (or absence) of ice during the winter months is known to have an effect on both regional climate and weather events (e.g., thermal moderation and lake-effect snow) [1]. Lake ice has also been shown to respond to climate variability, particularly changes in air temperature and snow accumulation [2]. Both long-term and short-term trends have been identified in ice phenology records and are typically associated with variations in air temperatures, while trends in ice thickness tend to be associated more with changes in snow cover [2].

Monitoring lake ice in vast areas such as Canada, particularly in the northern areas, is logistically challenging as many lakes are in remote regions with limited accessibility. Currently, there is a very limited surface-based network for monitoring lake ice in Canada. Volunteer-based lake ice observations have also been used for monitoring ice cover, for example, the “IceWatch” program, Canadian lake ice observations by citizens [3]; long-term monitoring of Crazy Lake, Nunavut by Nunavut Arctic College students [4]; in Alaska (USA) using data collected by teachers and students as part of the Alaska Lake Ice and Snow Observatory Network (ALISON) [5]. Digital camera imagery has been shown to be a useful tool for unattended lake ice monitoring [6] and while the imagery used for the aforementioned study was downloaded in situ, the technology is available for digital cameras to be connected to telemetry systems for real-time monitoring.

Remote sensing products can contribute towards monitoring lake ice cover where no in situ data exists. The Canadian Ice Service (CIS) currently monitors the fractional ice cover on 135 large lakes throughout Canada on a weekly basis using visual interpretation of satellite data. Products using a combination of several remotely sensed data sources such as the Interactive Multisensor Snow and Ice Mapping System (IMS, 4 km: [7, 8]) provide the opportunity for detection of daily ice cover extents across the entire Northern Hemisphere for larger lakes [9]. The use of optical imagery is limited in the northern regions due to frequent cloud cover and polar darkness (e.g., MODIS, AVHRR). However, microwave data provides the potential for obtaining ice phenology on a wide range of lake sizes [1012] as well as ice thickness on large lakes (using passive microwave) [13]. In latitudes where polar darkness does not occur, one product of particular interest for lake ice monitoring is the MODIS daily snow product [14, 15]. While primarily used for monitoring snow cover (e.g., [1618]), this product has a wide variety of applications including both the validation and incorporation into hydrological modelling (e.g., [1921]). MODIS surface reflectance products (MOD09A1) combined with eight-day snow cover composites (MOD10A2) [22] and calibrated radiance images/visible composites [23] have been used for identification of lake ice-on/off dates. Although lake ice and open water are identified in both IMS and the MODIS snow cover product, limited work has been done to examine the utility of extracting ice phenology from these products (e.g., [9, 24, 25]).

Satellite data presents the possibility of wide-scale monitoring of lake ice phenology, but ice thickness and composition are not readily obtainable through current techniques. Modelling can provide the opportunity to examine ice cover regimes on lakes, including the timing of breakup/freeze-up, ice thickness, and composition. Different types of models have been used to examine lake ice with varying degrees of complexity such as regression or empirical models (e.g., [26, 27]), energy balance models (e.g., [28, 29]), and thermodynamic lake ice models (e.g., [3033]). Additionally, lake ice models provide the ability to examine potential changes to the ice cover regimes through the use of future climate scenario data. Recent work using lake ice models suggests changes in the North American Arctic ice cover from the 1961–1990 mean to the 2041–2070 mean could result in a 10–60-day reduction in ice cover duration with the larger reductions found in coastal areas [34] and 15–50 days shorter ice cover duration for more southerly regions (40°N–75°N) across the Northern Hemisphere [35].

The most suitable way to undertake wide scale lake ice modeling is to force the models with climate model output or reanalysis data, rather than with in situ meteorological stations. In Canada, climate stations are typically confined to the relatively accessible areas of the mainland or coastal areas in the Arctic islands. Previous work has modelled ice phenology using grid-based climate data for hypothetical lakes of various depths (e.g., [34, 35]). In these studies, the one-dimensional models used produce the same ice phenology for the entire (reanalysis or climate model) grid cell for lakes of a given depth. However, in reality, a variety of different lake morphometric conditions could exist within a given grid cell leading to different durations of ice cover within the grid cell.

The objective of this paper is to determine how well the one-dimensional model CLIMo (Canadian Lake Ice Model) [32] is able to capture the spatial variability in lake ice-on/off dates present within a grid cell of forcing data. Specifically, this work aims to (1) examine the utility of both the IMS and MODIS snow products for identifying ice-on/off dates and determine the sub-grid cell variability using these products; (2) produce modelled ice cover which represents lake conditions within each grid cell of forcing data (rather than hypothetical lakes) and assess the model performance with respect to the sub-grid cell variability; (3) provide ice cover predictions based on future climate scenario data.

2. Study Area and Methodology

2.1. Study Area

The area of focus for this work is the province of Quebec, Canada (Figure 1). The climate in this region varies northwards from humid continental, to subarctic, to Arctic, with both continuous and discontinuous permafrost located in the far northern areas. Over 90% of the province is underlain by the Canadian Shield Precambrian rock formation and there are over 8000 lakes and reservoirs ranging in size up to 2500 km2 (93 > 100 km2, 20 > 400 km2).

Seasonal temperature increases throughout the province have been observed over the last 50 years, particularly in the summer, with the exception of two northern stations, Kuujjuaq and Schefferville, which showed nonsignificant cooling trends during the winter [36]. Rainfall also shows significant increasing trends across Canada with the most pronounced increases found in the spring and fall for Quebec [37]. Mixed trends were found for annual snowfall, although a notable decrease in winter snowfall was found in southern Quebec [37]. Trends in snow cover show a north-south gradient in maximum snow water equivalent with significant local decreases over southern Quebec and significant local increases over north-central Quebec [38]. Trends in the timing of lake ice formation and decay have also been identified in this region with varying degrees of strength. From 1961 to 1990, trends were predominantly towards later ice-on and earlier ice-off (though most were nonsignificant) and during the 1971–2000 period, trends showed earlier ice-off but only 2 stations were available for ice-on (one with an earlier trend, the other with a later trend) [39].

2.2. Subgrid Cell Lakes

The locations of lakes within the study area (including those smaller than can be resolved with satellite data used for this study) were identified using a modified version of the CanVec 1.1 files. CanVec is a digital cartographic product produced by Natural Resources Canada (obtained through GeoGratis), derived primarily from the National Topographic Data Base, with additional information from the GeoBase initiative, Landsat 7 or Spot imagery [40]. A 32 km resolution grid-based lake fraction for the province of Quebec, based on the North American Regional Reanalysis data grid, was determined including all sizes of lakes identified using the modified CanVec files (Figure 2).

2.3. Satellite Data

To examine the variability in both ice-on and off dates of the numerous lakes throughout Quebec, two satellite-derived products were used: the Interactive Multisensor Snow and Ice Mapping System (IMS, 4 km resolution Northern Hemisphere product, 2004–present [41]) and the MODIS snow product (MOD10A1, 500 m resolution global product, 2000–present [42]) both acquired through the National Snow and Ice Data Centre (NSIDC).

The IMS product provides categories of ice, water, land, and snow and is produced daily by analysts from numerous remote sensing products, mainly optical but with some passive microwave data, augmented by station data (see [7, 8]). The first date with open water for each lake pixel (change from ice to water) was used as the ice-off date, and vice versa for ice-on. IMS can provide ice phenology for lakes large enough to be resolved by the 4 km pixels [9] however, there are numerous smaller lakes throughout the province that are not resolved by the 4 km pixels.

The 500 m resolution MODIS snow product was evaluated for use in detecting lake ice phenology. The MODIS snow product (Collection 5) was obtained from the Terra sensor (MOD10A1) (2000–2011) with any missing days filled with the Aqua version of the snow cover product (MYD10A1). Snow cover is determined based on the algorithm Snow map, which uses a grouped criteria technique employing the Normalized Difference Snow Index (NDSI) and other threshold-based criterion to classify snow on a pixel-by-pixel basis [43] with an estimated accuracy of ~93%, varying with respect to snow conditions and land cover type [15]. Algorithm details can be found in the Algorithm Theoretical Basis Document (ATBD) [43, 44]. Inland water areas are also analysed for snow-covered ice conditions based on a land/water mask developed by Boston University, which includes water bodies greater than 1 km2 in dimension (removing isolated pixels). The 1 km land/water mask is applied to the four corresponding 500 m pixels in level 2 product used to generate the MOD10A1 product (MOD10_L2 swath product).

Cloud cover presents limitations in obtaining ice phenology from the daily MODIS snow product. In order to determine a “best-estimate” date for ice-off, the first day of open water was extracted for each pixel. The actual ice-off date could therefore be earlier depending on the persistence of the cloud cover. Ice-on was estimated as the day after the last day with open water being identified and could be later due to the precise ice-on date obscured by cloud cover. The extracted ice-on/off dates were then filtered using a 5 × 5 mode filter to reduce erroneous pixels found throughout the data set, likely as a result of cloud/ice confusion. Clouds are identified in the MODIS snow product using the MODIS cloud mask product and the misidentification of open water as ice may occur in situations where the water has a high turbidity, is shallow with a bright bottom, or is partially cloud-obscured resulting in the NSDI falling in the range considered to be snow or lake ice [44]. The IMS product does not exclude pixels with cloud; rather the IMS analyst provides the best estimate for conditions that day is based on previous imagery or, in the case of snow, carries the previous conditions forward until a cloud-free day is available [8].

The ranges of ice-on/off dates obtained from the MODIS snow product (2000–2010) were compared to those obtained by using the IMS data (2004–2010) as well as visual inspection of MODIS imagery in selected areas obtained from NASA LANCE Rapid Response (http://lance.nasa.gov/imagery/rapid-response/). The MODIS snow product will be referred to simply as MODIS for the purposes of this paper, and any other MODIS product used (e.g., reflectance imagery) will be identified as such. Both the lake-wide and grid-wide variability of ice-on/off dates were examined and mean difference in days was determined between them.

2.4. Lake Ice Model

The Canadian Lake Ice Model (CLIMo) is a one-dimensional thermodynamic model used for freshwater ice cover studies and is capable of simulating ice-on and ice-off dates, thickness, and composition of the ice cover (e.g., [6, 23, 34, 4547]). CLIMo has been shown to perform very well at simulating lake ice phenology when using input data that effectively represents the climate conditions at the lake, for example, from nearby meteorological towers [6, 32, 46, 48].

CLIMo has been modified from the one-dimensional sea ice model of Flato and Brown [49], which was based on the one-dimensional unsteady heat conduction equation with penetrating solar radiation presented by Maykut and Untersteiner [50], that is, where (kg m−3) is the density, (J kg-1 K−1) is the specific heat capacity, (K) is the temperature within the ice or snow, (s) is the time, (Wm-1 K−1) is the thermal conductivity, (m) is the vertical coordinate from the upper surface, positive downward, (Wm−2) is the downwelling short wave radiative energy flux, (Wm−2) is the fraction of short wave radiation flux that penetrates the top of the upper surface layer (snow, ice, or open water), is the surface albedo, and is the bulk extinction coefficient for penetrating shortwave radiation.

From this, the surface energy budget can then be calculated: where (Wm−2) is the net downward heat flux absorbed at the upper surface, is the surface emissivity, is the Stefan-Boltzmann constant (5.67 × 10−8 Wm−2  K−4), (Wm−2) is the downwelling long wave radiative energy flux, (Wm−2) is the downward latent heat flux, and (Wm−2) is the downward sensible heat flux [45, 46].

The snow layer (if present) is represented as a single layer and CLIMo has been shown to simulate the on-ice snowpack depth and melt effectively [6]. Snow ice, formed as a result of slushing on the ice surface, is created by the model if there is a sufficient amount of snow to depress the ice surface below the water level. The added mass of the water filled snow pores (slush) is added to the ice thickness as snow ice. The albedo parameterization in CLIMo is based mainly on surface type (ice, snow, or open water), surface temperatures (melting versus frozen), and ice thickness, with no distinction regarding ice composition.

CLIMo includes a fixed depth-mixed layer in the water column in order to represent an annual ice phenology cycle. When ice is present, the mixed layer temperature is fixed at the freezing point and when ice is absent, the mixed layer temperature is computed from the surface energy budget and hence represents a measure of the heat storage in the lake. The water column of shallow lakes is typically well-mixed isothermally from top to bottom during the ice-free period, permitting the mixed layer depth to be a good approximation of the effect of lake depth leading to autumn freeze-up [32]. A more detailed description of CLIMo can be found in Duguay et al. [32].

2.5. Forcing Data

CLIMo simulations were driven by National Oceanic and Atmospheric Administration (NOAA) North American Regional Reanalysis data (NARR, 32 km, 1979–2010) using air temperature, relative humidity, wind speed, cloud cover, and snow depth. Snow density was obtained from the Canadian Snow CD [51]. The snow densities were measured over land and as the density of on-ice snow has been found to be 120% higher than on-land snow nearby [52] the snow densities were adjusted accordingly.

CLIMo was also driven by future climate scenario results from the Canadian Regional Climate Model (CRCM 4.2.0, 45 km true at 60°N, provided by Consortium Ouranos). CRCM is a limited area model, originally developed at Université du Québec à Montréal (UQAM), driven at the boundaries by global climate models (GCMs) or reanalysis data. For a detailed description of CRCM see [53, 54]. Two 1961–2100 CRCM scenarios were used: driven at the boundaries with the Canadian Global Climate Model (CGCM 3.1/T47 member 4 (scenario 1) and member 5 (scenario 2)) following the SRES A2 greenhouse gas scenario. The two scenarios are the same except with a slight perturbation to the initial GCM conditions which allows the climate to evolve slightly differently, providing some insight into the interannual variability in the climate simulations. CGCM data is produced by the Canadian Centre for Modelling and Analysis (CCCma). The daily data from each scenario consisted of 2 m air temperature, humidity (specific humidity converted to relative humidity using a calculated saturated vapour pressure as a function of temperature [55] and a fixed air pressure of 1015 mb), wind speed, water equivalent of snow, snow density, and cloud cover fractions. A cold bias has been previously identified in the CRCM temperature data (e.g., [34, 56, 57]) and a bias correction was performed using the mean NARR data from 1979 to 2010 [34].

2.6. Lake Ice Simulations

The simulated ice cover was generated for both NARR (1981–2010) and CRCM (1961–2100) with and without snow cover on the ice surface in order to represent potential snow redistribution by wind across the ice surface. A range of lake depths (3 m, 5 m, 10 m, 20 m, 30 m) was simulated to represent a variety of potential lakes within each grid cell. When the depth of the lake is known and the climate conditions are well represented, ice-on has been simulated very well (e.g., [6, 32]). Lake depths are largely unknown through this region. In order to select the most appropriate depth for each grid cell from the range of depth simulations, the mean MODIS-detected ice-on for each grid cell was compared to each of the depth simulations and the most common best fit over the 11 years of satellite data was selected. Probabilities were determined for each grid cell to describe the model performance for the recent climate (2000–2010) with respect to the sub-grid cell variability from the satellite-detected lakes within the 32 km grid cells. The relative changes from the 1981–2010 CRCM-driven ice phenology means to the 2041–2070 mean and the 2071–2100 mean were extracted from the CRCM data at the central point of each NARR grid and applied to the 1981–2010 mean NARR climate data to simulate potential changes to the ice cover based on future climate scenarios.

2.7. Definitions of Terminology

Lake ice breakup and freeze-up are processes that occur over time and space as a lake may experience thaw or freeze in one area but not completely ice-free/ice-covered for several days later. Pixel-based observations are defined as ice-on/off, while lake-wide observations are two-dimensional and therefore spatial variations in terms of timing of ice-on/off may be detected. The earliest date of open water/ice detected would be defined as the beginning of the breakup/freeze-up process (i.e., a portion of the lake may still have ice if the lake includes several pixels), while water clear of ice (WCI)/complete-freeze over (CFO) are defined as the latest ice-off/ice-on date for the entire lake. Conversely, the one-dimensional model produces one value the same as resolution as the forcing data and will be referred to as ice-on/ice-off.

3. Remote Sensing of Ice Phenology

3.1. Lake Ice Phenology from IMS and MODIS Products

A comparison of the MODIS- and IMS-detected WCI/CFO dates was possible for 189 lakes in the study area (Figure 3). Spearman correlation values between the two products for the coinciding lakes were high, particularly for ice-off dates (Table 1). While ice-on dates were still significantly correlated, the lower values are indicative of not only the resolution differences between the products (with MODIS potentially identifying smaller areas remaining ice-free longer than IMS can detect) but also some suspected errors in the IMS ice-on dates (as discussed later in this section). As the earliest and latest dates detected across the lake will be more heavily influenced by resolution, the highest correlations are typically found for comparing the mean ice-on/off values for the lakes. The coefficient of variation (CV) for each lake was determined () for each year (Table 2). Little variation around the mean ice-on/off dates on the lakes was detected (~2% or less), particularly with IMS (<1%).

Thirteen water bodies (lakes or reservoirs) were selected for a more detailed comparison between the IMS/MODIS ice cover dates. The mean difference and overall mean absolute difference was determined (MODIS minus IMS) for both WCI and CFO (Table 3). WCI had an overall mean absolute difference for all 13 lakes of three days (with the largest mean difference of eight days, MODIS later than IMS). Variation is to be expected since the actual WCI date might be obscured by cloud cover for a few days. Also, the resolution differences would tend towards MODIS being slightly later than IMS as small patches of ice may remain. No pattern to the direction of the difference was found for WCI: six lakes had a mean difference indicating later WCI with MODIS, seven with IMS, though the individual years vary with no consistency. CFO had a higher average mean difference of nine days (eight lakes with late CFO for MODIS, five lakes later with IMS). The cloud cover in MODIS data may result in slightly earlier CFO dates than reality based on the current extraction technique; however, the spatial resolution differences between the satellite data products would tend towards later CFO as the MODIS pixels would be able to detect smaller areas of persistent open water. The resolution difference appears to be more influential for detecting CFO than WCI as more lakes had a tendency for later MODIS dates than with IMS.

Some of the differences between MODIS and IMS can be attributed to cloud cover and resolution. For example, at the Baskatong R., WCI was identified at April 29, 2009 with MODIS, while IMS identified WCI at May 4, 2009. Examining the MODIS reflectance imagery for the five days between the two WCI dates shows that the reservoir is obscured by cloud cover during this time. A similar scenario occurred in 2007 at L. à L’eau Claire, with IMS identifying WCI 10 days later than MODIS. MODIS reflectance data shows open water on the same day as the MODIS snow product (June 22, 2007), but cloudy days follow with only partially unobstructed views of the lake until July 2, 2007 when IMS identified ice-off. These two examples show that IMS missed the WCI on both lakes and was only able to assign the ice-free designation on the next cloud-free day (four and ten days later, resp.). It is likely that the finer resolution of MODIS provided a better “view” of the lake through partial cloud cover and was able to assign the open water category to the pixels while IMS could not.

In some cases, erroneous pixels remain in the MODIS ice cover dates after the mode filter was run, resulting in late WCI or CFO for the lake. For example, in 2005 (not shown), MODIS reflectance imagery shows what appears to be open water well before WCI is detected with the MODIS snow product for both the Baskatong and Gouin Reservoirs leading to 51 and 18 days difference between MODIS and IMS, both as a result of 1 erroneous pixel each.

For CFO, MODIS pixels are able to detect smaller areas of persistent open water. Two examples of this are as follows: L. Mistassini (Figure 4(a)) and the Manicouagan R. (Figure 4(b)). MODIS reflectance imagery shows a small area of open water remaining on L. Mistassini for several days after IMS detected CFO (IMS CFO = January 4, 2010; MODIS CFO = January 20, 2010). The Manicouagan R. has the largest difference for CFO between IMS and MODIS for all years, which occurred in 2009/2010 (IMS CFO = January 4, 2010; MODIS CFO = February 25, 2010), where imagery shows small clusters of MODIS pixels (that IMS is not able to resolve) remaining ice-free on the east side of the reservoir after CFO is indicated by IMS.

Of all 13 lakes examined in detail, L. à L’eau Claire had the highest mean difference of 32 days later for MODIS compared to IMS. In 2009, a widespread problem occurred with the IMS data where ice-on for the northern portion of the province occurred on October 29 or November 5, including L. à L’eau Claire (MODIS CFO = December 19). However, L. Bienville, just to the southwest, did not freeze until December 15 with IMS (MODIS = November 3), 40 days later. Visible imagery (Figure 5) shows open water on December 3, partial freeze-up on December 15, and December 25 appears frozen although somewhat obscured by cloud. A similar situation occurred in 2004 where a large portion of northern Quebec also experienced a uniform, early freeze-up, this time extending further southwest and including L. Bienville. Early freeze-up throughout the region again occurred in 2010/2011. While some of the discrepancies between CFO dates may be attributed to resolution, this area experiences frequent cloud cover, with more than 50% of the days in the year reported as cloud covered from 2000 to 2005 [58]. In areas with persistent cloud cover or low sun angle, IMS analysts incorporate microwave data rather than the more heavily used visible imagery [8]. If this was the case during the years with the ice-on discrepancies, the large satellite footprint for microwave data may be affecting the detection of ice cover on the lakes in this northern area.

Similar to findings regarding numerous snow cover studies comparing MODIS and IMS [16], neither product is clearly better than the other. The MODIS snow product does tend to produce better results in terms of spatial detail; however, even after a mode filter was run on the ice cover dates a few erroneous pixels remained. For studies interested in the break/freeze-up process rather than just WCI/CFO, the MODIS product outperforms the IMS product for spatial variability and detail on the lakes. In the Arctic region, however, MODIS may not be an option for freeze-up due to polar darkness and IMS is the best option available for daily lake ice cover. For ice-off, IMS performs comparably to MODIS in the Quebec region; however, ice-on was frequently too early. Early ice-on was also detected for Great Slave Lake (61°40′N, 114°W) and Great Bear Lake (66°N, 121°W) using IMS compared to ice phenology extracted from AMSR-E passive microwave data [13], suggesting that while IMS lake ice detection works well for ice-off, it experiences some limitations for the ice-on.

3.2. Grid-Wide Ice-On/Off from IMS and MODIS

While the previous section compared how well MODIS- and IMS-detected phenology dates compared to each other for entire lakes, the forcing data for the model is 32 km grid cells, so a grid-by-grid comparison was done as well for ice-on/off occurring within each grid cell. 788 of the possible 1452 grid cells covering the province contain lakes large enough to be resolved by the filtered MODIS product (≥1 km2) and 456 grid cells have lakes resolved by IMS (≥4 km2).

Comparison statistics were generated, similar to the previous lake-based data, for the coinciding grid cells with lakes identified by both products. Correlations between the products (Table 4) are lower than the lake-by-lake comparison as a result of resolution differences, with MODIS resolving smaller lakes, and more details in the larger lakes. This is particularly evident in the CFO values, though all values were still significantly correlated. The mean CVs for all years (Table 5) were fairly similar to those found in the lake comparison, with IMS having less variability in dates than MODIS. The average CV for each grid cell was determined for all available years (MODIS: 2000–2010, IMS: 2004–2010) to highlight the grid cells with the most deviation around the mean ice-on/off dates (Figure 6). The mean range in ice-off dates within the grids is 2 days for IMS (with a maximum of 97 days) and 14 days for MODIS (maximum of 116 days). The mean range in CFO dates within the grids is also 2 days for IMS (with a maximum of 76 days) and 7 days for MODIS (with a maximum of 81 days). Knowing the variability present at the sub-grid cell level provides insight into the representativeness of the grid-based simulations. The following section presents the simulated ice covers and examines the model performance with respect to the sub-grid cell variability.

4. Modelling Lake Ice Phenology

Figure 7 presents the mean ice-on/off dates (1981–2010) simulated by forcing CLIMo with NARR data across Quebec. Ice-off ranged from April 13 in the south to July 1 in the north along a northwest gradient (snow-covered data shown; snow-free simulations were an average of five days earlier). Ice-on ranged from October 7 to January 10, mainly on a northwest gradient, but dependent on the assigned depth of each grid cell, with ice-on occurring for deeper lakes later than shallow lakes under similar climate conditions.

The simulated ice-on/off dates were well correlated with those obtained from the satellite data, with all correlations significant at the 0.001 level (Table 6). Averaged over the entire region, the mean difference for all 11 years of MODIS data indicates that the simulated ice-off dates were only 1 day too early (5 days with no snow cover on the ice) and 1 day late for ice-on. Compared to IMS, the mean difference for ice-off was higher at 8 days too early (12 days too early with no snow on the ice) and 13 days early for ice-on. Mapping the average mean difference for each grid cell shows a mixture of early/late grid cells for both ice-on and ice-off with MODIS (Figure 8). The largest difference for ice-off was found in the northern portion of Quebec (near Kuujjuaq, 58° 6′N, 68°24′W). Simulations in this region range from 2 to 16 days too early.

An additional simulation was run using the Meteorological Service of Canada weather station data from Kuujjuaq to assess the NARR data for the grid cell encompassing this location. Simulations from both NARR and MSC vary similarly year to year, but MSC simulations were on average four days earlier. Both reflected the early ice-off for 9 of the 11 years simulated compared to satellite-detected ice-off for a lake 20 km away from the weather station (two years were represented well in the simulations). The amount of snow accumulation in the NARR input data is 1.6 times higher than that at the MSC station over the 30 year simulation span, although both underrepresent the normal annual snow amounts for the area (257.1 cm, 1971–2000 normal). As the MSC and NARR temperatures vary by less than 1°C average per season at this location, the difference between the simulations and the satellite-detected dates here is likely attributed to the misrepresentation of the actual snow conditions on the ice surface, as increasing the snow amount and density in the simulations delays the ice-off dates and approaches the satellite-detected dates. The area with a larger difference towards late simulations (east-central region of the province) is not likely a reflection of snow conditions, but rather more likely indicative of discrepancies between the NARR temperatures data compared to actual air temperature, as even with no snow cover on the ice, the simulations are still later than the observations.

The simulations are earlier than the IMS-detected ice cover dates for both ice-off and ice-on, more so for ice-on. The depth for each grid cell in the simulations was determined using MODIS, which can resolve the smaller shallower lakes that would freeze up earlier than the larger (and likely deeper) lakes that IMS can resolve, likely the source of the large differences for ice-on.

One lake in the study region (L. St. Jean) has a historical record of ice observations in the Canadian Ice Database [59] and provides the ability to compare all of the data sources. Figure 9 presents the long-term record of ice-on/off dates from observations (near-shore), simulations driven by NARR, and both IMS and MODIS satellite-detected dates (mean ice-on/off for the lake). The mean difference for ice-on between observations and simulations is 1 day (1979–1997), while ice-off is 4 days. No observations are available to compare to the satellite-detected dates, but simulations compare well for both ice-on and ice-off (MODIS: mean difference of 4 and 1 day, resp., IMS: mean difference of 1 and 5 days resp.).

Year-to-year variations in how well the NARR data is capturing the local climate vary. In order to represent how well the model is performing for the current climate (2000–2010: years with available satellite data), each year was tested to see if the simulations fell within ±1 standard deviations (σ) of the mean observations from all other years. A probability for each grid cell was then determined for the 11 years to represent the ability of the model to simulate the ice-on/off dates within the grid cell (Figure 10). Province-wide, the probability of the simulations for ice-off being within ±1σ of the mean observations is 62%, while ice-on is 80%. As little to no sub-grid cell variability was found with IMS, only MODIS was used to determine the probabilities. These probabilities can then be applied to the future ice cover scenarios, under the assumption that the variability captured in the 11 years of MODIS data remains the same in the future. The current data available only allows for an 11-year-time series; however, a long-time series (20 years or more) would provide greater confidence in the probabilities.

4.1. Future Ice Cover

Simulations were run using future climate scenario data from CRCM to investigate the potential changes for this region. Knowing how well the model performed for each grid cell for the contemporary climate provides some insight on the representativeness of the future predictions.

4.1.1. 1981–2010 to 2041–2070

Temperature changes between the 1981–2010 means and the 2041–2070 means predicted by CRCM are similar for spring, summer, and fall for both scenarios ranging from 1 to 3°C increase with greater warming in the south. The largest temperature increases are found in the winter, with more warming shown in the north (3-4°C increase for scenario 1 and 4-5°C for scenario 2) than the south (2-3°C increase) with the southeastern side of Hudson Bay experiencing the most warming in scenario 2 (up to 7°C). The largest changes in the area are found over Hudson Bay ranging from 7 to 9°C during the winter. Both scenarios suggest a decrease in snow water equivalent throughout most of the province, except for in the northern areas where an increase of up to 15% in the winter and spring is shown (with very small areas of < 5% increase in the fall). This is consistent with global climate model simulations of winter snow accumulation for this region, with increases in the north and decreases in the south [60]. Scenario 1 also indicates a region with a 5–10% increase around St. Lawrence River during the winter season.

The two scenarios show slightly different amounts of change, and in different regions, for both ice-on and ice-off (Figure 11). Scenario 1 shows changes in ice-off ranging from 10 to 15 days earlier, with a section of less than 5 days earlier in the central region of the province. Scenario 2 shows changes in ice-off from 15 to 25 earlier days in the south to less than 5 days in the northern half of the province (excluding an area of 5–10 days in the very north). Scenario 1 shows ice-on ranging mainly from 0 to 10 days later with a few grid cells of 10–20 days in the far southern areas. Scenario 2 shows much more pronounced change for ice-on, ranging from 20 to 25 days in the St. Lawrence region to a mixture of 5–15 days through the rest of the province. Combined, these changes to the ice-on/off dates result in reduced ICD of 5–25 days (with a few isolated pixels greater than 30 days) for scenario 1 and up to 50 days shorter for scenario 2.

Generally, scenario 2 shows more change than scenario 1. Using the same scenarios to simulate arctic ice cover changes [34] showed greater change with scenario 1 in the western Arctic and greater change in scenario 2 in the Baffin Island/Northern Quebec region. The areas of greater change for both ice-on and ice-off are the southern portion of the province in the St. Lawrence region. While previous research (not specifically in this region) has suggested greater changes to the ice-off dates [34, 35], in the case of scenario 2, greater change was found in ice-on dates. This is likely a reflection of the greater increase in fall temperatures experienced in this region for scenario 2 over scenario 1. The larger changes in both ice-on/off in scenario 2 lead to the larger reductions in ICD for this scenario.

4.1.2. 1981–2010 to 2071–2100

Extending the projected changes to the 2071–2100 means shows that temperature increases for the spring season are similar between the two scenarios, ranging from 4-5°C in the south to 2-3°C in the north. Summer is also similar with 4-5°C increase in the south, 3-4°C through most of central region, and 2-3°C in the north (with a section of 3-4°C increase at the far north for scenario 1). Scenario 1 has a fall temperature increase of 3-4°C throughout most of the province, while scenario 2 is again warming more in the fall with a 3-4°C increase in the south and 4-5°C in the north. Winter shows the greatest changes ranging from 1 to 6°C from south to north for scenario 1, and 3–7°C from south to north for scenario 2. Again, the greatest area of warming in the area is over Hudson Bay, with up to 11°C warming during the winter. Similar to the 2041–2070 means, snow water equivalent is shown to decrease throughout the province except for northern regions, which show an increase of 5−10% in the winter and up to 30% in the spring (up to 20% with scenario 1). Scenario 1 also shows a very small area of increase along St Lawrence River during the fall.

Figure 12 presents the ice cover changes for this time frame. Scenario 1 ranges from 10 to 30 days earlier for ice-off and mainly 5 to 30 days later for ice-on (with a maximum of 65 days in the southern area near St. Lawrence River). Scenario 2 ranges from 5 to more than 30 days earlier for ice-off with more pronounced changes in the St. Lawrence region than Scenario 1. Ice-on for scenario 2 also ranges from 10 to more than 30 days later (up to 42 days later), but differs spatially from scenario 1 with greater amounts of change occurring in the St. Lawrence region and in the eastern portion of the province. The resulting ICD ranges from 15 to nearly 100 days shorter, with most change occurring in the St. Lawrence region and decreasing northwards along a northwest gradient. For scenario 1, the southernmost region ranges from 60 to 94 days shorter ICD and in the St. Lawrence River region 40 to 60 days shorter. Scenario 2 has a larger area covered by large reductions in the ICD along St. Lawrence (60–72 days shorter), but less change in the northern areas than scenario 1.

5. Summary and Conclusions

Overall, the MODIS and IMS products are useful for monitoring daily ice cover changes on lakes, and while both can obtain comparable ice-off dates, MODIS outperforms IMS for ice-on due to finer resolution and the tendency towards early ice-on with IMS. Ice-off dates from both products are strongly correlated and while still significantly correlated, they were weaker for ice-on due predominantly to the resolution differences that allow for more spatial detail with the MODIS snow product. For the 13 relatively large lakes/reservoirs examined in detail, a mean absolute difference between MODIS and IMS of 3 days for ice-off and 9 days for ice-on was determined. The main differences between the products were a result of spatial resolution, cloud cover, and different methods for assigning the ice category. Whether or not lake ice is present within a pixel is at the IMS analyst’s discretion based on available data, while the MODIS algorithm is fully automated. While both MODIS and IMS have problems reported with overrepresented snow cover extent during snowmelt, attributed at least partially to the confusion of lake ice and snow [16], the change from ice-covered to ice-free for a lake (and to a lesser extent ice-free to ice-covered) is more readily distinguishable than snow versus lake ice. Some ice/cloud confusion does occur, leading to late WCI/CFO for some lakes, although filtering the MODIS dates removes the majority of these pixels. A new product has recently been developed to fill cloud gaps in the snow cover data [18] and this could be beneficial for future studies reducing the cloud interference for lake ice as well.

More variability was found in the ice-on/off dates at the sub-grid cell level using the MODIS snow product rather than IMS, as the resolution differences allow for more lakes and more detail on the larger lakes to be incorporated. The sub-grid cell variation was typically less than 2%, with a maximum of 10% of the mean value of the encompassing grid. The simulations for contemporary climate conditions were more comparable to the MODIS-detected ice phenology than IMS, again due to the incorporation of the smaller lakes in the MODIS data. An indication of whether or not the simulated ice phenology is within the 11-year sub-grid cell variability was determined. Averaged over the entire province the simulations were within the sub-grid cell variability 62% of the time for ice-off and 80% of the time for ice-on. Driving the model with the future climate scenarios predicts ice cover durations throughout the region will decrease by up to 50 days from the current 1981–2010 means to the 2041–2070 means and decrease from 15 to nearly 100 days shorter between the contemporary and 2071–2100 means.

The incorporation of remote sensing data into modelling can lead to advancements in the predictive abilities of lake/lake ice models by providing the ability to incorporate lakes at the sub grid cell scale. The inclusion of subgrid cell lakes in numerical weather predictions can contribute to improved forecast skill in some regions, particularly for spring and summer [24]. Unresolved lakes in the data used for this study (those less than 1 km2) are present throughout Quebec and are not accounted for in the grid-based analyses. The inclusion of the unresolved lakes could lead to an increase in the satellite-based variability in ice-on/off dates as the smaller lakes would likely become ice-covered sooner, and possibly ice-free earlier, than larger lakes in the same grid cell. The next generation of the CRCM model (version 5, version 4 was used in this study) incorporates a mosaic approach to parameterize the surface types using the CLASS 3.5 land surface scheme, with an additional surface-type added to represent lakes (all land water bodies). This provides, among other improvements, the ability to distinguish between open water and ice-covered parts of a lake for calculations such as latent and sensible heat fluxes [61].

One of the limitations in modelling ice cover is the lack of known lake depths, a problem particularly evident in Canada where numerous lakes are situated in northern or remote regions. Since lake depth is a controlling factor for ice-on, correctly representing the depth in the model is of key importance. Lake depths are often assumed in modelling at 10 m [62] or accounted for by numerous depth simulations [34]. In this study, we used an approach whereby the best depth for each grid cell was determined based on satellite data in the absence of measured depths, which could lead to some discrepancies between the simulations and satellite-detected dates for ice-on. Deviations in the NARR data from actual climate conditions, particularly with respect to temperatures in the summer and fall seasons, may have affected the selection of depths by advancing or delaying the ice formation. A global database and raster product of lake depths has been created for use in weather forecasting and climate modelling and is being implemented in several NWP and climate models [63]. A dataset of this nature could be very beneficial for lake ice modeling; however, with lake depths unknown throughout many parts of Canada there is still potential for the misrepresentations of depths.

Furthermore, an accurate representation of snow on the ice is essential for lake ice simulations. Reanalysis data or climate model output may not be able to represent the local snow conditions well within a grid cell and redistribution across the lake surface is difficult to capture. Using snowfall records to represent depth of snow on the ground is sometimes necessary, but weather station records have a tendency for snow undercatch [64]. The ability to obtain variables such as on-ice snow depths at fine spatial resolutions, such as that proposed by the Cold Regions Hydrology High-Resolution Observatory (CoReH20) candidate satellite mission of the European Space Agency [65], will be extremely beneficial for accurate lake ice modelling.

Acknowledgments

This research was supported by a NSERC Postgraduate Scholarship (CGS) to L. C. Brown and Discovery Grant to C. R. Duguay. The CRCM data was generated and supplied by Ouranos (with thanks to Ross Brown). The National Snow and Ice Data Center is acknowledged for providing access to the IMS and the MODIS snow cover products as well as NOAA for providing the NARR data (ftp://ftp.cdc.noaa.gov/).