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 15 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-day repeat cycle of Jason-2/3 in 88–90% of the cases analysed. The uncertainties in the river ice 20 thickness retrievals made via empirical relations between the satellite backscatter measurements and in situ observations, expressed as root mean square errors (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 Salekhard City can be predicted with an accuracy (expressed as RMSE) of 3–5 days. 25

stations, as well as construction and navigation activities. In arctic regions, frozen rivers provide a unique transportation infrastructure for the movement of merchandise and people via winter ice roads. The presence of river ice cover also provides local populations with access to fishing grounds and in some cases (e.g. Central Yakutia, Russia) to fresh water.
However, operational monitoring of ice in northern rivers is difficult because of low site accessibility. Moreover, icye conditions can be unsafe for people who are required to perform in-situ measurements, especially at the beginning or end of the ice season. Therefore, satellite remote sensing observation has been proposed as an alternative, or complement, to ground measurements, allowing for the characterisation of river ice at a temporal resolution suitable for addressing various climatic, scientific, and operational requirements.
Satellite-borne instruments provide observational capabilities for many river ice parameters. Optical sensors such as the Moderate Resolution Imaging Spectroradiometer (MODIS) and the Advanced Very High Resolution Radiometer (AVHRR) have been used to map river ice extent and phenology -freeze-up and breakup dates (Pavelsky and 50 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 ice monitoring river ice of rivers at high latitudes using optical sensors. Active sensors operating in the microwave spectrum are weather independent and provide a spatial resolution higher than that of the MODIS and AVHRR instruments. One such 55 sensor is synthetic aperture radar (SAR), a microwave sensor which has been used to monitor river ice phenology (Unterschultz et al., 2009;Mermoz et al., 2009;Sobiech et al., 2013;Sun and Trevor, 2018,), deformation (Unterschultz et al., 2009), and classification of ice types (Chu and Lindenschmidt, 2016). A comprehensive review of the SAR instrument application for remote sensing of lake and river ice can be found in Jeffries et al. (2005) and .

60
Ice thickness is another parameter which is of particular interest for operational purposes, such as public safety, ice road service, jam forecasts, and mitigation. Passive microwave and thermal satellite instruments for the retrieval of ice thickness have demonstrated capability for the retrieval of ice thickness for large lakes (Kang et al., 2014;Duguay et al., 2002Gunn et al., 2015;Kheyrollah Pour et al., 2017). However, the spatial dimension of rivers, notably the width of channels, limits the application of these instruments as they only provide coarse spatial resolution (km to tens of km). Several studies have used active microwave SAR images with high, 10-25 m, spatial resolution to retrieve river ice thickness (Unterschultz et al., 2009;Mermoz et al., 2014, Zhang et al., 2019, mainly, via establishing a statistical relation between backscatter and in situ ice thickness measurements. The SAR-based method of ice thickness long-term monitoring did not receive a large application because of limited asses to SAR data in years before launching of Sentinel 1 missions and because of sufficiently laborious processing of SAR 70 images for large spatial domain and for long period of time. Altimetric radars represent an excellent alternative or complement to SAR satellites. The satellite altimeters are also high resolution (200-400 m) weather-independent instruments with a long history of observations starting in the middle of 1980s. They have been largely used for monitoring the water level in narrow inland waterbodies and water courses starting at 100 m in width (Michailovsky et al., 2012). The altimeter data products are freely accessible via corresponding French and 75 European Space Agencys' portals and can be routinely and quickly processed over large areas for long period of time.
Radar signals incident upon the Earth's surface are modified according to the physical properties of surface materials. Similar to SAR systems, the signal recorded by radar altimeter can be interpreted as a function of changes in surface properties, and the backscatter coefficient (ratio between power of emitted and received signal) can be used to characterize surface state within a radar footprint. The value of radar backscatter over freshwater ice depends on radar configuration (radar type, frequency band) and on 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;Antonova et al., 2016;Gunn et al., 2018). Altimetric radars have already proven capability in monitoring ice phenology dates of large Eurasian lakes (Kouraev et al., 85 2007(Kouraev et al., 85 , 2015 and for measurements of ice thickness of large Canadian lakes (Beckers et al., 2017). However, the radar altimeters have never been used for monitoring ice regime of rivers. For the first time in this paper, we demonstrate the application of radar altimetry for interannual monitoring of river ice phenology and river ice thickness over long (400 km) section of the Ob River (Siberia, Russia).
The objective of the present study is to investigate the capacity of radar altimetry satellites 1) for retrieving river ice 90 onset/breakup dates and ice thickness, as well as 2) for providing local communities with operational information useful for ice bridge road construction and maintenance. Two altimetric satellite missions, Jason2 and Jason 3, were selected for this study for two reasons: 1) their best among altimetric missions temporal resolution (10 days); and 2) the long lifetime of this series of satellites providing observations on the same orbit, starting in 1992 with the Topex/Poseidon and continuing nowadays with the Sentinel-6 (on orbit since November 2020). 95 2 Regional setup and Data

Study Region
The study was conducted over an area that encompasses lower reaches of the Ob River (Siberia, Russia). The Ob 100 River drains the Western Siberian Plain and is the third largest river in the Arctic Ocean watershed, with a mean annual flow of 406 km 3 (Zakharova et al., 2020). The lower reaches of the Ob River extend for approximately 800 km and begin at the confluence of the Ob and the Irtysh rivers at 61.08°N (Figure 1). The reaches are characterised by a particularly wide floodplain (up to 50 km) and the presence of numerous branches. The easternmost channel is the main, largest branch, called the Big Ob. The second largest branch delineates the flood plain from the west. The

