<figure><img src="image_1.png" /></figure>

HAL
open science

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

Jawad Ziyad, Kalifa Goita, Ramata Magagi, Fabien Blarel, Frédéric Frappart

► To cite this version:

Jawad Ziyad, Kalifa Goita, Ramata Magagi, Fabien Blarel, Frédéric Frappart. Improving the Estimation of Water Level over Freshwater Ice Cover using Altimetry Satellite Active and Passive Observations. Remote Sensing, 2020, 10.3390/rs12060967. hal-04720566

HAL Id: hal-04720566

https://hal.science/hal-04720566v1

Submitted on 4 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Article

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

Jawad Ziyad <fnref n="1" />*, Kalifa Goïta <fnref n="1" />, Ramata Magagi <fnref n="1" />, Fabien Blarel <fnref n="2" /> and Frédéric Frappart <fnref n="2" /> ID

Received: 20 January 2020; Accepted: 12 March 2020; Published: 17 March 2020

<figure><img src="image_4.png" /></figure>

check for
updates

**Abstract:** 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.

**Keywords:** Inland water; radar altimetry; backscatter; brightness temperature; peakiness; ice cover; clustering; Jason-2

## 1. 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

<footnote n="1">¹ Centre d'applications et de Recherches en Télédétection (CARTEL), Université de Sherbrooke, Sherbrooke, QC J1K2R1, Canada; kalifa.goita@usherbrooke.ca (K.G.); ramata.magagi@usherbrooke.ca (R.M.)</footnote>

<footnote n="2">² LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS -14 avenue Edouard Belin, 831400 Toulouse, France; fabien.blarel@legos.obs-mip.fr (F.B.); frederic.frappart@legos.obs-mip.fr (F.F.)</footnote>

<footnote>* Correspondence: jawad.ziyad@usherbrooke.ca</footnote>
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–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–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–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.

## 2. Dataset Presentation and Methods

### 2.1. *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². 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²) and Lake Winnipeg (24 514 km²) 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²) with a short duration of ice cover, from December to May, and is covered by one ascending and one descending Jason-2 altimetry tracks.

<figure><img src="image_5.png" /><figcaption>**Figure 1.** Water bodies considered in the study: (a) Great Slave Lake, (b) Lake Athabasca, (c) Lake Winnipeg and (d) Lake of the Woods. Jason-2 ground-tracks are represented by red lines.</figcaption></figure>

<table><caption>**Table 1.** Parameters of selected Canadian water bodies and Jason-2 ground track numbers.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Latitude</th><th rowspan="2">Longitude</th><th rowspan="2">Surface Area (km²)</th><th rowspan="2">Frozen Period*</th><th colspan="2">Jason-2 (Track Number)</th></tr><tr><th>Ascending</th><th>Descending</th></tr></thead><tbody><tr><td>Great Slave Lake</td><td>61°40′N</td><td>114°W</td><td>27200</td><td>November to mid-June</td><td>#045</td><td>#076, #152, #178, #254</td></tr><tr><td>Lake Athabasca</td><td>59°16′N</td><td>109°27′W</td><td>7850</td><td>Late-November to mid-May</td><td>#095</td><td>#178, #254</td></tr><tr><td>Lake Winnipeg</td><td>52°7′N</td><td>97°15′W</td><td>24514</td><td>Mid-November to late-May</td><td>#195, #017</td><td>-</td></tr><tr><td>Lake of the Woods</td><td>49°14′59″N</td><td>94°45′03″W</td><td>4348</td><td>December to May</td><td>#169</td><td>#178</td></tr></tbody></table>

<footnote>* https://ccin.ca/ccw/lakeice/current/monitoring.</footnote>
## 2.2. Dataset

### 2.2.1. Altimetry Data, Jason-2 Mission

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–2002; Jason-1, 2002–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.legosobs-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).

### 2.2.2. 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.

<table><caption>**Table 2.** Characteristics of hydrometric stations for the study lakes.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Station Codes</th><th rowspan="2">Station Name</th><th colspan="2">Location</th></tr><tr><th>Latitude</th><th>Longitude</th></tr></thead><tbody><tr><td rowspan="2">Great Slave Lake</td><td>07OB002</td><td>Great Slave Lake at Hay River</td><td>60°51′35.5″N</td><td>115°44′3.1″ W</td></tr><tr><td>07SB001</td><td>Great Slave Lake at Yellowknife Bay</td><td>62°26′29.6″ N</td><td>114°20′59.3″ W</td></tr><tr><td>Lake Athabasca</td><td>07MC003</td><td>Lake Athabasca Near Crackingstone Point</td><td>59°23′4.0″ N</td><td>108°53′34.5″ W</td></tr><tr><td rowspan="3">Lake Winnipeg</td><td>05SA003</td><td>Lake Winnipeg at Victoria Beach</td><td>50°41′42.4″ N</td><td>96°33′43.8″ W</td></tr><tr><td>05SG001</td><td>Lake Winnipeg at Mission Point</td><td>53°11′28.0″ N</td><td>99°12′43.1″ W</td></tr><tr><td>05RE003</td><td>Lake Winnipeg at George Island</td><td>52°49′6.3″ N</td><td>97°37′10.4″ W</td></tr><tr><td>Lake of the Woods</td><td>05PD029</td><td>Lake of the Woods at Cyclone Island</td><td>49°19′42.2″ N</td><td>94°51′12.2″ W</td></tr></tbody></table>

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.
### 2.2.3. 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]:

