# River ice phenology and thickness from satellite altimetry: potential for ice bridge road operation and climate studies

Elena Zakharova, Svetlana Agafonova, Claude Duguay, Natalia Frolova, Alexei Kouraev

► To cite this version:

Elena Zakharova, Svetlana Agafonova, Claude Duguay, Natalia Frolova, Alexei Kouraev. River ice phenology and thickness from satellite altimetry: potential for ice bridge road operation and climate studies. The Cryosphere, 2021, 15, pp.5387-5407. 10.5194/tc-15-5387-2021. insu-03671316

HAL Id: insu-03671316

https://insu.hal.science/insu-03671316

Submitted on 18 May 2022

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_1.png" /></figure>

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

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

# River ice phenology and thickness from satellite altimetry: potential for ice bridge road operation and climate studies

Elena Zakharova<fnref n="1" />,<fnref n="2" />, Svetlana Agafonova<fnref n="3" />, Claude Duguay<fnref n="4" />,<fnref n="5" />, Natalia Frolova<fnref n="3" />, and Alexei Kouraev<fnref n="6" />

**Correspondence:** Elena Zakharova (zavocado@gmail.com)

Received: 2 November 2020 – Discussion started: 16 December 2020

Revised: 18 October 2021 – Accepted: 1 November 2021 – Published: 7 December 2021

**Abstract.** River ice is a key component of the cryosphere. Satellite monitoring of river ice is a rapidly developing area of scientific enquiry, which has wide-ranging implications for climate, environmental and socioeconomic applications. Spaceborne radar altimetry is widely used for monitoring river water regimes; however, its potential for the observation of river ice processes and properties has not been demonstrated yet. Using Ku-band backscatter measurements from the Jason-2 and Jason-3 satellite missions (2008–2019), we demonstrate the potential of radar altimetry for the retrieval of river ice phenology dates and ice thickness for the first time. The altimetric measurements were determined to be sensitive enough to detect the first appearance of ice and the beginning of thermal breakup on the lower Ob River (Western Siberia). The uncertainties in the retrieval of ice event timing were within the 10 d repeat cycle of Jason-2 and Jason-3 in 88 %–90 % of the cases analysed. The uncertainties in the river ice thickness retrievals made via empirical relations between the satellite backscatter measurements and in situ observations, expressed as the root mean square error (RMSE), were of 0.07–0.18 m. A novel application of radar altimetry is the prediction of ice bridge road operations, which is demonstrated herein. We established that the dates of ferry closing and ice road opening and closing in the city of Salekhard can be predicted with an accuracy (expressed as RMSE) of 3–5 d.

## 1 Introduction

River ice is a major component of the cryosphere. It affects an extensive portion of the hydrologic system in the Northern Hemisphere, where seasonal ice impacts 58 % of the total river length (Prowse et al., 2007). Hence, monitoring the interannual variability and changes in the river regime in response to climate change is critical for many environmental and societal applications. River ice plays a key role in the functioning of aquatic and riparian ecosystems (Prowse, 2001) and contributes to the erosion of channels and banks (Ettema, 2002) as well as the transport of sediments (Beltaos et al., 2018). River ice affects streamflow via the withdrawal (immobilisation) of part of the water during freeze-up and the consequent release during breakup. In addition, ice jams on rivers can cause catastrophic flooding (Beltaos, 2013). From socioeconomic and cultural perspectives, river ice affects the operation of hydropower stations as well as construction and navigation activities. In Arctic regions, frozen rivers provide a unique transportation infrastructure for the movement of goods and people via winter ice roads and provide local populations with access to fishing and hunting/trapping grounds and in some cases (e.g. central Yakutia, Russia) to freshwater. However, the operational monitoring of ice cover on northern rivers is difficult due to site accessibility. Moreover, ice conditions can be unsafe for people who perform in situ measurements, especially at the beginning or end of the ice season. Therefore, satellite remote sensing has been proposed

<footnote n="1">¹Water Problems Institute, Russian Academy of Science, Moscow, Russia</footnote>

<footnote n="2">²EOLA, Toulouse, France</footnote>

<footnote n="3">³Department of Land Hydrology, Moscow State University, Moscow, Russia</footnote>

<footnote n="4">⁴Department of Geography and Environmental Management, University of Waterloo, Waterloo, Canada</footnote>

<footnote n="5">⁵H₂O Geomatics, Waterloo, Canada</footnote>

<footnote n="6">⁶LEGOS, Université de Toulouse, CNES, CNRS, IRD, UPS, Toulouse, France</footnote>
as an alternate tool or as a complement to ground-based mea-
surements, allowing for the characterisation of river ice prop-
erties and processes at spatial and temporal resolutions suit-
able for addressing various scientific and operational require-
ments (Duguay et al., 2015).

Indeed, spaceborne instruments provide observational ca-
pabilities that can be used to map and monitor many river ice
parameters. For example, optical sensors, such as the Mod-
erate Resolution Imaging Spectroradiometer (MODIS) and
the Advanced Very High Resolution Radiometer (AVHRR),
have been used to map the river ice extent and phenology,
including freeze-up and breakup dates (Pavelsky and Smith,
2004; Chaouch et al., 2014; Chu and Lindenschmidt, 2016;
Muhammad et al., 2016; Cooley and Pavelsky, 2016; Beaton
et al., 2019). However, the presence of extensive cloud cover
for many months of the year and low solar illumination con-
ditions, particularly during the freeze-up period, are limiting
factors for monitoring river ice at high latitudes using opti-
cal sensors. Active sensors operating in the microwave re-
region of the electromagnetic spectrum are weather indepen-
dent and provide a spatial resolution higher than that pro-
vided by the MODIS and AVHRR instruments. Synthetic-
aperture radar (SAR) imaging systems, in particular, have
been used to monitor river ice phenology (e.g. Unterschultz
et al., 2009; Mermoz et al., 2009; Sobiech and Dierking,
2013; Sun and Trevor, 2018), deformation (Unterschultz et
al., 2009) and classification of ice types (e.g. Chu and Lin-
denschmidt, 2016). Comprehensive reviews on remote sens-
ing of lake and river ice can be found in Jeffries et al. (2005)
and Duguay et al. (2015).

Ice thickness (H ICE) is another parameter which is of
particular interest for operational purposes, such as public
safety, ice road service, ice jam forecasts and mitigation. Pas-
sive microwave and thermal satellite instruments have exhib-
ited potential for the retrieval of ice thickness for large lakes
(Kang et al., 2014; Duguay et al., 2002, 2015; Gunn et al.,
2015; Kheyrollah Pour et al., 2017). However, the spatial di-
mension of rivers, notably the width of the channels, lim-
its the application of these instruments, as they only provide
coarse spatial resolution (kilometres to tens of kilometres). A
few studies have used SAR images with high (10–25 m) spa-
tial resolutions to estimate river ice thickness (Unterschultz
et al., 2009; Mermoz et al., 2013; Zhang et al., 2019), mainly
through the establishment of the statistical relations between
backscatter and in situ ice thickness measurements. However,
these SAR-based methods have been limited to the estima-
tion of ice thickness for short time periods (one or a few ice
seasons) and small sections of rivers, primarily due to lim-
ited access to SAR imagery in years prior to the launch of
the Sentinel-1A and Sentinel-1B missions (2014) of the Eu-
ropean Space Agency (ESA). Altimetric radar represents an
excellent alternative or complement to SAR satellites. Satel-
lite altimeters also provide relatively high spatial resolution
(200–400 m) and are weather-independent instruments with a
long history of observations starting in the mid-1980s. Radar

altimeters have largely been used for monitoring water levels
of inland water bodies and water courses over 100 m in width
(Michailovsky et al., 2012). Radar altimeter data products are
freely accessible via portals, such as the ones operated by the
French space agency (CNES, Centre national d'études spa-
tiales) and ESA; data can be routinely and quickly processed
over large areas for relatively long time periods.

Radar signals incident upon Earth's surface are modified
according to the physical properties of the reflecting sur-
face. Similar to SAR imaging systems, the signal recorded
by radar altimeters can be interpreted as a function of the
changes in surface properties, and the backscatter coefficient
(ratio between the power of the emitted and received sig-
nals) can be used to characterise the surface state within a
radar footprint. Radar backscatter from freshwater ice de-
pends on the radar configuration (radar type, operating fre-
quency, viewing geometry and polarisation) and on the ma-
terial properties, such as ice thickness, layering, air bubble
inclusions, liquid water content, surface roughness and snow
on ice (Ulaby et al., 1986; Duguay et al., 2002; Leconte
et al., 2009; Atwood et al., 2015; Gunn et al., 2015, 2018;
Antonova et al., 2016). The potential of altimetric radars has
already been demonstrated for the monitoring ice phenology
dates of large Eurasian lakes (Kouraev et al., 2007, 2015) and
for the estimation of ice thickness on large northern lakes in
Canada (Beckers et al., 2017; Duguay et al., 2018). How-
ever, radar altimeters have not been used for monitoring the
ice regime of rivers. Therefore, for the first time, we demon-
strate the application of radar altimetry for monitoring the in-
terannual variability of river ice phenology and ice thickness
along a 400 km stretch of the Ob River (Siberia, Russia).

This paper describes the development and validation of the
algorithms for the retrieval of river ice phenology dates and
ice thickness and the use of these products for the predic-
tion of ice bridge road operations in the city of Salekhard
on the Ob River. For this purpose, data from two altimet-
ric satellite missions, Jason-2 and Jason-3, were selected for
two reasons: (1) they provide the best temporal resolution
(10 d) amongst the altimetric missions, and (2) the lifetime
of this series of satellites provides observations on the same
orbit, starting in 1992 with the TOPEX/Poseidon (Topog-
raphy Ocean Experiment) mission and continuing with the
Sentinel-6 mission (launched in November 2020), which is
offering an interesting perspective for a follow-up investiga-
tion on the response of river ice on the Ob to a changing
Arctic climate.

## 2 Study region

The study was conducted over an area that encompasses
the lower reaches of the Ob River (Siberia, Russia). The
Ob River drains the West Siberian Plain and is the third
largest river in the Arctic Ocean watershed, with a mean an-
nual flow of 406 km³ (Zakharova et al., 2020). The lower
<table><caption>**Table 1.** Gauging stations in the lower reaches of the Ob River.</caption><thead><tr><th>River - station</th><th>Distance from mouth (km)</th><th>Beginning of observation period</th><th>Data gap in observation record</th></tr></thead><tbody><tr><td>Ob - Polnovat</td><td>702</td><td>1970</td><td></td></tr><tr><td>Ob - Gorki</td><td>487</td><td>1935</td><td></td></tr><tr><td>Small Ob - Muzhi</td><td>463</td><td>1933</td><td></td></tr><tr><td>Ob - Kazym Mys</td><td>551</td><td>1979</td><td>1988-2003</td></tr><tr><td>Ob - Pitlar</td><td>386</td><td>1979</td><td>1990-2005</td></tr></tbody></table>

reaches of the Ob River extend approximately 800 km and begin at the confluence of the Ob and Irtysh rivers at 61.08° N (Fig. 1). The reaches are characterised by a particularly wide floodplain (up to 50 km) and numerous branches. The easternmost channel is the primary and largest branch and is called the “big Ob River”. The second largest branch delineates the flood plain from the west. The Ob River watershed is one of the largest peat bog systems in the world (Zakharova et al., 2014), and many settlements, which are located on the high terraces of the two main branches, have limited interconnection and access to supplies. The main branches of the river are navigable; however, they are covered by ice for 7 months of the year. In winter, when the bogs are frozen, local communities intensify their socioeconomic activities by constructing winter roads and ice bridges over the river crossings. River ice observations on the Ob are sparse and taken only at a few gauging stations dedicated to water level monitoring. In this study, we selected a section of the lower reaches of the Ob River located between two large administrative centres (Salekhard and Khanty-Mansiysk).

