Estimation of River Ice Thickness in the Shisifenzi Reach of the Yellow River With Remote Sensing and Air Temperature Data

River ice is an important part of the cryosphere, and effective monitoring of ice thickness information is essential for river ice research. The task of preventing and controlling winter ice damage in the Yellow River is severe, and it is significant to explore the applicability of different ice thickness estimation methods in this region. There is a lack of inversion studies on the thickness of the Yellow River ice. In this study, a typical reach of the Yellow River was selected as the research area. Based on the measured ice thickness data for many years, the applications of river ice thickness estimation models based on remote sensing and air temperature were explored respectively in this area. The results show that the VV polarized backscatter coefficient of Sentinel-1 SAR data has the highest correlation with the measured ice thickness, with a Pearson correlation coefficient of 0.702. SAR ice thickness inversion results are helpful to identify ice jam location. The simulation accuracy of the air temperature model for ice thickness is better than that of the SAR model. Long-term air temperature-based ice thickness estimates show the maximum ice thickness in the study area is thinning at a rate of –0.7 cm a decade.


I. INTRODUCTION
R IVER ice is an important part of the cryosphere [1], [2], The formation and break-up of river ice has an important impact on river morphology, bank stability, winter water supply, and hydraulic engineering [3]. River ice enriches human winter activities, such as ice-skating, ice fishing, and ice hockey [4], [5]. Ice floods caused by river ice seriously threaten the safety of people's lives and properties in some areas [6], [7], [8]. Whether it is ice use or ice disaster prevention, ice thickness is a key parameter and a difficult variable to monitor, and how to estimate it accurately and effectively has been the focus and difficulty in river ice research. The Inner Mongolia section of the Yellow River will occur in varying degrees of ice disaster every winter and spring [9], [10]. Since 1986, 13 dike breaches have occurred in the Inner Mongolia section of the Yellow River, resulting in huge economic losses and social impacts [9]. Therefore, effective monitoring of ice thickness information in this area is of great significance for the prevention and mitigation of disasters in the Yellow River during winter. At present, manual drilling in the field is the most reliable method to measure ice thickness [11]. However, it is timeconsuming and hardworking, and limited by the natural environment and transport conditions. Therefore, the measured data are often used to verify other technical methods [11]. Ice thickness variation is very sensitive to air temperature [12]. The degree-day method ice thickness estimation model has been widely used in the field of ice engineering since it is proposed by Stefan [13], [14], [15]. This method can reflect the trend of river ice formation and elimination well and refine the change of ice thickness. But it is difficult to reflect the spatial distribution of ice thickness.
The progress of remote sensing technology has promoted a large number of satellite-based earth observation studies. While optical remote sensing mostly uses the spectral characteristics of features to explore surface changes [16], [17], [18], the penetrating capacity of radar remote sensing on objects makes it increasingly used in the study of snow and ice monitoring [19], [20], [21]. With the development of radar technology, ground-penetrating radar [19], [22], airborne radar [6], [23], and satellite-based radar [20], [21], [24] all have been used for ice thickness monitoring. The first two are close-range measurements, mostly using nonimaging radar. the radar echo signal can distinguish air-ice and ice-water interfaces accurately, and calculate ice thickness through the propagation time of electromagnetic waves in different media. The satellite-based radar is long-range measurement, altimetry satellites such as ICEsat-2 [25] and CryoSat-2 [26] are mostly used for monitoring polar sea ice thickness. But the lower spatial resolution limits its application in river ice thickness monitoring. Many studies show that synthetic aperture radar (SAR) satellite with high spatial resolution imaging has great potential in river ice thickness inversion [21], [27], [28], [29]. Research on river ice thickness inversion has mainly focused on Canada [21], [28], [29], and a recent study used Sentinel-1 SAR data to invert the ice thickness of high-order rivers on the Tibet Plateau [24].While work on river ice thickness inversion of the Yellow River is extremely lacking. This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Generally speaking, the radar backscattering coefficient depends on the sensor properties, the dielectric properties, and the scattering properties of the imaging target [30]. The sensor characteristics of SAR are mainly determined by its wavelength, polarization mode, and incident angle. The first two are fixed for the sensor, while the backscattered signal from the same target varies with its microwave incidence angle [31]. The dielectric properties of a material are characterized by its dielectric constant. For freshwater ice, this constant hardly varies with ice temperature, unlike for sea ice [32]. The electromagnetic loss caused by the propagation process of electromagnetic wave in the medium is determined by the dielectric constant of the medium, low dielectric constant means low moisture content. For example, freshwater ice is a low-loss medium, and electromagnetic waves with frequencies in the range of 1-10 GHz have a penetration depth of 100-10 m for freshwater ice [33]. This gives the possibility of microwave monitoring of ice thickness changes. On the contrary, when the dielectric constant of the medium is large, the electromagnetic loss is high and the penetration depth is small. For example, liquid water is a high loss medium, when the water content in snow or ice increases, the penetration depth of electromagnetic waves will be reduced. For the 5 GHz C-band electromagnetic waves, only 5% of the volume of water content in the snow will reduce the penetration depth to about 10 cm [33]. At this point, backscattering records mainly snow information, not ice [34].
The scattering properties of river ice are mainly determined by two aspects: surface scattering and volume scattering. Surface scattering is dominated by specular reflection and diffuse reflection. Specular reflection occurs on smoother surfaces, which can be considered smooth within the incident wavelength scale of microwaves [30]. For example, a calm surface of water specular reflects the energy of an incident microwave away from the direction of microwave incidence, so the backscattering response is weaker and appears darker in radar images. While diffuse reflection occurs on relatively rough surfaces, the rough surface reflects incident microwaves in all directions and reflects some of that energy back to the sensor, so the rough ice surface will appear relatively bright on the radar image. Volume scattering occurs when incident microwaves penetrate the ice surface or dry thin snow and are repeatedly reflected at the dielectric discontinuity of the ice. The main factors affecting volume scattering in river ice are cracks, bubbles, and impurities. Volume scattering usually leads to an increase in radar backscatter [27], [35]. Recent modeling studies of electromagnetic radiative transfer from river ice have found that surface scattering at the ice-water interface makes a large contribution to backscattering [36], [37], [38]. However, the Yellow River is the most sediment-rich river in the world [39], and unlike other freshwater ice, it is unclear whether the body scattering caused by sediment impurities in its ice would result in stronger radar echoes. Therefore, the Yellow River should not be ignored in the research field of river ice thickness inversion based on spaceborne SAR.
Sentinel-1 remote sensing data, as the most widely used free satellite-based SAR data, is an effective data source for river ice monitoring. In this study, measured ice thickness data with spatial location information in 2019-2021 and ice thickness data published by hydrology stations at different times during the ice flood season of 2021-2022 were used to explore the applicability of ice thickness estimation models based on remote sensing and air temperature in the Inner Mongolia section of the Yellow River. First, the dual-polarized C-band Sentinel-1 data are used to explore the microwave inversion of the Yellow River ice thickness, which fills the gap of remote sensing inversion study of the Yellow River ice thickness. Then, the Yellow River ice thickness is estimated using the air temperature method and compared with the SAR method. Finally, the long-term trend of Yellow River ice thickness is predicted based on the historical temperature data. The objectives of this study are as follows.
1) Analyze the correlation between river ice thickness and SAR polarization parameters, and construct the SAR ice thickness inversion model to estimate river ice thickness. 2) Estimate river ice thickness using air temperature model. 3) Estimate the long-term trend of river ice thickness.