$$NDSI = \frac{\rho_{red} - \rho_{short wave infra-red}}{\rho_{red} + \rho_{short wave infra-red}} \qquad (1)$$

where $\rho_{red}$ and $\rho_{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].

### 2.3. 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.

<figure><img src="image_6.png" /><figcaption>**Figure 2.** Flowchart of processes defining the threshold of detection of ice and water classes to improve the estimation of water levels.</figcaption></figure>
### 2.3.1. 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]:

$$h = H - R - \sum \Delta R \qquad (2)$$

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].

$$\sum \Delta R = Corr_{wtropo} + Corr_{dtropo} + Corr_{iono} + Tide_{earth} + Tide_{pole} \qquad (3)$$

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]:

$$SWH = h - Geoid \qquad (4)$$

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).

### 2.3.2. Classification of Water and Ice

#### 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–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:

$$\overline{TB} = \frac{TB_{18.7} + TB_{37}}{2} \qquad (5)$$

Over each lake, the classification approaches exploited the averaged measurements of backscatter ($\sigma^o$), 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]:

$$ \overline{\sigma^0}(dB) = 10\log_{10}\left(10^{\frac{\sigma^0}{10}}\right) \qquad (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:

$$ s(i) = \frac{b(i) - a(i)}{\max\{a(i), b(i)\}} \qquad (7) $$

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:

$$ CH = \frac{SS_B}{SS_W} \times \frac{(N - K)}{(K - 1)} \qquad (8) $$

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.

#### 2.3.3. 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.

### 2.3.4. 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.

## 3. 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.

### 3.1. *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). Figures 3–5 represent the temporal profiles of backscatter, peakiness and average brightness temperature respectively. From January to mid-May, values of backscatter and peakiness are almost constant, i.e., about 22.2 ± 3.4 dB and 3.1 ± 0.56, respectively. From January to the end of May, average brightness temperatures are also quite stable over the lake surface, around 231.8 ± 52 K. From the end of May to mid-June, both backscatter and peakiness exhibit large variations. The amplitude of these changes can reach 40 dB for backscatter, and 35 for peakiness. For the month of June, average brightness temperature was observed to decrease from 265.2 ± 93 K to 200.7 ± 21 K. From the end of June to the end of November, all three parameters exhibit almost constant values of 17.1 ± 6.41 dB, 4.1 ± 4.17, and 167.1 ± 21.8 K for backscatter, peakiness and average brightness temperature, respectively. In contrast to the average brightness temperature, the backscatter and peakiness observed in December show an along-track variability like in June.

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–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.
<figure><img src="image_7.png" /><figcaption>Figure 3. Profiles of backscatter along Jason-2 track 045 over Great Slave Lake from January to December 2014.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Figure 4. Profiles of peakiness along Jason-2 track 045 over Great Slave Lake from January to December 2014.</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Figure 5.** Profiles of average brightness temperature along Jason-2 track 045 over Great Slave Lake from January to December 2014.</figcaption></figure>

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.
<figure><img src="image_10.png" /><figcaption>**Figure 6.** Monthly average (in red) and std (in grey) time series of backscatter, peakiness and average brightness temperature corresponding to the period 8 September 2008 to 21 September 2016 and for all Jason-2 tracks that have a cross-section with (a) Great Slave Lake, (b) Lake Athabasca, (c) Lake Winnipeg and (d) Lake of the Woods.</figcaption></figure>

## 3.2. 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–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><caption>**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.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">Cluster</th><th colspan="3">K-Means</th><th colspan="3">Hierarchical Clustering</th></tr><tr><th>σ° (dB)</th><th>peakiness</th><th>T&#773;B (K)</th><th>σ° (dB)</th><th>peakiness</th><th>T&#773;B (K)</th></tr></thead><tbody><tr><td rowspan="4">Great Slave Lake</td><td>1</td><td>17.1 (5.8)</td><td>2.7 (1.2)</td><td>165.6 (16.2)</td><td>17.3 (6.2)</td><td>2.8 (1.6)</td><td>165.6 (10.5)</td></tr><tr><td>2</td><td>21.8 (6.1)</td><td>3.2 (1.4)</td><td>195.6 (29.4)</td><td>22.5 (6.5)</td><td>3.4 (1.6)</td><td>190.5 (29.2)</td></tr><tr><td>3</td><td>32.6 (8.4)</td><td>9.2 (8.2)</td><td>255.2 (17.5)</td><td>32 (7.7)</td><td>7.3 (7.4)</td><td>251.5 (18.5)</td></tr><tr><td>4</td><td>35.9 (11.7)</td><td>13.8 (8.3)</td><td>213.9 (11.8)</td><td>43.1 (5.3)</td><td>18.1 (6.2)</td><td>215.5 (11.1)</td></tr><tr><td rowspan="4">Lake Athabasca</td><td>1</td><td>15.2 (6.6)</td><td>2.4 (1.8)</td><td>172.8 (15.1)</td><td>15.6 (3.4)</td><td>3.7 (1.9)</td><td>172.8 (15.3)</td></tr><tr><td>2</td><td>17.4 (8.1)</td><td>4.1 (4.1)</td><td>196.5 (18.5)</td><td>17.1 (3.7)</td><td>4.2 (2.1)</td><td>195.1 (18.2)</td></tr><tr><td>3</td><td>39.6 (9.3)</td><td>10.3 (7.6)</td><td>259.3 (26.1)</td><td>37.7 (8.8)</td><td>7.7 (2.6)</td><td>248.1 (25.1)</td></tr><tr><td>4</td><td>33.7 (8.2)</td><td>8.1 (5.5)</td><td>236.3 (27.5)</td><td>27.2 (8.5)</td><td>8.3 (2.3)</td><td>205.1 (25.2)</td></tr><tr><td rowspan="4">Lake Winnipeg</td><td>1</td><td>16.1 (5.6)</td><td>2.8 (1.8)</td><td>161.3 (26.1)</td><td>15.9 (5.2)</td><td>2.8 (1.8)</td><td>162.7 (27.1)</td></tr><tr><td>2</td><td>17.6 (7.3)</td><td>3.5 (2.8)</td><td>184.5 (16.6)</td><td>18.8 (5.5)</td><td>3.8 (1.3)</td><td>185.7 (17.1)</td></tr><tr><td>3</td><td>32.4 (7.4)</td><td>16.6 (6.1)</td><td>242.9 (28.8)</td><td>32.2 (7.5)</td><td>16.5 (2.9)</td><td>242.3 (29.2)</td></tr><tr><td>4</td><td>34.2 (9.6)</td><td>12.1 (6.7)</td><td>216.5 (25.9)</td><td>34.1 (5.4)</td><td>15.8 (4.9)</td><td>216.8 (25.2)</td></tr><tr><td rowspan="4">Lake of the Woods</td><td>1</td><td>16.2 (7.7)</td><td>2.9 (1.1)</td><td>176.6 (14.2)</td><td>15.1 (2.3)</td><td>2.8 (1.4)</td><td>192.7 (16.9)</td></tr><tr><td>2</td><td>15.4 (3.4)</td><td>2.9 (2.7)</td><td>196.7 (26.6)</td><td>17.1 (5.5)</td><td>2.9 (1.1)</td><td>174.1 (25.1)</td></tr><tr><td>3</td><td>21.1 (4.1)</td><td>3.7 (6.1)</td><td>226.6 (28.1)</td><td>21.2 (5.5)</td><td>7.5 (1.8)</td><td>228.8 (11.1)</td></tr><tr><td>4</td><td>32.6 (3.9)</td><td>7.3 (6.3)</td><td>248.5 (26.8)</td><td>31.8 (5.1)</td><td>4.8 (1.6)</td><td>234.4 (28.6)</td></tr></tbody></table>

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).