105
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 seven months of the year. In winter, when bogs are frozen, local communities intensify their socio-economic activities by constructing winter roads and ice bridges over river crossings. River ice observations on the Ob are sparse and taken 110 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-Mansyisk). 115 also provided at river crossings. The main map is produced using public The World Bank data (https://datacatalog.worldbank.org/dataset/major-rivers-world).

In-situ data and ice regime in studied river reach
The Russian Hydrometeorological Service monitors ice at all gauging stations measuring water level. There are five water level gauging stations in the studied Ob River reaches ( Figure 1, 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 on gauging stations includes: 1) daily visual observations of ice presence/absence and ice type/events; and 2) 3-6 times per month measurements of ice thickness and on-ice snow depth. Ice thickness was is measured by drilling one hole in ice using an ice auger. The snow depth corresponds to 130 the average value calculated from three snow-depth measurements located around the hole.
According to in-situ observations at the gauging stations, ice formation begins between 23and 27 October. For the last 20 years, the earliest and latest records for the last 20 years were 1 October and 18 November, respectively. Ice cover forms quickly on this section of the river, typically within just 2-3 days from appearance of first ice. However, in 15% of cases (for the last 20 years), the installation full freeze-up took up to 10 days. Ice grows rapidly during the According to observations on the gauging stations the temporal dynamics of ice growth on the large linear channel sections is similar throughout the studied reaches (Fig. 2a). However, in the south of the region, the ice thickness is 140 0.07-0.20 m less than in the north. Climate change has affected river ice in the Canadian Arctic (Prowse et al., 2011b) and the European part of Russia (Agafonova and Vasilenko, 2020) but has not yet resulted in a significant change in the ice regime in the lower Ob River. Mann-Kendall tests performed on ice time series with complete records (Polnovat, Gorki and Muzhi; see Table 1) reveal that long-term trends in ice onset dates are occurring significantly later (0.05 level) only at the Polnovat station, while maximum ice thicknesses have become

Altimetry data
The Jason-2 satellite is the third altimetric satellite of the Topex/Poseidon-Jason series. The satellite operated 155 between 2008-2016 and acquired data in a 10-day 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).
The altimetric radar aboard Jason-2 provided along-track measurements at Ku (13.6 GHz) and C (5.3 GHz) bands 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 160 of our study region (63-66°N), the cross-track oscillation band was approximately 400 m wide. Because the radar ground footprint in the Ku-band is lower than that in the C-band, the measurements in the Ku-band were used in this study. The theoretical footprint of the radar at the Ku-band is 10-12 km in diameter over the rough ocean surface.
However, this diameter decreases over smooth inland water and ice surfaces such that the main return signal can come from footprints of just a few kilometres in diameter (Legresy and Rémy, 1997).

