Improving the Estimation of Lake Ice Thickness with High Resolution Radar Altimetry Data

1. Introduction

Lake ice thickness (LIT) is a sensitive indicator of climate change recognized as a thematic variable of Lakes as an Essential Climate Variable (ECV) by the Global Climate Observing System (GCOS) of the World Meteorological Organization (WMO) [1]. LIT also plays a significant role on the economy of northern regions through its influence on transportation, travel, fishing, and recreation activities [2]. The thickness of lake ice has been decreasing with climate warming and is expected to further decrease in the foreseeable future [3]. Yet, continuous monitoring of LIT is limited as field measurements are sparse in both geographical coverage and time, and has been decreasing over the last few decades [4]. Remote sensing can help to fill this gap, ensuring global monitoring, but, to date, the number of studies exploring the potential of satellite-based LIT measurements has been limited and the GCOS requirements regarding spatial and temporal resolutions, measurement uncertainty and stability (change in bias over time) for LIT monitoring have not been fully met. According to GCOS [1], the minimum requirement (or threshold) LIT is a spatial resolution of 10 km, temporal resolution of 365 days (measurement at the time of maximum ice thickness), uncertainty of 15 cm (in units of 2 standard deviations), and stability of 10 cm. The GCOS ECVs requirements document indicates that a significant improvement (or breakthrough) would be achieved with a spatial resolution of 1 km, temporal resolution of 30 days, uncertainty of 10 cm (in units of 2 standard deviations), and an unspecified stability.

A recent review of radar remote sensing of lake ice that includes a review of algorithms and future directions for the retrieval of LIT is given in [4]. In previous studies, [5,6] demonstrated the high sensitivity of passive microwave brightness temperature satellite measurements to ice growth and the potential of AMSR-E 18.7 GHz (v-pol) channel to estimate LIT with an accuracy in the order of 20 cm on a 10-km grid (below the GCOS minimum requirement of 15 cm for uncertainty but in par with the 10 km requirement for spatial resolution). Radar altimetry currently has great potential to provide long LIT time series (20-30 years) needed for climate monitoring at the breakthrough level defined by GCOS as indicated above. [7] showed the potential of estimating LIT from the varying peak separation, linked to the ice growth, of Cryosat-2 Ku-band Synthetic Aperture Radar (SAR) waveforms, but achieving RMSE values ranging from 20 cm to 41 cm when compared to in situ measurements from Great Slave Lake and Great Bear Lake, Canada. More recently, [8] proposed to improve this approach by using ICESat-2 data to increase the accuracy of the Cryosat-2 peak selection, obtaining LIT estimates with a RMSE of 0.14 cm when compared to field measurements. It is worth noting here that assessing the accuracy of satellite-based LIT retrievals by using in situ measurements can be somewhat misleading as these are usually collected operationally near the shore of the lakes, while the satellite ground track segments used for the LIT estimates are taken over the lakes, the furthest away from the shore to avoid land contamination. Care should therefore be taken when evaluating along-track satellite altimetry retrievals of LIT against in situ measurements since the snow/ice properties and ice thickness can be quite different over distances that can be of the order of several tens of kilometers between the two sets of measurements.

Other recent studies [9,10,11,12] have exploited radar waveforms for LIT estimation from different missions (both Low Resolution Mode or LRM, and SAR mode), mainly focusing on the estimation of the bias introduced by the presence of ice in winter on lake water level (LWL) measurements, more than on obtaining an accurate LIT estimation. Additionally, the studies mentioned above, make use of empirical methods and/or retrackers that cannot correctly account for the specific waveform shape over ice covered lakes and also depend on threshold settings that are not accurate enough to follow the LIT evolution, in particular at the seasonal transitions when the LIT signature is less discernible in the radar waveforms, leading to biased and sub-optimal LIT estimation. Significant improvements in LIT estimation have recently been reported by [13] through the implementation of a physical based analytical retracker for the analysis of Ku-band LRM waveform data. The method provides a robust way to assess the accuracy of LIT retrievals that does not rely on in situ measurements. It provides estimates of LIT with an accuracy of the order of 10 cm. Also [13] demonstrated the capability of the analytical retracker to consistently capture the seasonal transitions of LIT during the freeze-up and breakup periods and ice growth over different winter seasons, making it a powerful and promising tool for monitoring the inter-annual variability and trends in LIT. The method has recently been implemented within the ESA Lakes_cci project 1 to generate long (20+ years) LIT time series from radar altimetry over different lake targets. The LIT time series are publicly available through the project portal according to the data release schedule. While [13] paved the way towards fulfilling the ambitious GCOS requirements for LIT, further improvements can still be achieved by using data at higher spatial resolution, that is SAR altimetry data.

In this paper, we present a novel LIT analytical retracker tailored for the analysis of Ku-band SAR waveforms from both Unfocused (UFSAR) [14] and Fully Focused (FFSAR) [15] modes that allows for a robust estimation of LIT and that significantly improves the accuracy of LIT retrievals with respect to previous studies, opening the possibility of precise LIT monitoring from current and future SAR altimetry missions. The paper focuses on the analysis of Sentinel-6 SAR data publicly available at nominal 20 Hz sampling. It also explores the benefit for LIT analysis of processed Sentinel-6 data at at higher posting rate (140 Hz) for both UFSAR and FFSAR modes. The SAR LIT analysis is evaluated through comparison with LIT thermodynamics model simulations [16] and LIT estimates obtained from Sentinel-6 and Jason-3 LRM data for the same lake targets and dates, during the tandem phase of the two satellites. Assessing the consistency between SAR based and LRM based LIT retrievals is important to ensure the continuity of the measurements issued from past, current and, eventually, future radar altimetry missions for climate monitoring.