### 3.3. 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.
<table><caption>**Table 4.** Results of the silhouette and Calinski-Harabasz objective criteria for every lake.</caption><thead><tr><th></th><th></th><th>K-Means</th><th>Hierarchical Clustering</th></tr></thead><tbody><tr><td rowspan="2">Great Slave Lake</td><td>silhouette</td><td>0.78</td><td>0.72</td></tr><tr><td>Calinski-Harabasz (10<sup>3</sup>)</td><td>3.63</td><td>3.32</td></tr><tr><td rowspan="2">Lake Athabasca</td><td>silhouette</td><td>0.78</td><td>0.73</td></tr><tr><td>Calinski-Harabasz (10<sup>3</sup>)</td><td>1.55</td><td>1.50</td></tr><tr><td rowspan="2">Lake Winnipeg</td><td>silhouette</td><td>0.78</td><td>0.71</td></tr><tr><td>Calinski-Harabasz (10<sup>3</sup>)</td><td>1.67</td><td>1.25</td></tr><tr><td rowspan="2">Lake of the Woods</td><td>silhouette</td><td>0.74</td><td>0.68</td></tr><tr><td>Calinski-Harabasz (10<sup>3</sup>)</td><td>1.53</td><td>1.34</td></tr></tbody></table>

<figure><img src="image_11.png" /><figcaption>**Figure 7.** Monthly mean distribution of the percentage of clusters occurring in each period for the different lakes that were considered (cluster1: open water class; cluster 2: ice cover class; cluster 3: ice freeze-up class; and cluster 4: ice break-up class).</figcaption></figure>