165
The satellite payload of Jason-2 also included a nadir-looking Advanced Microwave Radiometer (AMR), which provides measurements of brightness temperature in 18.7, 23.8, and 34.0 GHz bands at a nadir-looking angle with a sampling frequency of 1 Hz. Brightness temperature measurements acquired with by other passive microwave radiometers, such as SSM/I and AMSR-E, have demonstrated good performance for the retrieval of ice phenology dates and ice thickness onf large lakes in Russia and Canada (Kouraev et al.,2007, Kang et al., 2014. Because the 170 Jason AMR footprints are large, 42 km (18.7 GHz), 35 km (23.8 GHz), and 22 km (34.0 GHz) in diameter (Kouraev et al., 2007), the radiometric measurements over rivers are dominated by signals emitted mainly from land surfaces surrounding the river channels. That why in this study, we used Jason-2 and 3 AMR measurements only as auxiliary information to develop additional criteria for adjustment of the radar freeze-up/ breakup date retrieval algorithm.
In 2016, the successor to the Jason-2 satellite, the Jason-3, was sent into space with the same orbit and instruments as Jason-2. For 20 cycles, the two missions flew with an 80-second time lag, ensuring continuity of measurements.
During this period, the difference (bias) between Jason-2 and Jason-3 for Ku-band backscatter (Sig0) was within 1 dB. The difference in the brightness temperature measurements was within 3 °K.
For this study, the satellite radar measurements were extracted from a geophysical research data records product (GDR) distributed by AVISO+ data portal (avisoftp.cnes.fr, last access 20201/0306/2030). The GDR product 180 contains various parameters estimated from the radar return echo, which is represented as a waveform. In our study, we used the 20 Hz backscatter coefficient retrieved by an ICE1 algorithm (Bamberg, 1994). An altimeter waveform represents a histogram of energy backscattered by ground surface to the satellite with respect to time. The waveform ground-processing consists in its retracking aimed to obtain better range (distance from satellite to the Earth surface) estimates than those obtained with on-board algorithms. The ICE1retracking 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, 2020). The AMR measurements of the brightness temperature have a 1 Hz temporal frequency. In this study they were linearly interpolated to the coordinates of the 20 Hz radar measurements.
High-resolution (30 m) optical Landsat 8 images (https://earthexplorer.usgs.gov/) were used for the geographical 190 selection of altimetric radar and radiometer measurements over river channels using our own Python code that allowing the overlapping of along-track Jason measurements and Landsat images. The cross-section of an altimetric track with a river channel is called a virtual station (VS). The VS received 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 were extendedthe station names were extended with the corresponding 195 index "S_Ob" (if necessary).

Backscatter variability
An effect of different ice properties on return signal of nadir-looking altimeter instruments has been largely studied

200
for Antarctic and Greenland ice sheets (see for example Rémy et al., 2012;Nilsson, 2015;Larue et al., 2021), as well as for sea ice (see for example Kurtz et al.,2014;Guerreiro et al.,2016). However, the freshwater lake and river ice differs significantly from other types of ice either by its texture (crystal structure, density, layering, air content etc) or by chemical properties (salt content). At the same time, the freshwater ice properties affecting the satellite measurements in microwave range have been extensively studied in context of interpretation of measurements of 205 side-looking SAR instruments (Jeffries et al., 2005;Duguay et al., 2002;Unterschultz et al., 2009;Atwood et al., 2015 and many others). For a moment, there are no works synthesising the temporal variability of radar altimetry return signal over frozen rivers. Bellow, we present a short summary of variability of altimetric signal over river ice during different phases of ice formation, growth and breakup.
Previous studies (Kouraev et al., 2005;Duguay et al., 2018, Zakharova et al., 2019 showed that over Arctic 210 rivers and lakes, the altimetric backscatter has specific seasonal behaviour. This behaviour is strongly related to hydrological phases, especially, to presence of ice ( Figure 3a). We analysed a seasonal behaviour of the backscatter on the Ob River with respect to geophysical factors that can potentially affect the scattering properties of surface in microwave frequencies. Return signal of nadir-looking altimeter significantly differs from the return signal of sidelooking SAR radars. In nadir, smooth surface produces a higher return echo than rough surface (Ulaby et al.,1986).

215
The altimetric return signal is represented as a waveform ( Figure 3c). The shape of the waveform varies for different types of surfaces and can provide a variety of information useful for 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 reflected from lake ice. The authors interpreted this peak as reflection from the ice surface, while the main waveform peak was considered to originate from the ice bottom. The difference in ice surface -ice bottom 220 altimetric heights estimated from these two peaks was attributed to the ice thickness and used for its retrieval. On many of the radar waveforms extracted over river ice, we also detected this intermediate peak on the waveform leading edge ( Figure 3b). However, considering that the radar echoes over rivers come from very heterogeneous surfaces, we avoided referring this peak to any definitive reflecting boundary and suggested another approach for ice thickness retrievals. Figure 3b demonstrates that during ice growth, the main waveform evolution is related to a 225 decrease in power of waveform main peak. This can be explained by the volumetric scattering/absorption of signal within growing ice. 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 proposed to retrieve the ice thickness by establishing a simple statistical relation between radar backscatter and river ice thickness. This simple approach has already demonstrated good results for ice thickness retrievals in combination with SAR instruments (Unterschultz et al., 230 2009;Mermoz et al., 2014) .
High backscatter values of nadir-looking altimeters are observed when the instrument footprint contains large fraction of calm water (Fu and Cazenave, 1991). On rivers, over large flooded areas, the water surface exhibits a certain roughness owing to turbulent flow and wind impact and, thus, produces lower backscatter values comparing to calm water. On lakes, during the freezing, the presence of floating ice (floes or slush) prevents the development of 235 wind waves and increases the specularity of water (Kouraev et al., 2007). The freezing on rivers starts from formation of a fine skim ice along the banks. This ice can detach and drift (Hicks, 2009). The skim ice is characterised by a smooth surface and smooth bottom. We hypothesize that at the time of appearance of the first skim or drifting ice in late autumn, the peak on the backscatter time series can indicate the start of freezing ( Figure   3b).

240
The scattering and absorption of radar echo within ice depends on the penetration depth of emitted signal and on internal ice properties ( Rémy et al., 2012, Jeffries et al., 2012. In Ku-band the penetration depth into dry freshwater ice is in the order of 5 -12 m and depends on temperature and properties of the ice (Legrésy and Rémy, 1997) . It means that in the case of shallow lake and river ice, the penetration depth for waves emitted at altimetric radar frequencies (C and Ku-bands) is equivalent to the ice thickness. River ice mainly grows as water at the bottom of the 245 ice cover freezes . As ice grows, the volume scattering of radar signal within ice increases resulting in the altimetric backscatter decrease (Brown 1977;Rémy et al., 2012). On the Ob River backscatter time series, a distinct recession limb can be well observed in winter (see the red lines in Figure 3a). The Ob River ice gains approximately 30% of its total thickness during first freezing month ( Figure 2a). The highest temporal decrease in backscatter, ∆Sig0/∆t (where Δt is the time period between two consecutive observations), were 250 observed exactly during that period. During the winter, small peaks on backscatter time series can appear (see VS12 in Figure 3a). They can be explained by temporary changes of ice surface characteristics such as formation of polynya (open water), blowing away of snow, snow wetting during mechanical ice cracking, and occasional snow melt during warm sunny days in early spring. All these processes either change the roughness of the reflecting surface (polynya, snow crust) or decrease the radar penetration depth and, consequently, reduce the signal energy 255 loss within the ice.
River ice breakup is influenced by both thermodynamic and hydrodynamic processes, known as thermal and mechanical breakup, respectively (Ginzburg, 1977). Our studies used SARAL/AltiKa (Kouraev et al., 2015) and Jason-3 (Kouraev et al., 2021) altimeter data demonstrated that in early spring, when air temperatures are still mostly negative, lake ice undergoes to metamorphism under the influence of solar radiation. At that time, a drop in 260 backscatter in the order of 5-10 dB can be observed. Similar process can occur on the Ob River ice in pre-melting period (see, for example, 2012-2016 years in Figure 3a). When air temperatures become positive, snow on the ice surface melts, increasing the surface backscattering of the radar signal. Melting progressively affects the ice, and melt ponds with smooth water surface can appear on the ice surface. The presence of water on the ice surface and wet ice prevent the penetration of radar waves (Ulaby et al., 1986). During the melt, the surface scattering becomes prevailing in altimetric return signal. This, naturally, should result in increase of radar altimetry backscatter as the roughness of the surface of melting ice/melt ponds is low.
Mechanical breakup starts when the water level rises. At that moment, water can flood the ice due to earlier breakup in river headwaters or on tributaries or due to leakage through cracks in weakened/fractured ice (Ginzburg, 1977;Hicks, 2009). A high backscatter peak occurs at the beginning of the flood. The value of the peak varies in large range 25-50 dB. This blocks a development of a threshold algorithm for breakup date retrievals. As the river becomes free of ice, 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. The summer variability in backscatter depends on many factors, including, but not limited to, wind induced roughness, VS location (banks,

275
presence of islands, floodplain characteristics), the proportion of water within the footprint (intermittent summer rain flood inundation).
In our previous studies (Kouraev et al., 2005;Zakharova et al., 2019Zakharova et al., , 2020, we noted that over Arctic rivers, the returned altimetric radar signal (expressed as backscatter) has specific seasonal behaviour. This behaviour is strongly related to the hydrological phases, especially the presence of ice. In contrast to side-looking SAR 280 instruments, for nadir-looking altimetric radars, smooth surfaces produce a higher return echo than rough surfaces (Ulaby et al.,1986, Jeffries et al. ,200512). River ice mainly grows as water at the bottom of the ice cover (called congelation ice) freezes, and the latent heat of 295 crystallisation is conducted upwards through the ice and snow to the atmosphere. Growth can also occur on top of the ice cover when the snow load or hydrostatic pressure is high, and water seeps through cracks, wetting the snow.
The wet snow refreezes forming porous white ice, which is called snow ice. As ice grows and the volume scattering of the radar echo within ice increasesin the backscatter decreases and forms a well-detectable winter recession limb (see the red lines in Figure 3). On the Ob River, the ice gains approximately 30% of its total thickness during the

315
River ice break-up is influenced by both thermodynamic and hydrodynamic processes, known as thermal and mechanical break-up, respectively. First, when air temperatures are still mostly negative, ice undergoes metamorphism under the influence of solar radiation. At that time, a drop in backscatter in the order of 5-10 dB can be observed (see, for example, 2012-2016 years in Figure 3a). Similar phenomena have previously been observed on Lake Baikal ice using SARAL/AltiKa altimeter data (Kouraev et al., 2015). When air temperatures become positive, 320 snow on the ice surface melts, increasing the surface backscattering of the radar signal. Melting progressively affects the ice, and vast melt ponds can appear on the ice surface, leading to a sharp increase in the backscatter coefficient.
Mechanical breakup starts when the water level rises. Water can flood the ice surface due to earlier breakup on tributaries or due to leakage through cracks in weakened/fractured ice. A high (>25 dB) backscatter peak occurs at the beginning of the flood. The value of the peak ranges from to 25-50 dB, depending on the stage of breakup and 325 river morphology (channel width, banks, oxbow lakes). As the water becomes free of ice, the backscatter decreases

Waveform changes
The radar altimetry return signal is represented as a waveform (Figure 4). The shape of the waveform varies for different types of surfaces and can provide a variety of information useful for the interpretation of geophysical processes. When estimating lake ice thickness, Beckers et al. (2017) found the presence of an intermediate peak on

335
the leading edge of the waveform. The authors interpreted this peak asfrom the ice surface, while the main peak was considered to originate from the ice bottom. On many of the radar waveforms extracted over river ice, we also detected this intermediate peak on the waveform leading edge (Figure 4). Considering that the radar echoes over rivers come from very heterogeneous surfaces, we avoid referring this peak to any definitive reflecting boundary and suggest another approach. Figure 4 demonstrates that during ice growth, the main waveform evolution is related 340 to a decrease in power of waveform main peak. This can be explained by the signal volumetric scattering within the growing ice. We noted that the decrease is proportional to the value of the backscatter (which can be seen as the integral under the waveform). Based on this observation, we propose to retrieve the ice thickness by establishing a simple statistical relation between radar backscatter and river ice thickness.

350
Considering the behaviour of the backscatter over frozen rivers described in Section 4, we developed a dynamic algorithm for retrieving ice phenology dates based on the analysis of Sig0 time series specific for each VS. We suggested that the last annual peak of each year in the backscatter time series corresponded to the beginning of river ice formation. In the case of a multi-peaky recession limb, see for example VS12 in 2013 (Figure 3a), we selected Formatted: Font: Times New Roman, 10 pt the peak height of the order of spring-summer peak heights typical for this VS. If the selection of the peak was not 355 straightforward (for example, two high peaks within one month or the prominence of the peak was low), an additional criterion based on the brightness temperature difference between the 34.0 and 18.7 GHz frequencies (∆TB) was introduced. In this cases we selected the backscatter peak at a time t, if in a time frame of (t-1, t+2) satellite cycles at least three of the four ∆TB values are <2° K.
where -number of detected peaks, -number of dTB observation.
Radiometric brightness temperature measurements integrate emissions from larger surrounding areas than altimetric radar backscatter measurements. Freezing of small oxbow lakes on floodplain and soils and bogs on banks usually occurs earlier than in the big channels of the Ob river. By applying the (t-1, t+2) window, we ensured that the 365 freezing in the area of the VS progressed, and the backscatter peak was not caused by a synoptic-scale cooling episode or by calm weather conditions.
The beginning of the spring increase of the backscatter marks the The beginning of the ice cover decay (thermal breakup) leads to marks the beginning of the spring backscatter increase. The breakup detection algorithm searched for the spring peak in the backscatter time series. For multi-peak winters, the algorithm also useds the ∆TB 370 condition. In this case, the algorithm searched for the peak which is accompanied by a simultaneous increase in ∆TB in the order of mean summer ∆TB values for a given VS. 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 maximal increase in backscatter between two satellite cycles (∆Sig0/∆t) for the period from January to mid-June.

Ice thickness algorithm
Year  Hice_ alti= a×abs(Σ(∆Sig0/∆t)) b The coefficients a and b of the equation were estimated for each pair of gauging -VSs from the training set. Using the leave-one-year-out method (Picard and Cook, 1984) for each VS -gauge station pair, we obtained a set of a and b coefficients and estimated their mean values. These average mean values were, then, used for ice thickness estimation on the corresponding training VS. The accuracy of the ice thickness retrievals were evaluated using the correlation coefficient and RMSE calculated between retrieved and observed ice thickness for the 2008-2018 period.

420
Obtained coefficients a and b were applied to other 38 VSs of the main VS set. To do this, we constructed a Sig0 correlation matrix for all VS. Then, for each VS from the main set (VS i ), we used the coefficients a and b of the VS from the training set (VS jt ), which expressed the best correlation between Sig0 i and Sig0 jt . Finally, we used the estimated ice thickness at all 48 VSs to develop weekly ice maps covering the 400 km Low Ob River reach. A schematic representation of the processing steps is presented in Figure 5.

Elaboration of 2D spatio-temporal ice thickness product
Using the coefficients a and b of Equation (3) developed for the VS from the training set, we estimated the ice thickness at 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 at one of the 435 training VS t . Then, the coefficients a and b of VS t showing the highest correlation coefficient were applied to the VS considered.
Ice thickness retrieved at all 48 VSs were used for creation of weekly maps, which were generalised into a 2D spatio-temporal ice thickness product. For this, the altimeter-derived ice thicknesses were interpolated and smoothed in 2D spatio-temporal coordinates using a moving average filter. The size of smoothing window has to a)

Interaction nodes
An important effect of interpolation/smoothing procedure on Hice (i.e. rapid degradation of the R and RMSE statistics) occurred in 15-40 km spatial window range. In windows above 40-45 km the procedure resulted in little subsequent changes in smoothed product. In temporal domain, the application of windows bellow 15 days (for 0 km spatial window) and 40 days (for 60 km spatial window) also had minor effect (the correlation coefficient was above 0.98, RMSE and difference in maximal Hice were below 0.03 m).
The selected window of 40 km/30 days allowed for keeping the average Hice_max difference and RMSE below 0.04 and 0.03 m, respectively. The correlation between Hice time series extracted from smoothed product and retrieved at 20 VSs located on the main branch of the Ob River was 0.99.

Ice phenology algorithm verification
The ice onset was well detected by the Jason radar altimeter. Considering the 10-day repeat overpass of the satellite and the distance between the gauging stations and VS, we considered a 10-day time-step difference (i.e.. ±10 days)

460
as an acceptable accuracy for altimetry-derived ice phenology dates. For 90% of the manual routine retrievals using the visual backscatter time series interpretation (manual routine), the difference in the first ice events with in-situ observations was less than 10 days (Figure 6a). In 56% of the cases, the difference was 0 days. As the radar footprint over rivers is heterogeneous and is affected by signals from the frozen/unfrozen state of land/river/floodplain lakes, there are numerous variations in the behaviour of backscatter at the beginning of the freeze-up period. At this time,

465
the automated routine misses certain behaviour types, and detection is less accurate for the first ice appearance than in the case of the manual routine. Only 70% of the altimetric freeze-up dates fell within 10 days of the in-situ observations at gauges, and only 40% fell on the same day.
Breakup is a more complex process consisting of the thermal degradation of ice cover (breakup start) and its mechanical breakup and downstream movement (breakup end). Comparing the dates of altimetry-derived breakup 470 onset with the ice state flags provided by gauging stations, we can could conclude that the manual routine of our algorithm accurately detecteds the start of ice thermal degradation. In 88% of the cases, the difference between the manually retrieved breakup dates and in-situ observations of the first water appearance was less than ±10 days ( Figure 6b). The automatically derived breakup date estimations were less accurate for the detection of breakup start compared to the manually derived estimations. However, the automated routine was more adapted for the detection 475 of the breakup end than for detection of breakup start; the accuracy of less than ±10 days was achieved for 67% of cases (Figure 6 b, c).
The breakup algorithm was designed for detection of the breakup start, i.e. the wet ice, water-on-ice or open water events. Manual retrieving of breakup dates allowsed better control over the complex variability in of backscatter during the springand, than automated routine did. It is likely that in the complex cases the automated routine 480 detected the breakup end or even provided unrealistic early/late estimates of breakup dates. For example, for the whole set of 48 VSs, over the full 10-winter period of study , the automated detection failed, that is, detected unrealistic breakup dates before 10 April and after 10 June, in 10% of cases. As the manual routine demonstrated better accuracy than automated, it was selected for further analysis of results and for use with the ice thickness retrieval algorithm.

505
The accuracy of the ice thickness retrievals from altimetric measurements was assessed for the ten VSs from the training set located around the gauging stations (Figure 8). At the northern VSs (VSs 161, 138, 237S_Ob, and 249S_Ob), the correlation between retrieved and observed ice thickness was stronger than at other locations, and the errors in estimates of the altimeter-retrieved ice thickness were less than 0.12 m ( Table 2). For the southernmost VSs (VS 12 and 109), the error RMSE increased up to 0.14-0.18. For many years and many locations theThe relative uncertainty in ice thickness estimates was higher at the beginning (low Hice) of the ice period, compared to the Hice retrievals in the middle or end of winter (Figure 8). Inaccurate detection of ice onset can affect the accuracy of thickness estimates, especially at the beginning of the freeze-up. Another reason for the increased errors is the multipeak character of the backscatter winter recession curve and the residual noise that could remain in the backscatter time series after the application of the smoothing procedure. This occurred, for example, at VS12, where a polynya persisted until March in at least four years of the study period, producing a noisy backscatter time series in winter (see Figure 3b) and a high dispersion of points around 1:1 line on the Halti -Hinsitu scatterplot (Figure 8).
Except for VS109, the variability in the values of coefficients a and b in Equation (3) is 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 Hice_alti retrievals to fitting parameters consists 520 of in running a cross-validation experiment, i.g. retrieving the Hice using the coefficients obtained for an adjacent VS. Similar scores between retrieved and in situ ice thickness obtained in validation (  Table 3). fig.1

Ice thickness estimation for the entire studied river reaches
Using the coefficients a and b of Equation (3) developed for the VS 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 (Figure 9). The elaborated maps can be used for evaluation of ice thickness and ice 550 phenology dates in areas between virtual stations. For instance, two useful parameters could be extracted from the 2D product: the maximum ice thickness and ice thickness observed on December 1.
For this, the altimeter-derived ice thicknesses were interpolated and smoothed in 2D spatio-temporal coordinates using a moving average filter with a window size of 40 km/30 days. The size of the applied window allowed for a) preserving as much as possible the magnitudes and the spatial heterogeneity of ice thickness in the spatial domain, 555 and b) to reduce the residual noise in the temporal domain left after smoothing of the backscatter time series with the LOESS filter.
The along-river variability in ice thickness, controlled by local morphological factors, can be important. In the absence of validation data for the retrieved ice thickness for inter-station areas, we suggest examining the interannual dynamics of two parameters derived from altimetric and in-situ observations: the maximum ice 560 thickness and ice thickness observed on 1 December. From a practical standpoint, knowledge of the maximum river ice thickness is relevant for hydro-climate change monitoring, while the ice thickness determined on 1 December 1 is crucial for local and regional socioeconomic stakeholders, as this is the average date for the opening of the ice bridge road to the north of the study area at Salekhard. To assess the goodness of the 2D product for practical use, we compare the interannual dynamics of the mentioned parameters derived from this product and observed on 565 gauging stations.
For ~90% of the area of the studied river reaches, the interannual variability in maximum ice thickness retrieved from altimetric measurements at more than 90% of the VSs indicates indicated a clear decrease from 2008 to 2012 ( Figure 9). This tendency corresponds corresponded well to those observed at all gauging stations. (Figure 10a).
Since 2013, the maximum ice thickness has increased slowly. However, altimetric and in-situ observations after 570 2013 both exhibit spatio-temporal variability that is 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 station records, ridging events appear more frequently after 2012 at the northernmost gauging station, Pitlar, than at other gauging stations. We cannot fully exclude the effect of the simplicity of the retrieving algorithm, based on the empirical approach, as well 575 as the effect of spatio-temporal smoothing of the altimetric retrievals used in map production.
The difference between altimetric and in-situ ice thickness on December 1 observed in their interannual variability lay within algorithm uncertainties 0.07-0.18 m (Figure 10b). In addition to the geophysical reasons and algorithm simplicity noted above, the observed difference can be related to the degradation in the quality of the in-situ time series (gaps, unrealistic values) and the low representativeness of the one-hole sampling protocol.

Winter ice bridge roads 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 in summer are linked to supply centres only via expensive helicopters or boat transport, have the opportunity to directly access the main land transport arteries using frozenground and ice roads. The importance of ice roads is highest for the Arctic regions, where the construction of 595 permanent bridges is restrained by the presence of permafrost and its destabilisation. Using suggested criterion, we could predict the beginning of ice road traffic with an average accuracy of four days.
However, in half of the years, the predictions differed from observations by more than five days (11 days maximum). As noted in Section 6.1, at the beginning of freezing, the errors of Hice retrievals are quite high, and ice thickness is often overestimated (see Fig 10b). Based on this fact, we consider that at the moment the altimetric algorithm and the ice thickness product are not sufficiently accurate for the ice road opening forecast. Nevertheless, their accuracy is sufficient for the assessment of climatic perspectives as we capture the interannual variation of dates of ice road opening (Figure 11a). For the prediction of dates, when the ice road at Salekhard ceases its operation, the use of the northernmost VSs 630 iswas not possible as traffic on the ice road closes before the altimeter detects the breakup onset in this reach.
However, the satellite observations of the breakup onset at VSs located in the upper reaches can provide necessary information. Using 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 entire 400 km river reach.
This date servesd as a predictor of the date of ice-road closure at Salekhard. The correlation between AMO2 and 635 observations was significant (p-value = 0.025) at correlation coefficient 0.70 (Figure 12a). After applying a correction to AMO2 computed from the relationship in Figure 12a, one can obtained a forecast date of ice road closing at Salekhard and evaluate a residual difference between forecasted and observed dates. This difference, expressed as with an RMSE was equal toof 3 days (see Figure 11b). The leading time of the forecast can be determined from the plot shown in Figure 12b. It varied from 4 days (for late breakup start) to 22 days (for yearly early breakup start).

Factors affecting ice thickness retrievals from altimetry
Various factors can affect radar signal return echo and, consequently, the accuracy of river ice thickness retrieval.

650
One source of uncertainty could originate from the underestimation of the role of snow-on-ice in microwave signal scattering (King et al., 2013. However, Willatt et al. (2011) demonstrated that the Ku-band electromagnetic wave scattering by snow at nadir is low, and in our study, we neglected the presence of snow on ice. Using precipitation data from the nearest meteorological station, we noted that not all heavy snow accumulation episodes affected the backscatter over river ice. In several cases, snowfall resulted in backscatter changes in the order of 1.5 655 dB. The smoothing procedure applied to the cumulative ∆Sig0/∆t series helped to eliminate this effect. Moreover, after adding the in-situ snow depth to the ice thickness, we noted that the Hice -Σ(∆Sig0/∆t) power relationship becomes weaker.
Another factor potentially affecting the backscatter value over freshwater ice is ice roughness at the ice/water interface . The roughness of the ice bottom on the rivers is expected to be high at the beginning 660 of the freeze-up period in bridging areas, where floes juxtapose and accumulate underside. Any rough boundary dissipates the signal of nadir-looking radar instruments, resulting in a decrease in backscatter. Further congelation of the inter-floes volume as well as ice growth both lead to levelling of the ice low boundary. We suppose that the under-ice current may also contribute to the ice bottom levelling in 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 Figure 3a). Either the levelling effect is weak and is masked by high volumetric scattering of the radar echo within thickening ice, or at the location of our VSs, the ice juxtaposing (and consequently the ice bottom roughness) is unimportant. Future investigations with dedicated in-situ observations on river ice texture evolution are needed to understand the effect of bottom roughness on radar altimeter return signals. Slater et al., 2019). Under the climatic conditions of north of wWestern 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 portion of the study area and 1-2 weeks later to the north. During this time of the year, the ice was well developed and almost reached its maximum thickness. Layering can also occur after river water floods the ice surface through 675 cracks. According to in-situ observations at gauging stations, this phenomenon was observed in the last several years at the end of the ice season in the southern portion of the study area. Both warming episodes and flooding events lead to a backscatter increase in the order of 1-5 dB and render altimetric ice retrievals difficult by the end of the ice season. The highest underestimation of Hice (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., 680 2018). During ice formation, jamming and ridging, can occur in Arctic rivers, resulting in the formation of air pockets. Ridging is rare near gauging stations on the Ob River. However, there is no information about the state of the ice at other ungauged reaches, including the areas covered by VSs. We can only speculate that ridging/hummocking could be one of the reasons for the high difference in the coefficients of Equation (3) determined for VS109 and for lower accuracy of the Hice retrievals at this station comparing to northern virtual 685 stations. On Landsat-8 images acquired for two springs on 2019/04/25 April and 2015/05/01 (not shown), the irregular spatial ice structure in the area of VS109 indirectly confirms our hypothesis. More studies involving the simultaneous analysis of SAR imagery showing ice field irregularities (Unterschultz et al., 2009) and altimetric signals could help to clarify this issue.

