Assessing the Performance of Methods for Monitoring Ice Phenology of the World’s Largest High Arctic Lake Using High-Density Time Series Analysis of Sentinel-1 Data

: Lake ice is a dominant component of Canada’s landscape and can act as an indicator for how freshwater aquatic ecosystems are changing with warming climates. While lake ice monitoring through government networks has decreased in the last three decades, the increased availability of remote sensing images can help to provide consistent spatial and temporal coverage for areas with annual ice cover. Synthetic aperture radar (SAR) data are commonly used for lake ice monitoring, due to the acquisition of images in any condition (time of day or weather). Using Sentinel-1 A / B images, a high-density time series of SAR images was developed for Lake Hazen in Nunavut, Canada, from 2015–2018. These images were used to test two di ﬀ erent methods of monitoring lake ice phenology: one method using the ﬁrst di ﬀ erence between SAR images and another that applies the Otsu segmentation method. Ice phenology dates determined from the two methods were compared with visual interpretation of the Sentinel-1 images. Mean errors for the pixel comparison of the ﬁrst di ﬀ erence method ranged 3–10 days for ice-on and ice-o ﬀ , while average error values for the Otsu method ranged 2–10 days. Mean errors for comparisons of di ﬀ erent sections of the lake ranged 0–15 days for the ﬁrst di ﬀ erence method and 2–17 days for the Otsu method. This research demonstrates the value of temporally consistent image acquisition for improving the accuracy of lake ice monitoring.


Introduction
Lakes are a major component of the landscape at high latitudes. They can cover 15%-40% of the land in some regions [1]. Lakes and the ice cover they form are important in the local climate, as the presence of these features can have an impact on local energy budgets and precipitation amounts [2,3]. Recent analysis of lake ice phenology trends in the Canadian Arctic have shown that the end of break-up (BUE, when the lake is fully clear of ice) is occurring one day earlier per year (2000-2013, significant trends only) [4] and future predictions (2041-2070) of lake ice phenology estimate that the mean dates of when full ice cover occurs could shift to up to 15 days later, dates for the complete disappearance of ice up to 25 days earlier, and lead to an overall decrease in ice duration by up to 35 days [5,6]. Recently, observed changes in lake ice phenology have been related to recent, accelerated winter warming observed in the Arctic, resulting in shifts of >1 • C since 1990 [7]. These changes in lake ice dates have an impact on northern communities as the formation of lake ice establishes supply routes via ice roads, access to hunting camps and for travel between communities [8,9]. Recent research notes that a skating race in Sweden (the Vikingarännet, or The Viking Ride) had to be cancelled indefinitely in 2018 due method based on backscatter evolution and a two-class segmentation method to assess their accuracy in identifying lake ice phenology events for a large Arctic lake, Lake Hazen, using high-density time series of Sentinel-1 A/B data. The methods used in this study were developed to identify pixel-level lake ice phenology. For each pixel of the lake, ice-off is defined as the complete disappearance of ice and ice-on is defined as the appearance of consistent ice cover [14]. Lake ice phenology was also determined for different sections of the lake, these will be referred to as water clear of ice (WCI, the date when the section of the lake was free of ice) and complete freeze over (CFO, the date when the section is fully covered by ice) [14].

Study Area
This study focuses on Lake Hazen (71.05 • W, 81.78 • N), the largest freshwater lake in the High Arctic by volume, located in Quttinirpaaq National Park, on Ellesmere Island, Nunavut, Canada ( Figure 1) [33]. The nearest Environment and Climate Change Canada weather station is Alert (62.35 • W, 82.50 • N, 150 km north east of Lake Hazen), with climate normal data  showing maximum daily averages of 3.4 • C in July and minimum averages of −33.2 • C in January and February [34]. Lake Hazen experiences constant daylight from April 7 to September 4 and no daylight between October 17 and February 26 [35]. The lake has a surface area of 542 km 2 , a maximum depth of 267 m, and is located 158 m above sea level [33,36]. Lake ice begins to melt in late May or early June [37]. The lake may become ice free by late July or early August, before ice begins to form again in mid-September [18,37]. The average date for water clear of ice , determined by a combination of SAR and MODIS data [18]) is August 12 and the average date for initial ice appearance (2001-2018, determined from a combination of Sentinel-1 and MODIS data) is September 21. The lake is primarily fed by glacial runoff from the Northern Ellesmere Island Icefield, which flows between mid-June and the end of August [38]. The Lake Hazen watershed is important to study because the local climate has experienced notable increases in temperature, of 2.6 • C between 2000 and 2012, and increases in lake surface temperature of 0.1 • C per year during August, when the majority of the lake is ice free [37].