### 3.4. 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.

<figure><img src="image_12.png" /><figcaption>**Figure 8.** Examples of NDSI snow cover images over Great Slave Lake at different dates of the ice cycle: ice break-up (12 June 2014) and open water (11 August 2014), ground profiles of backscatter, peakiness and average brightness temperature that were observed on the same dates under the thresholds of discrimination between ice and water.</figcaption></figure>
<figure><img src="image_13.png" /><figcaption>**Figure 9.** Examples of NDSI snow cover images over Great Slave Lake at different dates of the ice freeze-up (28 November 2014), and ice cover (16 April 2015); ground profiles of backscatter, peakiness and average brightness temperature that were observed on the same dates under the thresholds of discrimination between ice and water.</figcaption></figure>

### 3.5. 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.
<figure><img src="image_14.png" /><figcaption>**Figure 10.** Temporal variation of Surface Water Height obtained from Jason-2 on selected Lakes: results without thresholding method (grey crosses) and results with thresholding method (blue circles) versus in situ water level (black lines).</figcaption></figure>

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><caption>**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.</caption><thead><tr><th>Lake</th><th>Station Identification</th><th>Number of Match-ups</th><th>Estimator</th><th>Open Water Cluster</th><th>All Clusters</th></tr></thead><tbody><tr><td rowspan="6">Great Slave Lake</td><td rowspan="3">07SB001</td><td rowspan="3">809</td><td>r</td><td>0.93</td><td>0.83</td></tr><tr><td>unRMSE (cm)</td><td>9.1</td><td>11.1</td></tr><tr><td>bias (cm)</td><td>45.8</td><td>49.1</td></tr><tr><td rowspan="3">07OB002</td><td rowspan="3">897</td><td>r</td><td>0.94</td><td>0.84</td></tr><tr><td>unRMSE (cm)</td><td>11.8</td><td>10.7</td></tr><tr><td>bias (cm)</td><td>37</td><td>42</td></tr><tr><td rowspan="6">Lake Athabasca</td><td rowspan="3">07MC003</td><td rowspan="3">559</td><td>r</td><td>0.97</td><td>0.93</td></tr><tr><td>unRMSE (cm)</td><td>11.3</td><td>18.7</td></tr><tr><td>bias (cm)</td><td>87.6</td><td>91.2</td></tr><tr><td rowspan="3">05RE003</td><td rowspan="3">595</td><td>r</td><td>0.89</td><td>0.78</td></tr><tr><td>unRMSE (cm)</td><td>16.4</td><td>22.1</td></tr><tr><td>bias (cm)</td><td>14.5</td><td>38.4</td></tr><tr><td rowspan="6">Lake Winnipeg</td><td rowspan="3">05SA003</td><td rowspan="3">535</td><td>r</td><td>0.83</td><td>0.74</td></tr><tr><td>unRMSE (cm)</td><td>20.6</td><td>24.3</td></tr><tr><td>bias (cm)</td><td>23.4</td><td>29.1</td></tr><tr><td rowspan="3">05SG001</td><td rowspan="3">545</td><td>r</td><td>0.81</td><td>0.74</td></tr><tr><td>unRMSE (cm)</td><td>20.8</td><td>24.7</td></tr><tr><td>bias (cm)</td><td>34</td><td>38</td></tr><tr><td rowspan="3">Lake of the Woods</td><td rowspan="3">05PD029</td><td rowspan="3">401</td><td>r</td><td>0.9</td><td>0.61</td></tr><tr><td>unRMSE (cm)</td><td>11.9</td><td>39.9</td></tr><tr><td>bias (cm)</td><td>11.8</td><td>12.3</td></tr></tbody></table>