Potential improvements of algorithms
The results obtained demonstrate that the altimetric retrievals of ice phenology dates and thickness are capable of representing the interannual variability of these parameters and can be potentially used in climate studies (see Figures 7 and 10a). In spite of potential for climate change monitoring, we cannot recommend the use of current first version of ice thickness altimetric product for winter road opening forecasts, as it seems that for many locations, 695 we overestimate the thickness at the beginning of the freeze-up period (see Figure 10b). Inaccurate detection of ice onset can affect the accuracy of thickness estimates, especially at the beginning of the freeze-up. Another reason of uncertainties is the multi-peak character of the backscatter winter recession curve and the residual noise that could remain in the backscatter time series after the application of the LOESS filter. This occurred, for example, at VS12, where a polynya persisted until March in at least four years of the study period, producing a noisy recession curve 700 (see Figure 3a) and a high dispersion of points around 1:1 line on the Halti -Hinsitu scatterplot (Figure 8).
Several improvements could be suggested for future studies. First, an improvement in the accuracy of detection of the freeze-up algorithm is envisaged. In our algorithm, the ice thickness estimation starts from the date of the first ice (bank ice or frazil floes) appearance. Usually, the river reach in the area of the VS at this moment is not fully frozen. The detection of the date of the first consolidated ice (i.e. fully frozen reach) may help to reduce the errors in 705 the retrievals in a low range of ice thickness.
Another improvement consists of in the use of other parameters of the altimetric radar waveform instead of (or in complement to) the backscatter coefficient. As shown in Figure 3b, the main peak of the radar waveform decreases as ice grows. We suppose that the amplitude of the main peak or the area under this peak may produce stronger variable part of the waveform. This peak corresponds to the return signals from the ice surface to the ice bottom (Backers et al., 2016). The other part of the waveform collects 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 surrounding areas that are different from the state in the near-nadir area (e.g. ice flooding or off-nadir polynya) lead to backscatter changes not related to ice growth. Consequently, the use of waveform parameters related to the main 715 peak could reduce the errors in the Hice estimates. Unfortunately, these parameters are not directly provided in the AVISO+ Jason GDR product, but they could potentially be estimated from the initial waveforms.