A. Study Area
The Inner Mongolia section of the Yellow River is located at the northernmost part of the Yellow River, the total length of about 830 km, between 39°∼ 41°N, 106°∼ 112°E, and an altitude of about 1000 m. The river flows from west to east, passing through four hydrological stations in turn: Bayangaole, Sanhuhekou, Baotou, and Toudaoguai (see Fig. 1). Due to its special geographical location and complex meteorological conditions, it leads to the order of freezing of the river from downstream to upstream during the ice period, and the opposite when the river is thawing, which makes it easy for ice flooding disasters to occur. The annual ice flood season usually starts from mid-late November and ends in mid-late March of the following year, lasting more than 100 days.
The Shisifenzi reach is located in the downstream of the Inner Mongolia section of the Yellow River in Tuoketuo County, with an "Ω" shaped bend and a river width of about 200-600 m. From November 1, 2021 to March 31, 2022, the average daily discharge in the study area was 602.8 m 3 /s, the average daily water level was 987.7 m, and the average daily temperature was -4.34°C. For many years, the bend has been a priority location for freezing of the Inner Mongolia section of the Yellow River (see Fig. 1). During the ice flood season, due to the narrow river channel at the bend, coupled with the growth of bank ice narrow river width, the river's ice transport ability significantly reduced. Ice flow and ice flowers are very easy to accumulate here to form ice jam, ice dam, thus causing high upstream water level, when serious can lead to the occurrence of ice flooding disaster. Therefore, this section of the river is a key section for ice disaster prevention and control, and the acquisition of ice thickness information is of great significance for the prevention and reduction of the Yellow River ice disaster.