## 4. 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–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–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].
<figure><img src="image_15.png" /><figcaption>Figure 11. *Cont.*</figcaption></figure>
<figure><img src="image_16.png" /><figcaption>**Figure 11.** Variations in backscatter (sigma 0) along the Jason-2 tracks over (a) Great Slave Lake during cycle 045 (4 August 2014), over (b) Lake Athabasca during cycle 095 (27 August 2014), over (c) Lake Winnipeg during cycle 195 (18 August 2014), and over (d) Lake of the Woods cycle 178 (27 August 2014).</figcaption></figure>

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.

## 5. 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.

**Author Contributions:** Conceptualization, J.Z., K.G., R.M., F.B., and F.F.; methodology, J.Z., K.G., F.B., and F.F.; software, F.B., and F.F.; validation, J.Z., K.G., F.B., and F.F.; formal analysis, J.Z., F.B., and F.F.; investigation, J.Z., F.B., and F.F.; resources, J.Z., F.B., and F.F.; data curation, J.Z., F.B., and F.F.; writing—original draft preparation,
J.Z.; writing—review and editing, J.Z., K.G., R.M., F.B., and F.F.; visualization, J.Z.; supervision, K.G., R.M., F.B., and F.F. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by the Natural Sciences and Engineering Research Council of Canada (NSERC), and the Canadian Mitacs Globalink Program.

**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 (COG), amplitude (A) and width (W):

$$COG = \frac{\sum_{n=1+\text{aln}}^{n=N-\text{aln}} ny^2(n)}{\sum_{n=1+\text{aln}}^{n=N-\text{aln}} y^2(n)} \qquad (A1)$$

$$A = \sqrt{\frac{\sum_{n=1+\text{aln}}^{n=N-\text{aln}} y^4(n)}{\sum_{n=1+\text{aln}}^{n=N-\text{aln}} y^2(n)}} \qquad (A2)$$

$$W = \frac{(\sum_{n=1+\text{aln}}^{n=N-\text{aln}} y^2(n))^2}{\sum_{n=1+\text{aln}}^{n=N-\text{aln}} y^4(n)} \qquad (A3)$$

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:

$$LEP = COG - 0.5 W \qquad (A4)$$

The backscattering coefficients (σ^°) and pulse peakiness (PP - [52]) are derived from the OCOG parameters as follows:

$$\sigma^° = 10\log_{10}(A) + K_{cal} \qquad (A5)$$

where Kcal 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:

$$PP = \frac{P_{\text{max}} \times N_{\text{right}}}{\sum_{i=1}^{N_{\text{WF}}} P_i} \qquad (A6)$$

where NWF is the total number of the waveform gates, and Nright is a number of points to the right from the tracking point (i.e., above the retracking point).

## References

1. Vincent, W.F.; Daniel, M.; Reinhart, P.; Isabelle, L.; Muir, D.C.G.; Young, K.L.; Bégin, Y. Chapter 4. Freshwater resources in a changing environment. In *Integrated Regional Impact Studies*; Allard, M., Lemay, M., Eds.; ArcticNet: Quebec City, QC, Canada, 2012; pp. 137–155.
2. Prowse, T.D.; Ommanney, C.S. *Northern hydrology: Canadian perspectives*; Institut national de recherche en hydrologie (Canada): Quebec City, QC, Canada, 1990; pp. 1–309.

3. Prowse, T.; Alfredsen, K.; Beltaos, S.; Bonsal, B.; Duguay, C.; Korhola, A.; McNamara, J.; Vincent, W.F.; Vluginsky, V.; Weyhenmeyer, G.A. Arctic Freshwater Ice and Its Climatic Role. *Ambio* **2011**, *40*, 46–52.
[CrossRef]

4. Lin, P.; Yang, Z.L.; Cai, X.; David, C.H. Development and evaluation of a physically-based lake level model for water resource management: A case study for Lake Buchanan, Texas. *J. Hydrol. Reg. Stud.* **2015**, *4*, 661–674.
[CrossRef]

5. Stammer, D.; Cazenave, A. *Satellite Altimetry over Oceans and Land Surfaces*; Stammer, D., Cazenave, A., Eds.; CRC Press: Boca Raton, FL, USA; Taylor & Francis: London, UK, 2017; p. 681.

6. Michailovsky, C.-I.; McEnnis, S.; Berry, P.; Smith, R.; Bauer-Gottwein, P. River monitoring from satellite radar altimetry in the Zambezi River basin. *Hydrol. Earth Syst. Sci.* **2012**, *16*, 2181–2192.
[CrossRef]

7. Baup, F.; Frappart, F.; Maubant, J. Combining high-resolution satellite images and altimetry to estimate the volume of small lakes. *Hydrol. Earth Syst. Sci.* **2014**, *18*, 2007–2020.
[CrossRef]