# 3 Data

## 3.1 In situ data and ice regime at studied river reaches

The Hydrometeorological Research Center monitors ice at all gauging stations that measure the water level. There are five water level gauging stations in the studied Ob River reach (Fig. 1 and Table 1). Four stations (Polnovat, Gorki, Kazym Mys and Pitlar) are located on the main branch of the Ob River, and one station (Muzhi) provides observations on the secondary channel called the “small Ob River”.

The standard protocol of river ice monitoring at the gauging stations involves daily visual observations of ice presence/absence and ice type/events and measurements of ice thickness and on-ice snow depth three to four times per month. Ice thickness is measured by drilling a single hole using an ice auger. The snow depth records are represented as the average value of three snow depth measurements located around the hole.

According to the in situ observations at the gauging stations, ice formation begins on average between 23 and 27 Oc-

<figure><img src="image_4.png" /><figcaption>**Figure 1.** The lower reaches of the Ob River and location of the virtual (circles) and gauging (red crosses) stations. Stations from the training set are shown as yellow circles. Jason-2 and Jason-3 satellite tracks numbers are also provided at the river crossings. The main map was produced using the public World Bank data (https://datacatalog.worldbank.org/dataset/major-rivers-world, last access: 2 December 2021).</figcaption></figure>

tober. In the last 20 years, the earliest and latest records are 1 October and 18 November, respectively. In the studied river reaches, ice cover forms quickly, typically within just 2–3 d from the first appearance of ice. However, according to the records provided by gauging stations, in 15% of the cases, full freeze-up (from ice onset to complete freeze-over) takes up to 10 d. Ice grows rapidly during the first month of the ice season. By the end of November, the river ice reaches a thickness of 0.23–0.30 m. Towards the end of the ice growth period (March–April), the average thickness of the ice cover is 0.80–1.0 m (maximum of 1.50 m). Snow depth on the ice surface varies from 0.09–0.13 m in November to 0.30–0.50 m in April. The temporal dynamics of ice growth on the large linear-channel sections of the Ob River is similar throughout the studied reaches (Fig. 2a). However, in the southern part of the region, the ice thickness is 0.07 to 0.20 m less than that in the north. Climate change has affected river ice in the Canadian Arctic (Prowse et al., 2011a, b) and the European part
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Examples of (a) seasonal evolution of ice thickness (1980–2017 mean) and (b) interannual variability of ice freeze-up and the end of breakup dates and maximum ice thickness at four stations (Gorki, Muzhi, Polnovat and Kazym Mys).</figcaption></figure>

of Russia (Agafonova and Vasilenko, 2020). Mann–Kendall tests performed on the ice time series with complete records (Polnovat, Gorki and Muzhi; see Table 1) have revealed that ice onset had occurred significantly later (0.05 level) only at the Polnovat station. Maximum ice thicknesses have become significantly thinner (0.05 level) only at the Gorki and Muzhi stations (Fig. 2b). The observed trends are of the order of 3 d every 10 years and 3 cm every 10 years for the ice onset and maximum ice thickness, respectively.

## 3.2 Altimetry and radiometry data

Data from the Jason-2 and Jason-3 satellite missions covering the period from July 2008 to August 2019 were used in this study. The Jason-2 satellite is the third altimetric satellite of the TOPEX/Poseidon-Jason series. The satellite operated between 2008 and 2016 and acquired data on a 10 d repeat orbit with an inclination of 66.08°. One satellite cycle consisted of 127 revolutions and 254 tracks (odd numbers for ascending and even numbers for descending orbits). In 2016, Jason-3 was sent into space with the same orbit and instruments as Jason-2. For 20 cycles, the two missions flew in a tandem mission with an 80 s time lag, ensuring continuity of measurements. The altimetric radar on board Jason-2 and Jason-3 provided along-track measurements at the Ku-band (13.6 GHz) and C-band (5.3 GHz) with a sampling frequency of 20 Hz, allowing for a 375 m distance between adjacent radar measurements. The ground track repeatability of the mission was maintained within a ±1 km cross-track at the Equator. At the latitudes of our study area (63–66° N), the cross-track oscillation band was approximately 400 m wide. In this study, only measurements obtained in Ku-band were used. The theoretical footprint of the radar at the Ku-band is 10–12 km in diameter over a rough ocean surface (Fu and Cazenave, 1991). However, this diameter decreases over smooth inland water and ice surfaces such that the main re-

turned signal can come from footprints of just a few kilometres in diameter (Legrésy and Rémy, 1997).

The satellite payload of Jason-2 also included the nadir-looking Advanced Microwave Radiometer (AMR), which provided measurements of brightness temperature at 18.7, 23.8 and 34.0 GHz frequencies with a sampling frequency of 1 Hz. Brightness temperature measurements acquired by other passive microwave radiometers, such as SSM/I (Special Sensor Microwave/Imager) and AMSR-E (Advanced Microwave Scanning Radiometer for EOS; Earth Observing System), have demonstrated good performance for the retrieval of dates associated with ice phenology (Kouraev et al., 2007; Kang et al., 2012; Du et al., 2017) as well as for the retrieval of the ice thickness (Kang et al., 2014) on large lakes in Russia and Canada. Because the Jason AMR footprints are large in diameter (42 km at 18.7 GHz, 35 km at 23.8 GHz and 22 km at 34.0 GHz) (Kouraev et al., 2007), the brightness temperature measurements over rivers are dominated by signals emitted predominantly from land surfaces surrounding the river channels. Hence, we used Jason-2 and Jason-3 AMR measurements only as a source of auxiliary information to develop additional criteria for the adjustment of the radar freeze-up/breakup date retrieval algorithm. During the tandem mission period, the differences in the Ku-band backscatter measurements (Sig0) between Jason-2 and Jason-3 were within 1 dB. The difference in the brightness temperature measurements was within 3 K.

In this study, the backscatter and brightness temperature measurements were extracted from a geophysical data record (GDR) product (AVISO+, 2021) distributed on the AVISO+ (Archiving, Validation and Interpretation of Satellite Oceanographic data) data portal (http://ftp-access.aviso.altimetry.fr/, last access: 2 December 2021). The GDR product contains various parameters estimated from the radar return echo, which is represented as a waveform. An altimeter waveform represents a histogram of energy backscattered by the ground surface to the satellite with respect to time.
The waveform ground processing (retracking) is aimed at obtaining a better estimate of range (distance from the satellite to Earth's surface) than those obtained with onboard algorithms. We used the 20 Hz backscatter coefficient retrieved from the ICE1 algorithm (Bamber, 1994). The ICE1 retracking algorithm calculates the centre of gravity, amplitude and width of a rectangular box using the maximum power of the waveform samples. The backscatter is estimated using an instrument link budget and the waveform amplitude (ESA, 2002). The AMR measurements of brightness temperature have a 1 Hz temporal frequency. Here, they were linearly interpolated to the coordinates of the 20 Hz radar measurements.

Landsat 8 images at a 30 m spatial resolution (https://earthexplorer.usgs.gov/, last access: 2 December 2021) were used for the geographical selection of altimetric-radar and radiometer measurements over the river channels using our own Python code, allowing us to overlap the along-track Jason-2 and Jason-3 measurements on the Landsat images. The cross section of an altimetric track with a river channel is called a virtual station (VS; see Fig. 1 for locations). Each VS was given a name containing the track number. To distinguish the VSs located on secondary branches from those located on the main Ob River channel, the names of the stations located on the secondary branch (i.e. small Ob River) were provided with the index “SOb”.

## 4 Temporal variability of radar altimetry signal over frozen rivers

Because there is an absence of work on the use of radar altimetry for river ice studies, we provide some elements of interpretation on the evolution of backscatter with freshwater ice properties. The effect of varying ice properties on the radar return from nadir-looking altimeter instruments has been reported in several investigations of the Antarctic and Greenlandic ice sheets (e.g. Rémy et al., 2012; Nilsson et al., 2015; Larue et al., 2021) as well as for sea ice (e.g. Kurtz et al., 2014; Guerreiro et al., 2016). However, freshwater ice differs from other types of ice in its textural/structural (e.g. crystal structure, density, layering and air bubble content) or chemical properties (e.g. salt content). The impact of freshwater ice properties on satellite measurements in the microwave region has been studied extensively in context of the interpretation of measurements from side-looking imaging SAR sensors (e.g. Duguay et al., 2002; Jeffries et al., 2005; Unterschultz et al., 2009; Atwood et al., 2015; van der Sanden et al., 2021). To date, no research has provided explanations on the temporal variability of radar altimetry measurements over frozen rivers. Below, we provide a summary of the variability of altimetric-radar signals over river ice during the freeze-up (ice formation), ice growth and breakup (ice decay) periods, based on our own interpretation of the altimetric data from this study.

Previous studies (e.g. Kouraev et al., 2005; Duguay et al., 2018; Zakharova et al., 2019, 2020) have demonstrated that over Arctic rivers and lakes, the backscatter from radar altimeters varies seasonally. This behaviour is strongly related to hydrological phases, especially to the presence of ice (Fig. 3a). We analysed the seasonal behaviour of backscatter from the Ob River with respect to geophysical parameters that can potentially affect the scattering properties of the surface at microwave frequencies. The radar return of nadir-looking altimeters significantly differs from that obtained with side-looking imaging SAR instruments (ca. 20–45° off nadir). At nadir, a smooth surface produces a higher return echo than a rough surface (Ulaby et al., 1986). The altimetric return signal is represented as a waveform (Fig. 3b). The shape of the waveform varies for different types of surfaces and can provide a variety of information that is useful for the interpretation of geophysical processes (Berry et al., 2005). For instance, Beckers et al. (2017) found the presence of an intermediate peak on the leading edge of waveforms that they attributed to the reflection of the radar signal from the ice surface (i.e. snow–ice interface). The main waveform peak was considered to originate from the ice bottom (i.e. ice–water interface). The difference in the ice surface–ice bottom altimetric heights estimated from these two peaks was found to be related to the ice thickness and therefore was used for its retrieval. This intermediate peak can also be detected on the leading edge of many waveforms extracted over river ice (Fig. 3b). However, considering that the radar echoes over rivers come from substantially heterogeneous surfaces, we avoided regarding this as a peak produced by the definitive reflecting boundary and suggest another approach for ice thickness retrievals. Figure 3b demonstrates that during ice growth, the primary waveform evolution is related to a decrease in power of the main peak. This can be explained by the volumetric scattering/absorption of the radar signal with ice growth. We noted that this decrease is proportional to the value of backscatter, which can be seen as the integral under the waveform. Based on this observation, we suggest that the retrieval of ice thickness is possible through the establishment of a simple statistical relation between the radar backscatter and in situ river ice thickness measurements. This simple approach has exhibited some potential for the retrieval of river ice thickness from imaging SAR instruments (Unterschultz et al., 2009; Mermoz et al., 2013). However, in contrast to nadir-looking radar altimeters, it has recently been suggested that backscatter from off-nadir looking imaging SAR during the ice growth season may be more related to increased roughness at the ice–water interface than ice thickness alone (Murfitt and Duguay, 2021).

High backscatter values of nadir-looking altimeters are observed when the instrument footprint contains a large fraction of open calm water with low roughness (Fu and Cazenave, 1991). In general, the surface roughness of inland water objects of small and medium size is low, and backscatter values are high, when compared with ocean conditions (Zakharova
<figure><img src="image_6.png" /></figure>

**Figure 3.** Variability of Jason-2 and Jason-3 measurements over sections of the Ob River: (a) interannual variability of backscatter and brightness temperature to the north (VS 161) and south (VS 12) of the study area (dark-red lines correspond to the ice cover period retrieved from altimetric measurements using the manual routine), (b) winter evolution of typical waveforms for track 88 over a river channel (coloured lines correspond to different dates during winter 2013–2014), and (c) typical seasonal variability of backscatter on the Ob River and related ice regime phases.

et al., 2014, 2018). However, on rivers, over large, flooded areas, the water surface can exhibit substantial roughness owing to turbulent flow or the impact of wind, thus resulting in lower backscatter values than those obtained over calm water. On lakes during freeze-up, the presence of floating ice (floes or slush) prevents the development of wind-induced waves and therefore increases the specular reflection off the water (Kouraev et al., 2007, 2021). Freeze-up on rivers starts with the formation of a fine skim ice layer along the banks (i.e. border ice). This ice can detach and drift (Hiks, 2009). The skim ice is characterised by a smooth surface and bottom. We suggest that the peak on the backscatter time series representing the time of the first appearance of the skim or drifting ice in late autumn can indicate the start of freezing (Fig. 3c).

The scattering and absorption of radar echoes within ice depends on the penetration depth of the incident radar signal and internal ice properties (Rémy et al., 2012; Jeffries

et al., 2012). At the Ku-band, the penetration depth into dry freshwater ice is of the order of 5–12 m and depends on the temperature and properties of the ice (Legrésy and Rémy, 1997). Therefore, in the case of thin lake or river ice, the penetration depth for waves emitted at altimetric-radar frequencies (C-band and Ku-band) allows for the estimation of ice thickness. River ice grows as water at the bottom of the ice cover freezes (Duguay et al., 2015). As ice grows, the volume scattering of the radar altimeter signal within the ice increases, resulting in a decrease in backscatter (Brown, 1977; Rémy et al., 2012). A distinct recession limb can be clearly observed on the Ob River backscatter time series in winter (see the red lines in Fig. 3a). The Ob River ice gains approximately 30% of its total thickness during the first month of ice growth (Fig. 2a). The highest temporal decrease in backscatter (Sig)), expressed as $\Delta Sig_0/\Delta t$ (where $\Delta t$ is the period between two consecutive observations), is observed during that period. During winter, small peaks can appear in
the backscatter time series (see VS 12 in Fig. 3a). The peaks can be explained as the temporary changes in the ice surface characteristics, such as the formation of polynya (open-water areas), redistribution of snow on the ice surface by winds, snow wetting during mechanical ice cracking (i.e. slushing) and occasional snowmelt during warm sunny days (and re-freezing at night) in early spring. All these processes either change the roughness of the reflecting surface (polynya and snow crust) or decrease the radar penetration depth and consequently reduce the signal energy loss within the ice.

River ice breakup is influenced by both thermodynamic and hydrodynamic processes known as thermal and mechanical breakup, respectively (Ginzburg, 1973). Studies using SARAL/AltiKa (Satellite with Argos and AltiKa; nadir altimeter at the Ka-band) (Kouraev et al., 2015) and Jason-3 (Kouraev et al., 2021) altimeter data have demonstrated that in early spring, when air temperatures are still mostly negative, lake ice undergoes metamorphism under the influence of solar radiation. At that time, a drop in backscatter of approximately 5–10 dB can be observed. A similar process can occur on the Ob River ice during the pre-melt period (Fig. 3c). When air temperatures become positive, snow on the ice surface melts, increasing the surface backscattering of the radar signal. Melt progressively affects the ice, and melt ponds with a smooth water surface can appear on the ice surface. The presence of water on the ice surface and wet snow/ice prevents the penetration of radar signals (Ulaby et al., 1986). During melt, surface scattering becomes dominant in the altimetric return signals. This leads to an increase in the radar altimetry backscatter because the roughness of the surface of the melting snow/ice and melt ponds is low.

Mechanical breakup begins when the water level rises. At this time, water can flood the ice surface due to the earlier breakup in the river headwaters or its tributaries or to water infiltration through cracks in the weakened/fractured ice (Ginzburg, 1973; Hiks, 2009). A high backscatter peak occurs at the beginning of flooding, which can vary significantly from 25 to 50 dB. Such a large range precludes the development of a threshold-based algorithm for the retrieval of breakup dates. As the river becomes ice-free, the backscatter decreases owing to the increased surface roughness induced by wind and turbulence.

During the open-water season, several backscatter peaks are frequently observed. Summer variability in backscatter depends on many factors, including wind-induced roughness, VS location (banks, presence of islands and floodplain characteristics) and the proportion of water within the footprint (intermittent summer inundation due to rain events), among others.

# 5 Methods

## 5.1 Ice onset and breakup algorithm

Considering the behaviour of radar backscatter over frozen rivers described in Sect. 4 and presented in Fig. 3c, we developed a dynamic algorithm for retrieving ice phenology dates based on the analysis of the Sig0 time series specific to each VS. We suggest that the last peak of each year in the backscatter time series corresponds to the beginning of river ice formation (ice onset). In the case of a multi-peak recession limb (see for example VS 12 in 2013; Fig. 3a), we selected the autumnal peak whose height is of the order of spring-summer peak heights typical for this VS. This helps to distinguish the peaks related to the first appearance of ice (higher peak) from the subsequent peaks related to the appearances of water on the ice (smaller peaks). If the peak selection is not straightforward (for example, two high peaks within 1 month or the prominence of the peak is low), an additional criterion based on the brightness temperature difference of the 34.0 and 18.7 GHz frequency channels (ΔTB) is introduced (Kouraev et al., 2005, 2007). Because the ΔTB values vary around zero and do not exceed 2 K (Fig. 3c), we selected the backscatter peak at time *t* if at least three of the four ΔTB values were < 2 K during the satellite cycle time frame of (*t* − 1, *t* + 2).

Brightness temperature measurements integrate emissions from larger surrounding areas than the altimetric-radar backscatter measurements. Small oxbow lakes on the floodplain as well as soils and bogs on banks usually freeze earlier than the big channels of the Ob River. By applying the (*t* − 1, *t* + 2) window, we ensured that the progression of freezing in the area of the VS was captured and that the backscatter peak was not caused by a synoptic-scale cooling episode or by calm-weather conditions.

The beginning of the ice cover decay (beginning of thermal ice melt) leads to an increase in spring backscatter. The melt detection algorithm searches for the spring peak in the backscatter time series. For multi-peak winters, the algorithm also uses the ΔTB condition. During the melt, the ΔTB is usually higher than 2 K (Fig. 3c). The algorithm searches for the peak, which is accompanied by a simultaneous increase in ΔTB in the same order of magnitude as the mean summer ΔTB value for a given VS. Among the spring peaks, the peak satisfying the condition min[ΔTB<sub>*t*</sub> − mean(ΔTB<sub>Jul-Aug</sub>)] (where mean(ΔTB<sub>Jul-Aug</sub>) is the mean summer brightness temperature difference) was selected. In a few instances, the spring peak was absent or could not be automatically detected because of a low prominence. In these cases, we used the date of the maximum increase in backscatter between two satellite cycles (ΔSig<sub>0</sub>/Δ*t*) for the period from January to mid-June.

A variety of combinations of different geomorphological (riverbanks, floodplain, river width and islands), meteorological (synoptic cooling/warming episodes) and ice
<figure><img src="image_7.png" /><figcaption>Figure 4. Processing scheme of retrieving ice phenology dates from altimetric measurements.</figcaption></figure>

cover/open-water (polynya and ridging) conditions can occur. Their complex impact on the variability in backscatter during freeze-up and breakup makes it difficult to address all the variations automatically. Therefore, we retrieved the ice phenology dates manually through a visual inspection of the backscatter (and ΔTB if necessary) time series for each VS to compare the performance of the automated and manual freeze/melt detection routines. Both the manual and automated routines use the same criteria. A summary of the processing scheme is illustrated in Fig. 4.

## 5.2 Ice thickness algorithm

Interannual variations in the backscatter values at the beginning of the freeze-up period (see Fig. 3a) may be caused by different land, water and ice proportions within the radar altimeter footprint, wind conditions, floating ice concentration or other factors. Assuming that the decrease in backscatter between two consecutive observations (ΔSig₀/Δt) is proportional to a gain in ice thickness (see Sect. 4), we used a relative backscatter decrease rather than the absolute backscatter values for the estimation of ice thickness. This allowed us to minimise the effect of the initial freezing conditions on the accuracy of ice thickness retrievals. For each year starting from the first date of freezing (i.e. ice onset determined manually), we estimated the backscatter cumulative difference Σ(ΔSig₀/Δt) and constructed the relationship between this parameter and in situ ice thickness (Hice) measured at the nearest gauging station (Fig. 5). The application of a locally estimated scatterplot smoothing filter (LOESS) on the Σ(ΔSig₀/Δt) parameter makes it possible to minimise the effect of secondary peaks on the winter backscatter curve.

Along the 400 km long lower Ob River reaches covered by the 20 northernmost Jason satellite tracks, 48 VSs were defined using the Jason-2 and Jason-3 tracks overlain on Landsat 8 images; 10 VSs nearest to the five gauging stations were chosen as a training set for calibrating the ice thickness retrieval algorithms and for estimating the uncertainties in the ice phenology dates and ice thickness retrievals. The results of ice onset and ice melt start date retrievals by both the automated and manual routines were compared with the ground

station records of ice types or water/ice cover state (ridging, polynya, water-on-ice and ice drift among others). Because the dates of the in situ measurements did not coincide exactly with the Jason overpass dates, the ice thickness values were linearly interpolated between the two adjacent dates of in situ observations for the dates of the satellite overpasses. The other 38 VSs (main set) were used to characterise the ice conditions within all the studied river reaches.

Among the tested fitting functions (linear, polynomial and power), the power equation (Eq. 1) produced the best fit between the cumulative backscatter difference and in situ ice thickness measurements. The selection of the fitting function was based on maximisation of the correlation between Σ(ΔSig₀/Δt) and in situ Hice and the minimisation of root mean square error (RMSE) calculated between the Hice retrieved from Eq. (1) and the observed values.

$$H_{\text{ice}} = a \times \text{abs}(\Sigma (\Delta \text{Sig}_0 / \Delta t))^b \quad (1)$$

Using the leave-1-year-out method (Picard and Cook, 1984) for each gauging station–VS pair in the training set, we obtained a set of a and b coefficients and estimated their mean values. These mean values were then used for ice thickness estimation on the corresponding training VS.

The algorithm was validated for the 10 virtual stations from the training set. For the ice phenology algorithm, the accuracy was evaluated using the percentage of good retrievals. Because the Jason satellite has a 10 d repeating orbit, a retrieval was considered good if the difference between it and the in situ observation was ≤ 10 d. The accuracy of the ice thickness retrievals was evaluated for the same 10 virtual stations using the correlation coefficient and RMSE calculated between the retrieved and observed ice thickness for the 2008–2018 period. A schematic representation of the processing steps is presented in Fig. 6.

Using coefficients a and b from Eq. (1) derived for the VSs from the training set, we estimated the ice thicknesses for the locations of the other 38 VSs. The following strategy was adopted. Firstly, for each VS, using the correlation matrix, we searched for the best correlation between its backscatter and the backscatter for one of the training VSs. Then, the coefficients a and b of VSs with the highest correlation coefficient were applied to the VS considered.

## 5.3 Generation of 2D spatiotemporal ice thickness product

The ice thickness retrieved at all 48 VSs was used for the creation of weekly maps, which were generalised into a 2D spatiotemporal ice thickness product. For this, the altimeter-derived ice thicknesses were interpolated and smoothed in 2D spatiotemporal coordinates using a moving average filter. The size of the smoothing window must (a) preserve the Hice magnitudes and the spatial heterogeneity of ice thickness in the spatial domain as much as possible and (b) reduce
<figure><img src="image_8.png" /><figcaption>**Figure 5.** Relationship between $\Sigma(\Delta\text{Sig0}/\Delta t)$ and in situ ice thickness at 10 virtual stations.</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Processing scheme of ice thickness retrieval from altimetric measurements.</figcaption></figure>

the residual noise in the temporal domain remaining after the backscatter time series was smoothed with the LOESS filter. The selection of the smoothing window was determined based on the correlation coefficient and RMSE estimated between the $H_{\text{ice}}$ time series retrieved at virtual stations and the $H_{\text{ice}}$ time series extracted from the smoothed product for the same VS locations. The effect of interpolation and smoothing on seasonal $H_{\text{ice}}$ magnitude was evaluated by estimating the differences between the maximum yearly values of $H_{\text{ice}}$.

Notably, a critical effect of the interpolation/smoothing procedure on $H_{\text{ice}}$ (i.e. rapid degradation of the correlation coefficient and RMSE statistics) occurred in the 15–40 km spatial-window range. In windows beyond 40–45 km,
<figure><img src="image_10.png" /><figcaption>**Figure 7.** Normalised probability density function (norm pdf) of difference in number of days between altimetric and in situ observed dates of (a) ice onset, (b) melt start and (c) melt end over the 2008–2019 period retrieved for 10 virtual stations from a training set of VSs using manual and automated routines.</figcaption></figure>

the procedure resulted in few subsequent changes in the smoothed product. In the temporal domain, the application of windows below 15 d (for 0 km spatial window) and 40 d (for 60 km spatial window) also had a minor effect (the correlation coefficient was above 0.98; RMSE and difference in maximum $H_{\text{ice}}$ ($H_{\text{icemax}}$) were below 0.03 m). The selected window of 40 km every 30 d allowed keeping the average $H_{\text{icemax}}$ difference and RMSE below 0.04 and 0.03 m, respectively. The correlation between the $H_{\text{ice}}$ time series extracted from the smoothed product and those retrieved at 20 VSs located on the main branch of the Ob River was 0.99.

# 6 Results

## 6.1 Ice phenology retrievals

Considering the 10 d repeat overpass of the Jason-2 and Jason-3 satellites and the distance between the gauging stations and VSs, we considered a 10 d time-step difference (i.e. ±10 d) as an acceptable accuracy for the altimetry-derived ice phenology dates. For 90 % of the manual-routine retrievals, the difference between the first ice events (ice onset) and in situ observations was less than 10 d (Fig. 7a). In 56 % of the cases, there was no difference in the number of days. The automated routine missed certain temporal Sig0 behaviour types, and detection was less accurate for ice onset than in the manual routine. Only 70 % of the altimetric ice onset dates were within 10 d of the in situ observations at the gauges, and only 40 % were on the same day.

Breakup is a more complex process consisting of the thermal degradation of ice cover (melt start) and its mechanical disintegration and downstream movement (melt end). Comparing the dates of the altimetry-derived melt onset (start) with those of the ice state flags provided by gauging stations, we could conclude that the manual routine of our algorithm accurately detected the start of ice thermal degradation. In 88 % of the cases, the difference between the manually retrieved melt dates and in situ observations of the first water appearance was less than ±10 d (Fig. 7b). The automatically derived melt date estimations were less accurate for the de-

tection of the melt start than the manually derived estimations. However, the automated routine was better adapted for the detection of the melt end than for the detection of melt start; an accuracy of less than ±10 d was achieved for 67 % of cases (Fig. 7b and c).

The melt algorithm was designed for the detection of the ice melt start (i.e. the wet-ice, water-on-ice or open-water events). The manual retrieval of breakup dates allowed for better control over the complex variability in the backscatter in spring than the automated routine. It is likely that in the complex cases the automated routine detects the melt end or even provides unrealistic early/late estimates of melt dates. For example, the automated routine detected unrealistic melting dates before 10 April and after 10 June in 10 % of the total cases for the entire set of 48 VSs over the 10 winter periods of study. Because the manual routine demonstrates better accuracy than the automated routine, it was therefore selected for further analysis of the results and for use with the ice thickness retrieval algorithm.

A comparison of the interannual variability of median dates estimated for VSs located on the main river branch (20 VSs in total) with the corresponding parameter estimated from the observations at four gauging stations (also located on the main river branch) demonstrated a good agreement between the satellite retrievals and in situ observations (Fig. 8). A significant difference between the gauges and VSs (approximately 20 d) was observed only for the melt start dates in 2014. This good agreement suggests that our approach is suitable for monitoring the interannual variability of ice events.

## 6.2 Ice thickness retrievals

The accuracy of the ice thickness retrievals from altimetric measurements was assessed for the 10 VSs from the training set located near the gauging stations (Fig. 9). At the northern VSs (VS 161, 138, 237 SOb and 240 SOb), the correlation between the retrieved and observed ice thickness was stronger than at other locations, and the RMSE in the estimates of the altimeter-retrieved ice thickness was less than 0.12 m (Table 2). For the southernmost VSs (VS 12 and 109),
<figure><img src="image_11.png" /><figcaption>**Figure 8.** Interannual variability of dates of (a) ice onset and (b) melt start for the main Ob River channel retrieved using the manual routine. The statistics are estimated using data from 4 gauging stations and 20 virtual stations located along the main Ob River channel.</figcaption></figure>

<table><caption>**Table 2.** Coefficients *a* and *b* for the built relations, correlation coefficient (*R*) and RMSE between the retrieved and in situ ice thickness for virtual stations from the training set. SOb refers to virtual stations located on the secondary (western) river branch of the Ob (see Fig. 1).</caption><thead><tr><th rowspan="2">Virtual stations</th><th rowspan="2">Corresponding gauging station</th><th rowspan="2">Distance to gauging station (km)</th><th colspan="2">Coefficients of Eq. (1)</th><th colspan="2">Validation scores</th></tr><tr><th>a</th><th>b</th><th>R</th><th>RMSE (m)</th></tr></thead><tbody><tr><td>161</td><td>Pitlar</td><td>5.6</td><td>8.69</td><td>0.39</td><td>0.94</td><td>0.07</td></tr><tr><td>138</td><td>Pitlar</td><td>5.1</td><td>6.54</td><td>0.42</td><td>0.94</td><td>0.07</td></tr><tr><td>237 SOb</td><td>Muzhi</td><td>15.9</td><td>7.64</td><td>0.42</td><td>0.90</td><td>0.10</td></tr><tr><td>240 SOb</td><td>Muzhi</td><td>7.7</td><td>7.96</td><td>0.41</td><td>0.90</td><td>0.10</td></tr><tr><td>240</td><td>Gorki</td><td>19.4</td><td>7.70</td><td>0.39</td><td>0.81</td><td>0.12</td></tr><tr><td>135</td><td>Gorki</td><td>28</td><td>6.88</td><td>0.42</td><td>0.87</td><td>0.11</td></tr><tr><td>164</td><td>Kazym Mys</td><td>22.4</td><td>8.83</td><td>0.35</td><td>0.84</td><td>0.10</td></tr><tr><td>211</td><td>Kazym Mys</td><td>20.1</td><td>10.7</td><td>0.31</td><td>0.76</td><td>0.12</td></tr><tr><td>12</td><td>Polnovat</td><td>10.9</td><td>8.23</td><td>0.41</td><td>0.77</td><td>0.18</td></tr><tr><td>109</td><td>Polnovat</td><td>33.4</td><td>2.92</td><td>0.55</td><td>0.84</td><td>0.14</td></tr></tbody></table>

the RMSE increased to 0.14–0.18 m. The relative uncertainty in the ice thickness estimates was found to be higher at the beginning (low $H_{ice}$) of the ice period than in the middle or at the end of winter. Except for VS 109, the variability in the values of coefficients *a* and *b* in Eq. (1) was low (Table 2), which indicates good stability in the established relationships and their potential validity for other VSs located far from the gauged reaches. One way to evaluate the sensitivity of the satellite $H_{ice}$ retrievals to the fitting parameters consists of running a cross-validation experiment (i.e. retrieving $H_{ice}$ using the coefficients obtained for an adjacent VS). Similar scores between retrieved and in situ ice thickness obtained in the validation (Table 2) and cross-validation (Table 3) experiments demonstrated the robustness of Eq. (1), resulting in low errors in the satellite $H_{ice}$ estimates for the northern VSs (VS 161, 138, 237 SOb and 240 SOb). Uncertainties in the $H_{ice}$ retrievals for the southernmost VSs were higher than for the northern VSs. However, the retrievals for the southern VSs in the cross-validation experiment can be improved by

selecting the equation parameters (*a* and *b*) not from the adjacent VS but from the training VSs, which provide the best Sig0 correlation. For example, when applying the coefficient derived for the VS 135–Gorki gauging station pair to VS 109 and VS 12 (backscatter of VS 135 demonstrated the highest correlation with the backscatter of VS 109 and VS 12), the RMSE of the retrieved $H_{ice}$ for these VSs decreased from 0.23 to 0.18–0.19 m (see scores in the denominator of corresponding lines in Table 3).

## 6.3 Ice thickness estimation for the entire studied river reaches

Using the *a* and *b* coefficients of Eq. (1) developed for the VSs from the training set, we estimated the ice thickness at locations of the other 38 VSs and produced weekly maps, which were generalised into a 2D spatiotemporal ice thickness product (Fig. 10). The elaborated maps can be used for evaluating the ice thickness and ice phenology dates in the
<figure><img src="image_12.png" /><figcaption>**Figure 9.** Ice thickness measured at gauging stations (x axis) and retrieved from the altimetric measurements (y axis) for the training VSs. The red dashed line corresponds to the 1 : 1 relation line.</figcaption></figure>

areas between the virtual stations. For instance, two useful
parameters can be extracted from the 2D product: the maxi-
mum ice thickness and the ice thickness observed on 1 De-
cember.

Generally, knowledge of the maximum ice thickness is rel-
evant for hydro-climate change monitoring (Vuglinsky and
Valatin, 2018); the ice thickness determined on 1 December
is crucial for the local and regional socioeconomic stakehold-
ers because this is the average date that the ice bridge road
to the north of the study area at Salekhard opens. To assess
the quality of the 2D product for operational use, we com-
pared the interannual dynamics of the mentioned parameters
derived from this product with those observed at the gauging
stations.

For approximately 90% of the area of the studied river reaches, the interannual variability in the maximum ice thickness retrieved from altimetric measurements indicated a clear decrease from 2008 to 2012 (Fig. 10). This tendency corresponds well to those observed at all the gauging stations (Fig. 11a). Since 2013, the maximum ice thickness has increased slowly. However, both altimetric and in situ observations after 2013 exhibited spatiotemporal variability that was not always in agreement. This disagreement may be related to environmental factors affecting ice growth, such as snow amount, autumn ice drift and accumulation, ridging/hummocking, and ice flooding (water on ice). For example, according to the station records, ridging events appear more frequently after 2012 at the northernmost gauging station, Pitlar, than at the other gauging stations. We could not fully exclude the effect of the simplicity of the retrieval algorithm, based on the empirical approach, as well as the effect of the spatiotemporal smoothing of the altimetric retrievals used in the map production.

Differences between altimetric and in situ ice thicknesses on 1 December lie within the algorithm uncertainties of 0.07–0.18 m (Fig. 11b). In addition to the geophysical reasons and simplicity of the algorithm noted above, the observed difference can be related to degradation in the quality of the in situ time series (gaps or unrealistic values) and the low representativeness of the one-hole sampling protocol for in situ ice thickness measurements.

6.4 Winter ice bridge road operation forecast

In many regions with seasonal ice cover, frozen rivers enhance the interconnection and supply of many small and even some large cities. Many remote villages that are linked in summer to supply centres only via expensive aircraft or boat transport can directly access the primary land transport arteries built on frozen ground and lake/river ice. The value of the ice roads is the highest for the Arctic regions where the construction of permanent bridges is restrained by the presence of permafrost and its destabilisation.

A good example is the city of Salekhard located on the Ob River near the polar circle in the zone of discontinuous permafrost. The city has 50 000 inhabitants and is supplied primarily via the Northern Railway, which connects the small town of Labytnangi on the left bank of the Ob River with the European part of Russia and the key supply centres. Merchandise from Labytnangi is delivered to Salekhard by ferry. Every winter, an ice road is constructed to ensure the transport of goods and people. Owing to security reasons, the ferry ceases operation after the appearance of the first ice. The ice road construction (ice artificially grown via pumping water onto the ice surface) begins when the ice thickness attains a minimum allowed value of 0.20–0.25 m (Instructions on Safety Organisation, 1969). The traffic of light vehicles is authorised when the ice thickness exceeds 0.30 m. According to information kindly provided by the State Traffic Service
<figure><img src="image_13.png" /><figcaption>**Figure 10.** Spatiotemporal ice thickness variability (m) for the main branch of all the lower Ob River reaches for the 2008–2018 period from the generalised weekly altimetric product. Distance (km) is indicated from the northernmost virtual station 187 in the upstream direction (south). Letters on the *x* axis correspond to the first letter of the month (December–April).</figcaption></figure>

<table><caption>**Table 3.** Correlation coefficient (*R*) and RMSE between the retrieved and in situ ice thickness in the cross-validation experiment.</caption><thead><tr><th>Virtual stations</th><th>Corresponding gauging station</th><th>VS whose coefficients <i>a</i> and <i>b</i> are used for cross-validation</th><th><i>R</i></th><th>RMSE (m)</th></tr></thead><tbody><tr><td>161</td><td>Pitlar</td><td>138</td><td>0.94</td><td>0.09</td></tr><tr><td>138</td><td>Pitlar</td><td>161</td><td>0.94</td><td>0.09</td></tr><tr><td>237 SOb</td><td>Muzhi</td><td>240 SOb</td><td>0.90</td><td>0.10</td></tr><tr><td>240 SOb</td><td>Muzhi</td><td>237 SOb</td><td>0.89</td><td>0.11</td></tr><tr><td>240</td><td>Gorki</td><td>135</td><td>0.81</td><td>0.13</td></tr><tr><td>135</td><td>Gorki</td><td>240</td><td>0.87</td><td>0.11</td></tr><tr><td>164</td><td>Kazym Mys</td><td>211</td><td>0.84</td><td>0.10</td></tr><tr><td>211</td><td>Kazym Mys</td><td>164</td><td>0.76</td><td>0.13</td></tr><tr><td>12</td><td>Polnovat</td><td>109/135*</td><td>0.76/0.76*</td><td>0.23/0.18*</td></tr><tr><td>109</td><td>Polnovat</td><td>12/135*</td><td>0.76/0.76*</td><td>0.23/0.19*</td></tr></tbody></table>

of the Yamalo-Nenetsky Autonomous District (Russia), the ice road operation usually starts 3–4 weeks after the beginning of freezing, with an average opening of 30 November to 1 December. The road operation closes gradually, starting by limiting the lorry load in the middle of April until fully

halted at the beginning of May. The ferry connection is restored approximately 3 weeks later. Between the ferry and the ice road operation, the connection is ensured via a hovercraft boat only for a limited number of passengers or for emergencies.

<footnote>* Two different pairs of *a* and *b* coefficients were used for cross-validation.</footnote>
<figure><img src="image_14.png" /><figcaption>**Figure 11.** Interannual variability in the ice thickness extracted from 2D ice thickness product at the Pitlar, Gorki and Polnovat stations and the values observed at these gauging stations: (a) maximum ice thickness and (b) ice thickness on 1 December.</figcaption></figure>

The dates of the autumnal halt of ferry operation for 2010–2018 agreed remarkably well with those of the first ice occurrence on the four northernmost tracks of the Jason satellite, located 65–75 km south of the city (Fig. 12a). For the short-term forecast, the satellite observations in the area of VS 112 and 187 are especially good and allow for the prediction of the end of ferry operations with a 5 d error expressed as RMSE (Table 4). We further assumed that the average dates when the satellite-derived $H_{ice}$ at the four northernmost VSs reaches 0.30 m may provide an estimate for the road-opening date. These dates were extracted from the spatiotemporally smoothed maps and compared with the dates that the ice road was opened, as provided by the State Traffic Service. The mean absolute difference between the satellite and observed dates was 3 d and the RMSE was 5 d. Despite the low mean error, the quality of the forecast was not satisfactory. The predicted dates were consistently earlier than the actual dates of the road opening (Fig. 12b). Moreover, for half of the years, the predictions differed from the observations by more than 5 d (11 d maximum), and the correlation between the predicted and observed dates was not significant.

For predicting the dates when the ice road at Salekhard ceases its operation, the ice road was closed before the predicted dates of the breakup onset detected by the altimeter; therefore, the use of the northernmost VSs was not possible. However, the satellite observations of the breakup onset at the VSs located in the upper reaches can provide the necessary information. Using the altimetric retrievals of the breakup start for the entire set of 48 VSs for each year, we searched the date when at least two altimetric breakup onsets (AMO2) were detected within the 400 km river reaches. This date served as a predictor of the date of the ice road closure at Salekhard. The correlation between the AMO2 and

<figure><img src="image_15.png" /><figcaption>**Figure 12.** Observed and predicted from the altimetric satellite observations dates of the Salekhard (a) ferry closing and ice road (b) opening and (c) closing.</figcaption></figure>

<table><caption>**Table 4.** Scores for the evaluation of accuracy of the prediction of ferry closing and ice road opening and closing dates in the city of Salekhard. DOY: day of the year.</caption><thead><tr><th>Statistic</th><th>Ferry closing</th><th>Ice road opening</th><th>Ice road closing</th></tr></thead><tbody><tr><td>Mean (DOY)*</td><td>309/305*</td><td>330/332*</td><td>123/122*</td></tr><tr><td>Standard deviation (days)*</td><td>8/8*</td><td>5/5*</td><td>4/6*</td></tr><tr><td>RMSE (days)</td><td>5</td><td>5</td><td>3</td></tr><tr><td>Mean absolute difference (days)</td><td>-4</td><td>3</td><td>0</td></tr><tr><td>Correlation coefficient</td><td>0.83</td><td>0.46</td><td>0.70</td></tr><tr><td><i>p</i> value</td><td>0.005</td><td>0.185</td><td>0.023</td></tr></tbody></table>

observations was significant ($p$ value of 0.025) with a correlation coefficient of 0.70 (Fig. 13). Using the linear equation computed from the relationship provided in Fig. 13, one can obtain a forecast date of the ice road closure at Salekhard (Fig. 12c). The uncertainty in the obtained forecast, expressed as RMSE, was 3 d (Table 4).

# 7 Discussion

## 7.1 Factors affecting the altimetric backscatter signal

Various factors can affect the radar return echoes and consequently the accuracy of the river ice thickness retrievals. One source of uncertainty could originate from the underestimation of the role of snow on the ice in the microwave signal scattering (King et al., 2013, 2015). However, Willatt et al. (2011) demonstrated that the Ku-band electromagnetic wave scattering caused by snow at nadir is low, and in our study, we neglected the presence of snow on ice. Using the

<footnote>* Mean and standard deviation are given for predicted (numerator) and observed (denominator) dates.</footnote>
<figure><img src="image_16.png" /><figcaption>**Figure 13.** Relation between the altimetric breakup onset dates (AMO2) and those observed for the Salekhard ice road closure dates expressed in the day of the year (DOY).</figcaption></figure>

precipitation data from the nearest meteorological station, we noted that not all heavy snow accumulation episodes affected the backscatter over the river ice. In several cases, snowfall resulted in backscatter changes of approximately 1.5 dB. The smoothing procedure applied to the cumulative $\Delta\text{Sig}_0/\Delta t$ series helped to eliminate this effect. Moreover, after adding the in situ snow depth to the ice thickness, we noted that the $H_{\text{ice}} - \Sigma(\Delta\text{Sig}_0/\Delta t)$ power relationship becomes weaker.

Another factor potentially affecting the backscatter value over freshwater ice is the ice roughness at the ice-water interface (Atwood et al., 2015). The roughness of the ice bottom on the rivers is expected to be high at the beginning of the freeze-up period in bridging areas, where floes juxtapose and accumulate underneath. Any rough boundary dissipates the signal of nadir-looking radar instruments, resulting in a decrease in backscatter. Further congelation of the inter-floe volume and ice growth both lead to the levelling of the ice low boundary. We suggest that the under-ice current may also contribute to the ice bottom levelling at the Ku-band radar wavelength scales (~1–1.5 cm). This levelling would increase the radar return power during winter. However, we did not find evidence of this process in the backscatter time series (see Fig. 3a). Either the levelling effect is weak and is masked by high volumetric scattering of the radar echo within the thickening ice, or at the location of our VSs, the ice juxtaposition (and consequently the ice-bottom roughness) is inconsequential. Future investigations with dedicated in situ observations on river ice texture evolution are necessary to understand the effect of bottom roughness on the radar altimeter return signals.

Ice internal layering is also key factor in the scattering of radar signals (Legrésy and Rémy, 1997; Nilsson et al., 2015; Slater et al., 2019). Under the climate conditions of north-western Siberia, ice layering (characterised by dense reflective icy surfaces) is rare, as the air temperature of winter warming episodes never approaches the melting point. Daily positive temperatures lasting several hours can occur starting from the end of March in the southern part of the study area

and 1–2 weeks later in the northern part. During this time of the year, the ice is well developed and almost reaches its maximum thickness. Layering can also occur after river water floods the ice surface through cracks. According to the in situ observations at the gauging stations, this phenomenon was observed in the last several years at the end of the ice season in the southern part of the study area. Both warming episodes and flooding events lead to a backscatter increase of approximately 1–5 dB and make the altimetric ice thickness retrievals difficult by the end of the ice season. The highest underestimation of $H_{\text{ice}}$ (0.15–0.20 m) was observed in such cases.

The internal ice structure can also affect the backscatter value, for example, via air bubble inclusion (Gunn et al., 2018). During ice formation, jamming and ridging can occur in Arctic rivers, resulting in the formation of air pockets. Ridging is rare near the gauging stations on the Ob River. However, there is no information about the state of the ice in other ungauged reaches, including the areas covered by the VSs. We could only speculate that ridging/hummocking might be one of the reasons for the high difference in the coefficients of Eq. (1) determined for VS 109 and for the lower accuracy of the $H_{\text{ice}}$ retrievals at this station compared to the northern virtual stations. On Landsat 8 images acquired for two springs in 25 April 2019 and 1 May 2015 (not shown), the irregular spatial ice structure in the area of VS 109 indirectly confirms our hypothesis. More studies involving the simultaneous analysis of SAR imagery exhibiting ice field irregularities (Unterschultz et al., 2009) and altimetric signals could help to clarify this issue.

## 7.2 Drawbacks and potential improvement of algorithms

Manual retrievals of ice phenology provided accurate estimates of the dates of the first appearance of ice (with 90% accuracy considering a 10 d difference) and the dates of the beginning of ice melt (with 88% accuracy considering a 10 d difference) on the Ob River. This accuracy is comparable with errors of 0–43 d and 66%–97% and 76%–99% accuracy obtained with other satellite sensors and reported in a summary review dedicated to lake ice (Murfitt and Duguay, 2021). The automated retrievals of the Ob River ice demonstrated lower accuracy than the manual retrievals with several possible explanations. Because the radar footprint over rivers is heterogeneous and is affected by signals from the frozen/unfrozen state of land, river and floodplain lakes, there are numerous variations in the behaviour of backscatter, especially at the beginning and end of winter. The automated routine misses some temporal Sig0 behaviour types, and detection is less accurate than that from the visual analysis of the Sig0 time series. The high difference between the manual and automated retrievals for the ice melt period is likely related to the approach used for adjusting the spring peak selection based on mean summer $\Delta$TB. In the man-
ual selection of peaks, the ΔTB criteria were used in fewer
cases than in the automated selection (see Sect. 5.1 for ex-
plications). It is likely that the ΔTB criteria were not valid
for all virtual stations because their configuration (presence
of oxbow lakes, fens, secondary branches and others) affects
the summer ΔTB values differently. A detailed study of the
reasons for the multi-peak backscatter and brightness tem-
perature variability in different VS configurations could help
to uncover universal criteria and improve the automated rou-
tine. The development of multi-satellite (combination of al-
timetry missions) and multi-sensor approaches (combination
of altimetry with SAR and optical sensors) could help to im-
prove the accuracy of the phenology date retrievals as well as
to enhance the spatiotemporal resolution of the product.

The 0.07–0.18 m accuracy of the altimetry-based ice thick-
ness estimates is also comparable with the accuracy reported
in the rare studies dedicated to satellite retrievals of freshwa-
ter ice thickness. The SAR-based retrievals for three Cana-
dian rivers achieved a 0.09 m accuracy (Mermoz et al., 2013).
For lake ice, the reported accuracy (expressed as RMSE)
reached 0.20 m for microwave retrievals (Kang et al., 2014)
and 0.12–0.40 m for altimetric retrievals (Beckers et al.,
2017).

The inaccurate detection of ice onset can affect the ac-
curacy of the ice thickness estimates, especially at the be-
ginning of freeze-up. Another reason for the uncertainties
is the multi-peak character of the winter backscatter reces-
sion curve and the residual noise that could remain in the
backscatter time series after applying the LOESS filter. For
example, this occurred at VS 12, where a polynya persisted
until March in at least 4 years of the study period, producing
a noisy recession curve (see Fig. 3a) and a high dispersion of
points around the 1 : 1 relation line on the H<sub>icealti</sub> − H<sub>icein situ</sub>
scatterplot (Fig. 9). Several improvements are proposed for
future studies. First, an improvement in the detection accu-
racy for the freeze-up is necessary. In our approach, the ice
thickness estimation starts from the date of the first ice (bank
ice or frazil ice floes) appearance. Usually, the river reaches
in the area of the VSs at this time of the year are not fully
frozen. The detection of the date of the first consolidated ice
(i.e. fully frozen reaches) using other satellite sensors may
help to reduce the errors in the retrievals in a low range of
ice thickness. A second improvement would involve the use
of other parameters of the altimetric-radar waveform instead
of (or as a complement to) the backscatter coefficient. As
shown in Fig. 3b, the main peak of the radar waveform de-
creases as the ice grows. We suggest that the amplitude of the
main peak or the area under this peak may produce stronger
relationships with the in situ ice observations than with the
backscatter coefficient alone. During winter, the main peak
is the most variable part of the waveform. This peak corre-
sponds to the return signals from the ice bottom (Beckers et
al., 2017). The other part of the waveform records the return
signal from the surrounding (off-nadir) areas and also con-
tributes to the total value of the backscatter coefficient. Any

temporal changes in the state of the surrounding areas that
are different from that in the near-nadir area (e.g. ice flood-
ing or off-nadir polynya) lead to backscatter changes that are
not related to ice growth. Consequently, the use of waveform
parameters related to the main peak could reduce the errors in
the ice thickness estimates. Unfortunately, these parameters
are not directly provided in the AVISO+ Jason GDR prod-
uct, but they could potentially be estimated from the initial
waveforms.

7.3 Ice road operation forecasting: problems and perspectives

As noted in Sect. 6.3, at the beginning of freezing, the errors
in the *H*<sub>ice</sub> retrievals are quite high, and the ice thickness is
often overestimated (see Fig. 11b). Therefore and consider-
ing the weak correlation between the observed and predicted
dates, we conclude that the predicting algorithm and the ice
thickness product are currently not sufficiently accurate for
the operational forecasting of the ice road opening. Never-
theless, their accuracy is sufficient for the assessment of cli-
mate impacts, as we captured an overall trend for an earlier
opening of the road demonstrated by both in situ and altimet-
ric time series (Fig. 12b). Moreover, the altimetric product is
accurate enough for the operational prediction of ferry and
ice road closing dates at Salekhard.

The example of the Salekhard ice road case study aimed
to demonstrate the capacity of the satellite product for appli-
cations in the socioeconomic domain. Only the most evident
parameters derived from the product (ice onset dates, date of
first or second consecutive upper reach melt onset, and date
when ice thickness reaches 0.30 m) were tested. It is likely
that other product-derived parameters may have higher pre-
dictive capacity. Their elaboration and testing will be the sub-
ject of future study.

The prediction algorithm for the Salekhard ice road op-
eration dates can benefit from the incorporation of other
satellites. The current altimetric product does not cover the
Salekhard river reach. The northernmost Jason tracks are lo-
cated 65 km south of the ice road. As the melting progresses
from the southern upper reaches, the use of the predictor for
an ice road closing related to these reaches was evident and
successful. However, the freezing on the Ob River progresses
from the northern lower reaches, which are not covered by
the Jason tracks. Our assumption that the ice growth dynam-
ics in area of the northernmost Jason tracks and those in the
area of the ice road are similar may not be fully correct and
may be another reason for the low accuracy in forecasting
the dates of the ice road opening. We expect that predictors
elaborated from satellites (altimetric or other) observing river
reaches north of the city of Salekhard may help to improve
the accuracy in forecasting the ice road opening dates. Other
improvements can be expected from an extension of the cur-
rent product with retrievals from prior altimetric satellites of
the same series: TOPEX/Poseidon (1992–2002) and Jason-
1 (2001–2012). This time-period extension will elongate the time series used for the ice road forecast test and will make it possible to create a robust forecasting algorithm with well-defined margins of applicability. Unfortunately, this was not possible to accomplish in this study due to the short time series of available simultaneous observations.

## 7.4 Potential application of radar altimetry to other rivers

One of the main advantages of satellite altimetry compared with imaging SAR is the relative ease in processing the measurements over the large hemispheric-scale domain. Because the developed algorithm for the ice phenology dates is independent of the availability of in situ observations, it can be directly applied to other rivers. The primary limitation of the algorithm would be the minimum width of the river, which, from our experience, is 100–200 m. The algorithm can be also applied to lakes and bogs because the temporal variability of the backscatter over these surfaces is similar to that described in Sect. 3. However, the algorithm would need an additional adjustment and verification for rivers with unstable winter ice cover, such as those at mid latitudes, as well as in areas of frequent rain-on-snow/ice episodes.

Unfortunately, similar to the proposed SAR-based methods (Unterschultz et al., 2009; Mermoz et al., 2013), the availability of ground truth data is indispensable for our ice thickness retrieval algorithm. Nevertheless, we demonstrated that for the long river reaches, the obtained empirical relations between the in situ data and satellite backscatter measurements are quite similar for different river segments. We suggest that the use of $\Delta\text{Sig}_0/\Delta t$, instead of the absolute value of $\text{Sig}_0$, allows for the minimisation of the perturbing effect of the surrounding banks (i.e. local morphological conditions). The characteristics of the surrounding banks affect the backscatter values at the beginning of the freeze-up period (i.e. ice onset). The effect of the initial ice conditions was thus reduced. The unification of the relations and the application to other reaches and other rivers is possible when using $\Delta\text{Sig}_0/\Delta t$. However, further investigations are necessary to assess the impact of ice textures and types on the parameters of the retrieval equation. Hence, the combination of data from imaging SAR (providing information on ice type) and altimetric satellite missions could be beneficial.

## 8 Conclusions

The low number of in situ observations and their general infeasibility to cover vast areas as well as the decreasing number of observational sites and degradation of the quality of their time series are good arguments for furthering the development of satellite remote sensing methods for freshwater ice monitoring. The present study demonstrates the potential of satellite radar altimetry for monitoring river ice parameters,

such as freeze-up and breakup dates and ice thickness, for a large Arctic river.

An algorithm based on the analysis of temporal changes of backscatter coefficients and the detection of the first peak preceding the winter backscatter recession from the Jason-2 and Jason-3 satellite altimeters provides an estimation of the river ice onset with an accuracy of $\pm 10$ d (corresponding to the 10 d satellite overpass frequency) in 90 % of the cases (manual routine). River ice breakup consists of two phases: thermal degradation and mechanical breakup. By identifying a peak in the backscatter in spring, the algorithm can detect the beginning of thermal degradation well, with the same accepted accuracy of $\pm 10$ d for 88 % of the cases. River ice thickness was retrieved from altimetric backscatter measurements via a simple empirical relationship with in situ observations, using 10 years of simultaneous in situ observations from five stations. The accuracy of the thickness retrievals (expressed as RMSE) ranges from 0.07 to 0.18 m (or 7 %–18 % considering the average ice thickness at the end of ice season).

The spatiotemporal interpolation and smoothing of satellite-derived river ice thickness retrieved for 48 VSs along the 400 km reaches of the lower Ob River allowed for the elaboration of weekly maps generalised in the form of an annual spatiotemporal product. The ice thickness time series can be extracted for any location and used for climate change monitoring and ice road operational purposes. Using the first version of the product, we demonstrated that the uncertainties in predicting the dates of ferry closing and ice road opening and closing near the city of Salekhard (expressed as RMSE) are 3–5 d. Despite these promising results, we perceived that the current (prototype) version of the product is not sufficiently mature for operational use because it overestimates the ice thickness at the beginning of the ice season. The accurate estimation of ice thickness is critical for safety. However, the algorithm and product could be significantly improved in the future through a multi-mission and multi-instrument (optical or SAR imagers) approach. With the use of the Copernicus satellite altimeters Sentinel-3A and Sentinel-3B, the retrieval of the ice thickness can be improved. These satellite missions carry more advanced altimetric SAR instruments with footprints representing a narrow band and return signals that are less contaminated by land than signals from the conventional Jason instrument. Although the nominal repeat frequency of the Sentinel-3 satellites (27 d) is not suitable for operational applications, they provide five overpasses within a 25 km distance around the Salekhard ice road, and thus, the temporal resolution of the observations may be significantly improved. The combination of the data from the Jason and Sentinel-3 missions could be beneficial.

The Salekhard ice road is exceptionally well instrumented, monitored and maintained by local authorities, due to the high demand for its use and high merchandise flow. In other regions, ice roads connecting small cities and villages are
less monitored, and access to operational information is poor. Moreover, many intermittent and unofficial river crossings are developed each year by the local people, and the lack of information on the state of ice often results in accidents and requires intervention by the emergency service. The demonstrated capability of the first version of the altimetric river ice product as a supporting tool for the ice road operations on the Ob River is quite promising. Further product improvements and the development of sophisticated prediction criteria for road operations that could be adapted to other reaches of the Ob River are planned.

**Code availability.** Free libraries used for data extraction and statistical analysis are available at https://www.lfd.uci.edu/~gohlke/pythonlibs/ (Gohlke, 2021).

**Data availability.** Data sets are available upon request by contacting the corresponding author.

**Author contributions.** All authors contributed to the data collection, algorithm development, analysis and presentation of results equally.

**Competing interests.** The contact author has declared that neither they nor their co-authors have any competing interests.

**Disclaimer.** Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Acknowledgements.** The authors express their gratitude to the staff of the State Traffic Service of the Yamalo-Nenetsky Autonomous District (Russia) for providing the valuable information about ice road operation. The authors would also like to thank Sofia Antonova and two anonymous referees for their thorough and constructive comments.

**Financial support.** This research was made possible with support from the RFBR (project no. 18-05-60021-Arctic), ESA EO Science for Society (LIAM project; contract no. 4000130930/20/I-DT) and ESA CCI Lakes+ (CCI+ Phase 1 new ECVs - Lakes; contract no. 4000125030/18/I-NB); the approbation of results was carried out in the framework of a government programme at the Water Problems Institute, Russian Academy of Sciences (no. 0147-2019-0004).

**Review statement.** This paper was edited by Moritz Langer and reviewed by Sofia Antonova and two anonymous referees.

## References

Agafonova, S. A. and Vasilenko, A. N.: Hazardous ice phenomena in rivers of the Russian arctic zone under current climate conditions and the safety of water use, Geogr. Environ. Sustain., 13, 43–51, https://doi.org/10.24057/2071-9388-2020-12, 2020.

Antonova, S., Duguay, C. R., Kääb, A., Heim, B., Langer, M., Westermann, S., and Boike, J.: Monitoring ice phenology and bedfast ice in lakes of the Lena River Delta using TerraSAR-X backscatter and coherence time series, Remote Sens., 8, 903, https://doi.org/10.3390/rs8110903, 2016.

Atwood, D. K., Gunn, G. E., Roussi, C., Wu, J., Duguay, C. R., and Sarabandi, K.: Microwave backscatter from Arctic lake ice and polarimetric implications, IEEE T. Geosci. Remote, 53, 5972–5982, https://doi.org/10.1109/TGRS.2015.2429917, 2015.

AVISO+: Jason-2 and Jason-3 geophysical data records (GDR), AVISO+ data portal [data set], https://www.aviso.altimetry.fr/es/data/data-access.html, last access: 2 December 2021).

Bamber, J. L.: Ice sheet altimeter processing scheme, Int. J. Remote Sens., 15, 925–938, https://doi.org/10.1080/01431169408954125, 1994.

Beaton, A., Whaley, R., Corston, K., and Kenny, F.: Identifying historic river ice breakup timing using MODIS and Google Earth Engine in support of operational flood monitoring in Northern Ontario, Remote Sens. Environ., 224, 352–364, https://doi.org/10.1016/j.rse.2019.02.011, 2019.

Beckers, J. F., Casey, J. A., and Haas, C.: Retrievals of lake ice thickness from Great Slave Lake and Great Bear Lake using CryoSat-2, IEEE T. Geosci. Remote, 55, 3708–3720, https://doi.org/10.1109/TGRS.2017.2677583, 2017.

Beltaos, S.: Hydrodynamic characteristics and effects of river waves caused by ice jam releases, Cold Reg. Sci. Technol., 85, 42–55, https://doi.org/10.1016/j.coldregions.2012.08.003, 2013.

Beltaos, S., Carter, T., Rowsell, R., and DePalmac, S. G. S.: Erosion potential of dynamic ice breakup in Lower Athabasca River. Part I: Field measurements and initial quantification, Cold Reg. Sci. Technol., 149, 16–28, https://doi.org/10.1016/j.coldregions.2018.01.013, 2018.

Berry, P. A. M., Garlick, J. D., Freeman, J. A., and Mathers, E. L.: Global inland water monitoring from multi-mission altimetry, Geophys. Res. Lett., 32, L16401, https://doi.org/10.1029/2005GL022814, 2005.

Brown, G. S.: The average impulse response of a rough surface and its applications, IEEE J. Ocean Eng., 2, 67–74, https://doi.org/10.1109/JOE.1977.1145328, 1977.

Chaouch, N., Temimi, M., Romanov, P., Cabrera, R., McKillop, G., and Khanbilvardi, R.: An automated algorithm for river ice monitoring over the Susquehanna River using the MODIS data, Hydrol. Process., 28, 62–73, https://doi.org/10.1002/hyp.9548, 2014.

Chu, T. and Lindenschmidt, K.-E.: Integration of space-borne and air-borne data in monitoring river ice processes in the Slave River, Canada, Remote Sens. Environ., 181, 65–81, https://doi.org/10.1016/j.rse.2016.03.041, 2016.

Cooley, S. W. and Pavelsky, T. M.: Spatial and temporal patterns in Arctic river ice breakup revealed by automated ice detection from MODIS imagery, Remote Sens. Environ., 175, 310–322, https://doi.org/10.1016/j.rse.2016.01.004, 2016.

Du, J., Kimball, J. S., Duguay, C. R., Kim, Y., and Watts, J.: Satellite microwave assessment of Northern Hemisphere lake
ice phenology from 2002 to 2015, The Cryosphere, 11, 47–63,
https://doi.org/10.5194/tc-11-47-2017, 2017.

Duguay, C. R., Pultz, T. J., Lafleur, P. M., and Drai, D.: RADARSAT backscatter characteristics of ice growing on shallow sub-arctic lakes, Churchill, Manitoba, Canada, Hydrol. Process., 16, 1631–1644, https://doi.org/10.1002/hyp.1026, 2002.

Duguay, C. R., Bernier, M., Gauthier, Y., and Kouraev, A.: Remote sensing of lake and river ice, in: Remote Sensing of the Cryosphere, edited by: Tedesco, M., Wiley-Blackwell, Oxford, UK, 273–306, ISBN 13:978-1118368855, 2015.

Duguay, C. R., Zakharova, E. A., Kouraev, A. V., Kheyrollah Pour, H., and Hoekstra, M.: Retrieval of ice thickness on large northern lakes from Jason-2 data, in: POLAR2018: Abstract Proceedings, Open Science Conference, Davos, 19–23 June 2018, Switzerland, Abstract No. 2638, p. 46, 2018.

ESA: ENVISAT RA2/MWR Product Handbook, RA2/MWR Products User Guide, available at: https://earth.esa.int/eogateway/documents/20142/37627/Envisat-RA-2-Level-2-Product-Handbook.pdf (last access: 2 December 2021), 2002.

Ettema, R.: Review of alluvial-channel responses to river ice, J. Cold Reg. Eng., 16, 191–217, https://doi.org/10.1061/(ASCE)0887-381X(2002)16:4(191), 2002.

Fu, L.-L. and Cazenave, A.: Satellite Altimetry and Earth Sciences: A Handbook of Techniques and Applications, Academic Press, San Diego, 1–495, ISBN 0080516580, 9780080516585, 1991.

Ginzburg, B. M.: Probabilistic characteristics of freeze-up and breakup dates on rivers and reservoirs of the Soviet Union, Hydrometeoizdat, Leningrad, p. 110, 1973.

Gohlke, C.: Unofficial Windows Binaries for Python Extension Packages [data set], https://www.lfd.uci.edu/~gohlke/pythonlibs/, last access: 2 December 2021.

Guerreiro K., Fleury, S., Zakharova, E., Rémy, F., and Kouraev, A.: Potential for estimation of snow depth on Arctic sea ice from CryoSat-2 and SARAL/AltiKa missions, Remote Sens. Environ., 186, 339–349, https://doi.org/10.1016/j.rse.2016.07.013, 2016.

Gunn G. E., Duguay, C. R., Brown, L. C., King, J. M. L., Atwood, D., and Kasurak, A.: Freshwater lake ice thickness derived using surface-based X- and Ku-band FMCW scatterometers, Cold Reg. Sci. Technol., 120, 115–126, https://doi.org/10.1016/j.coldregions.2015.09.012, 2015.

Gunn, G. E., Duguay, C. R., Atwood, D. K., King, J., and Toose, P.: Observing scattering mechanisms of bubbled freshwater lake ice using polarimetric RADARSAT-2 (C-band) and UW-Scat (X- and Ku-bands), IEEE T. Geosci. Remote, 56, 2887–2903, https://doi.org/10.1109/TGRS.2017.2786158, 2018.

Hiks, F.: An overview of river ice problems: CRIBE07 guest editorial, Cold Reg. Sci. Technol., 55, 175–185, https://doi.org/10.1016/j.coldregions.2008.09.006, 2009.

Instructions on Safety Organisation: Of rivers' and lakes' crossing, RD 34.03.221, INFORMENRGO, Moscow, 1969.

Jeffries, M. O., Morris, K., and Kozlenko, N.: Chapter 4 – ice characteristics and processes, and remote sensing of frozen Rivers and lakes, in: Remote Sensing in Northern Hydrology: Measuring Environmental Change, edited by: Duguay, C. R. and Pietroniro, A., American Geophysical Union, Washington, DC, 63–90, ISBN 13:978-0-87590-428-3, 2005.

Jeffries, M. O., Morris, K., and Duguay, C. R.: State of the earth’s cryosphere at the beginning of the 21st century: glaciers, global snow cover, floating ice, and permafrost and periglacial environments – floating ice: lake ice and river ice, in: Satellite Image Atlas of Glaciers of the World, edited by: Williams, R. S. and Ferrigno, J. G., USGS, Reston, VA, A381–A424, https://doi.org/10.3133/pp1386A, 2012.

Kang, K.-K., Duguay, C. R., and Howell, S. E. L.: Estimating ice phenology on large northern lakes from AMSR-E: Algorithm development and application to Great Bear Lake and Great Slave Lake, Canada, The Cryosphere, 6, 235–254, https://doi.org/10.5194/tc-6-235-2012, 2012.

Kang, K.-K., Duguay, C. R., Lemmetyinen, J., and Gel, Y.: Estimation of ice thickness on large northern lakes from AMSR-E brightness temperature measurements, Remote Sens. Environ., 150, 1–19, https://doi.org/10.1016/j.rse.2014.04.016, 2014.

Kheyrollah Pour, H., Duguay, C. R., Scott, A., and Kang, K.-K.: Improvement of lake ice thickness retrieval from MODIS satellite data using a thermodynamic model, IEEE T. Geosci. Remote, 55, 5956–5965, https://doi.org/10.1109/TGRS.2017.2718533, 2017.

King, J. M. L., Kelly, R., Kasurak, A., Duguay, C., Gunn, G., and Mead, J. B.: UW-Scat – ground-based dual frequency scatterometry for observation of snow processes, IEEE Geosci. Remote Sens. Lett., 10, 528–532, https://doi.org/10.1109/LGRS.2012.2212177, 2013.

King, J. M. L., Kelly, R., Kasurak, A., Duguay, C., Gunn, G., Rutter, N., Watts, T., and Derksen, C.: Spatio-temporal influence of tundra snow properties on Ku-band (17.2 GHz) backscatter, J. Glaciol., 61, 267–279, https://doi.org/10.3189/2015JoG14J020, 2015.

Kouraev, A. V., Zakharova, E. A., Samain, O., Mognard, N. M., and Cazenave, A.: Ob’ River discharge from TOPEX/Poseidon satellite altimetry (1992–2002), Remote Sens. Environ., 93, 238–245, https://doi.org/10.1016/j.rse.2004.07.007, 2005.

Kouraev, A. V., Semovski, S. V., Shimaraev, M., Mognard, N. M., Legresy, B., and Rémy F.: Observations of Lake Baikal ice from satellite altimetry and radiometry, Remote Sens. Environ., 108, 240–253, https://doi.org/10.1016/j.rse.2006.11.010, 2007.

Kouraev, A. V., Zakharova, E. A., Rémy, F., and Suknev, A. Y.: Study of Lake Baikal ice cover from radar altimetry and in situ observations, Mar. Geod., 38, 477–486, https://doi.org/10.1080/01490419.2015.1008155, 2015.

Kouraev, A. V., Zakharova, E. A., Kostianoy, A. G., Shimaraev, M. N., Desinov, L. V., Petrov, E. A., Hall, N. M. J., Rémy, F., and Suknev, A. Y.: Giant ice rings in southern Baikal: multi-satellite data help to study ice cover dynamics and eddies under ice, The Cryosphere, 15, 4501–4516, https://doi.org/10.5194/tc-15-4501-2021, 2021.

Kurtz, N. T., Galin, N., and Studinger, M.: An improved CryoSat-2 sea ice freeboard retrieval algorithm through the use of waveform fitting, The Cryosphere, 8, 1217–1237, https://doi.org/10.5194/tc-8-1217-2014, 2014.

Larue, F., Picard, G., Aublanc, J., Arnaud, L., Robledano-Perez, A., Le Meur, E., Favier, V., Jourdain, B., Savarino, J., and Thibaut, P.: Radar altimeter waveform simulations in Antarctica with the Snow Microwave Radiative Transfer Model (SMRT), Remote Sens. Environ., 263, 112534, https://doi.org/10.1016/j.rse.2021.112534, 2021.
Leconte, R., Daly, S., Gauthier, Y., Yankielun, N., Bérubé, F., and Bernier, M.: A controlled experiment to retrieve freshwater ice characteristics from an FM-CW radar system, Cold Reg. Sci. Technol., 55, 212–220, 2009.

Legrésy, B. and Rémy, F.: Surface characteristics of the Antarctic ice sheet and altimetric observations, J. Glaciol., 43, 265–275, https://doi.org/10.3189/S002214300000321X, 1997.

Mermoz, S., Allain, S., Bernier, M., and Pottier, E.: Investi-
gation of Radarsat-2 and TerraSAR-X data for river
ice classification, in: IEEE International Geoscience and
Remote Sensing Symposium, Cape Town, II-29-II-32,
https://doi.org/10.1109/IGARSS.2009.5417991, 2009.

Mermoz, S., Allain-Bailhache, S., Bernier, M., Pottier, E., van der Sanden, J. J., and Chokmani, K.: Retrieval of river ice thickness from C-band PolSAR data, IEEE T. Geosci. Remote, 52, 3052–3062, https://doi.org/10.1109/TGRS.2013.2269014, 2013.

Michailovsky, C. I., McEnnis, S., Berry, P. A. M., Smith, R., and Bauer-Gottwein, P.: River monitoring from satellite radar altimetry in the Zambezi river basin, Hydrol. Earth Syst. Sci., 16, 2181–2192, https://doi.org/10.5194/hess-16-2181-2012, 2012.

Muhammad, P., Duguay, C. R., and Kang, K.-K.: Monitoring ice break-up on the Mackenzie River using remote sensing, The Cryosphere, 10, 569–584, https://doi.org/10.5194/tc-10-569-2016, 2016.

Murfitt, J. and Duguay, C. R.: 50 years of lake ice re-
search from active microwave remote sensing: Progress
and prospects, Remote Sens. Environ., 264, 112616,
https://doi.org/10.1016/j.rse.2021.112616, 2021.

USGS: Landsat Science project, USGS [data set], https://earthexplorer.usgs.gov/, last access: 2 December 2021.

Nilsson, J., Vallelonga, P., Simonsen, S. B., Sorensen, L. S., Forsberg, R., Dahl-Jensen, D., Hirabayashi, M., Goto-Azuma, K., Hvidberg, C. S., Kjaer, H. A., and Satow, K.: Greenland 2012 melt event effects on CryoSat-2 radar altimetry, Geophys. Res. Lett., 42, 3919–3926, https://doi.org/10.1002/2015GL063296, 2015.

Pavelsky, T. M. and Smith, L. C.: Spatial and temporal pat-
terns in Arctic river ice breakup observed with MODIS and
AVHRR time series, Remote Sens. Environ., 93, 328–338,
https://doi.org/10.1016/j.rse.2004.07.018, 2004.

Picard, R. and Cook, R.: Cross-validation of regression models, J. Am. Stat. Assoc., 79, 575–583, https://doi.org/10.2307/2288403, 1984.

Prowse, T., Alfredsen, K., Deltaos, S., Bonsal, B., Duguay, C., Korhola, A., McNamara, J., Pienitz, R., Vincent, W. F., Vuglinsky, V., and Weyhenmeyer, G. A.: Past and future changes in Arctic lake and river ice, Ambio, 40, 53–62, https://doi.org/10.1007/s13280-011-0216-7, 2011a.

Prowse, T., Alfredsen, K., Deltaos, S., Bonsal, B., Bowden, W., Duguay, C., Korhola, A., McNamara, J., Vincent, W. F., Vuglinsky, V., Anthony, K., and Weyhenmeyer, G. A.: Effects of changes in Arctic lake and river ice, Ambio, 40, 63–74, https://doi.org/10.1007/s13280-011-0217-6, 2011b.

Prowse, T. D.: River-ice ecology: part B. Biological aspects, J. Cold Reg. Eng., 15, 17–33, https://doi.org/10.1061/(ASCE)0887-381X(2001)15:1(17), 2001.

Prowse, T. D., Bonsal, B. R., Duguay, C. R., and Lacroix, M. P.: River-ice break-up/freeze-up: A review of climatic drivers, historical trends, and future predictions, Ann. Glaciol., 46, 443–451, 2007.

Rémy, F., Flament, T., Blarel, F., and Benveniste, J.: Radar altimetry measurements over Antarctic ice sheet: a focus on antenna polarization and change in backscatter problems, Adv. Space Res., 50, 998–1006, https://doi.org/10.1016/j.asr.2012.04.003, 2012.

Slater, T., Shepherd, A., McMillan, M., Armitage, T. W. K., Oto-
saka, I., and Arthern, R. J.: Compensating changes in the
penetration depth of pulse-limited radar altimetry over the
Greenland ice sheet, IEEE T. Geosci. Remote, 57, 9633–9642,
https://doi.org/10.1109/TGRS.2019.2928232, 2019.

Sobiech, J. and Dierking, W.: Observing lake- and river-ice de-
cay with SAR: advantages and limitations of the unsuper-
vised k-means classification approach, Ann. Glaciol., 54, 65–72,
https://doi.org/10.3189/2013AoG62A037, 2013.

Sun, W. and Trevor, B.: A stacking ensemble learning framework for annual river ice breakup dates, J. Hydrol., 561, 636–650, https://doi.org/10.1016/j.jhydrol.2018.04.008, 2018.

Ulaby, F. T., Moore, R. K., and Fung, A. K.: Microwave Remote Sensing: Active and Passive, in: Radar Remote Sensing and Surface Scattering and Emission Theory, vol. 2, Addison-Wesley, Norwood, MA, USA, 1986.

Unterschultz, K., van der Sanden, J., and Hicks, F.: Poten-
tial of RADARSAT-1 for the monitoring of river ice: Re-
sults of a case study on the Athabasca River at Fort Mc-
Murray, Canada, Cold Reg. Sci. Technol., 55, 238–248,
https://doi.org/10.1016/j.coldregions.2008.02.003, 2009.

van der Sanden J., Drouin, H., and Geldsetzer, T.: An au-
tomated procedure to map breaking river ice with C-
band HH SAR data, Remote Sens. Environ., 252, 112119,
https://doi.org/10.1016/j.rse.2020.112119, 2021.

Vuglinsky, V. and Valatin, D.: Changes in ice cover du-
ration and maximum ice thickness for rivers and lakes
in the Asian part of Russia, Nat. Resour., 9, 73–87,
https://doi.org/10.4236/nr.2018.93006, 2018.

Willatt, R., Laxon, S., Giles, K., Cullen, R., Haas, C., and Helm, V.: Ku-band radar penetration into snow cover on Arctic sea ice using airborne data, Ann. Glaciol., 52, 197–205, https://doi.org/10.3189/172756411795931589, 2011.

Zakharova, E. A., Kouraev, A. V., Rémy, F., Zemtsov, V. A., and Kirpotin, S. N.: Seasonal variability of the Western Siberia wetlands from satellite radar altimetry, J. Hydrol., 512, 366–378, https://doi.org/10.1016/j.jhydrol.2014.03.002, 2014.

Zakharova, E. A., Kouraev, A. V., Guillaso, S., Garestier, F., Desyatkin, R. V., and Desyatkin, A. R.: Recent dynamics of hydro-ecosystems in thermokarst depressions in Central Siberia from satellite and in situ observations: importance for agriculture and human life, Sci. Total Environ., 615, 1290–1304, https://doi.org/10.1016/j.scitotenv.2017.09.059, 2018.

Zakharova, E. A., Krylenko, I. N., and Kouraev, A. V.: Use of non-polar orbiting satellite radar altimeters of the Jason series for estimation of river input to the Arctic Ocean, J. Hydrol., 568, 322–333, https://doi.org/10.1016/j.jhydrol.2018.10.068, 2019.
Zakharova, E. A., Nielsen, K., Kamenev, G., and Kouraev, A.: River discharge estimation from radar altimetry: Assessment of satellite performance, river scales and methods, J. Hydrol., 583, 124561, https://doi.org/10.1016/j.jhydrol.2020.124561, 2020.

Zhang, F., Li, Z., and Lindenschmidt, K.-E.: Potential of RADARSAT-2 to improve ice thickness calculations in remote, poorly accessible areas: A case study on the Slave River, Canada, Can. J. Remote Sens., 45, 234–245, https://doi.org/10.1080/07038992.2019.1567304, 2019.