Sentinel-1 Data
Sentinel-1 C-band (5.405 GHz) SAR data were acquired (downloaded from the Alaska Satellite Facility, ASF, https://www.asf.alaska.edu/ [39]) between November 29, 2014 and February, 2019; the images are Extra Wide (EW) mode Ground Range Detected Medium Resolution (GRDM) in either HH + HV or HH polarization with a pixel spacing of 40 m. This research used only HH-polarized images. For the best coverage of Lake Hazen, a combination of different acquisition paths and orbits (both ascending and descending) were selected, which resulted in a total of 1237 satellite images, providing 81% daily temporal coverage, with an average gap of 1.2 days.
The images were processed using the Sentinel Application Platform (SNAP). Pre-processing steps included thermal noise removal and radiometric calibration [40]. The Refined Lee filter was used to remove speckle; this filter uses the orientation of a pixel, which is determined within a 7 × 7 window using 3 × 3 subareas, and local mean and variation to reduce speckle within the images [41]. In order to ensure that images were aligned, terrain correction and projection to UTM Zone 18N was performed in SNAP and the 'warp' function from the Geospatial Data Abstraction Library (GDAL) was used to increase consistency [40,42]. The Arctic digital elevation model (DEM) was selected to correct the images. This high-resolution DEM is produced using high-resolution optical satellite images and has an approximate vertical and horizonal accuracy of 4 m [43]. The Arctic DEM is available in pixel sizes ranging from 2 m to 1 km [43]. For the purposes of this study, the 10 m DEM was selected.
Images acquired in EW mode display varying incidence angles (18.9 • -47 • ). Previous research conducted on sea and lake ice has shown that the varying incidence angle can affect the backscatter observed from different features [21,44,45]. Therefore, a linear correction was developed and applied to each image. This method of incidence angle correction has been applied in previous research conducted on Arctic sea ice monitoring and lake ice monitoring in Central Ontario [19,45]. Two hundred random points were generated over the surface of Lake Hazen and used to extract both HH and incidence angle values for 100 randomly selected images (n = 20,000). To ensure that points were independent of the shoreline and unique values were extracted, a buffer of 500 m was generated, and a minimum distance of 120 m was used between points. Zeroes and values less than −30 dB were removed from the sample to limit outliers in the data. The remaining data (n = 19,638) were split into a 66% training (n = 13,157) set and 33% validation (n = 6481) set. This resulted in a slope of −0.35 dB/ • (R2 = 0.26, p < 0.01, RMSE = 3.83). Results of the validation set showed a spearman's rho of 0.55 (p < 0.01) between predicted and observed HH values. The correction is larger compared to previously reported values (−0.23-0.26 dB/ • ) [19,45], likely a result of incorporating values from open water images into the model development. As previous research on sea ice has demonstrated, there is a steeper dependency between water and incidence angle compared to ice [46,47]. Using this slope value, the backscatter values of all images were normalized to the mean incidence angle of EW images, 32.95 • .

Validation for Ice Phenology Dates
Although optical imagery is available, cloud cover and high solar zenith angles during freeze-up limit the application of these images for observing patterns in ice phenology. Both MODIS MOD09GQ images (250 m surface reflectance of red (620-670 nm) and near-infrared (841-876 nm) bands from MODIS on the Terra satellite) and Sentinel-2 images (10 m) were investigated as sources of validation. Due to the quality of observations, manual determination of dates was made instead using SAR images. Sentinel-1 images were compared to the NIR band of the MOD09GQ [48] and Sentinel-2 NIR band (0.842 nm) during break-up to ensure that there was consistent identification of ice ( Figure 2). Sentinel-1 images provide a clear visual record for the evolution of ice growth and decay on Lake Hazen, allowing for visual comparison to be made between patterns in the images and the results of the proposed methodology for pixel-level ice-off and ice-on. Results of the lake ice detection methods were validated using 100 randomly selected points ( Figure 1). Both mean bias error (MBE, Equation (1)) and mean absolute error (MAE, Equation (2)) were calculated [49]. Additionally, in order to compare median WCI and CFO dates for Lake Hazen, the lake was split into three sections ( Figure 2). These three sections are approximately 22.3 km in length and better represent the variation in lake ice phenology for the large lake, as determined from visual inspection of ice cover freeze-up and break-up patterns on Sentinel-1 images for all ice seasons analyzed. Section 1 characterizes the southern area of the lake, Section 2 the middle portion of the lake, and Section 3 the northern section of the lake. The WCI and CFO dates used for validation of these three sections were determined through inspection of the Sentinel-1 imagery.

Physical Basis of Backscatter Difference
There is a common pattern in backscatter evolution for lakes in the Arctic and Subarctic from ice-free to ice-covered seasons (e.g., [21,24,25,50]). Prior to the formation of ice, backscatter values can be high for Arctic lakes due to Bragg scattering from waves generated by wind on the lake surface [20,24]. However, the initial formation of ice is generally associated with low values in backscatter, as the smooth ice-water interface of the thin ice acts as a specular reflector [24], but this can differ depending on initial ice formation [50]. Bright areas (high backscatter) are observed in SAR images across the ice surface during formation, due to the presence of cracks and deformation features [25]. Backscatter increases throughout the ice season due to a combination of the dielectric contrast between ice and water (3.2 and 80, respectively [22,51,52]) and surface scattering caused by a rough ice-water interface [53,54]. Indeed, models and field studies have recently shown that scattering at the ice-water interface is the dominant mechanism for radar return [53,54]. When temperatures increase later in the ice season, there is an increase in water content on the ice or snow-covered ice surface which can absorb the incoming signal or cause the smooth ice surface to act as a specular reflector, resulting in a decrease in backscatter [21,24]. When the lakes are clear of ice, depending on wind speeds, backscatter returns to the higher values associated with the open water season [23,24]. These patterns are also observed in the high-density time series of backscatter for Lake Hazen (Figure 3). During the ice-free season, typically for one month from mid-August to early September (based on Sentinel-1 observations in this study), there is a pattern of noisy backscatter associated with a mix of windy and non-windy days. Average backscatter during this time is −21.72 dB with a standard deviation of 5.30 dB. During initial ice formation, backscatter begins low before increasing, this is followed by a plateau and decrease when temperature rise. Average backscatter during full ice cover is higher compared to the ice-free months, −16.32 dB, with a smaller standard deviation of 2.51 dB. This typical pattern in backscatter evolution is used by thresholding methods to identify lake ice phenology events. These methods have been used to identify dates of freeze onset, melt onset, and WCI for large lakes in northern Canada (Great Bear and Great Slave Lake) as well as lakes small, shallow, lakes on the North Slope of Alaska [23,26]. These methods used thresholds that were based on observations of sampled backscatter evolution data [26]. Recently, a similar approach was used to identify multiple freeze and melt events for small/medium lakes in Central Ontario [19]. However, there are disadvantages to these methods: (1) human intervention or iterative processes are used to determine thresholds; and (2) restrictions for ice phenology dates are not included in these methods or must be determined by the user.

Identification of Lake Ice Phenology Using First Difference
The approach used in this research to extract lake ice phenology from backscatter profiles builds upon these methods and addresses the limitations identified above. The first difference is determined by subtracting the first image in the sequence from the next subsequent image in the time series on a pixel-by-pixel basis. This method has been used in the previous literature [19]. Additionally, to limit the influence of pixels below the minimum Sentinel-1 noise equivalent sigma zero (NESZ) (−22 dB) [55], pixels below −30 dB were set to a value of 100 and the associated pixel and pixels before and after are ignored when determining the occurrence of lake ice phenology. The threshold of −30 dB eliminates the significantly lower pixels but retains new ice and calm water pixels that are closer/below the noise floor. Large decreases in backscatter in the first difference raster are associated with pixel-level ice-off and ice-on events ( Figure 4). To extract the lake ice events from these rasters, the method iterates through each point to identify the latest and largest decrease in backscatter; this value is determined as ice-off or ice-on depending on the date ( Figure 4). The method does not use an assigned threshold, meaning the method could be applied to any lake/any year of data as the threshold will not change spatially or temporally [19,23,26]. Additionally, a 5 × 5 majority filter is used to smooth the final ice phenology raster and reduce the overall noise in the results [56]. The resulting raster is then clipped using a 300 m mask to remove pixels that may have been influenced by land contamination. In order to restrict the date range for determination of ice-off and ice-on, Sentinel-1 observations of initial ice decay, water clear of ice (WCI), and the first appearance of ice were made (Table 1). These date ranges can help to avoid the false detection of ice phenology events in SAR images due to errors, such as those introduced by wind [19]. However, this does mean some human intervention is needed in order to determine these dates. For the purposes of this research, a cut-off date associated with full ice cover for Lake Hazen (October 15) was used to improve processing time. For comparison with the manually observed Sentinel-1 sectional dates (WCI and CFO), the pixel results from the first difference method were extracted from each section and used to calculate the median date. Figure 5 shows the full workflow for determining ice phenology dates using the first difference method. The processing steps are performed in a Python environment using GDAL [42].

Physical Basis of Two Class Segmentation
The first difference method relies on changes in backscatter, these changes are also associated with an identified visual pattern in SAR images [57]. During break-up, bright tones represent lake ice while darker tones in the image represent open water due to the lake surface acting as a specular reflector away from the sensor [57]. During freeze-up, an opposite pattern is observed, with darker tones representing the newly formed ice and brighter tones representing rough, wind-blown water [24,57]. Although there can be a mix of these tones depending on conditions, these are the common patterns observed, and demonstrate that the images could be segmented into water and ice categories using backscatter values. Figure 6 provides visual support for these patterns on Lake Hazen during both freeze-up and break-up.

Identification of Lake Ice Phenology Using Otsu Segmentation
Although the Otsu method has been used to identify sea ice (both in simulated and natural conditions) and river ice concentrations, studies have relied on optical imagery acquired via field cameras and UAVs [58][59][60]. The Otsu method has seen limited application for obtaining sea ice concentration from SAR data [61] and has not been used for monitoring lake ice concentrations or phenology. The algorithm functions by evaluating the grey levels within an image and selecting an appropriate threshold to segment the image into foreground and background [62]. The first step of the algorithm is to identify the probabilities that a pixel will fall into one of these two classes. The calculation of these values is given by Equation (3): Where wf /wb indicate the two classes (foreground and background), the foreground ranges from the minimum pixel value to the threshold value, and the background is all other values [62]. pi denotes the probability distribution, determined by dividing the number of pixels in each class by the total number of pixels [62]. The next step is to determine the class level means using Equation (4) [62]: The overall computation of the Otsu threshold can be accelerated by ignoring within-class variances and using the between-class variance, which is given by Equation (5) [62]: The threshold associated with maximum between class variance is used to segment the grey scale image into two classes [62]. The distribution of values in the SAR images were altered to improve Otsu segmentation. This was done in a Python environment by converting the SAR float backscatter to absolute integers. The Otsu segmentation method was applied to converted SAR values using an iterative process that determined the range of available values and applied the Otsu algorithm assuming each value represented a possible threshold. The maximum between-class variance was determined from all possible thresholds and used to segment the image into water, a value of 1, and ice, a value of 3. Figure 7 shows the success of the Otsu method for image segmentation and ice identification compared to both Sentinel-1 and Landsat-8 Level-1 OLI/TIRS imagery. Additional high backscatter in the Sentinel-1 image (and Otsu result) may be the result of clear ice not visible in the Landsat image or wind causing a disturbance in this area of the lake. However, there are some limitations associated with this method. High returns caused by Bragg scattering result in misclassification of ice and water as the rough water surface is classified as foreground/ice (possibly observed in Figure 7). The occurrence of this issue increases with the proportion of open water. Furthermore, the Otsu method has difficulty identifying the difference between ice and water during freeze-up, due to the appearance of cracks and deformations. These features on the ice surface are known to cause higher backscatter, but the impact lessens throughout the ice season [25]. This means misclassification may occur early in the FU process.
In order to address these issues, a set of criteria were developed. To determine the ice-off dates, we iterated through each pixel and looked for the longest continuous coverage by ice that was followed immediately by water. This helped to avoid any false classifications before the date of ice-off and reduced the impact of Bragg scattering when ice-off was identified. To determine the ice-on date, we iterated through each pixel to determine the first sequence of continuous coverage by ice. This was defined as two sequential days. There were cases when two sequential days of ice cover did not occur, in these cases the pixel was assigned the date of initial ice formation observed in Sentinel-1 imagery. The occurrence of no sequential days occurred for less than 5% of pixels within the 300 m buffer for all years. The same observations used for the first difference method of initial ice decay, water clear of ice, and the initial appearance of ice were also used to establish a time frame for identifying pixel-level ice-on and ice-off dates (Table 1). Additionally, a 5 × 5 majority filter is used to smooth the final ice phenology raster and reduce the overall noise in the results [56]. Similar to the first difference method, the individual pixel results from the Otsu method were extracted and used to calculate the median WCI and CFO dates for the three sections of Lake Hazen. Figure 8 shows the full workflow and related steps for determining lake ice phenology dates using the Otsu segmentation method. The processing steps were implemented in a Python environment using GDAL [42].

Pixel-by-Pixel Comparison
The results of the first difference method were used to map the distribution of ice-off dates at the pixel scale for Lake Hazen (Figure 9). The comparison of dates visually extracted from Sentinel-1 images with the dates obtained from the first difference method shows a good agreement between the two datasets ( Table 2). When comparing the results on a pixel-by-pixel basis to the Sentinel-1 images, the overall MBE range is  Table 2. Mean bias and mean absolute errors for the pixel-by-pixel extraction of lake ice phenology events (ice-off and ice-on) using the first difference method. Similar to ice-off, the results of the ice-on analysis were used to map dates for individual pixels across the lake ( Figure 10). However, the analysis of ice-on shows higher error values compared to those from ice-off (Table 2). Dates determined by the first difference method are earlier compared to the Sentinel-1 observations; MBE ranges from −1 to −10 days. Additionally, MAE values for Sentinel-1 data are higher compared to ice-off, with values ranging 3-10 days.

Sectional Date Comparison
When comparing median WCI between Sentinel-1 observations and the first difference method for sections of Lake Hazen, higher differences are observed for individual years. Mean errors for WCI are all under eight days for 2015-2018 (Table 3). The year with the lowest error is 2016, with an MAE of three days. The largest mean error out of the four years is 2018, with an MAE of seven days. This is due to the large error of 12 days for Section 1. The dates determined by the first difference method are earlier compared to the actual observed dates. Table 3. Mean bias and mean absolute errors for the determination of overall lake ice phenology event dates (water clear of ice (WCI) and complete freeze over (CFO)) using the first difference method.

Water Clear of Ice (WCI)
Complete Freeze Over (CFO) The differences observed for complete freeze over (CFO) are larger compared to those for WCI (Table 3)

Pixel-by-Pixel Comparison
The results of the Otsu method were used to map ice-off dates for Lake Hazen from 2015-2018 ( Figure 11). Ice-off results for the Otsu method are earlier on a pixel-by-pixel basis compared to Sentinel-1 observations, MBE ranges from −2 to −7 days (Table 4). MAE error values vary for all years and have the same range as the MBE results of 2-7 days. The results from ice-on can also be used to represent the spatial pattern of events for Lake Hazen ( Figure 12). Similar to the first difference method, mean error values for ice-on are higher compared to ice-off (Table 4). Estimated ice-on dates are earlier compared to observed dates from Sentinel-1 data, MBE ranges from −6 to −10 days, and MAE for ice-on ranges 6-10 days.

Sectional Date Comparison
When evaluating the WCI dates determined for different sections, similar error values are observed compared to the pixel-by-pixel evaluation ( Table 5)  Larger errors are determined for CFO dates compared to the WCI dates (Table 5)

Sources of Error for the First Difference Method
The results obtained using the first difference method are similar to those obtained from previous threshold methods. Recent work published for Central Ontario showed MAE values for WCI ranging from 1.5 to 7.1 days for a series of 15 lakes [19], whereas our method obtains MAE values of 3-7 days for sectional WCI dates and 3-5 days for the ice-off pixel comparisons. Although the individual WCI sectional error values were larger, ranging from 0 days to 12 days, the values are comparable to research conducted on Great Bear Lake and on lakes on the North Slope of Alaska. The values for Great Bear Lake showed differences of 8-14 days from 2000 to 2006 when comparing QuickSCAT-determined dates to Canadian Ice Service (CIS) data [26] and errors obtained comparing dates determined using RADARSAT-2 and ASAR to the numerical lake ice model, CLIMo, showed differences of 0-9 days for lakes in northern Alaska [23]. The dates obtained by our method are also comparable to those determined by passive microwave data for lakes across the Northern Hemisphere [63]. Results from this previous analysis showed differences between algorithm determined dates and the CIS, ranging 2-12 days earlier [63]. Although pixel-by-pixel errors for ice-on were larger than ice-off, ranging 3-10, the values remain similar to those obtained in Central Ontario, which ranged from 2.5 to 10 days when comparing RADARSAT-2 to MODIS-determined dates [19]. The error values for individual lake section CFO dates, 3-15 days, are also similar to those reported in previous research. Overall errors for lakes in Barrow Alaska showed differences of 0-15 days in comparison to CLIMo [23] and differences of 4-18 days for Great Bear Lake when dates determined from QuikSCAT were compared to dates from the CIS [26]. Similar to WCI, passive microwave data used to determine CFO for Northern Hemisphere lakes show similar error ranges to our results, with differences of 0 to 11 days earlier [63].
Pixel-scale detection of ice-off using the first difference method is limited by the occurrence of multiple ice-off events. This occurs when the largest difference corresponds with the initial disappearance of ice, yet a shifting ice cover can delay the final date of ice-off. This results in an earlier estimate compared to the visual assessment of Sentinel-1 images. During break-up in 2018, a large portion of the ice cover melted on 22 July. However, on 23 July the same area was covered by ice ( Figure 13). The differencing method identified the first ice-off on 22 July and did not identify the final date that occurred on August 4, resulting in a difference of 13 days.
In addition to shifting ice cover, the other cause of earlier estimations is roughness caused by wind on the lake surface. For example, a pixel during break up in 2017 was estimated as having ice-off on 30 June. However, the actual observed date for that pixel is 2 August. Although the ice-off cut-off for this year was 3 August, wind resulted in increased backscatter values until 5 August, which was the closest and largest decrease (−8.2 dB) for this pixel; this would have resulted in a smaller error of two days.
Although the method uses the last and largest decrease in backscatter, the presence of smooth ice early in the freeze-up process can cause large differences in backscatter. This early ice may melt or shift, causing the method to miss the formation of consistent ice cover. Additionally, the wind can cause large increases and decreases in backscatter, which the method may incorrectly identify as the formation of ice cover.

Sources of Error for the Otsu Method
The results of pixel-by-pixel comparison for the Otsu method are comparable to accuracies reported in other research for ice-off that used backscatter evolution thresholds. The Otsu method accuracies range from two to seven days, and lake ice phenology dates from recent research in Central Ontario comparing RADARSAT-2 to MODIS images reported errors of 1.5 to 7.1 days [19]. The individual lake-section differences are less than 17 days, which are also comparable to results reported in the literature that use backscatter evolution thresholds. The threshold methods used to identify ice-off for Great Bear Lake reported differences of 14 days when comparing results to dates from the Canadian Ice Service [26], and similar values of 12 days were found when comparing passive microwave algorithms to CIS dates [63]. The results for CFO determination for different lake sections, ranging 5-16 days earlier, are also similar to the values reported in the previous literature using both passive and active microwave data, with differences ranging 0-18 days across the Northern Hemisphere, northern Canada, and Alaska [23,26,63].
Similar to the issues mentioned with the first difference method, the definition of WCI has an impact on the accuracy of the Otsu method as it identifies the first pixel following the largest number of continuous ice values. This results in the Otsu method missing the actual WCI date if wind or shifting ice covers a patch of open water following its initial disappearance. Although the criteria were used to limit the influence of wind on the Otsu method and provide the most accurate date for WCI, they also create a limitation for certain pixels. If there are consecutive dates of ice coverage followed by a misclassified water pixel early in the break-up period (between the initial disappearance of ice and WCI), this can limit the ability of the Otsu method to capture the actual WCI date later in the season. For example, in 2015, one pixel used in the validation process reported a WCI date that was 31 days earlier compared to Sentinel-1 observations, 5 July and 5 August. This was due to a group of consecutive ice dates followed by a misclassified open water pixel at the very beginning of the ice season, caused by misclassification of the Otsu method. This misclassification is due to a lack of contrast between ice and water, resulting in a noisy segmentation. These larger, earlier dates were observed for a small number of pixels in all years. The Otsu method shows consistently earlier dates for the determination of ice-on, due to confusion in the classification between open water and ice. During freeze-up, darker tones in the image are treated as ice features while brighter tones are treated as water. However, these bright tones can also be present on the ice surface due to deformations and cracks within the ice cover ( Figure 14) [25]. These cracks may persist throughout ice formation, resulting in later identifications of consistent ice cover as compared to areas where the initial ice has a darker tone. Additionally, a lack of waves on the surface of the lake can cause it to appear darker, resulting in further misclassification of ice and water. This confusion between calm water and ice results in the Otsu method making early identifications of consistent ice cover. This is also affected by the Sentinel-1 cut-offs, which are 13 to 20 days earlier compared to continuous ice cover for different sections of the lake. Additionally, on windy days there can be a mix of bright and dark tones in backscatter for sections of the lake. This mix of tones may result in the early identification of ice-on for certain pixels. For example, in 2017, pixels in the southern portion of the lake were 13 to 16 days earlier compared to Sentinel-1 observations because of a mix of bright and dark tones in this area of the lake. This was possibly the result of wind patterns, as the rest of the lake showed high backscatter values consistent with a wind-roughened surface.

Comparison of the Methods
Error values from the first difference and Otsu methods are comparable to the previously published literature on lake ice phenology event extraction [19,23,26,63]. Overall the first difference method demonstrates more consistent error values compared to the Otsu method for ice-off dates. All MAE values for the difference method are under five days, while the values for the Otsu method range from two to seven days. For ice-on, the first difference method shows the lowest error value for 2015 and 2016; other years were comparable to the two methods. For the sectional accuracies, the difference method outperforms the Otsu method. The WCI error values for the first difference method range from one to −12 days, while the Otsu errors range from −2 to −17 days. For CFO, both methods show similar differences between individual sections, yet, in 2016, the first difference method outperforms the Otsu method with error values ranging from −3 to−5 days compared to −9 to −15 days. Although both methods perform well and are comparable to previous research, recent literature notes that the date of minimum ice cover on Lake Hazen shifted earlier by only nine days between 1997 and 2011 [18]. This means that these methods may still introduce some errors when looking at trends in ice phenology within the last 10-20 years. Furthermore, although image frequency was almost daily during break-up and freeze-up processes, small gaps in the image record could cause errors.
The first difference and Otsu methods are also useful in tracking the progress of ice decay and formation. The Otsu method is easily able to capture the progression of ice decay due to the identification of ice and water. Figure 15 (see also Videos in Supplementary Materials) shows this during the end of break-up in 2018, as the Otsu method captures the disappearance of ice cover. The Otsu method shows a smoother ice cover compared to the first difference method, which shows inconsistent dates for certain areas. This relates back to the sources of error previously discussed, as the gaps that the first difference method captures are most likely caused by earlier melt for the associated pixels. However, the first difference method does accurately capture the disappearance of ice between August 9 and 10, whereas the Otsu method shows inconsistent dates for this area (Figure 15). Comparisons of ice cover fractions were derived from both methods to observations from the Canadian Ice Service (CIS) over 2017 and 2018 (limited by the availability of data from CIS for these two years only) to further demonstrate the capability of the methods to track ice cover ( Figure 16). Ice cover fraction (i.e., fraction of the lake covered by ice, in tenth) is determined on a weekly basis through the visual interpretation of SAR and optical imagery compiled over that week by ice analysts at CIS. The largest discrepancy for both methods is on 4 August 2017, when both methods report an ice cover fraction of 0/10 (i.e., no ice). However, the CIS reports 6/10, which is most likely due to the values from the CIS being based on weekly images and, therefore, may not provide an exact match to the Sentinel-1 images. Overall, both methods are in good agreement with fractions reported by CIS ( Figure 16); differences of 1/10 are likely caused by errors due to wind, shifting ice cover, or differences in image acquisition, as previously mentioned.
Tracking the progression of ice formation during freeze-up is more difficult compared to ice decay due to limitations such as deformations on the ice surface. This difficulty is to be expected based on the higher error values obtained for ice-on (pixel-scale) and CFO (lake sections). Figure 17 demonstrates that the first difference method is better for identifying the growth of ice as it is able to capture the appearance of larger areas of ice cover between 19 and 21 September. Additionally, between 22 and 23 September, it captures the appearance of ice at the southern end of the lake that is not identified by the Otsu method. The difficulty of capturing ice formation is also reflected when comparing the percentage of ice-on pixels determined by the first difference and Otsu method to the ice cover fraction reported by the CIS (Figure 18). However, the discrepancies in 28 September 2018 and 8 October 2018 are possibly due to the timing of observations from the CIS, as Sentinel-1 images show that there is almost full coverage or complete coverage on both dates.

Conclusions
This research presents two methods of extracting ice-on and ice-off dates from high-density time-series SAR imagery. The first difference method uses the evolution of backscatter that has been reported for lakes in the Arctic and temperate latitudes [20,21,24,25], and addresses one limitation (i.e., spatially dependent thresholds and dates captured after the true WCI date) of previous threshold-based methods [19,23,26]. This limitation was addressed by relying on the largest and last decrease to represent ice phenology dates and using lake-wide Sentinel-1 observations as cut-off dates for determining ice-on and ice-off. Pixel-by-pixel comparisons of this method show error values to be within 3-10 days for ice-off and ice-on and errors ranging from 0 to 15 days when comparing single dates for different sections of the lake for WCI and CFO. Errors for this method are primarily caused when pixels alternated between open water and ice due to wind or refreezing during ice-off. Errors in the determination of ice-on are caused by the Sentinel-1 dates used as a starting point for freeze-up. This method could be further improved by using lake ice models to provide the date ranges as this would eliminate the need to manually determine cut-off dates over longer timescales.
Additionally, this paper demonstrates the first application of the Otsu method [62] for the identification of lake ice and open water in SAR imagery. The method was successful in identifying two classes due to the contrast in tones between open water and ice (i.e., [57]). The classifications obtained from this method were used to identify ice-off and ice-on dates. This method shows error values ranging from two to 10 for the pixel-by-pixel comparisons of ice-off and ice-on, and errors of two to 17 days for single dates for different sections of the lake. The Otsu method is limited by surface conditions during different lake ice phenology events. Misclassification early in the break-up period caused by a lack of clear contrast between ice and water surfaces results in earlier identification of ice-off. During freeze-up, cracks and deformations on the surface of the ice also result in misclassifications that led to early or late identification of consistent ice cover.
Comparison of the two methods shows that errors are more consistent for the first difference method compared to the Otsu method. However, the first difference method does encounter issues tracking ice disappearance and is noisier compared to the Otsu results. Both methods experience similar difficulty identifying CFO dates, yet the first difference method identified the appearance of larger areas of ice cover. Ice cover fraction estimates from both methods are largely comparable to those determined through the visual interpretation of imagery by ice analysts at CIS. The application of these methods for ice cover fraction monitoring using consistent satellite data could provide longer term records of ice cover, which is an essential climate variable (ECV) according to the Global Climate Observing System (GCOS) [64].
Although both methods show that more temporally consistent (near daily) records of SAR imagery result in the identification of ice phenology dates generally close to observations, further research is needed to continue improving methods for studying changes in ice phenology in recent decades. Future explorations of identifying ice phenology dates from high-density SAR time series need to progress towards the integration of multiple data sources, such as lake ice models and visible remote sensing, through methods such as data assimilation and machine learning. The importance of temporally consistent images is especially relevant for the study of lakes across the Canadian Arctic with the recent launch of the RADARSAT Constellation Mission (RCM) by the Canadian Space Agency in June 2019. The launch of these satellites will improve the overall revisit time compared to the previous RADARSAT-2 system, allowing for multiple daily acquisitions and improved monitoring of lake ice in the Canadian North [65].