8. Tseng, K.-H.; Shum, C.K.; Yi, Y.; Fok, H.S.; Kuo, C.-Y.; Lee, H.; Cheng, X.; Wang, X. Envisat Altimetry Radar Waveform Retracking of Quasi-Specular Echoes over the Ice-Covered Qinghai Lake. *Terr. Atmos. Ocean. Sci.* **2013**, *24*, 615.
[CrossRef]

9. Crétaux, J.-F.; Nielsen, K.; Frappart, F.; Papa, F.; Calmant, S.; Benveniste, J. Hydrological applications of satellite altimetry: rivers, lakes, man-made reservoirs, inundated area. In *Satellite Altimetry Over Oceans and Land Surfaces*; Earth Observation of Global Changes; Stammer, D., Cazenave, A., Eds.; CRC Press: Boca Raton, FL, USA, 2017; pp. 459–504.

10. Birkett, C.M. The contribution of Topex/Poseidon to the global monitoring of climatically sensitive lakes. *J. Geophys. Res.* **1995**, *100*, 25179–25204.
[CrossRef]

11. Frappart, F.; Calmant, S.; Cauhope, M.; Seyler, F.; Cazenave, A. Preliminary results of ENVISAT RA-2-derived water levels validation over the Amazon basin. *Remote Sens. Environ.* **2006**, *100*, 252–264.
[CrossRef]

12. Kouraev, A.V.; Semovski, S.V.; Shimaraev, M.N.; Mognard, N.M.; Legresy, B.; Remy, F. The ice regime of Lake Baikal from historical and satellite data: Relationship to air temperature, dynamical, and other factors. *Limnol. Oceanogr.* **2007**, *52*, 1268–1286.
[CrossRef]

13. Kouraev, A.V.; Shimaraev, M.N.; Buharizin, P.I.; Naumenko, M.A.; Crétaux, J.-F.; Mognard, N.; Legrésy, B.; Rémy, F. Ice and Snow Cover of Continental Water Bodies from Simultaneous Radar Altimetry and Radiometry Observations. *Surv. Geophys.* **2008**, *29*, 271–295.
[CrossRef]

14. Kouraev, A.V.; Zakharova, E.A.; Rémy, F.; Suknev, A.Y. Study of Lake Baikal Ice Cover from Radar Altimetry and In-Situ Observations. *Mar. Geod.* **2015**, *38*, 477–486.
[CrossRef]

15. Kouraev, A.V.; Papa, F.; Mognard, N.M.; Buharizin, P.I.; Cazenave, A.; Cretaux, J.-F.; Dozortseva, J.; Remy, F. Sea ice cover in the Caspian and Aral Seas from historical and satellite data. *J. Mar. Syst.* **2004**, *47*, 89–100.
[CrossRef]

16. Kouraev, A.V.; Papa, F.; Mognard, N.M.; Buharizin, P.I.; Cazenave, A.; Cretaux, J.-F.; Dozortseva, J.; Remy, F. Synergy of active and passive satellite microwave data for the study of first-year sea ice in the Caspian and Aral seas. *IEEE Trans. Geosci. Remote Sens.* **2004**, *42*, 2170–2176.
[CrossRef]

17. Lambin, J.; Morrow, R.; Fu, L.-L.; Willis, J.K.; Bonekamp, H.; Lillibridge, J.; Perbos, J.; Zaouche, G.; Vaze, P.; Bannoura, W.; et al. The OSTM/Jason-2 Mission. *Mar. Geod.* **2010**, *33*, 4–25.
[CrossRef]

18. Dumont, J.; Rosmorduc, V.; Carrere, L.; Picot, N.; Bronner, E.; Couhert, A.; Desai, S.; Bonekamp, H.; Figa, J.; Lillibridge, J.; et al. *OSTM/Jason-2 Products Handbook*; CNES: Ramonville St Agne, France, 2015.

19. Hamilton, A.S.; Moore, R.D. Quantifying Uncertainty in Streamflow Records. *Can. Water Resour. J. Rev. Can. Ressour. Hydr.* **2012**, *37*, 3–21.
[CrossRef]

20. Hall, D.K.; Riggs, G.A.; Salomonson, V.V. Development of methods for mapping global snow cover using moderate resolution imaging spectroradiometer data. *Remote Sens. Environ.* **1995**, *54*, 127–140.
[CrossRef]

21. Riggs, G.A.; Hall, D.K.; Salomonson, V. V MODIS snow products user guide. *NASA Goddard Sp. Flight Cent Rep.* **2006**, *1–45*.

22. Riggs, G.A.; Hall, D.K.; Román, M.O. MODIS Snow Products Collection 6 User Guide. *NASA Goddard Sp. Flight Cent. Rep.* **2016**, *1–66*.