B. Filed Data
The field ice data collection site is located at the top of the bend in the Shisifenzi reach. The ice cover is thick during the  freezing period in the region, which was convenient for field investigation. The ice thickness is clearly different from the left bank to the right bank. River ice thickness information was collected here for three consecutive years, with collection dates from January 10 to 13, 2019, January 17 to 19, 2020, and January 12 to 15, 2021. According to the water resources industry standard of the People's Republic of China "Technical standard for observations of ice regime in open channels" for measurement sections layout. For safety reasons, the thickness of the ice near the open water surface was not measured manually. In this study, 259 effective ice thickness sampling points from 13 sections were selected (see Fig. 2). The sampling points were located using Global Positioning System Real-time kinematic equipment with centimeter-level positioning accuracy. The thickness of river ice at sampling sites was measured using an ice measuring ruler with millimeter-level accuracy. The ice thickness data for the 2021-2022 ice flood season came from the Toudaoguai hydrological station, where staff collected a total of 28 days of average ice thickness for the cross section (see Fig. 3).

C. Remote Sensing Data
Sentinel-1 is the first mission of the Copernicus program performed by European Space Agency (ESA), which is a constellation of two satellites, A and B, in the same orbital plane. A and B stars were launched in April 2014 and April 2016, respectively, with a single-star re-entry period of 12 days and a combined double-star re-entry period of 6 days, both carrying C-band  SAR with a wavelength of about 5.6 cm. Four data collection modes are supported: strip map mode, interferometric wide swath (IW) mode, extra wide swath mode, wave mode, Level 1 products include single-look complex (SLC) and multilook detected dual-polarization products. Only Sentinel-1A passes over the study area, so the re-entry period is 12 days, and these images cover the same area, with the same beam position and orbit.
Sentinel-2 is the second mission of the Copernicus program performed by ESA, which is also currently a constellation of two satellites, A and B, in the same orbital plane. A and B stars were launched in June 2015 and March 2017, respectively, with a single-star re-entry period of 10 days and a combined double-star re-entry period of 5 days, both carrying multispectral imagers. It can cover 13 spectral bands, with an image width of 290 km and a spatial resolution of up to 10 m.
In this study, 75 images of Sentinel-1 IW model SLC dualpolarization products were used to obtain the polarization information and analyze the relationship with river ice thickness. In addition, 6 images of Sentinel-2 products were used to extract the river channel boundary information. The satellite remote sensing data were downloaded from the official website of ESA (https://scihub.copernicus.eu/dhus/#/home/), and the detailed acquisition dates are shown in Fig. 4.

D. Air Temperature and Hydrological Data
The air temperature and hydrological data of the 2021-2022 ice flood season in the study area were obtained from Toudaoguai Hydrological Station (40.27°N, 111.06°E, about 4 km from the study area) (see Fig. 1). The information collected at this hydrological station includes daily maximum air temperature, daily minimum air temperature, daily average air temperature, water temperature, ice thickness, water level and discharge (see Figs. 3 and 5). In order to extend the temporal range of the study, day-by-day air temperature data from 1956 to 2022 were obtained for the meteorological station at Baita, Hohhot (40.85°N , 111.82°E, about 92 km from the study area) (see Fig. 1), where data for 1964-1972 were missing. Combined with a thermodynamic temperature model to estimate the long-term maximum ice thickness trends in the study area. Meteorological station data from the National Oceanic and Atmospheric Administration (https://www.ncei.noaa.gov/data/global-summary-ofthe-day/access/).

A. Ice Thickness Estimated With SAR
In order to explore the potential of Sentinel-1 SAR image parameters to invert the Yellow River ice thickness, three Sentinel-1 images with similar ice thickness collection times in the field were selected, specifically on January 8, 2019, January 15, 2020, and January 9, 2021 (see Fig. 4). The geographical locations of ice thickness sampling points were spatially matched with remote sensing images to extract the polarization parameters of the images, and analyze the correlation with ice thickness. Data processing is performed on Sentinel Application Platform, an open-source software developed by  ESA. The filter "Refined Lee" [40] is used to remove speckle noise from SAR images. Precision orbit data for orbit correction from ESA (https://scihub.copernicus.eu/gnss/#/home). Doppler topography correction was performed using shuttle radar topography mission digital elevation model data with a spatial resolution of 30 m published by the National Aeronautics and Space Administration (https://search.earthdata.nasa.gov). Five polarization parameters were extracted from Sentinel-1 image, namely VV polarization backscattering coefficient, VH polarization backscattering coefficient, entropy, scattering angle alpha, and anisotropy. The specific process is shown in Fig. 6. For the backscattering coefficient, VV polarization means vertical emission and vertical reception of electromagnetic wave, and VH polarization means vertical emission, and horizontal reception of electromagnetic wave. The backscattering coefficient is usually expressed by sigma (σ). VV polarization backscattering coefficient and VH polarization backscattering coefficient are expressed as σ VV , σ VH , respectively, in the following.
The H-Alpha dual-polarization decomposition algorithm is modified from the H-A-Alpha full-polarization decomposition algorithm [41], [42], and has been applied in different research areas [24], [43], [44]. With the scattering vector k, the coherence matrix T can be solved from the target scattering matrix S, then the eigenvalues and eigenvectors can be extracted, and the entropy, anisotropy and average scattering angle can be where L is the number of looks; k i is the scattering vector of the ith look; the superscript * indicates the complex conjugate transpose; q is the number of polarization channels, here q = 2; λi (i = 1, 2) is the eigenvalue of T , satisfying λ1 ≥ λ2 ; u i is the eigenvector of T ; α i and β i are the scattering angle of the eigenvector u i ; ϕ i and δi are the phases of the target. The entropy (H) characterizes the degree of randomness of the target and ranges from 0 to 1. The lower the entropy, the more certain the target is, and the higher the entropy, the greater the randomness of the target. It can be solved using the eigenvalues of the coherence matrix T where Pi denotes the contribution probability of the ith eigenvalue.
The average scattering angle (α) describes the main scattering mechanism of the target and ranges from 0°to 90°, which can well reflect the scattering characteristics of the target and has rotation invariance. α close to 0°indicates surface scattering, close to 45°indicates volume scattering, and close to 90°indicates dihedral angle scattering For fully polarized data, the randomness of the target is large when the entropy is very high (H > 0.7), making it difficult to determine the target type. The anisotropy (A) is then defined as a complement to the entropy such that the target remains distinguishable when H > 0.7. The anisotropy can be calculated using the second and third eigenvalues [42]. For dual-polarized data, the first and second eigenvalues are used to calculate the anisotropy In this study, the statistical regression method was used to establish the ice thickness inversion model. After analyzing the correlation between ice thickness and polarization parameters, the polarization parameters with the highest correlation with ice thickness were selected to establish the linear regression model. From the spatio-temporally matched ice thickness data, 70% were randomly selected as model construction data for determining the pending coefficients of the model, and the remaining 30% were used as model validation data. The coefficient of determination R 2 , root mean square error (RMSE), mean absolute error (MAE), and mean relative error (MRE) are used to characterize the model accuracy.
In order to eliminate the interference of land area on river ice information extraction, the river boundaries were extracted based on Sentinel-2 optical remote sensing data using modified normalized difference water index (MNDWI). The calculation formula is as follows: Where Green is the reflectance of the green band; and SWIR is the reflectance of the short-wave infrared band. After extracting the MNDWI, a suitable threshold was selected by visual interpretation so that the water body was clearly distinguished from the land. Then, binarization and vectorization were performed on the processed images, and the results were refined and modified to obtain the boundary information of the main channel of the Yellow River.

B. Ice Thickness Estimated With Air Temperature
In 1890, Stefan proposed the classic degree-day ice thickness estimation model [13], which simulates ice thickness growth with temperature. The model assumes that the heat released by freezing at the bottom of the ice is transmitted through the ice under a constant temperature gradient, which lays the foundation for the study of the thermodynamics of ice cover to some extent. Stefan's model does not take into account the thermal inertia, the internal heat source and the heat flux from water to ice, and replaces the temperature at the top of the ice layer with the air temperature, the heat transfer equation for ice can be simplified as follows [45]: the two ends of the equation represent the heat flux from ice to air, W·m -2 ; h is the thickness of ice, m; ρ i is the density of ice, kg·m -3 ; Li is the latent heat of melting of ice, J·kg -1 ; t is the time, s; ki is the thermal conductivity of ice, W·m -1 ·°C -1 ; Ttop is the upper surface temperature of the ice,°C, the air temperature is usually used instead; Tbottom is the lower surface temperature of the ice,°C, i.e., the freezing point. Note: 1W = 1J · s −1 .
Assuming t = 0 s, h = 0 m, T top < 0°C as the initial condition, the analytical solution of (8) is as follows: where a 0 is the ice thickness growth factor, m·°C -1/2 ·s -1/2 ; S is the freezing index,°C·s. In practical river ice or lake ice studies, cumulative freezing degree days (CFDD,°C·d) is commonly used instead of the freezing index S. CFDD is the accumulation of the difference between the daily average air temperature and the freezing temperature of river ice in the daily scale time range Then, the formula for ice thickness h can be simplified as follows: where for freshwater ice, 917×3.34×10 5 = 3.52 × 10 −2 (m·°C -1/2 ·d -1/2 ), the value is slightly larger than that under sea ice [45]. Stefan's model can only simulate the growth process of ice thickness, Bilello proposed to use cumulative thawing degree days (CTDD,°C·d) to simulate the decay process of ice thickness [14]. CTDD is the accumulation of the difference between the daily average air temperature and the melting temperature of river ice in the daily scale time range. Therefore, the formula for calculating the ice thickness during the melting period of river ice can be expressed as follows: where hmax is the maximum ice thickness during the ice period, m, which can be obtained from (13); a 2 is the ice thickness decay coefficient, m·°C -1 ·d -1 ; which can be solved by substituting the maximum ice thickness hmax, the river ice thickness of break-up date h0 and the corresponding CTDD into the above equation.

C. Evaluation Methods
Pearson correlation coefficient (R), coefficient of determination (R 2 ), RMSE, MAE, and MRE were chosen to evaluate the effectiveness of the two ice thickness estimation methods.
R is used to reflect the degree of linear correlation between two random variables, with R greater than 0 indicates a positive correlation and less than 0 indicates a negative correlation. The absolute value of R ranges from 0 to 1, and the larger the value, the stronger the correlation, as shown in Table I. The formula for calculating R is as follows: where n is the number of samples;x is the mean value of variable x; andȳ is the mean value of variable y. R 2 , also known as the degree of fit, is used to determine the degree of fit of the regression equation and has a value between 0 and 1. The greater the degree of fit, the higher the degree of explanation of the dependent variable by the independent variable, calculated as follows: whereŷ i is the predicted value; y i is the observed value, andȳ i is the average of the observed values.
RMSE is a measure of the error in the regression equation, which measures the deviation between the predicted and observed values and is calculated as follows: MAE is the arithmetic mean of the absolute error between the predicted and observed values, which can better reflect the actual situation of the error of the predicted values and is calculated as follows: MRE is the arithmetic mean of the ratio of the absolute error to the observed value, which better reflects the confidence level of the predicted value and is calculated as follows: IV. RESULTS

A. Estimation of Ice Thickness by SAR
First, a correlation analysis between ice thickness and polarization parameters was performed for all sampling sites in 2019-2021. It is found that both α and H are negatively correlated with ice thickness, and the correlation coefficient between α and ice thickness is −0.386, showing a weak correlation, while the correlation coefficient between H and ice thickness is −0.423, showing a moderate correlation. A, σ VV and σ VH are all positively correlated with the ice thickness, and the correlation coefficients with the ice thickness are 0.400, 0.483, and 0.411, respectively, showing moderate correlation.
In order to further analyze the correlation between ice thickness and polarization parameters, the sampling points of ice thickness with or without snow cover were distinguished. For the ice without snow cover (see blue dots in Fig. 7), α and H are negatively correlated with ice thickness, and the correlation coefficients with ice thickness are −0.579 and −0.591, respectively, showing a moderate correlation. A, σ VV and σ VH are all positively correlated with ice thickness, and the correlation coefficient between A and ice thickness is 0.591, showing a moderate correlation; the correlation coefficient between σ VV and ice thickness is 0.702, showing a strong correlation; the correlation coefficient between σ VH and ice thickness is 0.528, showing a moderate correlation. For the ice covered by snow (see red dots in Fig. 7), α and H are negatively correlated with ice thickness, and the correlation coefficients with ice thickness are −0.173 and −0.175, respectively, showing a very weak correlation. A, σ VV and σ VH are all positively correlated with ice thickness, and the correlation coefficient between A and ice thickness is 0.201, showing a weak correlation; the correlation coefficient between σ VV and ice thickness is 0.437, showing a moderate correlation; the correlation coefficient between σ VH and ice thickness is 0.123, showing a very weak correlation.
In total, the thickness of the ice with snow cover is poorly correlated with each polarization parameter, and the scatter distribution is quite scattered. It rained and snowed for two consecutive days from January 17 to 18, 2020. It should be due to the rainfall which caused the dry snow above the ice to become wet snow and limited the penetration ability of electromagnetic wave [33], thus affecting the correlation between ice thickness and polarization parameters. After excluding the influence of wet snow, the inversion model of ice thickness built is shown in Fig. 8(a), and the formula is as follows: where x denotes σ VV in dB and y denotes simulated ice thickness in cm. The model indicates that for every 1 dB increase in σ VV , the thickness of the ice increases by 4.591 cm. The R 2 of this linear regression model was 0.525. The RMSE, MAE, and MRE of the validation data were 11.75 cm, 9.70 cm, and 16.17%, respectively [see Fig. 8(b)]. The SAR ice thickness inversion model constructed in this article was applied to the radar images of the six ice flood  seasons from 2016 to 2022 to obtain the river ice thickness in the Shisifenzi river section. Fig. 9 shows the spatial distribution of river ice thickness during the 2021-2022 ice flood season. Ice thickness was very thin in November, increased since late November, reached a larger value by early January, was thickest in January, decreased in February but not much, and decreased rapidly in March until it melted completely. The bend was the first to freeze during the freezing of the river ice. On December 23, there was ice jam formation at the exit of the bend, which was due to the narrow river channel at the bend, the ice transport capacity of the river was significantly reduced, and the flow ice and ice flowers were very easy to accumulate here to form ice jam, and to develop gradually upstream to reach the mechanical equilibrium conditions when a stable ice jam was formed [46].  In November and March, some water bodies were mistaken for thin ice, and there was a significant overestimation of the ice thickness inversion results. Fig. 10 shows the results of the inversion of the average ice thickness in the study area for the six ice flood season from 2016 to 2022, where the changes in ice thickness basically follow a process of first increasing, then remaining essentially constant, and finally decreasing.

B. Estimation of Ice Thickness by Air Temperature
Based on the air temperature model estimating ice thickness, the daily average air temperature monitored by Toudaoguai hydrological station was used to simulate the river ice thickness in the winter of 2021 to 2022 in the study area. When simulating the growth process of ice thickness, CFDD was calculated from the 1st day when the daily average air temperature was less than 0°C for five consecutive days in order to exclude as much as possible the advancement of the simulation date due to the sudden change of air temperature within a short period of time. Similarly, in simulating the decay process of ice thickness, CTDD is calculated from the 1st day when the daily average air temperature is greater than −5°C for five consecutive days. Due to solar radiation and hydraulic factors, the thickness of river ice has begun to decay before the average daily air temperature rises to 0°C.
When the ice thickness growth coefficient a 1 is the ideal value of 3.52 × 10 −2 m·°C −1/2 ·d −1/2 , the ice thickness is significantly overestimated with an RMSE as high as 31.16 cm [see Fig. 11(a)]. The empirical ice thickness growth coefficient a 11 is 2.11 × 10 −2 m·°C −1/2 ·d −1/2 obtained by substituting the measured maximum ice thickness and the corresponding CFDD into (13). at which point the simulation is good [see Fig. 11(b)] and the RMSE decreases to 3.15 cm. December 19, 2021 is the freezing date of the river channel. At this time, the river ice has a certain thickness, and it has experienced ice run before, so it is appropriate that the simulated river ice starts to grow earlier than the freezing date of the Section. A low value of the measured ice thickness is not captured by the model when the CFDD is between 400°C and 450°C [see purple box in Fig. 11(b)], which is known from the Toudaoguai hydrological station data to have been collected on January 26, 2022. We found that the daily average air temperature was above −5°C on four of the five days before it, 0.6, −2.2, −3.8, and −2.6°C, respectively, and the higher air temperatures led to some decay in ice thickness, which resumed freezing on January 25 with a daily average air temperature of −10°C and a renewed increase in ice thickness. The simulated ice thickness reached a maximum of 56.44 cm on February 25, 2022, at which point the ice thickness stopped growing.
On March 15, 2022, the hydrological station observed the break-up of the river section, and the ice thickness of the section was 24 cm on March 12. According to the measured ice thickness, this article assumes that the ice thickness decreases to 20 cm when the river break-up. From the maximum ice thickness, the ice thickness of the break-up date and the corresponding CTDD to find the ice thickness decay coefficient a 21 is 2.4 × 10 −3 m·°C −1 ·d −1 , the end of the CTDD calculation is determined by the break-up date. Fig. 11(c) shows the simulation of the ice thickness decay process, which is overestimated throughout the simulation with an RMSE of 4.47 cm. Comparing the simulated values of ice thickness throughout the ice period with the measured values [see Fig. 11(d)], R 2 is 0.857 and RMSE is 3.77 cm, so the air temperature model can better reflect the process of river ice thickness.
In this study, only the winter air temperature data from 2021 to 2022 were obtained for the Toudaoguai hydrological station, and in order to analyze the multi-year trend of ice thickness, the day-by-day air temperature data from 1956 to 2022 were collected from the Hohhot Baita meteorological station, where the data from 1964 to 1972 were missing. Similar to the previous treatment, river ice thickness was simulated for the winter of 2021 to 2022 using meteorological station air temperature data (see Fig. 12). Simulated values of ice thickness throughout the ice period were compared with measured values [see Fig. 12(d)] with R 2 of 0.855 and RMSE of 4.52 cm. The empirical ice thickness growth coefficient a 12 of 1.90 × 10 −2 m·°C −1/2 ·d −1/2 is less than a 11 when simulating the ice thickness growth process, which is due to the proximity of the meteorological station to the urban area and the overall warmer air temperature than the field river air temperature. The RMSE is 4.13 cm, and the ice thickness is overestimated at the beginning of the simulation and underestimated at the end of the simulation [see Fig. 12(b)]. When simulating the ice thickness decay process, the ice thickness decay coefficient a 22 is 2.3 × 10 −3 m·°C −1 ·d −1 and the RMSE is 4.99 cm. In total, the simulation of ice thickness by air temperature at the Baita meteorological station is slightly worse than that at the Toudaoguai hydrological station, but it still reflects the process of river ice thickness.
Based on the multiyear air temperature data from the Baita meteorological station, the start date of CFDD and CTDD calculation and the simulated maximum ice thickness were obtained     Fig. 13). The simulated results may not be representative of the actual conditions in each year, but they can reflect the coldness of the ice period in each year because the determination criteria are the same. We found that the CFDD start date is delayed at a rate of 0.5 d per decade, the CTDD start date is advanced at a rate of 1.0 d per decade, and the simulated maximum ice thickness is thinned at a rate of −0.7 cm per decade. The estimates of maximum ice thickness are highly variable with a standard deviation of 4.90 cm, and the absence of air temperature data for the 10 ice periods from 1964 to 1972 also adds to the uncertainty of the estimates. Although a thinning trend in maximum ice thickness is indicated, these uncertainties should be kept in mind when interpreting the results. Combined with the multiyear break-up dates provided by the Toudaoguai hydrological station, we calculated the decay rate of ice thickness in different years (see Table II). Most of the break-up dates of the ice flood season from 2012 to 2022 are in mid to late November, and the decay rate of ice thickness in the last 5 years is less than that of 2012 to 2017.

A. Comparison of Ice Thickness Estimates Between SAR and Air Temperature Models
Compared with the measured ice thickness data, the RMSE of the constructed SAR ice thickness inversion model is 11.75 cm, and the RMSE of the empirical model of air temperature simulated ice thickness in this article is 3.77 cm, so the air temperature model is more accurate in terms of its effect on ice thickness estimation. The ice thickness of the Yellow River was estimated for the winter of 2021 to 2022 in the study area using two models, respectively. The estimation results of SAR model reflect the trend of ice thickness growth in the early stage of river ice growth, but the change of ice internal structure in the later stage of river ice growth leads to the underestimation of ice thickness. The study of lake ice thickness also found that the SAR method may have large errors when the ice thickness exceeds a certain value [47]. The advantage of this method is that the spatial distribution of ice thickness can be characterized to a certain extent, which is useful for the identification of ice jams. The ice thickness results estimated using air temperature are in better agreement with the measured ice thickness, and the trend of ice thickness is more refined.
Most of the current SAR ice thickness estimation methods use direct construction of the regression relationship between polarization parameters and river ice thickness [21], [24], [28]. The estimation results of river ice thickness depend on the parameters of regression equations, and the robustness is generally. SAR images contain complex information, and the contribution of changes in river ice type, structure, thickness, and other attributes to radar echoes is unclear, and the mechanism of change remains controversial [36], [48]. In addition, the physical properties of snow on the surface of river ice are affected by factors such as temperature and precipitation, and changes in surface snow change the propagation process of electromagnetic waves, which in turn affects the accuracy of river ice thickness estimation [49].
The method of air temperature estimation for ice thickness is divided into two processes, ice thickness growth and ice thickness decay, and the model for simulating ice thickness growth has a certain physical basis [13], [45]. However, the neglect of heat flux at the ice-water interface, snow and hydraulic conditions increases the uncertainty of ice thickness estimation. The main influence on the river ice-water heat flux is the frictional heat energy caused by the tangential stress of the water flow, which plays a suppressive role in the ice cover thickening. The heat retention effect of snow on the ice cover also inhibits the growth of ice thickness. The model for simulating ice thickness decay is empirical [14], the date of river break-up, and the corresponding ice thickness have a great influence on the ice thickness simulation results. In addition, solar radiation which is closely related to river ice melting, is not taken into account. In a study of lake ice thickness, the empirical value of CTDD corresponding to complete ice melting was determined to be 90°C·d based on satellite images [50]. In this article, it is more realistic to use the break-up date of the river provided by the hydrological station to determine the maximum value of CTDD.
Both ice thickness estimation methods have their advantages and disadvantages, and combining the ice thickness estimation results with hydrological data can explain the river ice processes in the study area during the ice flood season from 2021 to 2022. Combined with Figs. 3, 5, 9, and 14, it can be seen that in the early part of November 2021, both water level and discharge decreased, which was caused by the ice flood control operation of upstream reservoirs. On November 22, 2021, the daily average air temperature changed from positive to negative, the water temperature changed from positive to 0°C, and the ice run started to appear, and the water level and discharge increased. That day is also the start date of air temperature simulated ice thickness growth. The ice run continued until December 19, 2021, when the river surface froze, causing the water level to rise, at which time the air temperature simulated ice thickness grew to 22.51 cm. The period from December 19, 2021 to March 15, 2022 was the freezing period. The water level and discharge changed slightly in the early freezing period, and then tended to be stable. The ice thickness simulated by air temperature reflected the growth and decay process of ice thickness at this stage. From the SAR ice thickness estimation results, it can be seen that ice jam is formed at the exit of the bend and an open water surface exists downstream of the bend. On March 15, 2022, the water temperature turned from 0°C to positive, the river ice break-up, the water level dropped sharply, and the discharge increased and then decreased. Overall, the combination of ice thickness estimation results from different methods and hydrological data can help us understand the river ice processes better than using a single method.

B. Findings and Shortcomings of This Study
River ice thickness studies based on satellite-based SAR data are mainly on Canadian rivers [21], [27], [28], [51], and a recent study has been attempted on high-order rivers on the Tibetan Plateau [24]. The SAR satellites involved in these studies include Radarsat-1, Radarsat-2, and Sentinel-1 (see Table III). Most of the previous findings found that the backscattering coefficient is positively correlated with the ice thickness, which is consistent with the results of this article. However, in the correlation between α, H, and ice thickness, this article and Mermoz yielded different results [21], and the reason for this difference may be the difference in the physical structure of river ice. This article explores the potential of dual-polarized Sentinel-1 SAR data to invert the thickness of the Yellow River ice. The RMSE of the SAR ice thickness inversion model constructed in this article is greater than 10 cm, and there is uncertainty in the interpretation of the ice thickness inversion results.
In this article, the acquisition time of field measured ice thickness is not completely matched with that of SAR image, which has some influence on the research results. Although the number of river ice sampling points is large, the coverage area is small, which limits the large-scale application of the ice thickness inversion model, and the generalizability of the model needs further validation. The radiative transmission of radar electromagnetic waves in river ice is not well understood, which makes the method theoretically inadequate. In addition, the type of river ice, surface snow cover, sediment inside the ice, and precipitation all affect the scattering mechanism of river ice, which requires us to collect more information in the future study of river ice. SAR images are often expensive, and ESA's free distribution of Sentinel-1 data made this study possible. The current SAR satellite mostly selects C-band data, and the NISAR satellite planned to be launched in 2023 will distribute data containing both L-band and S-band for free, so that the application of electromagnetic waves of different wavelengths in river ice thickness inversion and ice jam identification can be further explored in the future.
In this article, when the air temperature model is used to simulate the growth process of ice thickness, CFDD is calculated from the 1st day when the daily mean air temperature is less than 0°C for five consecutive days. When the ice thickness decay process is simulated and the daily mean air temperature is greater than −5°C for five consecutive days, CTDD is calculated from the 1st day. This method can simulate the change process of river ice thickness in the study area in winter 2021 to 2022 better than the SAR model. During the ice thickness growth stage, the RMSE of ice thickness estimated by air temperature of Toudaoguai hydrological station is 3.15 cm, and that of ice thickness estimated by air temperature of Baita meteorological Station is 4.13 cm, indicating that the closer the air temperature collection site is to the study area, the more it can reflect the actual change of ice thickness. During the ice thickness decay phase, the accuracy of ice thickness estimation depends mainly on the break-up date of river ice and the corresponding ice thickness. This article assumes an ice thickness of 20 cm at the time of river ice break-up, and the applicability of this value to rivers in other regions in different years is not yet known. The time of river ice break-up in this article is from the Toudaoguai hydrographic station, which may not be applicable in the study of a large-scale river section. With the increase of available satellites, the study of river ice phenology based on remote sensing technology is expected to provide more reliable time points in the future. This article discusses the application of two ice thickness estimation methods on the Yellow River. In the future, if air temperature and remote sensing methods can be organically combined to establish an ice thickness estimation model suitable for different river ice types and different river ice development periods, and add parameters such as snow cover, velocity, precipitation, etc., it will be helpful to promote the in-depth development of river ice research.

VI. CONCLUSION
In this study, the river ice thickness in the Shisifenzi section of the Yellow River was estimated using dual-polarized Sentienl-1SAR data and air temperature data, and the main conclusions are summarized as follows.
1) The VV polarized backscatter coefficient of Sentinel-1 SAR data has the highest correlation with the measured ice thickness, with a Pearson correlation coefficient of 0.702. 2) The SAR ice thickness inversion results help to identify the location of ice jam.
3) The simulation accuracy of air temperature model is better than that of SAR model. 4) The estimates indicate that the maximum river ice thickness in the study area is thinning at a rate of -0.7 cm per ten decades. Our study shows that the uncertainty of SAR ice thickness inversion results is large, and the scattering mechanism of electromagnetic waves in river ice is unclear, but the inversion results help the identification of ice jams. The simulation of air temperature on ice thickness can reflect the variation of river ice thickness well and the physical mechanism is clear. In the future, we will consider the combination of the two methods, which will facilitate the further development of river ice thickness estimation methods to better monitor the ice thickness of rivers in cold regions.