Potential radar altimetry application for other rivers
The main advantages of the satellite altimetry comparing to SAR instruments consists in easy processing of satellite 720 measurements in large (hemispheric scale) spatio-temporal domain. As the developed ice phenology dates algorithm is independent of availability of in situ observations, it can be directly applied for other rivers. The main limitation of algorithm could be a width of river, which, from our experience, is limited by 100-200 m. The algorithm can be also applied for lakes and bogs, as the temporal variability of the backscatter over these objects is similar to that of described in the section 3. However, the algorithm needs an additional adjustment/verification for rivers with 725 unstable winter ice cover, as well as in areas of frequent occurrence of rain-on-snow episodes.
Unfortunately, similar to SAR-based methods (Unterschultz et al., 2009;Mermoz et al., 2014), for our ice thickness algorithm, the ground trust data are indispensable. Nevertheless, we demonstrated, that for the sufficiently long river reach, the obtained empirical relations between in situ data and satellite measurements are quite close each to other.
We suppose that the use of ∆Sig0/∆t instead of absolute value of Sig0, allowed for minimisation of perturbing effect 730 of surrounding banks (i.e. local morphological conditions). The character of surrounding banks affects the values of the backscatter in the beginning of freezing (i.e. the initial conditions). The effect of the initial conditions, was, thus, reduced. The unification of the relations and its application for other reaches and other rivers is likely possible, when using ∆Sig0/∆t. However, further investigations are needed for assessment of the effect of ice texture/type on equation parameters. Hence, the combination of SAR (providing the ice type) and altimetric instruments could be a 735 beneficial.

