Integrating intensity and context for improved supervised river ice classification from dual-pol Sentinel-1 SAR data

River ice is a major contributor to flood risk in cold regions due to the physical impediment of flow caused by ice jamming. Although a variety of classifiers have been developed to distinguish ice types using HH or VV intensity of SAR data, mostly based on data from RADARSAT-1 and -2, these classifiers still experience problems with breakup classification, because meltwater development causes overlap in co-polarization backscatter intensities of open water and sheet ice pixels. In this study, we develop a Random Forest classifier based on multiple features of Sentinel-1 data for three main classes generally present during breakup: rubble ice, sheet ice and open water, in a case study over the Athabasca River in Canada. For each ice stage, intensity of the VV and VH backscatter, pseudo-polarimetric decomposition parameters and Grey Level Co-occurrence Matrix texture features were computed for 70 verified sample areas. Several classifiers were developed, based on i) solely intensity features or on ii) a combination of intensity, pseudo-polarimetric and texture features and each classifier was evaluated based on Recursive Feature Elimination with Cross-Validation and pair-wise correlation of the studied features. Results show improved classifier performance when including GLCM mean of VV intensity, and VH intensity features instead of the conventional classifier based solely on intensity. This highlights the complementary nature of texture and intensity for the classification of breaking river ice. GLCM mean incorporates spatial patterns of the co-polarized intensity and sensitivity to context, while VH intensity introduces cross-polarized surface and volume scattering signals and is less sensitive to wind than the commonly used co-polarized intensity. We conclude that the proposed method based on the combination of texture and intensity features is suitable for and performs well in physically complex situations such as breakup, which are hard to classify otherwise. This method has a high potential for classifying river ice operationally, also for data from other SAR missions. Since it is a generic approach, it also has potential to classify river ice along other rivers globally.