The paper is organized as follows. sec:LIT signature describes the LIT signature in Ku-band SAR waveforms in the case of UFSAR at 20 Hz as well as UFSAR and FFSAR at 140 Hz, showing how the signal evolves during a typical ice season. sec:LIT rtks presents the mathematical formalism of the SAR_LITand FFSAR_LITretrackers, detailing the analytical modeling of the SAR radar waveforms with the LIT signature for UFSAR and FFSAR data, and describing the parameter estimation method. Details on the Sentinel-6 system and sensor settings used in the retrackers are given in the Appendix A. sec:dataandtargets presents the altimetry data used for the LIT analysis and the thermodynamic lake ice model simulations. It also describes the target lakes considered. sec:results describes the results and assessment of the LIT analysis. Finally, sec:conclusions summarizes the key findings of the study and provides avenues for future research.

2. Lake Ice Thickness Signature in SAR Radar Echograms and Waveforms

The presence of snow and ice modifies the propagation of the radar waves, and their penetration depends on the frequency of the incident signal and the properties of the different layers crossed. In the case of Ku-band (12-18 GHz), a specific signature is observed in radar altimetry waveforms consisting of a step-like break on the leading edge of LRM waveforms (see e.g. [13]) and a double peak in high resolution SAR waveforms in both unfocused and fully-focused modes (see e.g. [7,17] for UFSAR). The width of the step and the separation between the two peaks, respectively, is directly linked to LIT. This specific signature is related to the double backscattering of the radar waves at i) the snow-ice interface and ii) at the ice-water interface [7,18], which is particularly noticeable in the case of dry snow overlaying clear freshwater ice. In this paper, we focus on the analysis and modeling of high resolution radar waveforms data (UFSAR and FFSAR) over ice-covered lakes, with the aim of exploiting the double peak signature to obtain robust LIT estimates of higher accuracy compared to those retrieved from LRM data and in previous studies that make use of empirical methods (as e.g. [7,8,9,10]).

An illustrative example of the evolution of the double peak signature related to LIT in the case of 20 Hz Sentinel-6 UFSAR data from Great Slave Lake, Canada, is provided in Figure 1. The plots in the left column of the figure show Sentinel-6 UFSAR echograms, while those in the right column display the corresponding stacked waveforms, with, in black, the mean waveform at different periods: at the beginning of the ice season in December 2021 (top row), when the ice is well established on the lake in February 2021 (second row), when the ice is still present but the snow/ice interface is less discernible at the end of April 2021 (third row), and at the beginning of the melt period in May 2021 (bottom row). As can be seen in the figure, the LIT signature in radar waveforms evolves according to changes in the snow and ice conditions and thickness during the ice season. In December, the ice is thin and heterogeneous over the RoI, the peak separation is small with a slightly smaller amplitude of the first peak In February, when the ice is well established on the lake, the bimodal signature is clearly visible and the two peaks have a comparable amplitude. At the end of April, the double peak signature is again less perceptible, with a lower amplitude for the first peak. Finally, in mid-May, the bimodal feature is no longer present and this is likely due to the fact that surface (snow) melt has begun, which leads to a highly reflective surface where the two interfaces are no longer discernible. These features are in line with the ESA LIAM project 2 findings that support a lack of a peak at the start of the season, then a bigger step as the season progresses and snow is more prominent, then later in the season a loss of the peak. This evolution is primarily tied to the surface roughness but more exploratory work is needed to confirm this aspect of the waveform.

Figure 1. Illustration of the evolution of the bimodal lake ice thickness signature in Sentinel-6 UFSAR radargrams (left column) and the normalized waveforms (right column) at 20 Hz resolution at Great Slave Lake in December 2021 (top), February 2021 (second row), end of April 2021 (third row) and May 2021 (bottom). The black line in the plots of the right column corresponds to the mean waveform in the selected region of the lake.

Figure 1. Illustration of the evolution of the bimodal lake ice thickness signature in Sentinel-6 UFSAR radargrams (left column) and the normalized waveforms (right column) at 20 Hz resolution at Great Slave Lake in December 2021 (top), February 2021 (second row), end of April 2021 (third row) and May 2021 (bottom). The black line in the plots of the right column corresponds to the mean waveform in the selected region of the lake.

The Sentinel-6 FFSAR radargrams at 140 Hz resolution for the same lake target and dates are shown in Figure 2, revealing similar features as the UFSAR radargrams of Figure 1 but with increased resolution. The evolution of the LIT signature, being related to the properties of the lake surface, is visible in both UFSAR and FFSAR radargrams. The two peaks on the leading edge of the FFSAR waveforms appear, however, to be better resolved compared to those from UFSAR, suggesting possible improvements to be gained in the retrieval of LIT. Also, radar backscatter from off-nadir water surfaces is more clearly observed on the trailing edge of the Sentinel-6 FFSAR waveforms, better reflecting the complex nature of the surface within the radar footprint. This information can be very useful, for example, in gaining better insight into the evolution of LIT, especially in the seasonal transition, when the surface is not homogeneous.

Given the specific LIT signature in the high resolution Ku-band SAR waveforms described above, we developed a retracking approach based on the physical modeling of the signal and specific for the detection of LIT that allows to analytically model the LIT signal and its evolution in the case of UFSAR and FFSAR waveforms, as detailed in the following section.

Figure 2. Illustration of the evolution of the bimodal lake ice thickness signature in Sentinel-6 UFSAR (left column) and FFSAR radargrams (right column) at 140 Hz resolution at Great Slave Lake. From top to bottom the plots refer to December 2021, February 2021, end of April 2021 and May 2021 (same dates as in Figure 1).

Figure 2. Illustration of the evolution of the bimodal lake ice thickness signature in Sentinel-6 UFSAR (left column) and FFSAR radargrams (right column) at 140 Hz resolution at Great Slave Lake. From top to bottom the plots refer to December 2021, February 2021, end of April 2021 and May 2021 (same dates as in Figure 1).

3. Analytical Retrackers for the Estimation of LIT from High Resolution SAR Data

In this section, we provide details on the analytical retrackers developed to estimate LIT from UFSAR and FFSAR data that are referred to hereafter as the SAR_LITretracker and the FFSAR_LITretracker, respectively.

3.1. SAR_LITWaveform Model

In order to reproduce the double-peak signature related to LIT, the SAR waveform of ice covered lakes can be analytically modeled as the sum of two SAR waveforms. Making use of the formalism of [19,20], the resulting radar waveform for a single beam, ${\theta }_{\ell }$, and as a function of the range gate bin x and the parameter vector $\mathbf{p}$, takes the form3:

$$P_{SL}^{LIT}({\theta }_{\ell },x,\mathbf{p})\sim {\Gamma }_{k\ell }\left(0\right)[\frac{{\alpha }_1{f}_0\left(x_1\right)}{{\sigma }_c}+{\alpha }_2{f}_0\left(x_2\right)],$$

(1)

where ${\theta }_{\ell }$ is the look angle for the ℓ beam:

$${\theta }_{\ell }=\left(\frac{L_x}{h}\right)\ell ,$$

(2)

${\alpha }_1$ is the parameter related to the amplitude of the first peak, ${\alpha }_2$ is the amplitude of the second peak and:

$$\begin{array}{ccc}& & x_1=\frac{x-x_c}{2{\sigma }_c}\hfill \end{array}$$

(3)

$$\begin{array}{ccc}& & x_2=\frac{x-x_c-{\Delta }_{ICE}^{gates}}{2}\hfill \end{array}$$

(4)

where $x_c$ is the central gate, corresponding to the epoch of the first echo, ${\Delta }_{ICE}^{gates}$ the ice thickness parameter (in unit of range gates), that quantifies the peak separation. In the case of ice covered lakes, the term related to the significant wave height can be neglected, and the ${\sigma }_c$ term can be approximated by:

$${\sigma }_c\simeq {\sigma }_p\sqrt{1+{\left(\frac{{\theta }_{\ell }}{{\theta }_{LIM}}\right)}^2}$$

(5)

where ${\sigma }_p$ is the PTR bandwidth in range gate unit, and:

$${\theta }_{LIM}=\left(\frac{L_z}{{\alpha }_{Earth}{L}_x}\right)$$

(6)

with:

$${\alpha }_{Earth}=1+\frac{h}{R_{Earth}},$$

(7)

where h is the nominal orbit height, $R_{Earth}$ is the Earth’s radius, $L_x$ the along track resolution, $L_z$ the vertical resolution, respectively, defined as:

$$\begin{array}{ccc}& & L_x=\frac{ch{f}_p}{2v_t{f}_c{N}_b}\hfill \end{array}$$

(8)

$$\begin{array}{ccc}& & L_z=\frac{c}{4f_s}\hfill \end{array}$$

(9)

 

(10)

with $f_s$ the sampling frequency, c the speed of light, $f_p$ the Pulse Repetition Frequency, $v_t$ the nominal satellite velocity, $f_c$ the center frequency, and $N_b$ the number of pulses per burst. Following [20], the $f_0$ function can be expressed analytically as:

$$f_0\left(x_i\right)=\frac{\pi }{4}\sqrt{|x_i|}[{\mathcal{S}}_{-1/4}^{sc}\left({\xi }_i\right)+sign\left(x_i\right){\mathcal{S}}_{+1/4}^{sc}\left({\xi }_i\right)],$$

(11)

with ${\xi }_i=x_i^2,i=\{1,2\}$ and:

$${\mathcal{S}}_{\nu }^{sc}\left({\xi }_i\right)={\mathcal{I}}_{\nu }\left({\xi }_i\right)exp(-{\xi }_i)$$

(12)

with ${\mathcal{I}}_{\nu }\left({\xi }_i\right)$ the Modified Bessel function of the first kind and order $\nu$. In the case of no mispointing (${\theta }_{pitch}={\theta }_{roll}=0$), the ${\Gamma }_{k\ell }\left(0\right)$ function can be written as:

$${\Gamma }_{k\ell }\left(0\right)\sim e^{-{\theta }_{\ell }^2({\xi }_a+{\gamma }_x)}\left\{\begin{array}{cc}{e}^{[-({\gamma }_y+{\xi }_a)\frac{c(t-t_0)}{{\alpha }_{Earth}h}]}\hfill & \mathrm{if}x>x_c\hfill \\ 1\hfill & \mathrm{if}x\le x_c\hfill \end{array}\right.$$

(13)

where ${\xi }_a=1/s^2$ is the attenuation parameter defined as the inverse of the mean square slope parameter s as in [20], and:

$$\begin{array}{ccc}& & {\gamma }_x=\frac{8ln\left(2\right)}{{\theta }_x^2}\hfill \end{array}$$

(14)

$$\begin{array}{ccc}& & {\gamma }_y=\frac{8ln\left(2\right)}{{\theta }_y^2}\hfill \end{array}$$

(15)

with ${\theta }_x$ and ${\theta }_y$ the half-power beam widths, along track and across track, respectively.

Finally, the multi-look SAR waveform over ice covered lakes takes the form:

$$P^{SAR-LIT}(x,\mathbf{p})=\frac{1}{L}\sum _{\ell =-L/2}^{\ell =+L/2}{P}_{SL}^{LIT}({\theta }_{\ell },x,\mathbf{p}),$$

(16)

where L is the total number of Doppler beams and $\mathbf{p}$ is the parameters vector:

$$\mathbf{p}=\{{\Delta }_{ICE}^{gates},{\alpha }_1,{\alpha }_2,{\xi }_a,x_c\}.$$

(17)

Eq. 16 defines the SAR_LITretracker’s waveform model.

The values of the parameters related to the instrument in the SAR_LITretracker model for Sentinel-6 [21] used for the analysis of the SAR data from this mission are reported in Table A1. The ice thickness in unit of meters, ${\Delta }_{ICE}$, is defined by applying the following conversion from range gates to meters:

$${\Delta }_{ICE}={\left({\Delta }_{ICE}^{gates}\right)}_{SAR}\frac{c_{ice}}{4f_s}={\left({\Delta }_{ICE}^{gates}\right)}_{SAR}\frac{c}{n_{ice}4f_s}$$

(18)

where $f_s$ is the sampling frequency, $c_{ice}=\frac{c}{n_{ice}}$ is the speed of light in the ice, with c the speed of light in the vacuum and $n_{ice}$ the refractive index of ice (=1.7861). Note that a factor of four appears in the denominator of Eq. 18, accounting for the factor of two related to the travel time and the factor of two related the range gate oversampling of the SAR waveform.

3.2. FFSAR_LITWaveform Model

The FFSAR LIT waveform model shares the same formalism described in Section 3.1. However, for the Sentinel-6 mission, the FFSAR model is less complex and more computationally efficient than that for UFSAR. Unlike missions operating in closed-burst mode (or lacunar sampling) such as Cryosat-2 and Sentinel-3, the Sentinel-6 altimeter uses an interleaved measurement mode (open-burst mode) that significantly reduces the level of grating lobes in the Sentinel-6 FFSAR impulse response [15,21]. Consequently, the Sentinel-6 FFSAR waveforms are narrower than UFSAR waveforms and can be retracked with the zero-Doppler UFSAR beam model, as reported in [22]. The analytical formulation of the Sentinel-6 FFSAR waveform model is thus given by the single-looked echo of the central beam $L=\ell /2$, taking the form:

$$P^{FFSAR-LIT}(x,\mathbf{p})=P_{SL}^{LIT}({\theta }_{\ell },x,\mathbf{p}){|}_{\ell =L/2},$$

(19)

where $P_{SL}^{LIT}({\theta }_{\ell },x,\mathbf{p})$ is given in Eq. (1). This formulation is not valid for missions operating in closed-burst mode. Their grating lobes are strong and distort the altimeter’s impulse response, which makes the FFSAR waveforms resemble those of UFSAR, with thus similar resolution. In this case, the FFSAR waveform model can be directly emulated from the UFSAR model [22].

3.3. Parameter Estimation

The parameter estimation method is chosen with the goal of building consistent LIT time series with the associated uncertainty that can be exploited for LIT climatology. For this purpose, for each data cycle over a target lake, the LIT retracker analysis, for both UFSAR and FFSAR, is performed within a Region of Interest (RoI) defined as a segment of the satellite ground track with a latitude cut $L{W}_{LIT}=[la{t}_{min},la{t}_{max}]$ chosen to minimize land contamination, thus consisting in the part of the track furthest from the shore, and sufficiently large to include enough statistics to obtain robust and reliable estimates. The LIT analysis consists of two steps: 1) the optimization step, that is, the retracker analysis, involving the waveform fit with the corresponding LIT waveform model, and 2) the estimation of the parameters mean and standard deviation.

