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

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 as an alternate tool or as a complement to ground-based measurements, allowing for the characterisation of river ice properties and processes at spatial and temporal resolutions suitable for addressing various scientific and operational requirements (Duguay et al., 2015).

Indeed, spaceborne instruments provide observational capabilities that can be used to map and monitor many river ice parameters. For example, optical sensors, such as the Moderate 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 conditions, particularly during the freeze-up period, are limiting factors for monitoring river ice at high latitudes using optical sensors. Active sensors operating in the microwave region of the electromagnetic spectrum are weather independent and provide a spatial resolution higher than that provided 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 Lindenschmidt, 2016). Comprehensive reviews on remote sensing 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. Passive microwave and thermal satellite instruments have exhibited 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 dimension of rivers, notably the width of the channels, limits 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) spatial 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 estimation of ice thickness for short time periods (one or a few ice seasons) and small sections of rivers, primarily due to limited access to SAR imagery in years prior to the launch of the Sentinel-1A and Sentinel-1B missions (2014) of the European Space Agency (ESA). Altimetric radar represents an excellent alternative or complement to SAR satellites. Satellite 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 spatiales) 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 surface. 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 signals) can be used to characterise the surface state within a radar footprint. Radar backscatter from freshwater ice depends on the radar configuration (radar type, operating frequency, viewing geometry and polarisation) and on the material 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). However, radar altimeters have not been used for monitoring the ice regime of rivers. Therefore, for the first time, we demonstrate the application of radar altimetry for monitoring the interannual 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 prediction of ice bridge road operations in the city of Salekhard on the Ob River. For this purpose, data from two altimetric 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 (Topography Ocean Experiment) mission and continuing with the Sentinel-6 mission (launched in November 2020), which is offering an interesting perspective for a follow-up investigation 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 annual flow of 406 km³ (Zakharova et al., 2020). The lower 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).

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

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

Table 1 Gauging stations in the lower reaches of the 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.

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

According to the in situ observations at the gauging stations, ice formation begins on average between 23 and 27 October. 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 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 returned 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).

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.

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 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 $\mathrm{\Delta }{\mathrm{Sig}}_{\mathrm{0}}/\mathrm{\Delta }t$ (where Δ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 refreezing 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_t − mean(ΔTB_Jul–Aug)] (where mean(ΔTB_Jul–Aug) 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${}_{\mathrm{0}}/\mathrm{\Delta }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 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.

Figure 4 Processing scheme of retrieving ice phenology dates from altimetric measurements.

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${}_{\mathrm{0}}/\mathrm{\Delta }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${}_{\mathrm{0}}/\mathrm{\Delta }t)$ and constructed the relationship between this parameter and in situ ice thickness (H_ice) measured at the nearest gauging station (Fig. 5). The application of a locally estimated scatterplot smoothing filter (LOESS) on the Σ(ΔSig${}_{\mathrm{0}}/\mathrm{\Delta }t)$ parameter makes it possible to minimise the effect of secondary peaks on the winter backscatter curve.

Figure 5 Relationship between $\mathrm{\Sigma }(\mathrm{\Delta }\mathrm{Sig}\mathrm{0}/\mathrm{\Delta }t)$ and in situ ice thickness at 10 virtual stations.

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${}_{\mathrm{0}}/\mathrm{\Delta }t)$ and in situ H_ice and the minimisation of root mean square error (RMSE) calculated between the H_ice retrieved from Eq. (1) and the observed values.

$$\begin{array}{}\text{(1)}& H_{\mathrm{ice}}=a\times \mathrm{abs}{\left(\mathrm{\Sigma }\left(\mathrm{\Delta }{\mathrm{Sig}}_{\mathrm{0}}/\mathrm{\Delta }t\right)\right)}^b\end{array}$$

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.

Figure 6 Processing scheme of ice thickness retrieval from altimetric measurements.

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 VS_t. Then, the coefficients a and b of VS_t 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 H_ice magnitudes and the spatial heterogeneity of ice thickness in the spatial domain as much as possible and (b) reduce 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_ice time series retrieved at virtual stations and the H_ice time series extracted from the smoothed product for the same VS locations. The effect of interpolation and smoothing on seasonal H_ice magnitude was evaluated by estimating the differences between the maximum yearly values of H_ice.

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.

Notably, a critical effect of the interpolation/smoothing procedure on H_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, 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_ice ($H_{{\mathrm{ice}}_{\max }}$) were below 0.03 m). The selected window of 40 km every 30 d allowed keeping the average $H_{{\mathrm{ice}}_{\max }}$ difference and RMSE below 0.04 and 0.03 m, respectively. The correlation between the H_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 detection 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.

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.

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

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.

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), 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 areas between the virtual stations. For instance, two useful parameters can be extracted from the 2D product: the maximum ice thickness and the ice thickness observed on 1 December.

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

Generally, knowledge of the maximum ice thickness is relevant for hydro-climate change monitoring (Vuglinsky and Valatin, 2018); the ice thickness determined on 1 December is crucial for the local and regional socioeconomic stakeholders 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 compared 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.

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.

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.

Table 3 Correlation coefficient (R) and RMSE between the retrieved and in situ ice thickness in the cross-validation experiment.

^* Two different pairs of a and b coefficients were used for cross-validation.

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

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.

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.

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.

^* Mean and standard deviation are given for predicted (nominator) and observed (denominator) dates.

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

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

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 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 ΔSig${}_{\mathrm{0}}/\mathrm{\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_{\mathrm{ice}}-\mathrm{\Sigma }(\mathrm{\Delta }{\mathrm{Sig}}_{\mathrm{0}}/\mathrm{\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_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_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 ΔTB. In the manual selection of peaks, the ΔTB criteria were used in fewer cases than in the automated selection (see Sect. 5.1 for explications). 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 temperature variability in different VS configurations could help to uncover universal criteria and improve the automated routine. The development of multi-satellite (combination of altimetry missions) and multi-sensor approaches (combination of altimetry with SAR and optical sensors) could help to improve 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 thickness estimates is also comparable with the accuracy reported in the rare studies dedicated to satellite retrievals of freshwater ice thickness. The SAR-based retrievals for three Canadian 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 accuracy of the ice thickness estimates, especially at the beginning of freeze-up. Another reason for the uncertainties is the multi-peak character of the winter backscatter recession 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_{{\mathrm{ice}}_{\mathrm{alti}}}-H_{{\mathrm{ice}}_{\mathrm{in}\mathrm{situ}}}$ scatterplot (Fig. 9). Several improvements are proposed for future studies. First, an improvement in the detection accuracy 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 decreases 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 corresponds 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 contributes 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 flooding 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 product, 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_ice retrievals are quite high, and the ice thickness is often overestimated (see Fig. 11b). Therefore and considering 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. Nevertheless, their accuracy is sufficient for the assessment of climate impacts, as we captured an overall trend for an earlier opening of the road demonstrated by both in situ and altimetric 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 applications 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 predictive capacity. Their elaboration and testing will be the subject of future study.

The prediction algorithm for the Salekhard ice road operation 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 located 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 dynamics 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 current 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 ΔSig${}_{\mathrm{0}}/\mathrm{\Delta }t$, instead of the absolute value of Sig0, 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 ΔSig${}_{\mathrm{0}}/\mathrm{\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 ±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 ±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.