Introduction
River ice breakup is a major contributor to flooding in cold regions. The process may lead to the formation of ice jams that obstruct the flow of the river and thus cause flooding. Therefore, it is important to know when and where ice jams form so that emergency measures can be taken, e.g. the construction of (temporal) levees or the evacuation of populated areas. It is estimated that the annual direct costs, when keeping inflation in mind, are C$105 million in Canada (Gerard and Davar, 1995) and US$215 million in the United States (Carlson et al., 1989) alone, and that the indirect costs are even higher. Apart from the importance to monitor breakup occurrence, location and timing for local socio-economic effects, the timing is also an important factor in the context of global warming. For example, earlier breakup due to climate change has a positive feedback on Arctic sea ice melt (Park et al., 2020) and the intensifying hydrological cycle in polar regions (Déry et al., 2009, Peterson et al., 2002. Synthetic aperture radar (SAR) enables the monitoring of the ice cover during breakup with a high spatial and temporal resolution. The backscattered signal is a function of sensor and target characteristics. The melt of overlying snow marks the onset of the breakup season. This process causes the ice to become wet and causes a distinct shift in its backscatter behavior. In short, it changes from being a volume scatterer to a single bounce scatterer because its increased wetness obstructs the penetration of radar waves into the ice volume and so limits the ice backscattering to a process at the air-ice interface. Two principal types of breakup are distinguished, i.e. thermal breakups and dynamic breakups. Thermal breakups largely result from increasing air temperatures, which cause the ice to breakup apart in relatively large sheets, referred to as sheet ice, that remain stationary and thus melt in place. On the other hand, dynamic breakups are initiated by the sudden increase in the flow of water, which causes the ice to break apart in many small blocks, referred to as rubble ice. Depending on the situation, the flow moves the rubble ice downstream in the form of so-called ice run or the rubble ice may arrest at channel obstructions (e.g. sharps bends, islands, bridge piers) to form an ice jam. Rubble ice induces a much higher backscatter signal than sheet ice because its upper surface is much rougher. Van der Sanden et al. (2021) provide further details on the breakup process and its implication for the backscatter.
Our study aims to overcome three limitations of previously published approaches to classify the conditions of breaking river ice through application of radar satellite images. First, most existing approaches use the backscatter intensity comprised in single-polarized C-band SAR data-this holds true for the classification of river ice under freeze-up, winter, and breakup conditions (e.g. Jasek 2003;Puestow et al., 2004;Gauthier et al., 2010;Floyd et al., 2014;Van der Sanden et al., 2021). The most commonly applied polarization is Horizontal-Horizontal (HH). Relative to Vertical-Vertical (VV) polarized C-band (and X-band) signals, HH-signals are less sensitive to roughening effect of high winds on water (Long et al., 1996). This causes the backscatter of water to increase, which makes it more difficult to distinguish from ice cover. Also, the development of radar remote sensing for the characterization of river ice has largely been driven by Canadian scientists that used data from Canada's RADARSAT-1 and RADARSAT-2 satellites. RADARSAT-1 operated in the HH-polarization only. It appears, the HH-polarization remained the polarization of choice for later R&D with data available from the multi-polarization and full polarimetric RADARSAT-2 system. In this study, we assess and develop the utility of the European C-band Sentinel-1 satellite. Over our area of interest-the Athabasca River at Fort McMurray, Alberta, Canada (see section 2.1)-Sentinel-1 consistently operates in the VV-and VH-polarization. A classifier developed for either or both of these polarizations would improve the overall operational utility of radar remote sensing for river ice breakup monitoring because each additional applicable data source will augment the achievable observation frequency. As a rule, cross-polarized C-band measurements (e.g. VH and HV) have relatively poor signal-to-noise ratios and are therefore less suited for application to the discrimination of low backscatter targets such as wet sheet ice and calm water. Conversely, VH/HV signals maintain a better backscatter contrast between water and ice cover under high wind than both VV-and HHpolarized signals. Sentinel-1's standardized image acquisition schedule and cost-free data policy augment its potential benefit in support of operational river ice monitoring.
In addition to single-or dual-polarized backscatter intensities, other SAR image features may capture information that can facilitate the classification of river ice cover. For example, Gauthier et al. (2006) showed that image texture represents a good source of information for the classification of river ice types that form during the freeze-up process. Features that capture image texture provide information about the variation in the backscatter intensity over a certain image region--defined by the window size of the textural filter. As such, the inclusion of textural features in a classifier captures relationships between neighboring pixels, i.e. context. Similarly, Mermoz et al. (2009) and Lindenschmidt and Li (2018) demonstrated that polarimetric features can improve the understanding and classification of river ice processes during winter. Unlike backscatter intensity and image texture, features attained through the application of polarimetric decompositions capture information relating to scattering mechanisms, i.e. single (odd) bounce, volume, and/or double (even) bounce scattering. As such, polarimetric features hold promise for the detection of the transition from (dry) winter to (wet) spring conditions-a changeover that complicates the application of existing classification approaches. In theory, the backscatter of winter ice should result from volume and/or double bounce scattering within the ice mass, and/or single bounce scattering at the icewater surface. On the other hand, single bounce scattering at the air-ice surface governs the backscatter return signal of ice cover imaged under spring conditions; backscatter resulting from double bounce scattering is possible in the case of rubble ice. To this date, only Łoś and Pawłowski (2017) have published a paper addressing the use of Sentinel-1 image data for the classification of river ice cover. However, they dealt with ice under winter conditions.
Finally, most existing approaches to classify river ice involve socalled unsupervised algorithms. Especially, Fuzzy-K means clustering has been widely used for the identification of river ice types during freeze-up and breakup (e.g. Chu and Lindenschmidt, 2016;Gauthier et al., 2006;Sobiech and Dierking, 2013;Mermoz et al., 2009). Unsupervised algorithms operate by grouping image pixels in a pre-defined number of clusters regardless of their signal range. An analyst with access to ground reference data is required to associate each cluster with a thematic class of interest, e.g. an ice cover type. On the other hand, supervised algorithms allocate image pixels according to their signal range to a pre-defined number of thematic classes. The developer of the algorithm enables it to function by establishing the link between each class of interest and its signal range with the help of ground reference data, i.e. by training the algorithm. Without ground data and analyst intervention, unsupervised algorithms will not produce reliable classification results when one or more classes of interest are absent-unlike, well developed, supervised algorithms (Sobiech and Dierking, 2013;Van der Sanden et al., 2021).
Supervised classifiers can be developed in many different ways, but the Random Forest development approach has recently proven its success in studies concerned with the classification of lake ice and sea ice (e. g. Dabboor et al., 2018;Hoekstra et al., 2020;Shen et al., 2017). This approach offers several advantages. Firstly, multiple features are easily accommodated and so the complementary information comprised in backscatter intensity, image texture, and polarimetric features is fully exploited. Secondly, Random Forest provides insight into feature importance, which can be used to better understand how the SAR signals interact with river ice.
This study aims to provide insight into the utility of dual-pol Sentinel-1 SAR data for operational monitoring of river ice during breakup using Random Forest classification. Its main objective was to explore the potentials of intensity, pseudo-polarimetric and texture features for river ice classification. First, the Athabasca River test site, and the used data sets are presented (Section 2). Secondly, the method to select optimal features is discussed and the assessment of the accuracy used to compare the developed classifiers is explained (Section 3). Subsequently, the performance of different combinations of features is evaluated (Section 4). Then, implications of using this supervised classifier for other SAR sensors and other rivers are discussed, and recommendations to improve the current approach are provided (Section 5). Finally, the importance of supervised classification for river ice breakup is illustrated (Section 6).

Study area
The Athabasca River is one of Canada's longest rivers, which originates in a large icefield of the Rocky Mountains located along the borders of British Colombia and Alberta, and then travels 1538 km through the province of Alberta to eventually discharge in Lake Athabasca. Figure 1 shows the section of the Athabasca River near Fort McMurray that was selected for this study. The station numbers displayed on the river indicate distance in kilometers measured upstream from the mouth of the Athabasca River. The selected reach is about 160 km long and extends from the House River (km 445) to approximately 10 km downstream of the Clearwater River (km 285). The studied river stretch is particularly prone to ice jam flooding during the breakup season, with regular flooding events dating back to the 1870's when the area was first settled.
Breakup in this reach of the Athabasca River has a consistent pattern, even though the timing and the magnitude of the events vary from year to year. In spring, when temperatures rise and precipitation and snowmelt runoff increase, the deterioration of the ice cover starts upstream (southwest) and progresses downstream (northeast). Rubble ice accumulations develop mainly at the locations of the numerous rapids. The release of such accumulations causes waves, which may create more ice blocks that extend the existing small ice jams. Eventually, the weight of the ice jams located at the rapids is larger than the resisting force resulting in a sequence of ice runs and ice jams. Downstream of Fort McMurray, several small islands, a meander in the river and an abrupt decrease in river slope greatly increase the risk of ice jamming and flooding. The inflow from the Clearwater River provides an additional source of ice and water. The ice jams obstruct the water flow of the Clearwater and Athabasca River, leading to increasing water levels upstream and pose flood hazard to the city of Fort McMurray (She et al, 2009;Lindenschmidt and Li, 2019). Significant floodings of Fort McMurray have occurred as a result of this breakup pattern (for example in 1977, 1997 and 2020).

Synthetic Aperture Radar data
To develop a supervised classifier, SAR images during the 2018-2019 and 2019-2020 breakup of the Athabasca River were acquired (Table 1). SAR images were selected from start to end date of the breakup season, which was identified based on meltwater on the ice visible on reference data (Section 2.3). Sentinel-1 data are provided free of charge through the Copernicus programme of the European Space Agency.
The Sentinel-1 mission is a constellation of two polar-orbiting satellites launched in April 2014 -Sentinel-1A -and in April 2015 -Sentinel-1B. The Interferometric Wide swath (IW) acquisition mode was used, which captures the entire study area with its 250 km swath. The IW data products have a 5 m range by 20 m azimuth spatial resolution and dual-polarization capabilities (VV and VH) over the studied region. The used images have incidence angles ranging from 29 • to 46 • . We used the Level-1 Single-Look-Complex (SLC) data, as phase information is preserved (necessary for pseudo-polarimetric computations) and control over the entire preprocessing scheme is possible.

Reference data
Because of the ongoing risk of ice jam flooding in Fort McMurray, Alberta Environment and Parks (AEP) has an annual river ice monitoring and observation programme for the Athabasca River upstream of Fort McMurray (Sun et al., 2015). The on-site measurements by permanent cameras and observation flights operated and conducted by AEP were used as reference data for training and validation of sample areas. To ascertain or complement when AEP observations were unavailable, also Sentinel-2 data were used as reference data to differentiate between presence of either sheet or rubble ice cover and open water conditions.
The time difference between SAR images and reference data ranged from 3 hours to 1.5 days. The AEP ice progression maps and observation reports were mainly used as reference data when the time differences between the SAR acquisitions and flights were small, a few hours. However, through personal communication with the river ice experts of AEP, it could be confirmed that for some dates certain ice stages were stationary for a longer period of time. This enabled the use of AEP Figure 1. The study reach of the Athabasca River (km 285 -445). The river flows from south-west to north-east. Due to the complex channel conditions in this reach of the Athabasca River, it is highly susceptible to ice runs and ice jams during the breakup season. monitoring data, also when time differences were larger, for example one day. Figure 2 illustrates the overall image processing methodology of river ice during breakup. The Sentinel-1 SLC SAR images were prepared for classification using three different preprocessing schemes with ESA's Sentinel Application Platform (SNAP) software (version 7.0). This study assessed different realizations of Random Forest classifications, that combined intensity, pseudo-polarimetric, and texture features. The code, model and detailed workflow to perform river ice classifications (1) using the developed classifier and (2) by building a new Random Forest classification (for other SAR sensors or other rivers) are available on the GitHub repository at link: https://github.com /SdeRodaHusman/remotesensing-of-river-ice.

Image preprocessing and feature extraction
During preprocessing, Sentinel-1 Level-1 SLC products were converted to images representing a series of descriptive features that were subsequently applied to train a Random Forest classifier. The features computed relate to three different sources of information: backscatter, backscatter phase, and image texture. The backscatter intensity in the VV and VH polarization were computed according to the first preprocessing procedure. The second diagram provided H-A-α pseudo-polarimetric decomposition parameters (Cloude and Pottier, 1997). The texture scheme provided information about the spatial distribution of the backscatter, expressed as Grey Level Co-occurrence Matrix (GLCM) features (Haralick and Bryant, 1976).
The settings for Calibration, Speckle filtering, H-A-α decomposition and Image texture are discussed below. The default settings in SNAP were used for the other preprocessing steps.
In the intensity and pseudo-polarimetric schemes the digital pixel values were converted to the gamma nought (γ 0 ) calibrated product. Gamma nought is less sensitive to variations of the incidence angle, compared to radar backscattering coefficients like sigma nought and beta nought (Raney, 1998). In the polarimetric scheme a complex output was selected, because polarimetric phase information should be preserved for the computation of the decomposition parameters.
The equivalent number of looks (ENL) of Sentinel-1 SLC images corresponds to 1. After preprocessing an ENL around 10 is desired, because this results in a balanced trade-off between noise and resolution ( Van der Sanden et al., 2021). To reach an ENL of 10, a two-step approach was used, consisting of multilooking and speckle filtering. Multi-look of four at the range and one at azimuth was applied, resulting in near square-sized pixels with 15 meters pixel spacing. Speckle filtering was performed with the gamma-maximum-a-posteriori (MAP) filter (Lopes et al., 1993) with window size 3 x 3.
Polarimetric decomposition provides information about the scattering mechanisms from a target. Since the images used in this study are not quad-polarized, a modified version of the H-A-α decomposition was applied (Cloude and Pottier, 1997;Cloude, 2007). Based on visual comparisons, a window size of 7x7 was found to be optimal in discriminating the studied ice stages.
The texture features were calculated based on the Grey Level Cooccurrence Matrix (GLCM), which quantifies how frequent different combinations of pixel grey levels occur in the filter window (Haralick and Bryant, 1976). To allow for the comparison of texture over the entire breakup season, the data were normalized, from 0 to 255, by linear scaling between the lowest and highest backscatter values detected over the entire breakup season. The statistic group features (GLCM mean, GLCM variance and GLCM correlation) were computed on VV intensities. After some tests, a window size of 11 x 11 and displacement of 2 were selected.

Image classification
Random Forest classification was performed using the sklearn package in Python. Random Forest is an ensemble classifier consisting of a collection of tree-type classifiers. In the training process, the Random Forest creates multiple classification and regression trees, each of which is trained on a different bootstrap sample by randomly resampling the original training sample with replacement. During classification, each tree votes for the predicted class, and a pixel is labeled with the class having the most votes. In this research, the Random Forest algorithm was applied for different combinations of input features.
Based on the reference data described in Section 2.3, sample areas were created for rubble ice, sheet ice, and open water. The goal of sample area selection was twofold: i) to analyse the feature behavior per ice stage to gain a profound understanding of the river ice processed during breakup, and ii) to assess the accuracy of the classification results.
Sample areas of 100 pixels (or 1000 ENL) were selected manually over homogenous areas. A total of 210 sample regions were created, a third of them corresponding to rubble ice, a third to sheet ice and a third to open water. Per ice stage, 70% of the sample areas were used for training and 30% for validation of the developed classifiers.
All sample area pixels included values for intensity features VV and VH, polarimetric pseudo-decomposition parameters Alpha, Anisotropy and Entropy and texture features GLCM mean, GLCM variance and GLCM correlation. However, a Random Forest model becomes increasingly complex with an increasing number of features. In practice, not every feature carries information that is useful for discriminating ice stages. Some features are either redundant or irrelevant and hence can be removed. By only selecting important features, a model becomes easier to interpret, overfitting is suppressed, and the computational costs and time required to train the model are reduced.
Classifiers were developed for seven combinations of feature classes (Int, Pol, Tex, IntPol, IntTex, PolTex, IntPolTex) were Int stands for intensity features, Pol for pseudo-polarimetric features and Tex for texture features. For each classifier the optimal features were selected using two steps. First, the Recursive Feature Elimination (RFE) with Cross-Validation technique was applied to identify the importance of each feature. The first step of RFE was to divide the training data into 10 randomly divided subsets, also known as cross-validation folds. By using this approach, the validation data was a truly unseen data set for testing the final model (Cawley and Talbot, 2010). Next, for each subset the performance of a Random Forest model was evaluated, and the importance of each feature was computed. The least important features were removed. Then the Random Forest model was re-built and importance scores were computed again (Kuhn and Johnson, 2019). Second, the correlation between each feature was determined by computing a correlation matrix. To find the optimal combination of features for all the classifiers listed before, first the feature with the highest RFE importance score was selected, other features of classes of interest were only included when its correlation was higher than -0.85 or lower than 0.85.
After the extraction and selection of the most important, uncorrelated features, hyperparameters were tuned to optimize the performance of the Random Forest model. Four hyperparameters were adjusted one by one, that is (1) the number of trees in the forest of the model (n_estimators), (2) the maximum depth of each tree (max_depth), (3) minimum number of samples required to split an internal leaf node (min_samples_split) and (4) minimum number of samples required to be at a leaf node (min_sample_leaf). Validation curves were plotted to show the influence of a single hyperparameter on the training score and the validation score. A total number of 46 estimators, a maximum depth of 13, a minimum sample split of 2, and minimum sample at a leaf of 1 were found to be optimal. Then, seven Random Forest classifications were performed using the optimal hyperparameters and a combination of features.