The optimization step consists in performing a Least Square Levenberg-Marquardt weighted fit of each waveform in the chosen RoI with the waveform model described in sec:LIT rtks. As in [13], the weights are defined as the data variance within the RoI, ensuring an optimal parameter estimation 4. In practice, the function to be minimized is:

$${\chi }^2={\mathbf{r}}^T{\mathbb{C}}^{-1}\mathbf{r}\simeq \sum _i{\left(\frac{r_i}{{\sigma }_i}\right)}^2,$$

(20)

where, $\mathbf{r}=y\left(\mathbf{x}\right)-S_m(\mathbf{x},\mathbf{p})$ is the vector of residuals between the data, $y\left(\mathbf{x}\right)$, and the waveform model $S_m(\mathbf{x},\mathbf{p})$, and ${\mathbb{C}}_{ij}\equiv {\mathbb{C}}_{ii}={\sigma }_i^2$ is the diagonal noise covariance matrix, with ${\sigma }_i$ computed as the standard deviation of the echoes within the LIT analysis window. Weighting the fit by ${\sigma }_i$ allows to account for the waveform noise.

The second step of the retracking analysis is the parameters estimation, which provides, as the main output, the LIT estimate with the associated uncertainty. The estimation is done by computing the mean and standard deviation of the fitted parameters for each overpass (data cycle) in the RoI. To obtain the parameter constraints with the corresponding uncertainties, for each data cycle, the histograms of the five fitted parameters are computed and a Gaussian fit provides the mean and standard deviation for each parameter. In order to exclude possible outliers, a data editing is applied before computing the histograms, selecting the physical LIT values < 4 m and the LIT estimations within a window of 1 m centered on the mean, which is a conservative choice considering that, for a conservative LIT uncertainty around 15 cm, this corresponds to a window centered within ∼3-$\sigma$.

The LIT uncertainty, being calculated as the standard deviation of the LIT estimates in the LIT analysis window over the target lake, includes contributions related to the retracked radar waveform’s noise and the spatial variation of the LIT estimates within the RoI. In the case of homogeneous ice cover, the dispersion of the LIT measurements (i.e. the width of the histogram) is essentially due to the waveform noise, as is typically the case in the middle of the ice season when the ice is well established and the spatial variations in LIT are small. However, at the seasonal transitions, when the ice is initially forming or is melting, the lake surface is no longer homogeneous, leading to significant variability of the LIT signature in the radar waveforms. This translates in significant spatial variation of the retrieved ice thickness which dominates the error budget, resulting in an increased uncertainty in LIT with respect to the uncertainty expected when considering the waveform noise contribution only.

4. Target Lakes and Data

The lakes selected for the LIT analysis are described in Section 4.1. The data used for the LIT analysis are high resolution UFSAR and FFSAR mode data from the Sentinel-6 mission (Section 4.2.1) and LRM data from the Sentinel-6 and Jason-3 missions (Section 4.2.2). The analysis is performed over two ice seasons (2020-2021 and 2021-2022) during the tandem phase of the Sentinel-6 and Jason-3 satellites (December 2020 to April 2022) to ensure the consistency of the LIT estimates between the two missions, and also to assess the continuity between the LRM and SAR modes for the estimation of LIT. Thermodynamic ice model simulations of LIT time series [16] used for comparison with estimates from the retrackers are described in Section 4.3. It should be noted that in situ LIT measurements are not used for the validation of radar altimetry LIT estimates; no measurements are available for Great Slave Lake during the study period and the few from Baker Lake were collected close to the shore by local authorities at a distance of  30 km from the altimeter track (Figure 3). As previously indicated by [13], in situ LIT measurements gathered as part of the environmental monitoring programs of government agencies are usually taken near the lake shores, while the retrievals from altimetry are performed using data from a segment of the satellite’s ground track chosen the furthest away from shore to avoid land contamination, thus a very different location in terms of bathymetry, wind exposure, snow coverage and depth. These factors play a significant role on the timing of ice formation and the evolution of ice thickness, and this can lead to LIT differences on the order of tens of centimeters (e.g. [16,23]), making in-situ data of limited use for a precise validation of the satellite based LIT estimates.

