Evaluating RADARSAT-2 for the Monitoring of Lake Ice Phenology Events in Mid-Latitudes

: Lake ice is an important component in understanding the local climate as changes in temperature have an impact on the timing of key ice phenology events. In recent years, there has been a decline in the in-situ monitoring of lake ice events in Canada and microwave remote sensing imagery from synthetic aperture radar (SAR) is more widely used due to the high spatial resolution and response of backscatter to the freezing and melting of the ice surface. RADARSAT-2 imagery was used to develop a threshold-based method for determining lake ice events for mid-latitude lakes in Central Ontario from 2008 to 2017. Estimated lake ice phenology events are validated with ground-based observations and are compared against the Moderate Resolution Imaging Spectroradiometer (MODIS band 2). The threshold-based method was found to accurately identify 12 out of 17 freeze events and 13 out of 17 melt events from 2015–2017 when compared to ground-based observations. Mean absolute errors for freeze events ranged from 2.5 to 10.0 days when compared to MODIS imagery while the mean absolute error for water clear of ice (WCI) ranged from 1.5 to 7.1 days. The method is important for the study of mid-latitude lake ice due to its unique success in detecting multiple freeze and melting events throughout the ice season.


Introduction
Changes detected in the Northern Hemisphere (i.e., Central North America, Northern Europe, and Russia) lake ice conditions since 1855 have shown a delay in ice on dates by 10.7 days and earlier ice off dates by 8.8 days [1]. This extension of the open water season has implications for regional climate due to the impact lakes have on energy fluxes [2] and is also of interest for the study of under-ice biological community structure and the impact of warming temperatures on these communities [3,4]. Lake ice also plays an important role in recreational activities such as ice fishing and snowmobiling and serves economic purposes for northern communities that rely on the presence of ice roads, which makes the existence of a consistent record of lake ice important [5].
Although there have been observed changes in lake ice cover in recent years, there has been a decline in governmental in situ lake ice monitoring programs, which resulted in a paucity of lake ice observation records for Canada [6]. Currently, there are two main sources of lake ice data in Canada: a weekly product produced by the Canadian Ice Service (CIS) for~136 lakes [7] across the country with ice charts for Great Lakes (https://www.canada.ca/en/environment-climate-change/ services/ice-forecasts-observations/latest-conditions.html) and the IceWatch program where data are collected by volunteers [8,9]. While beneficial, both programs have spatial and temporal limitations with respect to the number of lakes monitored and the frequency of observations. The implementation variability of mid-latitude lake ice phenology as estimated by RADARSAT-2. To our knowledge, this is the first time that backscatter evolution has been reported for these lakes in Central Ontario and this is the first attempt to develop a threshold-based method for detecting ice phenology for lakes of this size in mid-latitudes that captures multiple mid-thaws experienced at this latitude.

Study Area
The geographic focus of this research is Central Ontario, Canada with a specific focus on lakes in the Haliburton region ( Figure 1). The 1981-2010 climatology for the region indicates an average temperature of −9.9 °C in January (coldest month) to an average temperature of 18.7 °C in August (warmest month) [28]. The Haliburton area receives an average of 1073.5 mm of precipitation a year with 26% (279.6 mm) falling as snow. The largest amounts of precipitation typically occur in the Fall (September to November) with the majority of the precipitation during this season falling as rain.
For this study, 15 lakes were selected in total ranging from 1.6 to 70.5 km 2 (80 to 2,000,050 meter pixels) with a latitudinal range from 45° to 45.7° (Figure 1, Table 1). The three primary lakes used for evaluating the accuracy of SAR phenology detection were Johnson Lake, MacDonald Lake, and Clear Lake, which was located in the Haliburton Forest and the Wildlife Reserve Ltd. (Figure 1) and was selected due to the availability of hourly in situ data.

RADARSAT-2
From October 2008 to April 2017, 382 RADARSAT-2 images were acquired at the C-band (5.405 GHz) HH polarization in ScanSAR Wide mode with a nominal spatial resolution of 100 m. Although some dual-polarized images were collected (HH + HV), only the HH polarization was used in order to establish the longest image record possible. Some images have different collection paths but were selected in order to provide the most complete record for the area around the Haliburton Forest and Wildlife Reserve. Images with both descending (acquisition time:~11:30 UTC,~07:30 EST) and ascending orbits (acquisition time:~23:10 UTC,~19:10 EST) were used in order to minimize the gaps in the image record as much as possible. The average gap between images was 3.5 days during the ice season excluding a large gap of 41 days between 21 December 2012 and 31 January 2013.
The RADARSAT-2 images were calibrated using the absolute radiometric calibration process in order to convert the raw digital numbers (DN) to sigma-naught (dB) [29,30], which is useful for detecting lake ice phenology events [12,19]. Additionally, the images were speckle-filtered by using a Refined Lee filter that determines the relation of a pixel to the orientation of an edge using a 7 × 7 window and 3 × 3 subarea window and uses the local measures of mean and variance to filter the images [31]. This process is necessary in order to reduce speckle that is present in all radar images [32]. Images were also terrain corrected by using a digital elevation model (DEM) to align the images and correct for elevation interference in order to ensure that the SAR images could be overlaid with other cartographic files [33,34]. The DEM used for this purpose was the Shuttle Radar Topography Mission (SRTM) 1 s grid (30 m), which was recorded by the space shuttle Endeavour from 11-22 February 2000 [35]. Seventy-three images from 2008-2010 had to be shifted due to the alignment issues between scenes. These shifts ranged from −0.013 to 0.003 decimal degrees and were identified visually by matching lake features such as peninsulas between images.
Images acquired by using the RADARSAT-2 ScanSAR wide mode contains a range of incidence angles between 20 and 49 • . Previous studies on both lake and sea ice have shown that incidence angle impacts the HH backscatter from different surfaces and the incidence angle also differs between images due to the combination of ascending and descending orbits used in this dataset [23,36,37]. Lower incidence angles result in increased interaction at the surface of the ice and more surface scattering while higher incidence angles are affected by volume scattering and characteristics of the ice column [23,[37][38][39]. In order to normalize the backscatter values for the images acquired in this study, 1094 points were randomly generated for lakes in Central Ontario between 2008 to 2016. The backscatter was extracted from each image for every point (n = 308,468). There was a minimum distance of 500 meters between each point to ensure that each point has a unique value for each image. Linear regression models were used to examine the impact that the incidence angle had on the HH backscatter. Sixty random locations (out of the 1094) were selected in order to provide samples for developing the linear model, which resulted in over 15,000 distinct data values from the different images. These data were split into 70% training (n = 11,827) and 30% validation (n = 4991). This resulted in a slope of 0.258 dB/ • (R 2 = 0.16, p < 0.01, RMSE = 4.2), which is slightly higher than values used Remote Sens. 2018, 10, 1641 5 of 26 in previous research [37]. This slope was used to normalize backscatter values for all images in the dataset to an incidence angle of 39 • , which was the median of the extracted incidence angle values.

In Situ Data
A meteorological tower was constructed on the shore of MacDonald Lake (Figure 1) for collecting temperature, relative humidity, shortwave/longwave radiation, precipitation, snow depth, and air pressure since the fall of 2015 (averaged or sampled on the hour). Additional meteorological data from the nearest government weather towers were acquired from the Environment and Climate Change Canada's historical dataset for ice seasons between 2009 to 2015 and for any gaps in the record during the 2016 and 2017 ice seasons (available at http://climate.weather.gc.ca/historical_ data/search_historic_data_e.html). The majority of data were acquired for the Haliburton 3 station (78.53 • W, 45.03 • N, 22.6 km south of MacDonald Lake tower), which was the closest to all of the lakes selected for study ( Figure 1). The Bancroft Auto station (77.88 • W, 45.07 • N, 57.7 km south-east of MacDonald Lake tower) was used to fill in any gaps that existed in the Haliburton 3 station data ( Figure 1). These temperature data were also used to create seasonal means (fall-winter means used September, October, November, and December (SOND) daily mean temperatures and spring means used March, April, and May (MAM) daily mean temperatures) for correlation analysis with freeze onset dates and water clear of ice (WCI) dates.
A network of 9 Reconyx Hyperfire cameras with a resolution of 3.1 megapixels were set up surrounding MacDonald, Clear, and Johnson Lakes in the fall of 2015 ( Figure 1). Cameras were lashed to trees and were pointed towards the lakes to record the full season of ice formation and decay. Each camera has a 40-degree field of view and took one picture every hour from 06:00 EST to 18:00 EST by automatically switching from RGB images to infrared images in low/no light conditions. In addition to the shoreline cameras, a Campbell Scientific CC5MPX camera with a 5-megapixels resolution was mounted on the weather station pointed towards MacDonald Lake, which provides additional in situ data. All of this data were used to validate the results of the ice phenology detection.

MODIS Data
MODIS imagery was used to provide validation data in lieu of observations for the years between 2008 to 2015 prior to the installation of the in situ camera network. The MODIS images were viewed by using the MODIS worldview web application (https://worldview.earthdata.nasa.gov) and further visual inspection was conducted using the MOD09GQ product (ground surface reflectance, primarily band 2 (841-876 nm)) with a spatial resolution of 250 m [40]. Band 2 was selected due to the near infrared band allowed for the clearest discrimination of ice and water, which makes it optimal to provide a reference for ice that was first visible at the beginning of the ice season (freeze onset) and when there was no ice cover remaining at the end of the season (water clear of ice). This data source is the best continuous record of image coverage during the ice season. However, the MODIS images lack the confidence of the in situ camera network and there can be temporal gaps due to instances of cloud cover. Mean bias error and mean absolute error (Equations (1) and (2) [41][42][43]) were used to compare the results of the ice phenology detection and the 'observations' obtained from the MODIS imagery where n is the number of lakes used in the analysis.

Mean Absolute Error
(1) Freeze onset dates in 2009 and dates detected for Bernard Lake and Lake Opeongo were excluded from both MODIS comparisons due to temporal limitations in the RADARSAT-2 image record.

Mid-Latitude Lake Ice Phenology Detection
The formation of ice results in an increase in backscatter from the lake ice surface, which results in lighter tones in SAR imagery [15,18,20,23] (Figure 2). The melt process results in a wetter surface and the formation of melt ponds, which causes backscatter to decrease and results in darker tones in SAR imagery [23] (Figure 2). WCI has a similar response to melt as the open water surface results in lower backscatter [12]. In order to establish representative phenology thresholds for mid-latitude lake ice, we first examine the temporal evolution of the RADARSAT-2 backscatter for 15 different lakes from 2015-2017 when the temporal resolution was the highest over the study period. The mean backscatter was calculated by using a 150 m buffered mask from each lake shoreline.

Mid-Latitude Lake Ice Phenology Detection
The formation of ice results in an increase in backscatter from the lake ice surface, which results in lighter tones in SAR imagery [15,18,20,23] (Figure 2). The melt process results in a wetter surface and the formation of melt ponds, which causes backscatter to decrease and results in darker tones in SAR imagery [23] (Figure 2). WCI has a similar response to melt as the open water surface results in lower backscatter [12]. In order to establish representative phenology thresholds for mid-latitude lake ice, we first examine the temporal evolution of the RADARSAT-2 backscatter for 15 different lakes from 2015-2017 when the temporal resolution was the highest over the study period. The mean backscatter was calculated by using a 150 m buffered mask from each lake shoreline.  Figure 3 illustrates the evolution of average backscatter from freeze onset and the transition from open water to the ice cover to WCI for a subset of six lakes in our study area with the remaining lakes shown in Appendix A. In December, the backscatter is relatively low with values typically ranging from −20 to −25 dB. There are some spikes observed in the backscatter during late November and early December that were likely the result in brief freezing events (e.g., 19 December 2015, seen in shoreline cameras) or windy conditions (e.g., 3 December 2016, seen in shoreline cameras) ( Figure 3). Overall, the backscatter patterns observed for the Central Ontario lake ice during freeze up to complete ice cover correspond well with previous SAR observations conducted in the Arctic and midlatitudes (Northern Montana, 48° N) [15,18], but a consideration for Central Ontario mid-latitude lakes was required to account for, and subsequently detect, the frequently observed freeze-thaw events.

Freeze Onset Threshold Physical Basis
Once the ice cover has completely formed, the backscatter begins to increase, which is attributed to the continued thickening of ice and the subsequent return of more energy back to the sensor and eventually backscatter remains relatively stable for the rest of the season (Figure 3). Values remained  Figure 3 illustrates the evolution of average backscatter from freeze onset and the transition from open water to the ice cover to WCI for a subset of six lakes in our study area with the remaining lakes shown in Appendix A. In December, the backscatter is relatively low with values typically ranging from −20 to −25 dB. There are some spikes observed in the backscatter during late November and early December that were likely the result in brief freezing events (e.g., 19 December 2015, seen in shoreline cameras) or windy conditions (e.g., 3 December 2016, seen in shoreline cameras) ( Figure 3). Overall, the backscatter patterns observed for the Central Ontario lake ice during freeze up to complete ice cover correspond well with previous SAR observations conducted in the Arctic and mid-latitudes (Northern Montana, 48 • N) [15,18], but a consideration for Central Ontario mid-latitude lakes was required to account for, and subsequently detect, the frequently observed freeze-thaw events.

Freeze Onset Threshold Physical Basis
Once the ice cover has completely formed, the backscatter begins to increase, which is attributed to the continued thickening of ice and the subsequent return of more energy back to the sensor and eventually backscatter remains relatively stable for the rest of the season (Figure 3). Values remained low (−20 to −25 dB) during late December until the start of ice formation, which corresponds with Remote Sens. 2018, 10, 1641 7 of 26 the increase seen in the dB values ( Figure 3). The backscatter ranged from −19.2 to −21.8 dB during these early stages of ice formation. The observed backscatter increased until full ice cover was established and a consistent average backscatter is reached, which was found to range from −14.9 to −18.4 dB. For Arctic lakes, research from Reference [20] reported a decrease in backscatter during late December/early January for Tazlina Lake in Southern Alaska, which is attributed to the formation of stable ice cover during this time. Backscatter values during the beginning of the ice season are different between Arctic/Subarctic and Central Ontario lakes because no decrease in backscatter during the initial formation of ice was observed most likely due to there being less wind action on lakes at mid-latitudes. Therefore, the threshold method for small and medium lakes in mid-latitudes uses an initial increase in backscatter to indicate the formation of ice while, in the Arctic/Subarctic region, a decrease may be more appropriate. low (−20 to −25 dB) during late December until the start of ice formation, which corresponds with the increase seen in the dB values ( Figure 3). The backscatter ranged from −19.2 to −21.8 dB during these early stages of ice formation. The observed backscatter increased until full ice cover was established and a consistent average backscatter is reached, which was found to range from −14.9 to −18.4 dB. For Arctic lakes, research from Reference [20] reported a decrease in backscatter during late December/early January for Tazlina Lake in Southern Alaska, which is attributed to the formation of stable ice cover during this time. Backscatter values during the beginning of the ice season are different between Arctic/Subarctic and Central Ontario lakes because no decrease in backscatter during the initial formation of ice was observed most likely due to there being less wind action on lakes at mid-latitudes. Therefore, the threshold method for small and medium lakes in mid-latitudes uses an initial increase in backscatter to indicate the formation of ice while, in the Arctic/Subarctic region, a decrease may be more appropriate.

Melt Onset Threshold Physical Basis
The mean backscatter range observed for Central Ontario during the melting period was −18.6 to −21.1 dB. This was lower than the mean backscatter observed during the full ice cover varying

Melt Onset Threshold Physical Basis
The mean backscatter range observed for Central Ontario during the melting period was −18.6 to −21.1 dB. This was lower than the mean backscatter observed during the full ice cover varying from −14.9 to −18.4 dB (Figure 3), which is also reflected in the transition between brighter and darker tones shown in Figure 2. This decrease in backscatter values from full ice cover to the melting period occurs between the end of February and the end of March (depending on the year) and can be used to identify melting events ( Figure 3). The decrease in backscatter captures the increased water content in the melting snow and ice cover. However, fluctuations between lower backscatter and higher backscatter values are observed throughout the ice season, which indicates the occurrence of multiple thaw and refreeze events as shown in previous research [18]. These freeze-thaw events are associated with backscatter decreases and increases and follow similar transition values as the primary freeze and melting events.

Water Clear of Ice Threshold Determination
For Central Ontario lake ice, the backscatter decreases are similar during melt onset and (WCI), but WCI backscatter values return to the low open water values that are typically observed prior to the beginning of the ice season ( Figure 3). Large decreases in backscatter are not evident when WCI occurs (since the melt process and decreases in backscatter have already begun). Since the backscatter response for melting and WCI are similar and the evolution in backscatter is nearly indistinguishable (except for the return to a consistent backscatter after WCI), it is appropriate to use the last decrease in backscatter as the WCI date for mid-latitude lakes. The difference in backscatter evolution between mid-latitudes and Arctic/Subarctic latitudes is that, for small and medium lakes in mid-latitudes, backscatter does not return to high values after WCI, which is noted in Reference [20]. In Central Ontario, the values remain low throughout the year, which is most likely due to less wind and a smaller fetch on the lakes ( Figure 3). Figure 4 illustrates the application of the threshold-based method applied to RADARSAT-2 imagery acquired over MacDonald Lake for the 2015/2016 ice season and Figure 5 provides the processing chain used to identify freezing events, melting events, and WCI events. However, it does not discriminate between ice and water and instead focuses on the identification of the listed events. The first step determined the change in backscatter between consecutive images on a pixel-by-pixel basis. In order to detect if there is a phenology event, threshold values were determined based on our previous analysis of evolution of the mean backscatter from the 15-different lakes and applied the equations below.

Ice Phenology Detection
where σ • Th(Melt) and σ • Th(Freeze) represent the change thresholds for melt and freeze, respectively. σ • freeze , σ • melt , and σ • ice represent the respective mean backscatter values during that stage of the ice season and σ • St Dev. represents the standard deviation in backscatter for the respective time. Mean backscatter values were calculated by using all lake pixels within a 150 m buffer that included both ice and open water pixels (depending on the image). To determine the thresholds, the midpoint (transition) value from full ice cover to the melting period (Equation (3)) or from the initial stages of freezing (including instances of open water and some initial ice formation along the shoreline) to full ice cover (Equation (4)) is calculated and the mean backscatter during either full ice cover (Equation (3)) or initial stages of freezing (Equation (4)) is subtracted to determine the decrease (melt) or increase (freeze) in backscatter that will define the transition. This type of threshold approach has been successfully applied to Arctic sea ice by the study from Reference [37]. For the optimal calculation of thresholds using the equations above, in situ data is needed in order to accurately determine the timing of the difference stages (freeze, full ice cover, and melt) during the ice season. These thresholds were developed to ensure the ability to capture multiple melt-freeze events (which corresponding to distinct slushing events) during the ice season. These mid-season thaws are uncommon in the Arctic and the ability to capture multiple melt-freeze events is not factored into previous methods that have been applied at higher latitudes [12].
Remote Sens. 2018, 10, x FOR PEER REVIEW 9 of 26 and the ability to capture multiple melt-freeze events is not factored into previous methods that have been applied at higher latitudes [12]. To detect freeze events, every pixel in the backscatter difference rasters (created in the previous processing step, see Figure 5) was examined and it was noted if there was an increase greater than or equal to the freeze threshold (σ°Th(Freeze) = 2.4 dB). This was established by using the mean backscatter values during the initial freeze and full ice cover in the 2015/2016 and 2016/2017 ice seasons. In Figure  4, it is apparent that increases in backscatter greater than the threshold correspond to lower temperatures, which indicates possible freeze events. The increase of 2.  To detect freeze events, every pixel in the backscatter difference rasters (created in the previous processing step, see Figure 5) was examined and it was noted if there was an increase greater than or equal to the freeze threshold (σ • Th(Freeze) = 2.4 dB). This was established by using the mean backscatter values during the initial freeze and full ice cover in the 2015/2016 and 2016/2017 ice seasons. In Figure 4, it is apparent that increases in backscatter greater than the threshold correspond to lower temperatures, which indicates possible freeze events. The increase of 2.4 dB captures the transition from water to ice for 6 of 8 years from the 2009/2010 to 2016/2017 ice seasons (the 2008/2009 ice season is excluded due to only three images being available in December) and the remaining year thresholds range from 2.1 to 2.3 dB. This method is similar to that used for the detection of lake ice freeze and melt events in northern Alaska using SAR and for large lakes in northern Canada (Great Bear and Great Slave Lake) using QuikSCAT data where thresholds were established for freeze onset, melt onset, and WCI [12,44].
To detect melt events and estimate WCI, every pixel in the backscatter difference rasters was examined and it was noted if there was a decrease less than or equal to the melt threshold (σ • Th(Melt) = −1.9 dB). As evident in Figure 4, decreases in backscatter that are less than −1.9 dB align with warmer temperatures, which indicates melt events. Similar patterns are seen for the lakes that are not shown in Figure 4. This threshold accurately captures the transition for five out of nine seasons from 2008/2009 to 2016/2017 and the remaining yearly thresholds range from −1.5 to −1.8 dB. If the increase (decrease) in backscatter is detected as being above (below) the appropriate threshold, the pixel is assigned a value equal to the date of the image acquisition date.
The final step of the threshold-based method was to determine the occurrence of events throughout the ice season. Using a 150 m buffer to limit shoreline interference, the total number of pixels in each image that reported an event were tallied and used to determine the percentage of pixels for each study lake that reported an event. The percentages were filtered by each individual lake so that events were only considered if ≥40% of pixels reported an event. This value was selected as it was determined to be the third quartile of all percentages calculated for freeze and melt events between 2009 to 2017. For the 2015/2016 and 2016/2017 ice seasons, the analysis was focused on events that occurred on two or more of the main study lakes (MacDonald, Clear, and Johnson Lakes). These events were classified as major events. Minor events (events that occurred on only one of the main study lakes) are also mentioned briefly in the results section. When the threshold-based method is used without the addition of shoreline observations, minor events are determined for all lakes. However, the addition of shoreline observations allows for further refinement to major events.
Remote Sens. 2018, 10, x FOR PEER REVIEW 10 of 26 range from 2.1 to 2.3 dB. This method is similar to that used for the detection of lake ice freeze and melt events in northern Alaska using SAR and for large lakes in northern Canada (Great Bear and Great Slave Lake) using QuikSCAT data where thresholds were established for freeze onset, melt onset, and WCI [12,44].
To detect melt events and estimate WCI, every pixel in the backscatter difference rasters was examined and it was noted if there was a decrease less than or equal to the melt threshold (σ°Th(Melt) = −1.9 dB). As evident in Figure 4, decreases in backscatter that are less than −1.9 dB align with warmer temperatures, which indicates melt events. Similar patterns are seen for the lakes that are not shown in Figure 4. This threshold accurately captures the transition for five out of nine seasons from 2008/2009 to 2016/2017 and the remaining yearly thresholds range from −1.5 to −1.8 dB. If the increase (decrease) in backscatter is detected as being above (below) the appropriate threshold, the pixel is assigned a value equal to the date of the image acquisition date.
The final step of the threshold-based method was to determine the occurrence of events throughout the ice season. Using a 150 m buffer to limit shoreline interference, the total number of pixels in each image that reported an event were tallied and used to determine the percentage of pixels for each study lake that reported an event. The percentages were filtered by each individual lake so that events were only considered if ≥40% of pixels reported an event. This value was selected as it was determined to be the third quartile of all percentages calculated for freeze and melt events between 2009 to 2017. For the 2015/2016 and 2016/2017 ice seasons, the analysis was focused on events that occurred on two or more of the main study lakes (MacDonald, Clear, and Johnson Lakes). These events were classified as major events. Minor events (events that occurred on only one of the main study lakes) are also mentioned briefly in the results section. When the threshold-based method is used without the addition of shoreline observations, minor events are determined for all lakes. However, the addition of shoreline observations allows for further refinement to major events.

Validation of 2016 and 2017 Freeze Events for MacDonald, Clear, and Johnson Lakes
Potential major freeze events during the 2015/2016 ice season were identified ( Figure 6) and five of these events were validated by using shoreline observations and meteorological data ( Table 2).
The freeze event on 19 December 2015 was captured and the shoreline cameras support the occurrence of the event as frost flowers and a rough ice surface was visible (Figure 7). Freeze onset events were detected on 1 January 2016 and 2 January 2016 for all lakes. These events coincide with the formation of ice along the shorelines that is visible in the shoreline images ( Figure 6). Events for all lakes were captured on 25 February 2016 ( Figure 6) by following a brief increase in temperature to 2.7 • C after the previous image (19 February). The warmer temperatures caused visible snow melting in the cameras and may have resulted in changes in the surface characteristics of ice and higher backscatter values.
Remote Sens. 2018, 10, x FOR PEER REVIEW 11 of 26 Potential major freeze events during the 2015/2016 ice season were identified ( Figure 6) and five of these events were validated by using shoreline observations and meteorological data ( Table 2). The freeze event on 19 December 2015 was captured and the shoreline cameras support the occurrence of the event as frost flowers and a rough ice surface was visible (Figure 7). Freeze onset events were detected on 1 January 2016 and 2 January 2016 for all lakes. These events coincide with the formation of ice along the shorelines that is visible in the shoreline images ( Figure 6). Events for all lakes were captured on 25 February 2016 ( Figure 6) by following a brief increase in temperature to 2.7 °C after the previous image (19 February). The warmer temperatures caused visible snow melting in the cameras and may have resulted in changes in the surface characteristics of ice and higher backscatter values.  Potential major freeze events during the 2015/2016 ice season were identified ( Figure 6) and five of these events were validated by using shoreline observations and meteorological data ( Table 2). The freeze event on 19 December 2015 was captured and the shoreline cameras support the occurrence of the event as frost flowers and a rough ice surface was visible (Figure 7). Freeze onset events were detected on 1 January 2016 and 2 January 2016 for all lakes. These events coincide with the formation of ice along the shorelines that is visible in the shoreline images ( Figure 6). Events for all lakes were captured on 25 February 2016 (Figure 6) by following a brief increase in temperature to 2.7 °C after the previous image (19 February). The warmer temperatures caused visible snow melting in the cameras and may have resulted in changes in the surface characteristics of ice and higher backscatter values.   shows the presence of frost flowers and a rough lake surface, which possibly results in a higher backscatter value. This is compared to the later image (right), which shows a relatively smooth surface. Both images were acquired from the same camera on Johnson Lake. Black squares represent images that were acquired in the PM and grey squares represent images acquired in the AM. The blue lines indicate when a major freezing event was determined and the solid lines are the dates of the sample images. The 13 December image illustrates the initial conditions captured by RADARSAT-2 and the 2 January images illustrates full ice cover on the lake following melting on 27 December. The 13 December photo was acquired from MacDonald Lake and the 2 January photo was acquired from Clear Lake.
Freeze events were also detected in the 2015/2016 ice season on 6 and 14 April 2016. The freeze event that was detected on April 6 temperatures prior to image acquisition were 2.7 °C and the event follows a period where maximum daily temperature ranged from −4.4 to 9.1 °C and melt ponds were visible on the surface of the ice for MacDonald, Johnson, and Clear Lakes. The event on 14 April occurred at the beginning of a period of warmer temperatures. However, shoreline images showed that the surface of the lake may have become rough following the melt and refreeze from the previous day ( Figure 6). This event showed the impact that diurnal fluctuations in temperature can have on the detection of events. In addition to the events detected for multiple lakes, a freeze event was detected for only Johnson Lake on 5 February following periods of warmer temperatures (0.3 to 5.7 °C) between 1 February and 4 February.
Eight potential major events were identified during the 2016/2017 ice season for the three primary study lakes ( Figure 8) and six events were validated by using shoreline observations and meteorological data ( Table 2). The first identified freeze on 13 December 2016 matches the rougher surface conditions of the ice seen in the shoreline cameras. This would have resulted in a higher mean backscatter of −14.9 dB (Figure 8). This event followed the initial formation of skim ice on 10 December (−24.7 dB, surface conditions not shown). Another freeze event was detected on 2 January 2017 and was confirmed by the presence of ice cover in shoreline images (Figure 8) following melting that occurred on 26/27 December where maximum temperatures reached 3.9 °C. Additional freeze events were detected on 6 January and 26 January in 2017. Both of these events are confirmed as daily Black squares represent images that were acquired in the PM and grey squares represent images acquired in the AM. The blue lines indicate when a major freezing event was determined and the solid lines are the dates of the sample images. The 13 December image illustrates the initial conditions captured by RADARSAT-2 and the 2 January images illustrates full ice cover on the lake following melting on 27 December. The 13 December photo was acquired from MacDonald Lake and the 2 January photo was acquired from Clear Lake.
Freeze events were also detected in the 2015/2016 ice season on 6 and 14 April 2016. The freeze event that was detected on April 6 temperatures prior to image acquisition were 2.7 • C and the event follows a period where maximum daily temperature ranged from −4.4 to 9.1 • C and melt ponds were visible on the surface of the ice for MacDonald, Johnson, and Clear Lakes. The event on 14 April occurred at the beginning of a period of warmer temperatures. However, shoreline images showed that the surface of the lake may have become rough following the melt and refreeze from the previous day ( Figure 6). This event showed the impact that diurnal fluctuations in temperature can have on the detection of events. In addition to the events detected for multiple lakes, a freeze event was detected for only Johnson Lake on 5 February following periods of warmer temperatures (0.3 to 5.7 • C) between 1 February and 4 February.
Eight potential major events were identified during the 2016/2017 ice season for the three primary study lakes ( Figure 8) and six events were validated by using shoreline observations and meteorological data ( Table 2). The first identified freeze on 13 December 2016 matches the rougher surface conditions of the ice seen in the shoreline cameras. This would have resulted in a higher mean backscatter of −14.9 dB (Figure 8). This event followed the initial formation of skim ice on 10 December (−24.7 dB, surface conditions not shown). Another freeze event was detected on 2 January 2017 and was confirmed by the presence of ice cover in shoreline images (Figure 8) following melting that occurred on 26/27 December where maximum temperatures reached 3.9 • C. Additional freeze events were detected on 6 January and 26 January in 2017. Both of these events are confirmed as daily temperatures returned to below 0 • C following temperatures above 0 • C (−1.0 to 0.6 • C on 3 January and maximum hourly temperatures of 1.1 • C between 18 January and 20 January) during previous acquisition dates. Temperatures are −15.2 • C and −1.2 • C prior to image acquisition on 6 January and 26 January, respectively. Diurnal fluctuations resulted in a freeze event on 2 April in 2017. The event was validated by using meteorological data, which showed warmer temperatures on 1 April (2.7 • C) before temperatures dropped to lows of −2.7 • C overnight. There was an additional freeze event that occurred on 9 April that is also confirmed with temperature data, which shows an average temperature of −2.9 • C prior to the acquisition of the image. In addition to the events detected for multiple lakes, a minor freeze event on MacDonald Lake was detected on 13 January following a return to colder temperatures (−3.3 to −10.4 • C) in the evening and morning of 12 January and 13 January. Table 2. Summary of the validation of major freeze events for the 2015/2016 and 2016/2017 ice seasons. "Validated" events were freeze events confirmed by shoreline observations and climate data. "Not Validated" events were identified as freeze by the threshold method but did not match observations in shoreline cameras or local climate data. "Not Captured" events indicate that no freeze event was reported for that lake on the specific date.

Event Date
Johnson Lake MacDonald Lake Clear Lake Freeze events were also identified that were likely erroneous, according to shoreline observations and local climate data. There was a false freeze event recorded on 9 December 2015 since no ice was present on the lake and the shoreline cameras suggest detection was likely the result of a wind action given the roughness created by waves. On 9 January 2016 and 26 December 2016 (Table 2), false freeze events were detected and likely attributed to a rougher ice surface created during the early stages of a mid-winter melt-freeze event. Errors on 23 January 2016 and 15 March 2017 were reported when there was no discernible difference in surface conditions between these dates and previous dates with RADARSAT-2 data. However, the error on 23 January 2016 is most likely due to a four-day gap in the imagery, which results in an immediate increase in backscatter that would be indicative of increasing ice thickness ( Table 2). In addition, three minor freeze events in 2015/2016 and four events in 2016/2017 could not be validated by using shoreline observations or climate data.

Validation of 2016 and 2017 Melt Events for MacDonald, Clear, and Johnson Lakes
Eight distinct major events for the melting period were determined in 2015/2016 for the three primary study lakes (Figure 9) in which six were confirmed by using the available meteorological data and shoreline observations ( Table 3). The first event identified was on 22 December due to the high backscatter values noted on 19 December attributed to the short-lived frost flowers mentioned previously (Figures 6 and 7). The event on 13 March 2016 corresponded to warmer temperatures of −1.4 to 9.2 • C during the day and a lack of snow cover on the ice. This resulted in a reflective surface, which lowered backscatter. The next event was detected on 31 March due to the appearance of melting ponds along the shore caused by warmer temperatures (3.0 • C on 31 March prior to image acquisition) and precipitation evident in the shoreline images. The next event was detected on 17 April. This corresponded to melt following the previously mentioned freeze events on 14 April and corresponded to high maximum temperatures of 14.3 and 17.5 • C on 15/16 April. One minor melting event was observed on 20 March 2016 for Clear Lake following warm afternoon temperatures of 2.7 • C.
Remote Sens. 2018, 10, x FOR PEER REVIEW 14 of 26 data and shoreline observations ( Table 3). The first event identified was on 22 December due to the high backscatter values noted on 19 December attributed to the short-lived frost flowers mentioned previously (Figures 6 and 7). The event on 13 March 2016 corresponded to warmer temperatures of −1.4 to 9.2 °C during the day and a lack of snow cover on the ice. This resulted in a reflective surface, which lowered backscatter. The next event was detected on 31 March due to the appearance of melting ponds along the shore caused by warmer temperatures (3.0 °C on 31 March prior to image acquisition) and precipitation evident in the shoreline images. The next event was detected on 17 April. This corresponded to melt following the previously mentioned freeze events on 14 April and corresponded to high maximum temperatures of 14.3 and 17.5 °C on 15/16 April. One minor melting event was observed on 20 March 2016 for Clear Lake following warm afternoon temperatures of 2.7 °C. Nine distinct major events were reported in 2017 ( Figure 10) in which six were confirmed by using the available meteorological data and shoreline observations (Table 3). 27 December ( Figure  10) was confirmed as the first melting event and corresponded to increases in temperature (0.1 to 3.9 °C) that occurred during the overnight period prior to image acquisition between 26/27 December. The following event on 12 January was associated with positive temperatures that resulted in a snowmelt visible in the shoreline cameras. The spring deterioration of the ice cover was first detected on 26 March 2017 (Figure 10), which was shown by the extensive pooling that is visible on the surface of the lakes. The melt event on 12 April corresponded with warmer temperatures of 2.4 °C compared to the previous image on 9 April (−2.9 °C) and shoreline cameras showed noticeable melting between images captured during the previously confirmed freeze event (9 April) and 12 April. The 19 April melting event corresponded with a large decrease in ice cover from the previous observations in Figure 9. Backscatter and temperature evolution for the 2015/2016 ice season. Black squares represent images that were acquired in the PM and grey squares representing images acquired in the AM. The orange lines indicate when a major melt event was determined. The red lines indicate WCI and the solid lines are the dates of the sample images. The two images are from the lake ice monitoring network and display the surface conditions corresponding to the respective event. The 13 March image was acquired from MacDonald Lake and the 30 April image was acquired from Johnson Lake.
Nine distinct major events were reported in 2017 ( Figure 10) in which six were confirmed by using the available meteorological data and shoreline observations (Table 3). 27 December (Figure 10) was confirmed as the first melting event and corresponded to increases in temperature (0.1 to 3.9 • C) that occurred during the overnight period prior to image acquisition between 26/27 December. The following event on 12 January was associated with positive temperatures that resulted in a snowmelt visible in the shoreline cameras. The spring deterioration of the ice cover was first detected on 26 March 2017 (Figure 10), which was shown by the extensive pooling that is visible on the surface of the lakes. The melt event on 12 April corresponded with warmer temperatures of 2.4 • C compared to the previous image on 9 April (−2.9 • C) and shoreline cameras showed noticeable melting between images captured during the previously confirmed freeze event (9 April) and 12 April. The 19 April melting event corresponded with a large decrease in ice cover from the previous observations in shoreline imagery from 12 April. Two minor melting events were captured by the method on 19 January when temperatures were around 0 • C and there was some surface melting visible on 19 February for MacDonald Lake when warmer temperatures (maximum 7.4 • C) were observed in the afternoon and evening prior to image acquisition.  The threshold-based method was over-sensitive to melting events and identified events that did not correspond to a freeze event or noticeable change in ice conditions. On 12 January 2016, an event was incorrectly reported in response to the high backscatter values seen at the beginning of a melting event on 9 January (Table 3). Errors on 12 February 2016 and 3 December 2016 were most likely due to the mix of orbits used (ascending and descending) ( Table 3). Events detected on 9 and 20 December 2016 and 1 April 2017 could not be validated by using shoreline observations or available climate data. On 30 April 2016 a decrease in the backscatter was detected for Clear Lake six days after the in situ observed WCI event (24 April  The threshold-based method was over-sensitive to melting events and identified events that did not correspond to a freeze event or noticeable change in ice conditions. On 12 January 2016, an event was incorrectly reported in response to the high backscatter values seen at the beginning of a melting event on 9 January (Table 3). Errors on 12 February 2016 and 3 December 2016 were most likely due to the mix of orbits used (ascending and descending) ( Table 3). Events detected on 9 and 20 December 2016 and 1 April 2017 could not be validated by using shoreline observations or available climate data. On 30 April 2016 a decrease in the backscatter was detected for Clear Lake six days after the in situ observed WCI event (24 April) that matched with the RADARSAT event observed on 24 April. However, only 47% of pixels reported the event on 30 April vs 82% on 24 April, which shows that the event had more extensive coverage on 24 April. Events identified in the shoreline imagery on 30 January 2016 and 6 March 2016 were missed due to gaps in the RADARSAT-2 imagery. Six minor melt events in 2015/2016 and six events in 2016/2017 could not be validated by using shoreline observations or local climate data. Table 3. Summary of the validation of major melt/WCI events for the 2015/2016 and 2016/2017 ice seasons. "Validated" events were melted or WCI events confirmed by shoreline observations and climate data. "Not Validated" events were melted events identified by the threshold method but did not match observations in shoreline cameras or local climate data and "Not Captured" events indicate that no melt or WCI event was reported for that lake on the specific date.

Event Date
Johnson Lake MacDonald Lake Clear Lake

Validation of 2016 and 2017 WCI Events for MacDonald, Clear, and Johnson Lakes
In 2016, the RADARSAT-2 detected WCI date was 24 April for both MacDonald and Clear Lakes (Figure 9), which is only one day difference from the WCI date observed in the shoreline imagery for MacDonald lake and coincides directly with WCI from the shoreline imagery for the Clear Lake. The RADARSAT-2 WCI date determined for Johnson Lake was on 30 April (Figure 9), which was the next radar image available after the observed WCI in the shoreline cameras on 27 April.
In 2017, the WCI event for MacDonald and Clear Lakes were successfully detected by RADARSAT-2 as 25 April, which was the subsequent image following the WCI event recorded in the shoreline cameras on the evening of 19 April (the previous RADARSAT-2 image was from the morning of 19 April prior to WCI).
WCI was not accurately determined for Johnson Lake in 2017 since the last event was on 19 April, but shoreline observations still showed extensive ice cover on this date and did not observe WCI until 21 April. No event was detected in the later RADARSAT-2 images acquired on 25 April or 26 April for Johnson Lake due to the similarity in backscatter between these dates and 19 April (Figure 10). This was caused by a wet or possibly smooth ice surface on 19 April for Johnson Lake, which was visible in the shoreline cameras and may have resulted in a similar backscatter value to open water ( Figure 11). Remote Sens. 2018, 10, x FOR PEER REVIEW 17 of 26 Figure 11. Comparison of ice surface conditions between the three primary study lakes on 19 April 2017.

Validation of 2009-2017 Ice Seasons using MODIS Imagery
Between RADARSAT-2 and MODIS visible imagery, the freeze onset mean bias error (MBE) ranged from 6.1 days earlier to 10 days later for the initial appearance of ice and 2.5 to 10.0 days in difference for the mean absolute error (MAE) ( Table 4). The error values correspond well with previous research conducted [16] in which MAE values of 9.8 days were found for the equivalent of freeze onset determined from AVHRR data and compared dates from the Canadian Ice Database (CID) and Global Lake and River Ice Phenology Database (GLRID). The error values observed are also similar to other threshold-based methods used in northern Alaska that showed a difference of 0-15 days between mean SAR and CLIMo observations [12]. Additionally, there are similar differences for freeze onset for Great Bear Lake in northern Canada between QuickSCAT and the Canadian Ice Database (CID) of 4 to 18 days [45]. The MODIS visible images also corresponded with the in situ ground observations of freeze onset, which range from 3 to 6 days earlier for MacDonald, Johnson, and Clear Lakes for 2016 and 2017. Additionally, visual inspection of freeze up in 2010/2011 shows that spatial patterns of lake ice in the RADARSAT-2 images are similar to that observed in the MODIS visible imagery (Figure 12). For example, both datasets show extensive ice cover across the Lake of Bays except for the center of the lake (the dark tones in both images represent open water conditions while the light grey tones in RADARSAT-2 and bright tones in MODIS visible imagery represent ice) ( Figure 12). However, RADARSAT-2 identified freeze onset dates for this year that were later than dates showing the first appearance of ice in the MODIS visible imagery (Table 4) due to the lack of consecutive images until 10 December and 18 December, which caused RADARSAT-2 to miss freeze onset.
During the 2011/2012 ice season freeze, events were identified on 16 December, which is 7 days earlier than ice was first seen in the MODIS visible imagery. The lakes were obscured by clouds in the MODIS imagery from 3 December to 17 December. However, maximum daily temperatures between 1 December and 16 December ranged from −3.5 to 12.5 °C, which helped confirm that the early detection by RADARSAT-2 on 16 December is likely false and may have been related to a rough lake surface rather than ice formation.

Validation of 2009-2017 Ice Seasons using MODIS Imagery
Between RADARSAT-2 and MODIS visible imagery, the freeze onset mean bias error (MBE) ranged from 6.1 days earlier to 10 days later for the initial appearance of ice and 2.5 to 10.0 days in difference for the mean absolute error (MAE) ( Table 4). The error values correspond well with previous research conducted [16] in which MAE values of 9.8 days were found for the equivalent of freeze onset determined from AVHRR data and compared dates from the Canadian Ice Database (CID) and Global Lake and River Ice Phenology Database (GLRID). The error values observed are also similar to other threshold-based methods used in northern Alaska that showed a difference of 0-15 days between mean SAR and CLIMo observations [12]. Additionally, there are similar differences for freeze onset for Great Bear Lake in northern Canada between QuickSCAT and the Canadian Ice Database (CID) of 4 to 18 days [45]. The MODIS visible images also corresponded with the in situ ground observations of freeze onset, which range from 3 to 6 days earlier for MacDonald, Johnson, and Clear Lakes for 2016 and 2017. Additionally, visual inspection of freeze up in 2010/2011 shows that spatial patterns of lake ice in the RADARSAT-2 images are similar to that observed in the MODIS visible imagery ( Figure 12). For example, both datasets show extensive ice cover across the Lake of Bays except for the center of the lake (the dark tones in both images represent open water conditions while the light grey tones in RADARSAT-2 and bright tones in MODIS visible imagery represent ice) ( Figure 12). However, RADARSAT-2 identified freeze onset dates for this year that were later than dates showing the first appearance of ice in the MODIS visible imagery (Table 4) due to the lack of consecutive images until 10 December and 18 December, which caused RADARSAT-2 to miss freeze onset.
During the 2011/2012 ice season freeze, events were identified on 16 December, which is 7 days earlier than ice was first seen in the MODIS visible imagery. The lakes were obscured by clouds in the MODIS imagery from 3 December to 17 December. However, maximum daily temperatures between 1 December and 16 December ranged from −3.5 to 12.5 • C, which helped confirm that the early detection by RADARSAT-2 on 16 December is likely false and may have been related to a rough lake surface rather than ice formation.  3.2 Similar to freeze, visual inspection of spatial patterns in the distribution of ice are similar between MODIS and RADARSAT-2 ( Figure 12). For WCI RADARSAT-2, dates ranged from 5.1 days earlier to 2.7 days later for MBE and 1.5 to 7.1 days in difference for MAE from 2009 to 2017 when compared to MODIS visible imagery (Table 4). This range of errors also compares well with work conducted using AVHRR, which reports the error for the equivalent of WCI at 4.2 days [16]. These results are also similar to studies conducted at higher latitudes by using threshold methods, which show differences between 0 and 9 days when comparing mean SAR detected dates and mean dates from the Canadian Lake Ice Model (CLIMo) [12]. Additionally, similar differences were noted between the QuickSCAT threshold methods and the CID 8 to14 days for Great Bear Lake in northern Canada [44]. The MODIS visible images also corresponded with the in situ ground observations of WCI, which ranged from 0 to 3 days earlier for MacDonald, Johnson, and Clear Lakes in 2016 and  Similar to freeze, visual inspection of spatial patterns in the distribution of ice are similar between MODIS and RADARSAT-2 ( Figure 12). For WCI RADARSAT-2, dates ranged from 5.1 days earlier to 2.7 days later for MBE and 1.5 to 7.1 days in difference for MAE from 2009 to 2017 when compared to MODIS visible imagery (Table 4). This range of errors also compares well with work conducted using AVHRR, which reports the error for the equivalent of WCI at 4.2 days [16]. These results are also similar to studies conducted at higher latitudes by using threshold methods, which show differences between 0 and 9 days when comparing mean SAR detected dates and mean dates from the Canadian Lake Ice Model (CLIMo) [12]. Additionally, similar differences were noted between the QuickSCAT threshold methods and the CID 8 to14 days for Great Bear Lake in northern Canada [44]. The MODIS visible images also corresponded with the in situ ground observations of WCI, which ranged from 0 to 3 days earlier for MacDonald, Johnson, and Clear Lakes in 2016 and 2017. RADARSAT-2 dates occasionally detected events that occurred after the WCI date in visible MODIS images. In 2015, the WCI date is observed as 28 April. However, RADARSAT-2 also determined two dates on 7 May and 14 May. These dates are unlikely to be representative of WCI as temperatures had been above 0 • C since early April. The threshold-based method was applied to all images provided regardless of the weather in a given season since some ice seasons are longer than others. Our analysis attempted to cover the beginning of November to the end of May depending on available imagery. For example, in 2015, the images extended to 23 May while, in 2012, the images only extended to 4 April. Therefore, these erroneous events can be removed as errors in the radar record, which resulted in the lower MAE value of 1.9 that is reported in Table 4 (Table 5). While fall/early winter 2014 was not especially cold (4.0 • C) 2014 was the second coldest year overall (4.1 • C) with the latest WCI timing (Table 5). This, in combination with a slightly cooler than average summer, would have limited the possible heat storage in the lakes during the open water season, which leads to earlier freezing in the fall. The latest date of freeze onset was observed during the 2015/2016 ice season on 3 January (Table 5), which coincides with the warmest seasonal average of 6.8 • C. There was a strong correlation (p = 0.07), r = 0.67, observed between the date of freeze onset and seasonal temperatures. As lake depth/volume is one of the main controls for the timing of freeze onset, some variability in the temperature and phenology correlation is to be expected when using the average freeze onset dates for all of the study lakes since they encompass a variety of lake morphometries. Given the short record of our estimated RADARSAT-2 freeze onset, long term trends cannot be assessed. However, from Table 5, interannual variability is considerable. A variety of trends have been reported for the region, which varies by the length of the time span examined. Reference [6] reported a lack of a significant freeze up trends for this region from 1961-1990 and Reference [14] report a very slight cooling during the fall months with delayed freeze up (non-significant) from the brief span of 2001-2014. However, this is a short snapshot superimposed on the longer term temperature trends, which show increasing fall and winter temperatures across Canada between 1950-2010 and 1900-2010 [46]. Reference [1] report a shift to later freeze for a small sample of North American lakes (8) from 1975-1976 to 2004-2005 and across the whole Northern Hemisphere report a shift of 10.7 days later per century, which is in agreement with the reported temperature increases reported by Reference [46]. Similar shifts in freeze timing are identified for lakes in the Arctic between 1950-2011 with later freeze onset occurring 5.9 days later based on modelled data [25].
The average date of WCI for the 2008-2017 ice seasons was 20 April. The earliest average WCI date was during the 2011/2012 ice season on 26 March (Table 5). This ice season also experienced warmer spring temperatures (7.1 • C) compared to other years and the warmest full-season temperature overall (7.4 • C) resulted in the shortest ice cover duration (94 days). The latest WCI date was during the 2013/2014 ice season on May 6 ( Table 5) and occurred during the year with the coldest spring and full-season temperatures, which resulted in the longest ice cover duration of 140 days. A strong significant negative correlation was found between WCI dates and spring temperatures (r = −0.91, p < 0.001). This result is in agreement with previous findings that have also shown strong correlations between ice break up dates and the mean temperature for lakes in Central Ontario and r values ranging from −0.47 to −0.74 by using a record extending over 90 years [47]. Interannual variability in RADARSAT-2 detected WCI in Central Ontario from 2009-2017 is considerable and similar to ice freeze up, which is a variety of trends that have been reported for break up in the region. Reference [6] describes a finding that, across Canada, break up dates generally trend earlier, but a mix of non-significant trends in the Central Ontario region were found during the time spans examined (1950-1980, 1961-1990, and 1971-2000). Reference [14] describes a finding that no significant changes in spring temperature or break up dates in the region from 2001-2014. Looking at the long-term record for nine lakes across south-Central Ontario, Reference [47] reports a mix of trend direction from 1905-1991 with stronger trends towards earlier break up identified in the more recent years. For the northern hemisphere as a whole, Reference [1] describes a finding that a shift towards 8.8 days per 100 years earlier for break up, which is in alignment with the long-term temperature trends that show warming temperatures across Canada during the spring [46]. In the Arctic region, however, clear trends toward earlier break up have been identified [48,49].

Method Limitations
The lake ice detection method developed in this paper is unique from those developed in previous studies [12,44] since it is able to identify multiple freezing and melting events in a single ice season, which are common occurrences in the mid-latitudes. The RADARSAT-2 method still requires the use of a > 40% filter for the percentage of pixels that detect an event in order to help remove erroneous dates. Additional manual inspection of dates was needed in order to remove errors in the 2015/2016 and 2016/2017 ice seasons for the primary study lakes. For these years during the freeze, filtering helped to remove seven events that would have been considered errors and removed 12 erroneous events for melting.
The sensitivity of the RADARSAT-2 method is affected by the variability in surface meteorological conditions because of the strong association between backscatter and surface characteristics. Previous research has shown that detection of lake ice using HH backscatter is affected by wave action/Bragg scattering (water-wavelength resonates with the SAR wavelength causing increased backscatter) on the lake surface [18,26]. This can cause artificial increases and decreases in backscatter that can be mistaken for freezing and melting events. The formation of skim ice may act to reflect emitted radiation away from the sensor, which results in lower backscatter. During the beginning of ice formation, we would expect to see an increase. This is seen at the beginning of the 2017 ice season on 10 December 2016 for MacDonald, Johnson, and Clear Lakes when backscatter reached an average minimum of −24.7 dB. Additionally, backscatter can be affected by other meteorological events such as wet snow, which results in a lower backscatter as the water content increases the radar signal absorption compared to dry snow, which the radar signal penetrates [50,51]. The increases in water content and signal absorption could result in confusion between snowmelt events and ice melt events. Additionally, in some cases, the method is not sensitive enough and misses the event. For example, on 19 February 2017, warmer temperatures were reported but only an event for MacDonald Lake was identified and only reported as occurring for 36% of pixels on Johnson and Clear Lake, which means it was not classified as a major event ( Figure 13). The lack of reporting for all lakes may be caused by local patterns in the melt and the buffer excluding shoreline pixels where the melting period may have occurred.  13). The lack of reporting for all lakes may be caused by local patterns in the melt and the buffer excluding shoreline pixels where the melting period may have occurred. Uncertainty in the exact timing of freeze and melt events is introduced by the limited number of available RADARSAT-2 acquisitions since there was an inconsistent (ideally daily) record for backscatter, as previously discussed. In 2013 and 2014, the end of the melt season was missed due to the limited RADARSAT-2 acquisitions. The use of ascending and descending images can also contribute to temporal uncertainty since they are collected at different times of the day, which could result in the detection of freeze or melt events when there is no change in ice conditions [12]. For example, in 2017, a melt event was successfully validated on 26 March. However, on 1 April, a large decrease in backscatter was detected, but there was no clear change in ice conditions between the two dates. This change in backscatter was most likely due to the different orbits between the two images, which resulted in a difference of 2.72 dB while the threshold was 1.9 dB. The use of ascending and descending orbits may also cause the detection of additional events due to overnight refreeze between images taken in the evening and those acquired the following morning. With the upcoming launch of the RADARSAT Constellation Mission (RCM), which will provide improved temporal resolution, the previously mentioned limitations can be improved in future research [52].
The spatial resolution of RADARSAT-2 ScanSAR images limit the size of the lake that can be studied and may cause detection issues on smaller lakes especially with respect to specific ice features such as cracks and melt ponds that are more difficult to detect on smaller lakes. Deformation of the ice surface such as cracks have been known to have an impact on backscatter and can cause a difference of 10 dB between surface deformation features and smooth ice [20]. Additionally, the spatial coverage of smaller lakes such as MacDonald Lake is reduced due to the shoreline 150 m buffer that was used to eliminate the impact of shoreline pixels. The method could be tested by using finer resolution images that may provide a better indication of the spatial patterns in ice formation and decay on lakes in Central Ontario. However, the 50 m pixel spacing has some improvement over the use of other remote sensing products (i.e., MODIS snow and ice product at 500 m resolution) since the smaller pixel size can reduce resolution-related errors.

Conclusions
We developed a RADARSAT-2 lake ice phenology threshold-based method that was successful at detecting multiple freezing and melting events for lakes in Central Ontario when compared to both Uncertainty in the exact timing of freeze and melt events is introduced by the limited number of available RADARSAT-2 acquisitions since there was an inconsistent (ideally daily) record for backscatter, as previously discussed. In 2013 and 2014, the end of the melt season was missed due to the limited RADARSAT-2 acquisitions. The use of ascending and descending images can also contribute to temporal uncertainty since they are collected at different times of the day, which could result in the detection of freeze or melt events when there is no change in ice conditions [12]. For example, in 2017, a melt event was successfully validated on 26 March. However, on 1 April, a large decrease in backscatter was detected, but there was no clear change in ice conditions between the two dates. This change in backscatter was most likely due to the different orbits between the two images, which resulted in a difference of 2.72 dB while the threshold was 1.9 dB. The use of ascending and descending orbits may also cause the detection of additional events due to overnight refreeze between images taken in the evening and those acquired the following morning. With the upcoming launch of the RADARSAT Constellation Mission (RCM), which will provide improved temporal resolution, the previously mentioned limitations can be improved in future research [52].
The spatial resolution of RADARSAT-2 ScanSAR images limit the size of the lake that can be studied and may cause detection issues on smaller lakes especially with respect to specific ice features such as cracks and melt ponds that are more difficult to detect on smaller lakes. Deformation of the ice surface such as cracks have been known to have an impact on backscatter and can cause a difference of 10 dB between surface deformation features and smooth ice [20]. Additionally, the spatial coverage of smaller lakes such as MacDonald Lake is reduced due to the shoreline 150 m buffer that was used to eliminate the impact of shoreline pixels. The method could be tested by using finer resolution images that may provide a better indication of the spatial patterns in ice formation and decay on lakes in Central Ontario. However, the 50 m pixel spacing has some improvement over the use of other remote sensing products (i.e., MODIS snow and ice product at 500 m resolution) since the smaller pixel size can reduce resolution-related errors.

Conclusions
We developed a RADARSAT-2 lake ice phenology threshold-based method that was successful at detecting multiple freezing and melting events for lakes in Central Ontario when compared to both shoreline imagery and visible remote sensing images. For the 2015/2016 ice season, five freeze events and six melting events for the three primary study lakes were confirmed by using a combination of shoreline imagery and meteorological data and, for 2016/2017, six freeze events and six melting events were successfully validated. The multiple freeze-thaw events included the detection of a brief freeze-up event early in the 2015/2016 ice season and a freeze event visible in shoreline cameras following an early season melt event. This is important since the occurrence of mid-winter thaw is common in mid-latitudes and not considered with previous ice detection using SAR in northern latitudes. Water clear of ice dates were detected within one to six days of the date observed in the shoreline imagery with the six-day difference due to a gap in the RADARSAT-2 imagery. Error measurements ranged from 2.5 to 10.0 days for freeze events and 1.5 to 7.1 days for melting events for 2009-2017 when compared to MODIS visible imagery. Moving forward, combining optical and RADARSAT-2 images could help improve the overall accuracy of lake ice phenology detection by providing accurate information on timing as well as providing detailed information on spatial patterns in ice phenology with high-resolution optical images (e.g., Landsat-8 Operational Land Imager (OLI)). In addition, the variability observed in RADARSAT-2 retrieved dates corresponded to variability observed in the temperature, which is important information for communities and recreational users throughout Central Ontario who rely on lake ice for activities such as snowmobiling that may be affected by variability in temperature in future years.
Our mid-latitude lake ice detection method did encounter issues due to environmental and meteorological conditions that affected the surface cover. Depending on the water content of the snow cover, backscatter was seen to increase (dry snow) and decrease (wet snow), which created artificial events in some instances. Wave action was also shown to affect backscatter and create artificial events. Additionally, the combination of ascending and descending images resulted in some artificial events due to the transition between melting in the evening and possible re-freeze overnight due to decreasing temperatures since diurnal fluctuations around the freezing points are common during the shoulder seasons in the region.
The automation process of the method could be improved through the addition of temperature criteria or manual cut off to the ice season that would force the method to not consider any dates before (freeze) or after (melt) a certain point. The temporal resolution of the RADARSAT-2 image acquisitions resulted in gaps in the image record ranging from four to seven days, which resulted in some missed freezing and/or melting events. Furthermore, it may cause increases or decreases in backscatter due to gaps in the images. This problem will be alleviated in future work, which is made possible through the use of the RADARSAT Constellation Mission launch that is planned for 2018. This will have an effective revisit time of one day provided a complete record of imagery is available [52].  Black squares represent images that were acquired in the PM and grey squares represent images acquired in the AM. The black line represents daily mean temperatures and the solid line is 0 • C. The dotted red lines represent a range of backscatter between −20 and −25 dB. Data for other ice seasons was obtained but is not shown.