23. Pavelsky, T.M.; Smith, L.C. Spatial and temporal patterns in Arctic river ice breakup observed with MODIS and AVHRR time series. *Remote Sens. Environ.* **2004**, *93*, 328–338.
[CrossRef]
24. Ault, T.W.; Czajkowski, K.P.; Benko, T.; Coss, J.; Struble, J.; Spongberg, A.; Templin, M.; Gross, C. Validation of the MODIS snow product and cloud mask using student and NWS cooperative station observations in the Lower Great Lakes Region. *Remote Sens. Environ.* **2006**, *105*, 341–353. [CrossRef]

25. Frappart, F.; Fatras, C.; Mougin, E.; Marieu, V.; Diepkilé, A.T.; Blarel, F.; Borderies, P. Radar altimetry backscattering signatures at Ka, Ku, C, and S bands over West Africa. *Phys. Chem. Earth Parts A/B/C* **2015**, *83–84*, 96–110. [CrossRef]

26. Normandin, C.; Frappart, F.; Diepkilé, A.T.; Marieu, V.; Mougin, E.; Blarel, F.; Lubac, B.; Braquet, N.; Ba, A. Evolution of the Performances of Radar Altimetry Missions from ERS-2 to Sentinel-3A over the Inner Niger Delta. *Remote Sens.* **2018**, *10*, 833. [CrossRef]

27. Fu, L.-L.; Cazenave, A. *Satellite Altimetry and Earth Sciences: A Handbook of Techniques and Applications*; Academic Press: Cambridge, MA, USA, 2001; pp. 1–463.

28. Wingham, D.J.; Rapley, C.G.; Griffiths, H. New Techniques in Satellite Altimeter Tracking Systems. In *Proceedings of the Digest—International Geoscience and Remote Sensing Symposium (IGARSS)*, Zurich, Switzerland, 8 September 1986; pp. 1339–1344.

29. Sulistioadi, Y.B.; Tseng, K.-H.; Shum, C.K.; Hidayat, H.; Sumaryono, M.; Suhardiman, A.; Setiawan, F.; Sunarso, S. Satellite radar altimetry for monitoring small rivers and lakes in Indonesia. *Hydrol. Earth Syst. Sci.* **2015**, *19*, 341–359. [CrossRef]

30. Yi, Y.; Kouraev, A.V.; Shum, C.K.; Vuglinsky, V.S.; Crétaux, J.-F.; Calmant, S. The performance of altimeter waveform retrackers at Lake Baikal. *TAO Terr. Atmos. Ocean. Sci.* **2013**, *24*, 513. [CrossRef]

31. Chelton, D.B.; Ries, J.C.; Haines, B.J.; Fu, L.-L.; Callahan, P.S. Chapter 1 Satellite Altimetry. In *Satellite Altimetry*; International Geophysics; Academic Press: Cambridge, MA, USA, 2001; pp. 1–131.

32. Pavlis, N.K.; Holmes, S.A.; Kenyon, S.C.; Factor, J.K. The development and evaluation of the Earth Gravitational Model 2008 (EGM2008). *J. Geophys. Res. Solid Earth* **2012**, *117*, 1–38. [CrossRef]

33. Mercier, F.; Rosmorduc, V.; Carrere, L.; Thibaut, P. *Coastal and Hydrology Altimetry Product (PISTACH) Handbook*; CNES: Toulouse, France, 2010; pp. 1–64.

34. Hartigan, J.A.; Wong, M.A. *Algorithm AS 136: A K-Means Clustering Algorithm*. *Appl. Stat.* **1979**, *28*, 100. [CrossRef]

35. Rokach, L.; Maimon, O. Clustering Methods. In *Data Mining and Knowledge Discovery Handbook*; Springer: New York, NY, USA, 2005; pp. 321–352.

36. Rousseeuw, P.J. *Silhouettes: A graphical aid to the interpretation and validation of cluster analysis*. *J. Comput. Appl. Math.* **1987**, *20*, 53–65. [CrossRef]

37. Calinski, T.; Harabasz, J. A dendrite method for cluster analysis. *Commun. Stat. Theory Methods* **1974**, *3*, 1–27. [CrossRef]

38. Saitta, S.; Raphael, B.; Smith, I.F.C. A comprehensive validity index for clustering. *Intell. Data Anal.* **2008**, *12*, 529–548. [CrossRef]

39. Sarmiento, S.E.; Khan, S.D. Spatial–Temporal Variability of Great Slave Lake Levels From Satellite Altimetry. *IEEE Geosci. Remote Sens. Lett.* **2010**, *7*, 426–429. [CrossRef]

40. Brown, G. The average impulse response of a rough surface and its applications. *IEEE Trans. Antennas Propag.* **1977**, *25*, 67–74. [CrossRef]