4.1. Target Lakes

The two target lakes selected for the LIT analysis are Great Slave Lake and Baker lake, both located in the subarctic continental climate zone of Canada (Figure 3). Great Slave Lake is the second largest lake in Canada, covering an area of ≃ 28 000 km², with mean and maximum depths of 41 and 614 m, respectively. Baker Lake is located at a higher latitude (64°09’ N) than Great Slave (61°40’ N), it is also smaller, covering an area of ≃ 1887 km², thus making it a more challenging target for LIT estimation with altimetry data, in particular in LRM, because of the unavoidable land contamination given the size of the footprint which is in the order of 10 km for LRM data. The two lakes are characterized by different environmental conditions that lead to varied snow and ice properties and evolution, making them interesting case studies for the LIT analysis. It is important to note that there are missing data in the SAR HR mode for several cycles over Baker Lake during the first ice season considered for the LIT analysis (2020-2021). This is due to the fact that at the beginning of the mission, the operational acquisition mode was set to LRM with specific masks that included the Baker Lake.

4.2. Altimetry Data

4.2.1. High Resolution Data

The Sentinel-6 UFSAR waveform data at 20 Hz are provided by the Payload Data Acquisition and Processing (PDAP) component of the mission ground segment and are publicly available. Additional Sentinel-6 UFSAR and FFSAR waveform data were generated at a higher posting rate (140 Hz) using an updated version of the SMAP code [24], to assess the potential improvements brought by these new Doppler processing configurations on the LIT estimation. By increasing the posting rate of the altimeter observations beyond the standard 20 Hz, which corresponds on-ground to the UFSAR altimeter resolution, more measurements are made along the satellite track to possibly extract more information. The along-track ground spacing is reduced from 350 m to 50 m. However, in contrast to FFSAR, the spatial resolution of the UFSAR data remains unchanged and is equal to that of the 20 Hz UFSAR whatever the posting rate used.

4.2.2. Low Resolution Data

The data used in the analysis are the 20 Hz LRM waveform data from the Jason-3 mission (SGDR L2 products available from the Aviso+ portal5) and the Sentinel-6 PDAP LRM waveform data provided by the Sentinel-6 ground segment, which are also publicly available.

4.3. Lake Ice Model Simulations

Lake ice models provide an alternative to in situ measurements that are scarce in space and in time. Such models are particularly valuable for the development of new remote sensing based LIT retrieval algorithms and their evaluation since in situ measurements that match the timing and locations of overpasses are most often not available or limited in terms of coverage (i.e. area covered and temporal frequency of measurements). Using lake ice models, ice thickness can be simulated on a daily basis for areas of differing water depth, snow density, and snow depth. The Canadian Lake Ice Model (CLIMo, [16]) has been used regularly in development of remote sensing algorithms as well as interpreting remote sensing signals ([6,18,25,26,27]). CLIMo has also been used to investigate the impact of contemporary and future climate conditions on ice phenology (freeze-up, break-up, and ice duration), LIT, and the types of ice (snow ice and clear ice) (e.g., [3,23,28]). CLIMo is used in this study to simulate LIT for comparisons with retrievals from the altimetry data for the areas of Great Slave Lake and Baker Lake highlighted in Figure 3. CLIMo is a 1-D thermodynamic ice model and is forced using meteorological variables that can be provided from weather stations, atmospheric reanalysis, or climate model gridded products. Input data for CLIMo is provided daily and consists of 2-m air temperature, relatively humidity, cloud cover, wind speed, and snowfall (or snow depth from a nearby land site when available). Additionally, a fixed mixed layer depth and average snow density is required for the lakes of interest. The outputs from CLIMo include energy/radiation balance components, on-ice snow depth (both accumulation and melt), the temperature profile at specified depths within the snow/ice, and ice thickness. The model also provides the freeze-up and break-up dates for the lake of interest.

For this study, 2-m air temperature, wind speed, relative humidity, and cloud cover were retrieved from the ERA5 atmospheric reanalysis product from the European Centre for Medium-Range Weather Forecasts (ECMWF). Snow depth data were not used from ERA5 and were instead acquired from Environment and Climate Change Canada weather stations on the shores of GSL and Bake Lake. For GSL, the Hay River weather station was used (115.78° W, 60.84° N) and for Baker Lake the Baker Lake station was used (96.08° W, 64.3° N). Simulations were performed for the years corresponding with the altimetry acquisitions with a mixing depth of 30 m for GSL and 20 m for Baker Lake. Average snow density for GSL was set to 300 kg m^-3, which has been used for previous simulations of this lake ([6]). For Baker Lake, simulations were run at varying snow densities (200, 250, 300, and 350 kg m^-3) as fewer past simulations have been conducted for this lake and this will account for variations in the snow cover. To account for variations in snow drifting on the lake ice surface CLIMo was run with four sets of snow depth scenarios: no snow, 25, 50, and 75%. From the simulation results on-ice snow depth, ice thickness and melt onset were retained to evaluate the estimates obtained from the altimetry data. From past simulations of LIT using CLIMo, RMSE values have been reported between 5-17 cm when compared to in situ measurements collected near Yellowknife ([23,27]), similar errors are expected in this study.

5. Results and Discussion

This section presents the results of the LIT analysis performed with the data described in Section 4. The results focus on the LIT retrievals obtained with the SAR_LITretracker applied to the Sentinel-6 UFSAR data at 20 Hz, which are compared to the LIT estimates obtained with Sentinel-6 data at higher resolution, that is, the UFSAR and FFSAR data processed at 140 Hz (Section 5.1). The assessment of the SAR_LITanalysis, consisting of the comparison with the LIT analysis performed on LRM data from both Sentinel-6 and Jason-3 missions and with the output of the thermodynamic CLIMo simulations of the same targets and periods, is presented in Section 5.2. Consistency checks of the LIT retrievals from radar altimetry with satellite images are also shown in this section.