Conclusion
The decreasing number of in-situ observations and degradation of the quality of their time series is a good argument 740 for boosting the development of satellite methods for freshwater ice monitoring. The present study demonstrates the potential of satellite radar altimetry for monitoring river ice parameters such as freeze-up, breakup, and ice thickness for the large Arctic Ob River.
An algorithm based on the analysis of backscatter coefficients from the Jason-2 and -3 satellite altimeters provides an estimation of river ice onset with an accuracy of ±10 days (corresponding to the 10-day satellite overpass frequency) in 90% of the cases.
River ice breakup consists of two phases: thermal degradation and mechanical breakup. The algorithm detects the beginning of thermal degradation well with the same accepted accuracy of ±10 days for 88% of the cases. River ice thickness was retrieved from altimetric backscatter measurements via simple empirical relations with in-situ observations. The accuracy of the thickness retrievals (expressed as RMSE) ranges from 0.07 to 0.18 m.

750
The spatio-temporal interpolation and smoothing of satellite-derived river ice thickness retrieved for 48 VSs along the 400 km reach of the Lower Ob River allowed for the elaboration of weekly maps generalised in the form of an annual spatio-temporal product. The ice thickness time series can be extracted for any location and used for climate change monitoring and ice road operational purposes.
Using this first version of the product, we demonstrated that the dates of opening of ice road near Salekhard City can 755 be predicted from altimetric ice onset retrievals with an accuracy of 4 days. Errors in the prediction of dates of ice road closure were within 3 days. Despite these promising results, we consider that the current first version of the product is not sufficiently mature for operational use as it overestimates 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 and/or SAR imager) 760 approach. We are hopeful that with the use of the Copernicus satellite altimeters Sentinel-3A and 3B, an improvement in the retrieval of ice thickness can be made. 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 days) is not suitable for operational applications, they provide five overpasses within a 25 km distance around the Salekhard city ice road and, thus, the temporal resolution of observations may be significantly improved.
The combination of data from the Jason and Sentinel-3 missions could be fruitful.
The Salekhard ice road is very well instrumented, monitored, and maintained by local authorities, thanks 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 770 crossings are developed each year by local people. Often, the lack of information on the state of ice results in accidents and requires intervention by the emergency service. The demonstrated capacity of the first version of the altimetric river ice product as a supporting tool for the ice road operation to the north ofon the Ob River is quite promising. Further product improvements and the development of sophisticated prediction criteria for road operation that could be adapted to other reaches of the Ob River are planned.
Author contribution. All authors contributed to the data collection, algorithm development, analysis and presentation of results equally.