41. Wang, H.; Chu, Y.; Huang, Z.; Hwang, C.; Chao, N. Robust, Long-term Lake Level Change from Multiple Satellite Altimeters in Tibet: Observing the Rapid Rise of Ngangzi Co over a New Wetland. *Remote Sens.* **2019**, *11*, 558. [CrossRef]

42. Santos da Silva, J.; Calmant, S.; Seyler, F.; Rotunno Filho, O.C.; Cochonneau, G.; Mansur, W.J. Water levels in the Amazon basin derived from the ERS 2 and ENVISAT radar altimetry missions. *Remote Sens. Environ.* **2010**, *114*, 2160–2181. [CrossRef]

43. Frappart, F.; Biancamaria, S.; Normandin, C.; Blarel, F.; Bourrel, L.; Aumont, M.; Azemar, P.; Vu, P.-L.; Le Toan, T.; Lubac, B.; et al. Influence of recent climatic events on the surface water storage of the Tonle Sap Lake. *Sci. Total Environ.* **2018**, *636*, 1520–1533. [CrossRef] [PubMed]

44. Pfeffer, J.; Seyler, F.; Bonnet, M.P.; Calmant, S.; Frappart, F.; Papa, F.; Paiva, R.C.D.; Satgé, F.; Silva, J.S. Da Low-water maps of the groundwater table in the central Amazon by satellite altimetry. *Geophys. Res. Lett.* **2014**, *41*, 1981–1987. [CrossRef]
45. Medina, C.E.; Gomez-Enri, J.; Alonso, J.J.; Villares, P. Water level fluctuations derived from ENVISAT Radar Altimeter (RA-2) and in-situ measurements in a subtropical waterbody: Lake Izabal (Guatemala). *Remote Sens. Environ.* **2008**, *112*, 3604–3617. [CrossRef]

46. Zakharova, E.A.; Kouraev, A.V.; Rémy, F.; Zemtsov, V.A.; Kirpotin, S.N. Seasonal variability of the Western Siberia wetlands from satellite radar altimetry. *J. Hydrol.* **2014**, *512*, 366–378. [CrossRef]

47. Kim, J.W.; Lu, Z.; Lee, H.; Shum, C.K.; Swarzenski, C.M.; Doyle, T.W.; Baek, S.H. Integrated analysis of PALSAR/Radarsat-1 InSAR and ENVISAT altimeter data for mapping of absolute water level changes in Louisiana wetlands. *Remote Sens. Environ.* **2009**, *113*, 2356–2365. [CrossRef]

48. Okeowo, M.A.; Lee, H.; Hossain, F.; Getirana, A. Automated Generation of Lakes and Reservoirs Water Elevation Changes From Satellite Radar Altimetry. *IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens.* **2017**, *10*, 3465–3481. [CrossRef]

49. Shu, S.; Liu, H.; Beck, R.A.; Frappart, F.; Korhonen, J.; Xu, M.; Yang, B.; Hinkel, K.M.; Huang, Y.; Yu, B. Analysis of Sentinel-3 SAR altimetry waveform retracking algorithms for deriving temporally consistent water levels over ice-covered lakes. *Remote Sens. Environ.* **2020**, *239*, 111643. [CrossRef]

50. Rybushkina, G.; Troitskaya, Y.; Soustova, I. Ice cover determination of the Volga and the Don River reservoirs on the base of Jason-2 sattelite observations. In Proceedings of the 2014 IEEE Geoscience and Remote Sensing Symposium, IEEE, Quebec City, QC, Canada, 13–18 July 2014; pp. 149–152.

51. Bonnefond, P.; Verron, J.; Aublanc, J.; Babu, K.; Bergé-Nguyen, M.; Cancet, M.; Chaudhary, A.; Crétaux, J.-F.; Frappart, F.; Haines, B.; et al. The benefits of the Ka-Band as evidenced from the SARAL/AltiKa altimetric mission: quality assessment and unique characteristics of Altika data. *Remote Sens.* **2018**, *10*, 83. [CrossRef]

52. Laxon, S.W.; Rapley, C.G. Radar altimeter data quality flagging. *Adv. Sp. Res.* **1987**, *7*, 315–318. [CrossRef]

53. ESA, ENVISAT RA2/MWR Product Handbook; ESA (European Space Agency): Paris, France, 2007; pp. 1–204.

54. Peacock, N.R.; Laxon, S.W. Sea surface height determination in the Arctic Ocean from ERS altimetry. *J. Geophys. Res. C Ocean.* **2004**, *109*, C07001. [CrossRef]

<figure><img src="image_17.png" /></figure>

© 2020 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (http://creativecommons.org/licenses/by/4.0/).