5.1. LIT Analysis of Sentinel-6 High Resolution Data

An illustrative example of the SAR_LITanalysis of high resolution Sentinel-6 SAR data at 20 Hz for Great Slave Lake during the 2020-2021 ice season is given in Figure 4. Examples of Sentinel-6 waveform data (in black) with, in blue, the SAR_LITfit are shown in the left column. From top to bottom the plots correspond to December 2020, February 2021, April 2021 and mid-May 2021, respectively. The evolution of the bimodal LIT signature is clearly visible, starting with a small peak separation during the ice formation period, which increases through the winter, indicating thickening of the ice, and then disappears at the beginning of the melt period. Plots in the right column show the histograms of the LIT measurements obtained with the SAR_LITretracker in the lake RoI on the same dates. The blue curve corresponds to the Gaussian fit of the histograms. The changes in the histograms distribution and mean reflect the evolution of the ice thickness within the considered region and during the ice season. The LIT standard deviation ranges from 3.6 cm to 5.9 cm during the ice season, except during the melt period, where it is bigger, in the order of 15 cm, due the limited penetration depth of the signal and higher spatial variation of the ice surface within the RoI. In fact, because of the melting snow/ice, most of the waveforms do not contain the LIT signature (as in the bottom-left plot example), yet, some waveforms still have it, as the ice is still present and the melting is not homogeneous over the considered region. In this case, the spatial variation of the LIT is higher, thus resulting in a higher standard deviation of the LIT measurements within the RoI.

Similar performances and behavior of the LIT evolution are obtained for Baker Lake. A representative example showing the spatial evolution of LIT at Great Slave and Baker lakes in February 2021 is given in the left and right upper panels plots of Figure 5, respectively. The bottom panels show the evolution of the reduced ${\chi }^2$ statistics (Eq. (19)) at the best fit, providing a goodness of fit metric. The fit performances are good for both lake targets, showing an overall reduced ${\chi }^2$ of one, meaning that the model and the data are in agreement within the variance. It can be noted that for Baker Lake there are less data points in contract to Great Slave Lake, as it is a smaller lake and thus the ground track segment considered for the LIT analysis is shorter (see Figure 3), including only ∼ 40 waveforms, instead of ∼ 120 waveforms in the RoI of Great Slave. This shows the limitation of building reliable LIT time series from small target lakes by using 20 Hz data. A possible solution to overcome this limitation is to perform the LIT analysis on the mean waveform in the RoI, instead of fitting the individual waveforms. Yet, if the number of exploitable waveform data is too small (< 10/20 waveforms) this could also lead to less robust LIT estimates. Another possibility would be to process the data at a posting rate higher than 20 Hz to increase the statistics, as shown for instance in [20,29]. With respect to this, as detailed below, the LIT analysis has been performed on Sentinel-6 SAR data processed at 140 Hz, which has the advantage of increasing the number of waveforms in the RoI, and thus the spatial sampling of the surface, by a factor of seven in comparison to the 20 Hz data.

Figure 6 shows representative results of the LIT estimation obtained with high resolution Sentinel-6 data for a complete ice season at Great Slave (left plot) and Baker (right plot) lakes. The red circles correspond to the LIT estimates obtained with the Sentinel-6 UFSAR data at 20 Hz, the purple circles to the UFSAR data at 140 Hz, and the cyan circles to the FFSAR data at 140 Hz. The shaded regions in the upper panels and the lines in the bottom panels correspond to the LIT 1-$\sigma$ uncertainties. Overall there is a good consistency in the LIT estimates obtained with the three datasets for both lake targets. As expected, the LIT retracker analysis performed on data at higher posting rate (140 Hz) shows an increased performance with respect to the 20 Hz data, especially at the melt transition, because of the increased statistics (a factor of seven more data points and, therefore, a higher spatial sampling of the surface). Also, at an equivalent posting rate, the FFSAR data allow for a better accuracy of the LIT estimates, overall yielding ∼20% smaller uncertainties, in particular in the middle of the ice season, once the ice is well established on the lakes. This is related to the increased spatial resolution of the FFSAR data resulting in a smaller and more homogeneous footprint from which the LIT estimates are retrieved. Finally, noteworthy is that a drop in LIT is noticeable in the three datasets at Baker Lake, at the beginning of May 2022, before the final transition into the melt period, which is due to the temporary melt followed by refreezing of the snow on the ice surface, as discussed in more details in the next section.

Overall, the results show a significant improvement in the accuracy of LIT estimates in comparison to previous analysis, in particular when using conventional altimetry data, thanks to the higher resolution of the Sentinel-6 SAR data.

5.2. Evaluation and Consistency

The evaluation of the SAR_LITretracker analysis is done by quantitatively comparing the Sentinel-6 UFSAR 20 Hz LIT estimates to those obtained from the thermodynamic simulations (CLIMo), described in sec:CLIMo, and to the LIT estimates derived from the LRM_LITanalysis [13] of Low Resolution Mode (LRM) data of both Sentinel-6 and Jason-3 during their tandem phase. Assessing the consistency between the SAR and LRM estimates is particularly important to ensure the stability of the LIT measurements between current and future altimetry missions. The LIT indicators and metrics used for the comparison are, as summarized in Table 1 for both lake targets, the LIT maximum and the corresponding date, the mean LIT in the middle of the ice season when the ice is well established on the lakes (1st February to mid-April), with the associated uncertainty, the Mean Bias Error (MBE), and the Root Mean Square Error (RMSE) between the Sentinel-6 UFSAR 20Hz and the other datasets. It is worth noting here that in the Table 1 the LIT indicators are reported for the 2020-2021 season for Great Slave and for the 2021-2022 season for Baker due to a fair amount of missing waveform data in the other ice seasons (missing input data for 5 cycles for Great Slave in the 2021-2022 season and missing input data for 10 cycles for Baker in the 2020-2021 season). As shown in the Figure 7 and Figure 8, despite the missing data, the LIT evolution is consistently captured for both seasons and lakes, and the LIT indicators reported are, therefore, representative of a typical ice season.