Accuracy assessment
Accuracy assessment was conducted for the validation sample areas. We examined three commonly used evaluation indices, confusion matrix, overall accuracy, and Cohen's Kappa coefficient. The confusion matrices were obtained by comparing the classification results with the actual labels based on the reference data. The diagonal elements of the confusion matrix provide overall accuracy, where the Kappa coefficient takes all the elements in the confusion matrix into account. The Kappa coefficient (Cohen 1960), and its estimated standard deviation (Fleiss et al. 1969), were calculated for each confusion matrix to evaluate the agreement between the classification results and the reference data.
In this study, intensity, pseudo-polarimetric and texture features were combined to find the optimal Random Forest classifier for river ice breakup. Hypothesis testing was used to examine whether one classifier performed significantly better than another one. The Kappa coefficient and its standard deviation were used to construct a hypothesis test to identify significant differences between classifiers. The test statistic ΔK is given by where K is the Kappa coefficient, σ 2 its corresponding variance and A and B represent the two classifiers that are compared (Bishop et al., 2007;Congalton et al., 1983). In this study a confidence level of 95% was used, meaning two classifiers may be considered significantly different when ΔK > 1.96 (Benson and DeGloria, 1985).

Results
Figure 3 presents the normalized feature importance of the Random Forest classification and the overall accuracy depending on the number of features selected. The Random Forest classification with all features was selected to investigate the relative importance of each feature for river ice classification. Figure 3 clearly demonstrates that features from all feature classes (intensity, pseudo-polarimetric, texture) are important for river ice classification. GLCM mean, pseudo-anisotropy, and VH intensity were found to have the highest feature importance. Pairwise correlation between features had been considered, as RFE with crossvalidation identifies the optimal features by eliminating the redundant features during the cross-validation phase. This explains the relatively low feature importance of the VV intensity, which was 92.0% correlated to GLCM mean. Figure 4 shows the boxplots of four of the studied features, divided into the three ice stages. As expected, for rubble ice the largest VV and VH intensities are found. Due to the roughness of the surface, both the co-and cross-polarized channels show a high backscatter. Open water has a low VV and VH backscatter. The smooth water surface behaves like a mirror, and since no small incidence angles were used in this study, this results in specular reflection directed away from the sensor. Sheet ice has the largest range of VV and VH backscatter values. The sample areas at the beginning of the breakup season show high values. Volume scattering is enhanced, because the dielectric constant is low, leading to penetration through the ice cover. When the melt starts, the volume scattering is impeded and specular reflection dominates, resulting in low VV and VH intensities.
The highest GLCM mean values are found for rubble ice, due to the homogeneity in high backscatter over these rough surfaces. Open water shows spatially consistent low backscatter, resulting in the low GLCM mean values. Again, the texture of sheet ice differs throughout the breakup season. Sheet ice in the early phase of the breakup season has a uniform, high backscatter, which results in a high GLCM mean value. Later in the season, when it starts to melt, lower GLCM mean values are found. However, the backscatter pattern is less uniform than for open water, with larger variability in backscatter over the melting area. Hence slightly higher GLCM mean values are found for sheet ice under melting conditions than for open water.
Anisotropy provides information on the relative importance of secondary scattering mechanisms, which can be single bounce, volume, or double bounce scattering. However, since only dual-polarized signals were studied, only pseudo-polarimetric decomposition parameters were computed. Therefore, the pseudo parameters have no explicit physical interpretation. However, the pseudo-anisotropy feature was found to help distinguishing open water from the other ice stages.
Given the results obtained from the feature importance analysis and the pairwise correlation between features, the most important and noncorrelated features were extracted. Table 2 shows the different classifiers that were built, where Int represents the intensity features, Pol the polarimetric features, and Tex the texture features. The overall accuracy of the Int-classifier was 83.3%, and the overall accuracy of the classifier that performed best, the IntTex-classifier, was 86.8%. Table 3 shows the confusion matrices of the Int-and IntTex-classifiers. Even though the IntTex-classifier using Sentinel-1 data lead to a significantly higher classification accuracy (Kappa = 0.805) than the Intclassifier (Kappa = 0.750), confusion between sheet ice and open water (omission error of 10%) and sheet ice and rubble ice (omission error of 16%) remain a source of uncertainty for the ice classification. Most  confusion between sheet ice and rubble ice took place during the morning observations, when refreezing of meltwater is probable due to colder temperatures overnight. On the other hand, under melting conditions sheet ice and open water can be confused, due to low backscattering values from both ice stages. Figure 5 (a and b) illustrates the VH intensity and GLCM mean feature for three SAR images, just before, during and after an ice jam event during breakup season 2019-2020. Figure 5 (c) shows the classification maps of the IntTex-classifier during that period, combined with an ice progression map in Figure 5 (d) and a visual assessment in Figure 5 (e). Most classification results are comparable to reference data. However, some of the pixels in the upstream part of the river of the Sentinel-1 image of 5 May 2020 may have been misclassified as sheet ice. This might be caused by wind speeds that result in a rough surface in this east-west part of the river. Classified Sentinel-1 images later in May result in open water pixels in this stretch of the river, which indicates that also for this smaller east-west oriented part of the river the classifier is able to generate accurate classification maps. Figure 5 clearly indicates the added value of GLCM mean, when distinguishing water and ice. The sheet ice patches are captured by the texture features. The preprocessed Sentinel-1 images have a pixel spacing of 15 m x 15 m, resulting in a window of 165 meters for the computation of the texture features.

