Improving the Estimation of Water Level over Freshwater Ice Cover using Altimetry Satellite Active and Passive Observations

Owing to its temporal resolution of 10-day and its polar orbit allowing several crossings over large lakes, the US National Aeronautics and Space Administration (NASA) and the French Centre National d’Etudes Spatiales (CNES) missions including Topex/Poseidon, Jason-1/2/3 demonstrated strong capabilities for the continuous and long-term monitoring (starting in 1992) of large and medium-sized water bodies. However, the presence of heterogeneous targets in the altimeter footprint, such as ice cover in boreal areas, remains a major issue to obtain estimates of water level over subarctic lakes of similar accuracy as over other inland water bodies using satellite altimetry (i.e., R ≥ 0.9 and RMSE ≤ 10 to 20 cm when compared to in-situ water stages). In this study, we aim to automatically identify the Jason-2 altimetry measurements corresponding to open water, ice and transition (water-ice) to improve the estimations of water level during freeze and thaw periods using only the point measurements of open water. Four Canadian lakes were selected to analyze active (waveform parameters) and passive (brightness temperature) microwave data acquired by the Jason-2 radar altimetry mission: Great Slave Lake, Lake Athabasca, Lake Winnipeg, and Lake of the Woods. To determine lake surface states, backscattering coefficient and peakiness at Ku-band derived from the radar altimeter waveform and brightness temperature at 18.7 and 37 GHz measured by the microwave radiometer contained in the geophysical data records (GDR) of Jason-2 were used in two different unsupervised classification techniques to define the thresholds of discrimination between open water and ice measurements. K-means technique provided better results than hierarchical clustering based upon silhouette criteria and the Calinski-Harabz index. Thresholds of discrimination between ice and water were validated with the Normalized Difference Snow Index (NDSI) snow cover products of the MODIS satellite. The use of open water threshold resulted in improved water level estimation compared to in situ water stages, especially in the presence of ice. For the four lakes, the Pearson coefficient (r) increased on average from about 0.8 without the use of the thresholds to more than 0.90. The unbiased RMSE were generally lower than 20 cm when the threshold of open water was used and more than 22 cm over smaller lakes, without using the thresholds.


Introduction
Lakes are important components of subarctic regions. They have an enormous significance to physical, ecological and socio-economic systems. In many places, these freshwaters provide diverse habitats for plants, animals and microbes. Further, they offer important potential hunting and fishing grounds, together with drinking water supplies, while constitute a key resource for certain industries such as hydroelectricity, transport and mining [1].
The hydrological regime of subarctic lakes is strongly impacted by the seasonal ice cover during several months [2]. Several factors (climatic, hydro-morphologic, and physiographic) lead to the modification of ice cover regime, including magnitude, timing, location, and duration of ice cover. These directly or indirectly control water storage, either at local or at regional scales [2,3]. For example, when combined with higher temperatures during the increasingly long open-water period, decreases in lake ice duration will lead to increased evaporation and lowering of lake levels [1]. Lake water level is a crucial parameter for flood management and dam operations [4].
Monitoring of the dynamics of freshwater ice cover systems is conducted through hydrological networks. Given its temporal coverage, repeat period, spatial coverage and frequency, satellite radar altimetry could be useful for filling the gaps in hydrometric data over un-gauged water bodies [5]. Yet, the presence of heterogeneous reflective surfaces within the altimetry footprint strongly affects the accuracy of altimetry measurements [6,7]. The presence of ice and mix of water and ice modifies the amplitude and the shape of the radar echo (waveform) acquired by the altimetry satellite affecting the height retrieval [8]. This leads to the following research question: During the seasonal cycle of ice, when is the altimetry satellite able to provide reliable estimates of water level over freshwater ice cover?
Time series of water levels are obtained by averaging several altimeter measurements along the satellite tracks using every cycle at the cross-section of the altimeter track and the water body (see [9][10][11] for a review). Discrimination of water and ice along the altimeter tracks is expected to improve the estimation of water levels from radar altimetry measurements. An extensive analysis of simultaneous acquisitions of radar altimetry backscattering coefficients at Ku-band and brightness temperatures from the radiometer on-board TOPEX/Poseidon, Jason-1, ENVISAT, Geosat, and AltiKa as a function of the different surface states of subarctic lakes was achieved [12][13][14][15][16]. These previous studies were able to discriminate between two surface states corresponding to open water and ice-cover and to identify important dates related to changes of the lake surface states: first appearance of ice, formation of stable ice cover, first appearance of open water, complete disappearance of ice [12][13][14][15][16].
Similarly, using the complementarity between active and passive microwave observations from the radar altimeter, operating in Low Resolution Mode (LRM), and the microwave radiometer on-board Jason-2, the goal of our study is to automatically identify the altimetry data acquired over open water to improve water level estimates during the periods characterized with the presence of ice, using only the point measurements of open water. The Jason-2 mission was chosen as it is representative of the majority of the missions which has been operating since the beginning of the high precision radar altimetry era that started with the launches of ERS-1 and Topex/Poseidon, respectively in 1991 and 1992, and as its data have been continuously acquired since 1992. For this study, four subarctic water bodies within Canada that are influenced by ice were selected: Great Slave Lake, Lake Athabasca, Lake Winnipeg, and Lake of the Woods. These water bodies were selected not only according to particular characteristics in terms of geography and climate conditions (localization, surface area, and freezing period), but also the ready availability of Jason-2 altimetry mission adequate coverage.