Figure 7 and Figure 8 summarize the comparison of LIT among the different datasets for Great Slave Lake and Baker Lake, respectively. In both figures, the top panel shows LIT estimates obtained with the Sentinel-6 UFSAR at 20 Hz data (red), the Sentinel-6 LRM data (green) and the Jason-3 data (blue) for the 2020-2021 and 2021-2022 ice seasons during the Sentinel-6 and Jason-3 tandem phase (December 2020-April 2022). The shaded regions of the corresponding colors refer to the LIT error envelops at 1-$\sigma$ for each case. The orange shaded area shows the evolution of LIT obtained with the CLIMo thermodynamic simulations with different snow-on-ice scenarios: the upper bound refers to the simulations without snow on the ice surface, while the lower bound to the simulations with 75% of snow with respect to snow depth measured at nearby meteorological stations on land (see Section 4.3 for details). The middle panel in both figures shows the evolution of the mean 2-m daily air temperature (black line) with the grey shaded area delineating the minimum and maximum daily temperatures extracted from ERA5 hourly data. Plotting air temperature measurements along with the LIT estimates is useful since it provides some information about the status of the lake surface, revealing snow-on-ice melt when the air temperature is near or greater than 0°C, and refreeze when temperatures drop again below 0°C during the spring transition period. Finally, the bottom panel highlights the evolution of the LIT standard deviation derived from the three datasets (same legend as in the top panel).

Overall, there is a very good consistency among the SAR and LRM LIT results. The MBE and RMSE of Sentinel-6 SAR and LRM LIT estimates are 1.7 cm and 3.3 cm, respectively, for Great Slave Lake and 2.3 cm and 3.2 cm, respectively, for Baker Lake. The agreement with Jason-3 is also very appreciable, with MBE and RMSE of 4 cm for Great Slave and 7 cm for Baker. The maximum LIT and the LIT mean are also estimated consistently within a few centimeters difference among the different datasets. The dates of maximum LIT are also fully consistent between the SAR and LRM. The LIT evolution is in agreement with the output of CLIMo simulations, within the uncertainty envelop given by the different snow-on-ice scenarios used as input. For both lake targets and for the seasons considered for the quantitative analysis, the satellite based LIT estimates seem to be more in agreement with a small amount of snow (25 %) or no snow on ice in terms of LIT maximum (values and date) and mean, and with MBE and RMSE values that are the smallest for the simulations with no snow on-ice, being -7 cm and 8 cm, respectively, for Great Slave and -3 cm and 5 cm, respectively, for Baker Lake.

As shown in the figures, LIT drops are detected right before the maximum is reached, at the end of the ice season, during the 2021-2022 season for Great Slave and during both seasons for Baker Lake. It can happen that the snow/ice on the lake surface starts to melt due of an increase of temperatures, as discussed above, in which case the LIT signature is no longer present in the radar waveforms, so no or a very small LIT is retrieved. Then, with refreezing, the LIT signature is again detectable. As shown in the middle panels of Figure 7 and Figure 8, a clear correlation can in fact be seen between the detected LIT drops and the rise of the air temperature near or above 0°C, supporting this explanation. As pointed out in [13], it is worth reiterating here that radar altimetry LIT retrackers can indeed capture the seasonal transitions of ice forming and melting but cannot precisely follow the ice evolution at the transitions because of the difficulty of retracking heterogeneous surfaces when the ice is too thin (fall freeze-up period) and when the snow on the ice surface begins to melt (spring break-up period). For this reason, as can be seen in the figures, the satellite based LIT estimates drop at the beginning of melt, thus earlier with respect to thermodynamic CLIMo simulations that realistically generate LIT estimates throughout the whole melt period.

It can be noted that the LIT evolution is different at Great Slave and Baker, which is expected since of the lakes are located at different latitudes and have distinct characteristics as indicated earlier. For instance, the LIT maximum measured with Sentinel-6 20 Hz UFSAR occurs on April 15 and is $1.486\pm 0.025$ m at Great Slave Lake, while for Baker Lake it generally occurs roughly one month later, on May 17, and is 50 cm greater reaching $2.005\pm 0.0027$ m. Also, the mean LIT in the middle of the ice season (as measured from February to April) is lower at Great Slave ($1.22\pm 0.047$ m) compared to Baker ($1.709\pm 0.05$ m). These results show that the SAR_LITretracker can precisely follow the LIT evolution on different target lakes.

As expected, and as shown in the bottom panel of Figure 7 and Figure 8, the LIT estimates obtained from SAR data, because of the improved spatial resolution, show a significantly smaller dispersion with respect to the LIT estimates obtained with LRM data, thus yielding to a LIT uncertainty that is a factor of ∼ 2, up to 3, times smaller, depending on the season and target, with respect to the LRM estimates. For instance, as detailed in Table 1, the accuracy of the maximum LIT estimation during the 2020-2021 season at Great Slave Lake with Sentinel-6 SAR data is 2.5 cm, compared to 5.1 cm with Sentinel-6 LRM data and 6.4 cm with Jason-3 LRM data, thus a factor of 2 and 2.6 better, respectively.

Finally, we performed consistency checks by superimposing the LIT estimates obtained with the SAR_LITretracker applied to Sentinel-6 20 Hz UFSAR data onto optical (Landsat-8) or radar (Sentinel-1) images, according to the images availability on the Sentinel-6 flyover dates over the two lakes. This is illustrated in Figure 9 for Great Slave Lake (left) and Baker Lake (right), where the Sentinel-6 LIT estimates along the ROI track segment are overlaid on the images and color coded, ranging from 0 to 1.6 m for Great Slave and from 0 to 2 m for Baker. The images and LIT estimates are shown for representative dates following the typical duration of the ice season for the two lake targets, from January to May for the Great Slave lake and from December to June for the Baker lake. Overall, there is a very good consistency between the LIT estimates and the lake surface conditions for both target lakes, with thinner ice detection at the beginning of the ice season when the images show an heterogeneous surface of initial ice formation and open water (top panels), consistently growing LIT in the middle of the ice season when the snow/ice cover uniform over the lake surface (second and third panels), and drop of the LIT estimates to small values at the end of the season with snow melt, as shown in the bottom panels.