Discussion
Our results reveal that the developed classifier provides a promising means to distinguish river ice during breakup. This analysis overcomes many of the current limitations that exist in river ice classification by creating a supervised classifier based on multiple features from Sentinel-1 data.
The results showed that including features other than solely intensity, significantly increases classification results. We analyzed the performance of three classes of features, i.e. intensity, polarimetric and texture features, for river ice classification. Intensity and texture features result in a highly positive impact on classification accuracies, whereas the impact of polarimetric features is limited. Especially the GLCM mean texture feature is beneficial during breakup, because it detects the patches that are formed during breakup. Polarimetric decomposition parameters are more promising for ice classification during winter (Mermoz et al., 2009;Lindenschmidt and Li, 2018) than during breakup. This can be explained by a transition in scattering mechanisms. When spring arrives, mostly odd bounce scattering -at the air-ice surface -dominates the backscatter signal, instead of a combination of volume scattering -within the ice mass -and odd bounce scattering -at the ice-water surface -during winter. We hoped that polarimetry would allow us to detect the transition from winter to spring, but we did not find evidence for this. More research is needed to study this for fully polarized data, which might contain more insight than the dual-polarized data used in this study.
The developed classifier with best performance used a combination of the GLCM mean texture feature and VH intensity. Rubble ice has high GLCM mean and VH intensity values, because of the homogenously high VV backscatter and a large amount of depolarization. The GLCM mean of sheet ice is consistently lower than for rubble ice. Depending on the advancement of the melting process, the VH intensity of sheet ice can be high or low. Open water has a low GLCM mean and low VH intensity signal, because much of the signal is reflected of its smooth surface in the form of specular reflection. A great advantage of GLCM mean compared to the commonly used VV intensity, is the sensitivity to spatial context. Texture features consider the spatial relationship between VV intensity pixels, which is beneficial for classification since river ice has specific spatial patterns. However, further research should focus on the optimal settings of GLCM features, which depend on the spatial resolution and river dimensions. This paper introduces the first river ice breakup classifier based on Sentinel-1 data. The developed classifier is highly beneficial for monitoring purposes, because of the good classification results and Sentinel-1's great accessibility. Although developed for dual-polarized Sentinel-1 data, this method could also be applied to other rivers and SAR sensors. On the GitHub repository, a roadmap following two paths is provided. Following the first path, our developed classifier can be implemented for classification of river ice breakup using Sentinel-1 data. However, for SAR sensors that are differently polarized or rivers with other breakup patterns (for example due to a different climate, orientation, surrounding) it is recommended to create a separate training data set. Consequently, the second path leads to the development of a new Random Forest classifier. For this approach, new sample areas should be collected for each ice stage based on reference data. This was tested for breakup observations of the Athabasca River from RADARSAT-2 (2018RADARSAT-2 ( -2019 and RADARSAT Constellation Mission (2019. Random Forest classification of intensity, polarimetric and texture features resulted in a 91.2% (Kappa = 0.867) classification accuracy from RADARSAT-2 data and a 91.0% (Kappa = 0.865) accuracy from RCM data, showing the potential of the implementation of the proposed method on a larger scale.
Even though the developed supervised classifier results in high classification accuracy, different ice stages (e.g., open water and sheet ice covered by wet snow) with similar feature characteristics remain a challenge, albeit less of a challenge after inclusion of texture features.