Study Area
The study area is composed of four bodies of freshwater that are located within Canada and which are seasonally ice-covered ( Figure 1). The ground-tracks of Jason-2 are superimposed on the four lakes, which are shown in the figure. Table 1 presents the characteristics of each selected freshwater body and the number of Jason-2 ground tracks covering them. Great Slave Lake is ranked as the largest lake that is entirely located inside Canada, with an area over 27 200 km 2 . It is covered by five Jason-2 altimetry tracks (1 ascending and 4 descending). It has a long duration of freeze cover, from November to June. Lake Athabasca (7850 km 2 ) and Lake Winnipeg (24 514 km 2 ) are covered respectively by three (one ascending and two descending) and two ascending Jason-2 altimetry tracks. They are characterized by a shorter frozen period, generally from mid-November to the end of May. Lake of the Woods has the smallest area (4348 km 2 ) with a short duration of ice cover, from December to May, and is covered by one ascending and one descending Jason-2 altimetry tracks.  The Ocean Surface Topography Mission (OSTM)/Jason-2 radar altimetry mission was launched in June 2008, and is the follow-up to a series of altimetry satellites (Topex/Poseidon, 1992-2002Jason-1, 2002Jason-1, -2008. The satellite orbits at a 1336 km altitude on a quasi-polar orbit covering Earth's surface between ± 66 • of latitude during its 10-day orbital cycle. Its equatorial cross-track is 315 km-wide. The Jason-2 payload includes two main nadir-looking instruments: Poseidon-3 altimeter, operating at two frequency bands (13.57 GHz, Ku band; 5.3 GHz, C band); and Advanced Microwave Radiometer (AMR), operating at three frequency channels (18.7, 23.8 and 37 GHz) [17].
Jason-2 data are made available through the Geophysical Data Records (GDR), which contain the time and location of the measurements, together with all of the parameters that are necessary to estimate the altimeter height (orbit, range, and corrections). A detailed description can be found in [5] for altimeter height estimates over land surfaces, the parameters related to the radar echoes, and the brightness temperatures, among other data.
In this study, the Jason-2 GDR data are obtained from each cycle of each ground track passing over lake. Every one track contained~300 repeat cycles, corresponding to the nominal ground track, which were observed from July 2008 to August 2016 [18]. GDR data at high frequency (20 Hz) that were used were made available by the French Centre de Topographie des Océans et de l'Hydrosphère (CTOH -http://ctoh.legos.obs-mip.fr). In addition to altimeter data, backscattering coefficients (σ • ) and peakiness (i.e., the ratio between the maximum value of the altimeter waveform and its average) in the Ku-band were used, together with brightness temperatures at two frequencies (18.7 and 37 GHz).

In situ Water Levels
Hydrometric data in Canada are collected by the Water Survey of Canada, and compiled under the HYDAT database (National Hydrological Data Archive). HYDAT provides daily and monthly mean of water levels collected at hydrometric stations for all seasons. Water level data that were related to the four water bodies under consideration were extracted from the archive. It should be noted that for ice conditions, erroneous or altered measurements are indicated in the gauge records [19]. Table 2 shows the 7 selected hydrometric stations over the sampled water bodies. The time series of in situ water level covers the acquisition periods of the Jason-2 altimeter in this study, from 8 September 2008 to 21 September 2016. In situ hydrometric station data are referenced to the Geodetic Survey of Canada Datum and to assumed datum or arbitrary datum, while Jason-2 altimetry data are referenced to the Topex/Poseidon ellipsoid. For purposes of comparison, the Jason-2 altimetry data were converted to the hydrometric station datum, by defining the systematic bias between datum. The bias is obtained from the mean difference between the water level time series derived from hydrometric data and from the Jason-2 altimetry data.

NDSI from MODIS
The Moderate Resolution Imaging Spectroradiometer (MODIS) is a sensor that was included in the payload of the Terra (launched in 1999) and Aqua (launched in 2002) satellites. MODIS is a multispectral instrument measuring radiances in 36 bands. Its spatial resolution varies from 250 m to 1-km depending upon the spectral band. MODIS has permitted the development of a wide range of standard products for applications that are related to land, ocean, and atmosphere. Among these, a snow cover product was defined, based upon the normalized difference snow index (NDSI), using MODIS reflectance acquired in the red (0.545-0.565 µm=band 4) and short wave infra-red (1.628-1.652 µm = band 6), as follows [20,21]: where ρ red and ρ short wave infra-red are spectral reflectances in the red and shortwave infra-red bands. Daily NDSI snow cover data consist of 1200 km by 1200 km tiles, at a spatial resolution of 500 m [19]. The NDSI snow cover product can take the following values: 0 (land), 1-100 (NDSI snow cover in %), 200 (missing data), 201 (no decision), 211 (night), 237 (inland water), 239 (ocean), 250 (cloud), 254 (detector saturated), and 255 (fill) [22]. NDSI can also be used to discriminate between ice and open water [23,24].

Method
The following flowchart ( Figure 2) presents the different steps of the approach that was developed to discriminate between open water, ice and the transitions, and to improve the estimates of water level during freeze and thaw periods over the four lakes that were considered in this study.

Preliminary Processing Using MAPS
The Multi-mission Altimetry Processing Software (MAPS) [25,26] was used to pre-process the altimetry data. Along-track altimeter heights (h, km) were computed from 20 Hz Jason-2 data that were contained in the GDR, as follows (equation 2) [27]: where H is the altimeter orbit, R is the range, (distance between the satellite and the surface, km), and ∆R represents the corrections to the range, which must be applied to obtain accurate estimates of height. Range values that are derived from the Offset Center Of Gravity (OCOG) retracker [28] are used to estimate altimeter height. Ranges derived from this retracker, rather than from Ocean retracker, also present in the GDR, were chosen as it was found to be robust and to provide accurate values of water levels over inland water bodies [9,11,29,30]. As presented in equation (3), conventional corrections use Wet tropospheric corrections Corr wtropo (atmospheric refraction range delay due to water vapour and cloud liquid water content of the troposphere), Dry tropospheric correction Corr dtropo (atmospheric refraction range delay due to the dry gas component of the troposphere), Ionospheric correction Corr iono (atmospheric refraction range delay due to the free electron content associated with dielectric properties of the ionosphere), and the corrections that respectively account for crustal vertical motions that are due to solid Earth Tide earth and polar Tide pole tides [31].
Lake surface level, defined by surface water height (SWH), are obtained by removing a geoid model (EGM 2008 in this study [32]) from the ellipsoidal heights [33]: The selection of valid altimetry measurements corresponding to the lake surfaces then is performed manually for each altimetry track, using the MAPS graphical user interface (see [26] for more details).

Identification of Thresholds by Clustering Methods
A clustering method is applied to classify the measurements acquired along the altimetry tracks for each cycle (or 10-day period) into different subgroups, or clusters, during the ice covered period. The objective of clustering algorithms is to partition a set of unlabelled data into distinct homogenous groups. Unlabelled means that the algorithm is not previously trained to recognize classes (i.e., unsupervised learning). K-means [34] and hierarchical clustering [35] are commonly used for unsupervised classification. These algorithms were tested to obtain the clusters of the most distinctive classes along the tracks per cycle. In this study, the Euclidean distance estimates (L1 distance) was used in both clustering algorithms to calculate straight lines between two points. Three different variables were used in the clustering process in each case (backscattering coefficient; peakiness; brightness temperature) to discriminate water and ice classes. Previous works [12][13][14][15][16] have suggested using the average value of the brightness temperatures (TB) at two different frequencies (18.7 GHz and 37 GHz). The average value of the brightness temperatures at two different frequencies is defined as follows: Over each lake, the classification approaches exploited the averaged measurements of backscatter (σ • ), peakiness and average brightness temperature along ascending and descending tracks. For the backscatter, the average value along the track is calculated in dB (decibels) using the following equation [25]: σ 0 (dB) = 10log 10 10 σ 0 10 (6) For peakiness, the formulation can be found in Appendix A. The main objective of clustering methods is to allocate more homogenous sets of data to different classes. The cluster centre of each group of data that are classified and its associated standard deviation define the thresholds of variation for each class. For K-means clustering, the cluster centre is defined by the centroid; for hierarchical clustering, the cluster centre is defined by the barycentre of all points in each cluster.

•
Optimal clustering determination To determine the better technique of clustering between K-means and hierarchical clustering, two objective criteria were selected, i.e., Silhouette [36] and Calinski-Harabasz Index [37]. These criteria measure the tightness of sample grouping within each cluster and the separation between clusters. The Silhouette method determines the quality measure of each object within the cluster by calculating the similarity (cohesion) of an object in its cluster, and its separation from other clusters (separation). The silhouette coefficient for object i is defined as: where a(i) is the average distance between object i and the other objects within cluster A; b(i) is the average distance between object i and objects that are in the "second closest" cluster B; max is the maximum operator. Silhouette values range from -1 to +1. A value of s(i) that is close to +1 means a good assignment of object i to cluster A. If s(i) is close to 0, the object is likely between clusters A and B. If s(i) is close to -1, the assignment of object i to cluster A reveals that it might be closer to members of cluster B [36]. The Calinski-Harabasz Index (or Variance Ratio Criterion) calculates the relationship between the distance from a point within a cluster to its centroid (cohesion), and the distance from the group centroid to the global centroid (separation). Calinski-Harabasz (CH) is expressed as: where SS B is the sums-of-squared distances between the centre of each cluster and the centroid of the data set, weighted by the size of the cluster; SS W is the sums-of-squared distances between the centre of each cluster and every point within the cluster; K is the number of clusters; and N is the number of observations. To obtain well-separated and compact clusters, SS B is maximized and SS W minimized. Therefore, the maximum value for CH indicates a suitable partition for the data set [38]. •

Validation of thresholds using NDSI from MODIS
To validate the classification results, MODIS-based NDSI data acquired at dates that were representative of surface status of the lakes (ice melting, ice cover, ice break, and open water) were used. The results of the two clustering approaches that were applied to the geophysical parameters (backscatter, peakiness, and average brightness temperature) were compared to snow cover extent that were obtained using the MODIS NDSI product.

Estimation of Water Level over Lake
Lake water levels or Surface Water Heights (SWH) are determined by computing the median (and its associated deviation) of the altimetry measurements for each cycle and each ground-track (see Section 2.3.1). Two different time series are produced. The first is obtained by using all measurements over a given lake based upon MAPS selection. The second is based only on the measurements corresponding to open water using the thresholds that were obtained from the cluster analysis.

Evaluation of Altimetry-based Water Levels from in-situ Water Levels
Altimetry-based water levels of the lakes are obtained using either all altimetry measurements or only those corresponding to open water, based upon thresholds that were obtained as results of the clustering process. Both of these are compared against in-situ water level obtained from hydrometric stations of the different lakes. The Pearson correlations (r), unbiased root mean-square error (unRMSE) and bias are used to evaluate the radar altimetry estimated water levels (or Surface Water Heights) in comparison with real in-situ water levels. Comparisons between in-situ and altimetry-based water levels are performed if the data were acquired the same day.

Results
In the following sections, we will analyze variation in backscatter, peakiness, and average brightness temperature over the lakes under different surface conditions, i.e., in presence of freshwater, ice, or a mixture of both. Results of the clustering process, and the estimation and validation of water levels are then presented.

Temporal Evolution of Jason-2 Parameters
For illustrative purposes, the results of temporal variation in the three parameters are shown for the Great Slave Lake, for the Jason-2 track number 045. The period shown extends from 6 January 2014 (cycle 203) to 08 January 2015 (cycle 240). Temporal variation in backscatter, peakiness and average brightness temperature are closely related to the annual cycle of the lake surface states: presence of ice during winter, break up, open water, and freeze up. This agrees with previous work based on backscatter and average brightness temperature [12][13][14][15][16]. To better illustrate and understand the behaviour of the different parameters, their monthly averages and standard deviations (std) are computed from the mean value of parameters along ascending and descending tracks covering each water body. Results of the annual time series of these parameters are shown in Figure 6.   Behaviour similar to that depicted in Figures 3-5 can be observed for each lake. These temporal variations are related to changes of state for the lake surfaces. The presence of ice corresponds to decreased values of backscatter, low values of peakiness, and high average brightness temperature. This is observed from January to the end of April for Great Slave Lake, from January to the beginning of April for Lake Athabasca, from January to mid-March for Lake Winnipeg, and from January to the beginning of March for Lake of the Woods. As expected, the decrease in the icy period moves southward. In this period, the values of peakiness are particularly stable, i.e, around the minimum.
During ice break up and ice melting phases, values of backscatter, peakiness and average brightness temperature reach their maxima, but they decrease quickly towards the minima after melting. This latter phase is observed from the beginning of May to mid-June for Great Slave Lake, from the middle of April to the middle June for Lake Athabasca, from the end of March to the beginning of June for Lake Winnipeg, and from mid-March to the end of May for Lake of the Woods.
In the presence of open water, the backscatter, peakiness and average brightness temperature are relatively stable around their minima until the date of ice formation. The open water phase lasts from the end of June to the end of October for Great Slave Lake, from mid-June to the end of October for Lake Athabasca, from the beginning of June to mid-November for Lake Winnipeg, and from the end of May to mid-November for Lake of the Woods.
The ice formation phase exhibits a behaviour that is opposite to the ice break-up phase for the three parameters that were considered in this study. Freeze-up of ice is observed from the beginning of November to the end of December for Great Slave Lake and Lake Athabasca, and from mid-November to the beginning of December for Lake Winnipeg and Lake of the Woods.

Along-Track Clustering of Surface States of the Lakes
The results that are presented above suggest that it could be possible to determine the dates of ice formation and ice break-up automatically and, consequently, the presence of open water and ice over the lake surfaces using Jason-2 data. To do so, K-means and hierarchical clustering algorithms (see Section 2.3.2) were applied using the backscatter, peakiness and average brightness temperature parameters. Table 3 summarizes the results of the two clustering methods that were used to classify lake surface conditions over the four Canadian lakes that were considered. Following the classification that was proposed in [9][10][11][12][13] for subarctic lakes, we considered a total of four clusters to discriminate, which would classify the chosen lake surfaces in this study (ice cover, ice break up, ice formation, and open water). Table 3. Thresholds defined for each class depending upon the clustering technique. The centre of the class and associated deviation are indicated. The clusters 1, 2, 3, or 4 correspond to one of these classes: open water, ice cover, ice break-up, or ice freeze-up. Based upon the cluster centres of each class, given by each technique, thresholds were defined for the different surface states (see Section 2.3.2).

Clustering Analysis
The objective criteria that had been chosen, i.e., Silhouette and Calinski-Harabasz, were applied to evaluate the results obtained from the K-means and hierarchical clustering algorithms. Both criteria try to measure the separation between clusters with different perspectives (see Section 2.3.2). The results that are presented in Table 4 show that the K-means algorithm performed slightly better than the hierarchical clustering technique according to both Silhouette and Calinski-Harabasz criteria. Silhouette and Calinski-Harabasz values are higher for K-means compared to hierarchical clustering. This means that separation between clusters is greater using K-means for all selected water bodies.
Based upon its performance, the clustering results that were produced by the K-means algorithm were used in the rest of the study. Hence, for each lake, monthly percentages of each cluster were produced using the measurements that had been acquired on every Jason-2 track over 9 years. The results are compiled in Figure 7. Cluster 1, representing between 90 and 100% of the measurements acquired in summer period that is not visible during winter, can be associated with open water. In contrast, Cluster 2, which had exhibited the opposite distribution, can be associated with ice cover. Clusters 3 and 4, which were only present during spring and autumn, can be associated with ice freeze-up and break-up, respectively.

Evaluation of the Thresholds using NDSI from MODIS
To verify the accuracy of the thresholds that were used to discriminate between ice and water obtained by K-means, we used in situ condition of lake surfaces that were derived from the NDSI snow cover images. The along track measurements of backscatter, peakiness and average brightness temperature parameters were compared with the same date images of NDSI snow cover. Figures 8 and 9 shows an example of the variations of ice cover that were determined using NDSI snow cover over Great Slave Lake and the corresponding along track profiles of backscatter, peakiness and average brightness temperature parameters that were obtained on the same date. For open water and ice cover, the comparison to NDSI snow cover confirmed the great capacity of the defined thresholds (Table 3) to classify the ground profiles of the three parameters provided by Jason-2. During the date of ice break-up (lake opening) on 12 June 2014 and the date of ice freeze-up (lake close-up) on 28 November 2014, the along-track variations of radar altimetry parameters over melting ice are close to thresholds of ice freeze-up (cluster 3) and ice break-up (cluster 4). On these dates, we observed a high degree of correspondence among the variations in backscatter, peakiness and average brightness temperature (delimited using the thresholds in Table 3) with the transitions (water-ice) seen along tracks. The different thresholds that were defined in the study proved to be efficient in discriminating between lake surface conditions. Therefore, they could be used to classify the measurements of the Jason-2 radar altimeter in terms of ice freeze-up, ice cover, ice melting, and open water for each lake.

Evaluation of Water Levels Estimated under Thresholds of Open Water
Following clustering, the time series of lake levels along each altimetry track were obtained using (i) all data over the selected lakes using MAPS, and (ii) constraining the former selection to the open water cluster only. The resulting time series are presented in Figure 10, along with the temporal variation in water levels that were acquired by in situ gauge stations. The time series of water levels that were obtained using only measurements of the "open water" cluster exhibit temporal variations that were close to the in situ ones. This greater agreement is especially remarkable during freeze and thaw periods in autumn and spring, when the lake surface is partly frozen. Statistical comparisons (Pearson product-moment correlations, unbiased RMSE and bias) against in situ measurements are presented in Table 5. The number of match-ups between in situ and altimetry-based water stages ranges from more than 400 to a little bit less than 900 depending on the in-situ gauge station. Greater agreement between altimetry-based water levels and in situ measurements is found when only considering altimetry data from the open water cluster for all gauge stations that were considered over the different lakes. In this case, altimetry-based water levels exhibit very good agreement with in situ gauge records, with correlations r higher than 0.9 in most cases, unbiased RMSEs generally lower than 20 cm, and bias varying from 11.8 cm to 87.6 cm. These results are clearly better than those obtained when all clusters are considered together, where half of the lakes exhibit unRMSE errors that are well above 20 cm. Biases are also reduced in all cases when using only data acquired over open water. Table 5. Comparisons between altimetry-based water levels and in situ gauge records using Pearson product-moment correlations (r), unbiased root mean-square error terms, and bias.

Discussion
This research showed how changes in lake surface states affected the measurements of both passive and active sensors onboard the Jason-2 satellite altimetry mission. This agrees with comprehensive analyses that were provided by previous studies based upon temporal variation in backscatter and average brightness temperature, respectively [12][13][14][15][16]. The present study goes further by automatically discriminating between the different states of lake surfaces. To do so, an unsupervised K-means clustering algorithm was proposed to discriminate between open water, ice freeze-up, ice and ice break-up, which are the four main cycling conditions of the lake surface in subarctic regions. More importantly, the algorithm uses not only the backscatter and the average brightness temperature, but also the information on the shape of the altimetry waveform through the peakiness parameter.
Each of the three parameters (backscatter , peakiness, and average brightness temperature) shows a potential for distinguishing between the different lake surface states, according to Figures 3-6. However, larger deviations can be observed for backscatter and average brightness temperature than for peakiness. These variations are most likely due to: i) Differences in the lake surface covered by the radar and the microwave sensor from one track to another, which affects the average backscatter and average brightness temperature; or ii) Changes in lake surface conditions (e.g., changes in roughness that are induced by winds or rain from one cycle to another [27] or due to the presence of currents procuring slope of the lake surface [39]), which modulate the backscattering value. Over the Great Slave Lake, up to 23.7 cm variability was observed in lake surface height along the altimeter tracks [39].
In contrast, the peakiness parameter that is related to the shape of the waveform is less sensitive to the aforementioned factors. It is particularly stable in open water and ice cover conditions. The goal here was to identify the altimetry measurements corresponding to open water, and to improve the quality of water level estimates in the presence of ice (during freeze-thaw periods). As shown in Figures 3-6, the radar altimetry data that were acquired during these periods are contaminated by the presence of ice in the altimeter footprint. This situation causes changes in the waveform shape over large lakes, from an ocean-like waveform close to the Brown model [40] in the summer, to a peaky waveform during the winter [41]. The changes in shape and amplitude strongly affect backscatter, peakiness and the altimetry-based water level estimate. If the Ocean retracker is likely to provide good estimates of the water level during ice-free periods because the waveform shape is close to the ocean-one, it is likely to provide results with lower accuracy for peaky waveforms (see, for instance, [11,29,30]). On the contrary, OCOG/Ice retracker was found to be very robust and is commonly used in the estimate of surface water over inland water bodies (i.e., lakes, rivers and wetlands) [11,26,29,[42][43][44][45][46][47]. Moreover, water strongly absorbs microwave radiation, thereby producing a substantial decrease in average brightness temperature under open water conditions. K-means clustering that was applied to the three parameters appears to be sufficient for monitoring changes in surface states of the subarctic lakes. Figure 7 summarizes the mean annual distribution of surface states for each lake over the nine years that were considered. It clearly showed that open water cluster dominates during summer months (July, August, September, and October) as expected, and is also detectable during the dates of break-up (May or June) and freeze-up (November or December). Figure 10 and Table 5 demonstrate how the estimation of water level is improved following the use of the proposed approach. This is ascertained by the strong correlations that were obtained with comparisons with the in situ station measurements (r varying from 0.81 to 0.97). Developing the capacity to determine accurately water levels using satellite data, especially during ice break-up in the spring period, is of great importance for the prediction and management of eventual flood events.
Despite good performance of the clustering process, Figure 7 shows some misclassifications during the summer period (July to September) for all lakes. However, the presence of the other classes in the summer is less than 15%. Misclassifications might be attributed to variability in the Ku-band backscatter along the altimeter tracks over lake surfaces, as observed from Jason-2 data in Figure 6 and in previous studies [48]. Figure 11 shows examples of variations of backscatter along the altimeter tracks for the four selected water bodies during the summer period. Over Great Slave Lake, backscatter exceeded the threshold of ice cover at 20.23 dB on 04 August 2014. On 27 August 2014, backscatter reached the threshold of ice break-up at 32.45 dB over Lake Athabasca. Over Lake Winnipeg on 18 August 2014, the value of backscatter approached the threshold of ice melting and break-up at 32.16 dB. Finally, over Lake of the Woods, backscatter exceeded the value of ice break-up at 31.45 dB on 27 August 2014. These variations are above 6 dB over distances of several tens of km. As the measurements over the lake surface are acquired in less than a minute by the altimeter, only spatial variation in lake surface roughness can cause such large differences. Changes of roughness may result from the effects of wind, rainfall on lake surfaces (presence of waves and attenuation due to rain, [27]), and lake circulation inducing slopes [39].  To improve the classification results, other approaches could be tested in future studies. The application of supervised classification approaches, using NDSI from MODIS to discriminate between altimetry measurements over ice and open water to create the training samples, is likely to increase the efficiency of the clustering technique. Another possibility is the analysis of the radar altimetry waveform to detect the signature of open water and ice in the radar altimeter footprint. In a recent study, the signature of ice was detected in the waveforms from Sentinel-3A, a recently launched altimetry mission operating in Synthetic Aperture Radar (SAR) mode, with the presence of a secondary maximum [49].
Even though Jason-2 can be considered as a good candidate for the monitoring of lakes and rivers due to its temporal repeat period of 10 days [50], it suffers from two major limitations: it does not cover the latitudes above 66 • N, and the footprint diameter in the Ku-band is about 10 km when it operates in LRM. Other radar altimetry missions covering higher latitudes may offer better capabilities of providing continuous monitoring of continental water bodies in subarctic regions. For instance, SARAL, the first altimeter that has been operating in the Ka-band since 2013, has a much smaller footprint (~4 km) and a larger bandwidth (480 MHz) [51]. All of these characteristics lead to greater sensitivity to different surface conditions and provide a clearer identification of the different lake surface states [14]. Sentinel-3A and 3B have been operating in SAR mode with greater along-track spatial resolution. They already demonstrated very strong capabilities for better discriminating the surface states of the lakes [49]. Due to the coarse spatial coverage of altimetry ground-tracks, even at high latitudes, the monitoring of the subarctic lakes water levels would benefit the combination of the data from the different altimetry missions operating simultaneously, not only from the 10 (Topex/Poseidon, Jason-1/2/3), 27 (Sentinel-3A and B) and 35 (ERS-1/2, ENVISAT, SARAL) -day repeat orbits, but also from Cryosat-2, first mission to operate in LRM, SAR, and SAR interferometry (SARin) modes, long-term (369 days) repeat orbit and from the lidar altimetry missions ICESat and ICESat-2. For these two latter missions, operating in the visible domain and not in the microwave ones, our approach would need to be adapted to the specific parameters describing the waveform. In a close future, with the launch of the wide-swath altimetry missions as Surface Water and Ocean Topography (SWOT), both the spatial and the temporal samplings of the subarctic lakes will increase.

Conclusions
This study presents an efficient approach for automatically identifying the surface states (open water, freezing, ice presence and break-up) of subarctic lakes, based upon active and passive microwave measurements from the Jason-2 radar altimetry mission. The method is based upon unsupervised clustering using the backscattering coefficient, the peakiness, and the brightness temperature. The approach was applied to Jason-2 data that were acquired during the nominal phase of the mission (June 2008 -September 2016) over four Canadian freshwater lakes (Great Slave Lake, Lake Athabasca, Lake Winnipeg, Lake of the Woods), which experience ice coverage during colder months.
Strong discrimination of the four lake surface states (open water, freezing, ice presence and break-up) was found using the Silhouette coefficient and Calinski-Harabaszi index. Better results were obtained with the K-means technique compared to hierarchical clustering. The K-means algorithm was then used to classify the lake surface conditions along the altimeter tracks during the study period. Overall, the approach is efficient for determining the dates of freezing and ice break of the subarctic lakes, despite some minor misclassifications. Comparisons of the estimated water level for the different lakes to in situ measurements show improved results when the approach is applied (increase in r up to 0.29, and a decrease in unbiased RMSE down to 28 cm). The proposed approach significantly improved water level estimates from radar altimetry, especially during the period of ice break-up in the spring. Its implementation in global databases of radar altimetry-based water levels of lakes would improve the quality of the data made available to the hydrological community and may increase the interest for the use of these datasets for flood event prediction and management even at high latitudes. Acknowledgments: The authors thank the Natural Sciences and Engineering Research Council of Canada (NSERC), and the Canadian Mitacs Globalink Program for the financial support of this project, as well as the French space agency (Centre National d'Etudes Spatiales -CNES) through the TOSCA Hydroweb and OSTST PRIAM grants. The authors also thank CNES for supporting Jason-2 radar altimetry data, as well as Environment and Climate Change Canada for supporting hydrometric HYDAT data. We thank three anonymous Reviewers for their helpful comments on our manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest.

Appendix A
The OCOG retracking algorithm approach consists of approximating the radar altimetry waveform with a rectangle which has the same center of gravity as the waveform [28]. The rectangle is characterized by the three following parameters -center of gravity ( where N is the total gate number, aln is the number of gates in the beginning and the ending of the waveform which are not considered, and y(n) is the power of the n th gate. The leading edge point (LEP) is given by: The backscattering coefficients (σ • ) and pulse peakiness (PP - [52]) are derived from the OCOG parameters as follows: σ • = 10log 10 (A) + K cal (A5) where K cal is the power scaling factor obtained taking into account the Automatic Gain Calibration (AGC), the antenna gain and a correction for Earth curvature [53]. The peakiness (PP), which provides information on the waveform shape, is computed as the ratio of the maximum power (highest bin value) to the accumulated echo power above the retracking point [54]. It is computed as follows: where NWF is the total number of the waveform gates, and N right is a number of points to the right from the tracking point (i.e., above the retracking point).