6. Conclusions and Outlook

In this paper, we presented a novel, efficient, analytically based retracking approach for the estimation of LIT from high-resolution Ku-band SAR (FFSAR and UFSAR) altimetry data. We carried out the LIT analysis with Sentinel-6 UFSAR data at 20 Hz and with UFSAR and FFSAR data at 140 Hz. We evaluated the SAR_LITand FFSAR_LITretracker analysis with LIT obtained from numerical lake ice model simulations and estimates from Sentinel-6 and Jason-3 conventional LRM altimetry data at 20 Hz for two ice seasons during the tandem phase of the two satellites. We performed consistency checks by superimposing the along-track Sentinel-6 LIT estimates on optical and radar images. We performed the LIT analysis and assessment for Great Slave Lake and Baker Lake (Canada), two lakes that differ in terms of their size, bathymetry, and snow/ice properties. The results demonstrate the ability of the SAR_LITand FFSAR_LITretrackers to retrieve robust LIT estimates over different ice seasons and for lakes of varied characteristics. The estimates are in high agreement with those obtained with LRM data, ensuring the continuity of LIT measurements and time series from conventional altimetry missions as well as current and future SAR altimetry missions. The typical LIT uncertainty obtained with the SAR_LITretracker is of the order of 5 cm, when the ice is well established on the lake, representing a significant improvement with respect to previous studies (as e.g. [7,8,9,10,11]). As expected, LIT estimates are of significantly higher accuracy, by a factor of 2 to 3 better, when using Sentinel-6 UFSAR data at 20 Hz in contrast to LRM data due to increased resolution. The FFSAR processing at 140 Hz provides even better LIT estimates with  20% smaller uncertainties. The LIT retracker analysis performed on data at the higher posting rate (140 Hz) shows increased performance in comparison to the 20 Hz data, especially during the melt transition period, due to the increased statistics (a factor of seven times more data points and, therefore, a higher spatial sampling of the lake surface).

As previously reported by [13], the main limitation of LIT retrackers based on Ku-band waveform data is that they work well if the ice-related signature, that is, the double peak feature in the Ku-band radar waveforms, is present. Indeed, the freshwater ice signature depends on the properties and thickness of the snowpack and the ice cover. The double peak feature in the waveform may be absent if some conditions are not met, for instance in the case of melting snow on the ice surface or snow-free lake ice. For this reason, the SAR_LITand FFSAR_LITretrackers can capture snow/ice melt onset on a lake, but they cannot precisely follow the evolution of LIT during melt periods since the radar waves are reflected by the surface (i.e., they do not penetrate through wet snow and ice). Also, there may be lakes for which the ice signature is not clearly visible or not present, as is the case for snow-free areas on lakes. In such a situation, consistent LIT estimates and time series cannot be generated for these lake targets. Overall, as long as the LIT signature is present in the radar waveforms, the SAR_LITand FFSAR_LITretrackers can precisely capture the seasonal evolution and the inter-annual variability in LIT, making them powerful tools for obtaining robust LIT estimates required for climate monitoring.

This research has identified several topics that merit further investigation. First, Ku-band SAR LIT time series should be generated for different lake targets and from different SAR altimetry missions, for instance, Sentinel-6, Cryosat-2, and Sentinel-3, enriching the LIT time series datasets currently being generated within the ESA Lakes_cci project, which consists of LIT long time series generated from LRM data (the LIT product over one target, the Great Slave lake, has been publicly released at the end of 2023, LIT products for more lakes will be released in 2025). Sentinel-3 is of particular interest since it allows for the inclusion of lakes located at latitudes higher (above 66 deg N) than is possible from Jason/Sentinel-6 altimetry missions, enabling a broader spatial coverage of LIT monitoring. Another key topic is understanding the propagation of radar waves through snow and ice of various properties at different frequencies (e.g., Ka and Ku-band). Within this context, the analysis of the LIT signature and its evolution in SAR waveforms over well-characterized lake targets, such as Great Slave Lake, could be used as a reference to better understand the properties of the different layers the radar waves pass through and the interfaces leading to the propagation and/or reflection of the radar waves. Radiative transfer modeling experiments on the impact of ice and overlying snow properties on altimeter waveforms (Ka/Ku-band) at several lake sites are needed in parallel, building on the recent ESA LIAM project [30] that established a modeling framework that integrates the Snow Microwave Radiative Transfer (SMRT) model and the lake ice model CLIMo with simulations that can currently be performed in LRM, but that should be extended to SAR mode. A multi-sensor/frequency approach is also needed to assess the possibility of retrieving on-ice snow depth. This topic is particularly relevant in view of the future Copernicus Polar Ice and Snow Topography Altimeter (CRISTAL) mission6, planned for launch in 2027, that will carry a dual-frequency Ku/Ka-band radar altimeter [31]. In the case of freshwater lake ice covered by dry snow, Ka and Ku-band radar waves may be reflected at the air/snow (or at other interfaces within the snowpack) and snow/ice interfaces, and the snow/ice and ice/water interfaces, respectively. However, the propagation of the radar waves at different frequencies depends on the properties of ice and overlying snow properties (e.g., ice type, snow density and surface roughness) and is not fully understood yet. The analysis of Ka/Ku-band radar waveforms on target lakes of known ice/snow properties observed during field campaigns would offer a robust way to characterize the signals. Additionally, future work should investigate the impact of changing ice roughness on the retrieval by retrackers. The roughness of interfaces within the ice and snow column has been an area of interest in recent years and further investigation could provide better insights into the impact of these properties on radar waves and could be further supported by radiative transfer models such as SMRT. The lakes could be considered as cal/val sites for regular monitoring and assessing the performance of the dual-band measurements of the CRISTAL upcoming mission. A SAR_LITretracker, based on the work presented in this paper and tailored for the CRISTAL mission, is currently being developed and validated within the ESA CLE2VER project and will be implemented in the CRISTAL prototype processor.