Table 2
Hypothesis test comparing the performance of Random Forest classifications using different combinations of features to the Int-classifier. Note that for the classification results that are significantly different (accuracies at 95% confidence level, shown in bold), ↓ and ↑ indicate that the classifier performance is worse or better than the Int-classifier, respectively.  Temporal patterns could further improve classification results, as was proposed by Van der Sanden et al. (2021) as well. Figure 6 shows the typical backscatter behavior over one breakup season, subdivided into i) thermal breakup and ii) dynamic breakup patterns. Depending on the nature of the breakup, two or three peaks in backscatter can be expected, for thermal and dynamic breakups, respectively. For both breakup processes, the pre-breakup phase starts with high backscatter signals in winter conditions (Figure 6-1), which reduces due to melt of snow and ice ( Figure 6-2 and 6-3). Subsequently, the backscatter increases again, due to enhanced roughness of the ice cover's surface because of meltwater drainage (Figure 6-4). Then the thermal and dynamic breakup paths separate. A thermal meltout decreases the backscatter in case of a thermal breakup (Figure 6-6a and 6-7), whereas ice jam formation leads to an increased backscatter in a dynamic breakup ( Figure 6-5), followed by an ice run (Figure 6-6b) and/or melt of the grounded rubble ice ( Figure 6-6a and Figure 6-7). The presented typical temporal backscatter patterns could potentially be exploited to improve classification accuracies. As an example, Markov chains have already proven their value in the application of SAR images for flood mapping (e.g. Martinis and Twele, 2010) and land cover mapping (e.g. Hagensieker et al., 2017). This temporal pattern recognition technique exploits observations relating to past events to predict future events and their observational characteristics. A study to further assess and develop the utility of Markov chains for the classification of breaking river ice by means of SAR images is recommended.

Conclusion
In this article, a classification approach was proposed which addresses (1) the potential of a multi-feature approach for river ice classification, (2) the lack of a Sentinel-1 classifier, especially during breakup, and (3) the problem of poor accuracies when not all water and ice stages are present, which applies to unsupervised classification. This study therefore focused on the use of different features for a supervised classification using Sentinel-1 data. The results demonstrate that multiple features enhanced river ice cover information and, accordingly, the classification accuracy. Based on the integration of feature importance analysis and pair-wise correlation, a high classification accuracy of 86.83% was attained using a combination of GLCM mean of VV, and VH intensity features. This represents a significant overall improvement of 3.35% for classification using only intensity features. The results show that the inclusion of texture features augment sheet ice versus water discrimination, providing additional information unavailable from backscatter intensity.
In this article, Sentinel-1 data was used to classify the breakup season of 2018-2019 and 2019-2020 of the Athabasca River. The Random Forest approach was also tested on two other SAR missions, which resulted in overall accuracies over 85% for both missions. The developed method is flexible, and allows for adding different features or using other classifiers than Random Forest. Therefore, it is expected that this supervised approach will get further attention in future studies.

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Figure 6. Illustration of VV backscatter over one breakup season, divided into thermal breakup (red) and dynamic breakup (blue). The main backscatter signatures both start with high backscatters due to volume scattering in winter conditions (1). When temperatures start to rise the backscatter typically decreases (2, 3), followed by a jump in the backscatter induced by the rougher surface resulting from the penetrated snow melt (4). In case of a dynamic breakup, ice jams lead to a high backscatter signal (5), while during a thermal breakup the ice slowly decays leading to a smaller returned signal (6). When the river is free of ice, a low backscatter is observed (7).