Sea ice and water classification on dual-polarized Sentinel-1 imagery during melting season

We provide a new sea ice and water classification product with high spatial and high temporal coverage using 10 Sentinel-1 Synthetic Aperture Radar (SAR) data. The classification is applied in the Fram Strait region in the Arctic during melting seasons, when the contrast between backscatter intensities of different ice types observed by SAR is reduced due to the melted ice surface and wet snow on sea ice. The wet or melted snow strongly reduces the SAR penetration depth and thus suppresses the volume scattering contribution of sea ice. Furthermore, within the marginal sea ice zone (MIZ) ambiguities between ice and water can result from the effects of winds and ocean currents on the ocean SAR backscatter. 15 On the other hand, under calm conditions the contrast between thin ice and flat open water can be reduced, and thus decrease the separability of some ice. In summary, the melting season represents the most challenging time of the year for reliable ice-water classification from SAR data. We propose here a new approach to overcome these problems by using a mixture statistical distribution based conditional random fields (MSTA-CRF) model. To obtain reliable ice-water classification whilst maintaining a fast computation time suitable for operational applications, the MSTA-CRF adopts a 20 superpixel approach in the fully connected CRF model. The MSTA-CRF is a semantic model, which integrates statistical distributions (Gamma, Weibull, Alpha-Stable, etc.) to model the backscatters of ice and water and overcome the effects of speckle noise and wind-roughened water. Dual-polarization Extended Wide (EW) mode Sentinel-1A/1B SAR data with 40 m spatial resolution is available several times per day within the Fram Strait region. Observations from June to September during the six years 2015-2020 are collected and classified into ice and water categories. The classification 25 performance of algorithm is evaluated using ice charts from the Ice Service at the Norwegian Meteorological Institute (MET Norway). The methods of training sample selection, and their application to processing large data volumes and automatic classification of ice-water are discussed. In the experiment part, we demonstrate that the MSTA-CRF can provide a good performance with about 90% accuracy for ice-water classification, which is better than most of other state-of-the art algorithms. Compared with the 89 GHz microwave radiometer ASI sea ice concentration product, the sea ice extent in 30 Fram Strait derived from MSTA-CRF algorithm is lower during melting seasons from 2015 to 2020, and the monthly June to September sea ice area does not change so much in 2015-2017 and 2019-2020, but it has a significant decrease in 2018.


Introduction
Synthetic Aperture Radar (SAR) is widely used in sea ice monitoring of sea ice concentration, area, leads, ice floe, ice edge, ice classification and deformation (Dierking, 2010). With its all-weather day-and-night, high spatial resolution and subsurface imaging capabilities, SAR has become the most important observation technique for operational sea ice monitoring related to maritime navigation and search-and rescue (Liu et al., 2015;Zakhvatkina et al., 2017). From the 5 early 1990s, spaceborne C-band SAR such as ERS-1/-2, ENVISAT, RADARSAT-1/-2(RS1/RS2) and Sentinel-1A/-1B have been used as the primary data sources for the operational ice services since it provides good contrast between open water (OW) and ice categories (Johannessen et al., 2007). Sentinel-1 is the major data source for operational sea ice and iceberg monitoring in the polar regions from the European Commission's Copernicus programmer, and is freely available to registered users. 10 Ice charts are important for operational applications like ship navigation, and scientific studies such as the mass balance between atmosphere-ocean-ice systems. Earlier generations of SAR sensor provided only single-polarization HH or VV (horizontal or vertical transmit and receive polarization) where it can be difficult to distinguish sea ice from windroughened open water (Scheuchl et al., 2005). Newer SAR sensors introduced a cross polarization HV (horizontal transmit and vertical receive polarization), or the opposite VH, channel where the backscattering coefficients are relatively 15 independent of the water surface roughness conditions caused by high wind speed (Dierking, 2013). This can be utilized to improve water detection as it provides good contrast between wind-roughened water and ice during freeze-up periods.
However, the HV channel usually has a lower signal to noise ratio. In some cases, polarimetric features for improved distinguishing of different ice types are used when the full polarimetric data is available (Ressel et al., 2016;Johansson et al. 2017). The co/cross polarization (i.e. HH/HV or VV/VH; Komarov et al., 2021) is useful in sea ice classification since 20 it not only provides a better contrast between ice and water, but also decreases the instability of backscatter value in open water areas caused by high speed wind. (Liu et al., 2015;Zhu et al., 2016;Zakhvatkina et al., 2017). The Map-Guide Sea Ice Classification System (MAGIC) integrates Iterative Region Growing using Semantic (IRGC) classification (Yu and Clausi, 2008) with pixel-based Supported Vector Machine (SVM) classifier. A classification accuracy of 96.5% with respect to manually drawn ice charts was achieved using a limited number (20) of RS-2 scenes in Beaufort Sea area (Clausi 25 et al., 2010;Ochilov and Clausi, 2012;Leigh et al., 2014). Moreover, the operational Ice Service at the Norwegian Meteorological Institute (MET Norway) use high-resolution SAR data from Sentinel-1, RS-2 and COSMO SkyMed, in combination with optical imaging from Sentonel-2 and -3, NASA Suomi NPP VIIRS, NOAA AVHRR for manual interpretation, with low resolution passive microwave form low resolution microwave from AMSR2 used as a last resort if no other data is available. Backscattering characteristics determined by sea ice surface roughness and its dielectric properties can be fully explored to separate different sea ice types and water. As the backscattering is usually affected by ocean waves propagating into the ice area, it is not enough to only rely on backscattering intensity in identifying the different ice types and water.
Additional information such as statistical distribution of sea clutter performs well in distinguishing water and ice. For https://doi.org /10.5194/tc-2021-85 Preprint. Discussion started: 30 March 2021 c Author(s) 2021. CC BY 4.0 License. modeling the backscatter signal in SAR imagery, the K-distribution is a good model for four-look HH-polarized SAR data from Arctic sea ice (Roberts and Furui, 2000;Joughin et al., 1993). Since scatters of the sea clutter cannot be completely uniformly distributed (Ward et al., 2006), the mixture of several statistical distributions is needed. Statistical distributions combined with advanced classifiers have been widely used to improve ice and water classification (Barber andLedrew, 1991, Nystuen et al., 1992). Besides the statistical distribution-based methods, Artificial Neural Network (ANN) 5 characterized by modeling nonlinear relationships between the input coherency matrix features and ice categories (Ressel et al., 2016) shows promising potential in ice water analysis. SVM realize ice water classification by training the kernel with the transformation into high dimensional space, which solves the problem of undistinguishable features in low dimensional space. Textual feature based neural network methods also shows good performance in ice-water classification, especially in the marginal ice zone (Aldenhoff et al., 2018). Murashkin et al. (2018) use textural features and a random 10 forest classifier to detect leads within the ice pack from Sentinel-1 data but do not focus on th MIZ. Markov random field (MRF) has been widely used in sea ice segmentation and classification since MRF framework modeling the contextual information by utilizing the Gaussian statistics to exploit backscatter characteristics also shows its efficiency in sea ice classification (Clausi, 2001). However, MRF usually assumes that the observation vectors are conditionally independent, which ignores the interactions between ice and water, thus it is difficult to describe the dynamics of sea ice scenes especially 15 in melting seasons. Conditional random fields (CRF) can capture the contextual information in both labels and observed data, and directly model the posterior as a Gibbs distribution and thus allow capturing the dependence of the observed data (Wang, 2016). The CRF model solve the maximum posterior inference over the defined superpixels obtained by the segmentation, in order to optimize the heterogeneity of the superpixels, and the pixel-wise ice maps obtained by SVM are utilized as the category prior to updating the label in superpixel (Zhu et al., 2016). 20 In this paper, the Mixture Statistical distribution based fully connected CRF (MSTA-CRF) framework is proposed for accurate and efficient ice-water classification. Our objective is to use this statistical information in a Bayesian classification scheme for an operational ice-water classification algorithm. The higher order potentials method not only considers the neighboring information but also the long-range interactions between pixels. Statistical distribution such as Rayleigh, Weibull, Gamma, Log-Normal and Alpha-Stable are integrated in the CRF framework, parameter estimation of mixture 25 statistical distribution is investigated in this work. Considering the heterogeneity of the superpixel, a sub-superpixel based FCNN (Fully connected neural network) is used to represent the hidden information corresponding to the spatial and semantic relationship for each superpixel. We propose a mixture of statistical distributions based on fully connected CRF for ice-water classification using Sentinel-1 dual polarized SAR data with the following goals: 1. Roughened ice-free water caused by winds and ocean currents make it difficult to distinguish water from ice in 30 the co-polarization (HH or VV) channel of SAR data. Statistical distribution is proposed for modeling the SAR backscatter signal to deal with the wind-roughened ice-free pixels. Then, we incorporate the statistical distributions into the proposed CRF framework which considers the spatial and semantic information.
2. During the melt seasons, sea ice deformation features shows much more textual that decrease the severability between flat thin ice and calm water. Most sea ice classification methods rely on multiple features to input, e.g., 35 https://doi.org /10.5194/tc-2021-85 Preprint.  into SVM. To deal with the large computational cost, these approaches usually use a down-sampling strategy, which results in losing a lot of thin structures such as leads. We propose sub-superpixel, based MSTA-CRF models, to reduce processing time as well as improve the accuracy of ice-water classification.
The statistical distribution is used to reduce the influence of the wind roughened water. The use of superpixels can improve ice-water classification results by preserving leads and ice edges. An operational ice-water classification algorithm 5 by MSTA-CRF model is proposed for addressing problems of ice-water classification during melt seasons in this paper. This paper is organized as follows: In section 2, we give a brief description of the research area and data used. In section 3, we introduce the framework of the proposed algorithm, and give details of the methodology including preprocessing, training data selection, MSTA-CRF modeling and classification. We evaluate the ice-water classification results with the help of MET Norway ice charts to manually delineate the testing and training samples of open water and 10 ice for MSTA-CRF in Section 4. In section 5, we discuss the dynamics of sea ice during melting seasons. After that, we give our conclusion in Section 6.

Research area
To consider the spatial contextual information and preserve the spatial details of each pixel in SAR imagery, the 15 energy function based maximum a posteriori (MAP) estimation in MSTA-CRF framework is proposed for operational ice water classification during melting seasons in Fram Strait. Fram Strait is between the northeast of Greenland and northwest of Svalbard, and is the main outflow passage for sea ice export of the Arctic Ocean thereby controlling the mass balance of the Arctic Ocean sea ice cover (Maslanik et al., 2009;Spreen et al., 2009, Ludwig, et al., 2020. Increasing ice exports of Fram Strait in summer will contribute directly to open water further north (Metfies et al., 2017;Spreen et al., 2020). 20 Figure 1 shows an overview of the research area and some satellite scenes used in this manuscript.  Table 1.

b) MET sea ice chart
The operational Ice Service at Norwegian Meteorological Institute (MET Norway) uses high-resolution SAR data 10 and optical data from Sentinel-2 and 3, NASA Suomi VIIRS, NOAA AVHRR for manual interpretation, with low resolution passive microwave from AMSR-2 used as last resort if no other data is available for manual interpretation by its ice analysts.
The sea ice concentration classes, including whether the ice is drifting or (land) fast, together with reference colors are shown in Table 2. For the ice-water classification in this study, the Ice Free, Open Water and Very Open Drift Ice are defined as water categories with the sea ice concentration (SIC) less than 1 (tenths level), and Open Drift Ice, Close Drift 15 Ice, Very Close Drift Ice and Fast Ice are defined as the ice categories. The accuracy of different algorithms is calculated pixel-by-pixel with the MET Norway ice charts. The validation of sea ice classification using MET Norway ice chart will be discussed in Section 4.

c) ASI sea ice concentration 20
The daily ASI sea ice concentration (Spreen et al., 2008) product from AMSR2 is also used to validate the performance of MSTA-CRF algorithm. The ASI product used in the manuscript is retrieved using the 89 GHz passive microwave dataset with the spatial resolution of 3.125km, and it is downloaded from the following website: www.seaice.uni-bremen.de. We use the following framework for supervised sea ice classification from Sentinel-1 SAR data. Here we give a brief summary of the three involved steps (Figure 2), while the details will be explained in the following sections.
(I) Preprocessing: we apply a thermal noise removal and incident angle normalization on the dual polarization Sentinel-1 SAR images, before calculating the SPAN of the HH and HV channels (√ 2 + 2 ). The SPAN takes the 5 advantage of the two channels and is used in all following steps.
(II) Prior to classification we model the SAR data by combining a Mixture of three STAtistical (MSTA) distributions (Log-Normal, Rayleigh, and Alpha-Stable) with conditional random field (CRF) theory (MSTA-CRF). To reduce computation time we follow a hierarchical approach: the SAR images are first divided into several roughly homogeneous segments called superpixels based on mean shift theory (>5000 pixels at 40 × 40 m 2 ). Then we use another, more 10 restricted, mean shift procedure to segment the superpixel into even more homogeneous sub-superpixels (the mean size of the sub-superpixels contains 24 pixels, 0.038km 2 ). These sub-superpixels are the smallest entity our ice-water classification is based on.
(III) The sub-superpixel classification uses three weighted decision criteria referred to as potentials within the CRF context: The first potential uses a support vector machine (SVM) classification on all pixel in a sub-superpixel based on 15 their backscatter. As second potential, the similarity of the previously fitted MSTA backscatter distributions for the ice and water class to the backscatter distribution found in the sub-superpixel is used. The third potential describes the heterogeneity of the superpixel by using a fully connected network on the relationship of different sub-superpixels which accounts for their location and backscatter similarity. During training, we use square training samples as superpixels for each category (Figure 2 top-middle) and segment them into sub-superpixel. Then the parameters of the three potentials described above as well as their corresponding weights will be trained to generate the MSTA-CRF classifier for sea ice and water classification.

HH HV
Finally, we use a graph cut classifier based on these three potentials to maximize the probability of the MSTA-CRF 5 model. As a result each sub-superpixel is labeled with either the ice or water category. The details of MSTA-CRF including the above mentioned steps will be described in the following sections.  For C-band dual-polarization Sentinel-1 data, the backscatter coefficients at different incidence angles can vary significantly from the near to far range. For the HH channel of sample image in Figure 3 (a), backscattering coefficients at 10 near range is much larger than that at far range. In the HV channel, the dependence of backscatter on the incidence angle is not that obvious as it also contains bias of backscattering coefficient for different sub-swath. In addition to the angular dependence in the HH channel, contrast between open water and ice is much lower compared to the HV channel ( Figure   3(b)). Ice-water classification along marginal ice zone can be significantly influenced by winds and ocean waves due to wind-roughened water surface leading to much larger backscattering values. However, the thermal noise including the 15 scalloping noise and noise floor in the HV data limited its application for automatic analysis algorithm (Park et al., 2018).

Pre-processing of SAR data: incidence angle normalization and noise reduction
To reduce the effect of incidence angle, radiation calibration is required especially for the HH channel. Noise reduction is required for HV data (Figure 3 (b)). Preprocessing methods including incidence angle normalization and noise reduction methods are performed as described below.
The backscattering power from sea ice and open water is determined by several factors including orientation, slope, 20 roughness, surface cover and permittivity. The combination of these factors may lead to additional heterogeneity in vertical or horizontal directions, which can further impact the intensity and amplitude of the backscattering signal. The radar angular configuration, polarization, and frequency are the major instrument-related parameters. The radiometric angulardependence is corrected by a cosine law (Gauthier et al., 1998;Mladenova et al., 2013) using the local incidence angle.
This approach assumes that the amount of power that is backscattered towards the sensor direction is proportional to the 5 path length and thus to the cosine of the incidence angle, Since the radiation variance as a function of the observed area is also cosine function, the measured backscatter is described as: Where the measured backscatter 0 is related to the sea ice surface roughness dependence parameter n and it can be derived from the slope of the measure backscatter 0 and the local incidence angle by the linear regression model. and 0 0 are the reference incidence angle and nadir backscatter respectively. The normalized backscatter ref 0 is 10 defined by calculating the ratio of measured backscatter 0 at the reference angle of and incidence angle with a cosine function. This normalization method provides the observation at a constant incidence angle according to the application. To deal with the SAR images acquired in different conditions and compare the backscatter values of sea ice in the reference angles, absolute radiation calibration is carried out to convert the amplitude value of HH channel to 0 . Then the incidence angle normalization is utilized to transform the backscattering value into the reference angle of 23°, the 15 details of reference angle selection will be discussed in the experiment part in Section 4.1.
Due to the sparse ground control points in the XML file (121 points for nearly 400 km swath) of EW mode Sentinel-1 data, the noise floor of the neighboring sub-swaths in HV channel is discontinuous. The noise level of Sentinel-1 data before September 2016 and afterwards are different. The noise reduction mainly relies on the XML file provided by Sentinel-1 head file data and derived by the interpolation of the noise tie points for each pixel. Then the corrected SAR 20 backscatter is calculated by the subtraction of the backscatter values of the SAR and the noise pattern. 0 is the sigma naught in the look up table.
The subtraction of the noise of the entire image for each pixel from the digital number (DN) is divided by the backscatter coefficients provided in the XML file. With respect to the first sub-swath of each SAR scene, FFT transform is utilized to remove the remaining scalloping noise with high frequency (the amplitude of FFT transform of high frequency 25 is set to zero) (Choi et al., 2019). To prevent the yielding of noise floor of the neighboring sub-swath, a 10-pixel wide stripe of data along the edge of each sub-swath is not taken into consideration for the further analysis. Figure 4 illustrates the incidence angle correction and thermal noise removal results on two sample images. For better illustration, the contrast of the image is enhanced. During the melting seasons, the ice-covered areas, and the rough OW areas appear both bright in the HH channel and are therefore difficult to distinguish, which can be seen in Figure 4  also result in misclassification. Both HH and HV channels have been corrected as shown in Figure 4 (b) and (d). Angular dependence is corrected using the reference angle of 23°. The incidence angle correction will inevitably lead to the loss of some information and the contrast of the corrected image is decreased compared with the original image. Nevertheless, it shows a high distinguishability between sea ice and open water. For HV channel, the thermal noise caused by the discontinuity of noise floor of the neighboring sub-swath at far ranges have been removed. For the near range (first and 5 the second sub-swath), the overall image quality is greatly improved, although there is a small amount of noises remaining.
In the corrected HH channel, backscattering coefficients over wind-roughened water surface make it difficult to be distinguished from thin ice along marginal ice zone. In corrected HV, backscatter coefficients of thin ice are greatly reduced.
The Span image takes advantage of both HH and HV channel, with SPAN being defined as the total power calculated by the square root of HH and HV. The equation for SPAN is given as below, and are backscatter value of HH and 10 HV channel respectively.
It is clear to see that the rough area in the near range caused by the wind and ocean current is removed as shown in Figure 5 (e). Span data as combination of the HH and HV channels reduces the amount of data and thus is helpful for icewater classification. We will used the Span in the following for all classifications.  and (d) are the incidence angle normalization of HH channel and noise removal for HV channel respectively, (e) is the span image calculated by √ 2 + 2 .

Training samples selection 15
The training samples for sea ice and water classification is selected from the Sentinel-1 data acquired from 2015 to 2018 for the Sentinel-1 SAR imagery in the Fram Strait area. We first select one image for each day from 2015 to 2018 during June to September (488 images in total). Then 10 patches for each category (i.e., ice and water) are selected from the image dataset randomly using the MET Norway ice chart, and the size for each patch is 64*64 pixel. As a result, we https://doi.org/10.5194/tc-2021-85 Preprint. Discussion started: 30 March 2021 c Author(s) 2021. CC BY 4.0 License. got 9760 patches (4880 for ice and 4880for water), which is called the training dataset.
During the training step, we first select 100 patches for each categories to generate the MSTA-CRF model, each patch of the training samples is labeled pixel by pixel. This, e.g., means we know where the ice sub-superpixel is for an otheraise roughly homogenous water training sample. For the statistical distributions training the contribution of the ice subsuperpixel for the probability of the corresponding water superpixel is set to 0, while for the fully pairwise potential, the 5 ice sub-superpixel will be used to calculate the pairwise potential with the other sub-superpixels. The performance of the classifier is validated on all the training datesets. If the overall accuracy (OA) is lower than 99%, we add 100 patches (50 for ice and 50 for water) from the rest of the training dataset to train the revised model, and it will also be tested on the whole training dataset. For the sea ice classification in this paper, we repeat the training procedure for ten times and find that when the training samples reaches 100, the OA is over 99%, which means that 1000 patches are enough to get a 10 satisfied MSAT-CRF model, and it only occupies 10% of the full training dataset.

Mixture statistical distribution based fully connected CRFs
CRF (Conditional Random Fields) is a framework for constructing probabilistic models, which have been widely used in segmentation and classification (e.g., Tuia et al., 2018). For the operational ice-water classification, we define CRF on superpixels instead of performing a pixelwise classification. A superpixel is a group of homogeneous pixels that rendered 15 with uniform backscatter coefficients. These superpixels are obtained by the unsupervised mean-shift algorithm (Comaniciu and Meer, 2002). Generation of accurate superpixels is difficult, especially thin structures, such as leads are misclassified. We divide the superpixels into smaller sets as sub-superpixels using mean shift method. The subsuperpixels are now the smallest units we consider for the statistical distribution calculation. Then, we have the mixture statistical distribution based fully connected CRF as the form: 20 where − is the partition function. and are weight parameters of the MSTA-CRF model.
The unary potential is defined on the individual sub-superpixels at site . accounts for SAR backscatter coefficients pixel and denotes its category (ice or water).
( | ) deals with the single pixels individually and is represented as where δ( , ) is the Dirac function, where δ( , ) = 1 for = for the ice type, and δ( , l) = 0 for ≠ 25 with open water category. The local conditional distribution ( | ) is obtained via Support Vector Machine (SVM).
( | ) is calculated by statistical distribution models and the fully connected potential, which will be given in following section.

Mixture statistical distribution
The conditions of the scatters in a resolution cell can be treated as the sum of elementary scattering areas on a rough 30 surface due to its random walk characteristics in the complex plane (Goodman, 1976). The statistical model for speckle noise is on the assumption that each pixel contains a great number of scatters of radiation with a wavelength of C-band https://doi.org /10.5194/tc-2021-85 Preprint.  that is comparable to the roughness of the sea ice surface. For the lack of detailed information on the microscopic structure of the surface, the statistical attributions of the speckle patterns seem to be a better way to understand the characteristics of the radar signal. Backscatter signals of SAR image originally forms by the random scatters in the radar backscattering, and can be modeled by statistical distribution models, e.g. Rayleigh, Gamma, Log-normal, Weibull and Alpha-stable, which have been exploited to model the heavy-tailed and sharp-peaked statistical properties of SAR imagery under 5 complicated coherent noise condition (Yijing, et al., 2013).
We define the statistical potential ( | ) on sub-superpixels and model each superpixel with different statistical distributions. The parameter estimation methods for each statistical distribution are presented in Table3. Table 3. Single statistical distribution model and parameters estimation method used in this paper.

Statistical distribution model Parameters Method of Logarithmic Cumulants
Rayleigh  For each single statistical distribution, logarithmic cumulants k can be used to estimate the corresponding distribution parameters as shown in Table 3. The PDF reconstruction for each distribution is shown in Section 4.3.
With respect to the alpha-stable distribution without analytical expression, we present it as follows: where sign is the sign function. {α, β, γ, μ} are the parameters of the alpha-stable model. α is the characteristic exponent and β is the skewed parameter. γ is the dispersion parameter and μ stands for the location parameter. The 15 distribution parameters of the alpha-stable distribution are estimated by using the simulated annealing method (Grosse, 2007;Geman 1984).
With the increasing spatial resolution of SAR, image scene becomes more and more complex, to deal the SAR image with small inter-class differences, a mixture statistical distribution is integrated into the CRF framework to discriminate the ice category from open water using Sentinel-1 SAR images. After estimating parameters of a single statistical 20 distribution, we can formulate these statistical distributions in a weighted way as The weighting parameter θ~( , ), i = 1, … , M is estimated using an expectation-maximization (EM) algorithm (Saldju, 2000).

Fully connected CRFs
The unary potential of the CRF model can capture the dependence of the observed data. For each sub-superpixels, the fully connected network is constructed to represent the hidden information corresponding to the spatial and semantic relationship. Figure 5. Illustration of fully connected network of superpixel. The span image is first divided into superpixel using mean-shift method. For each superpixel, it is treated as the heterogeneous area and will be divided into sub-superpixel under more restrict mean-shift procedure. The relationship of different sub-superpixels which accounts for their location and backscatter similarity will be calculated in the fully connected network. Then we give the labels for each sub-superpixel and finally obtained the classification results.
For maintaining the accuracy of statistical distribution parameter estimation in the fully connected CRFs, the size of each superpixel is larger than 5,000 pixels. Each superpixel has a label to represent the dominating corresponding category of ice or OW. Since, however, there are still heterogeneous area within a certain superpixel, it will lead to a significant 5 misclassification in the results. Thus, each superpixel is divided into several sub-superpixel using a random number (less than 50), then we can get a much more homogeneous area within a certain superpixel by integrating the Gaussian potential in the fully connected framework. It may have different labels to represent details with in the superpixel, but for the calculation of the pairwise potential between different superpixels, they are treated as a homogeneous area as usual.
Local fully connected potential of CRF models ( , | ) imposes the spatial interaction and is defined as 10 and are the positions and the backscatter coefficients for each sub-superpixel respectively. The Gaussian kernel ( ) ( , ) includes an appearance kernel that favors similar backscatter coefficients pixels and a smoothness kernel that removes small isolated ice floes. When nearby pixels are assigned different labels, a penalty term ( , ) with 15 smoothness kernel is introduced and they can be calculated from the Potts model (Krahenbuhl and Koltun, 2012). The parameters and in the appearance kernel control the similarity and they can be learned from the training dataset.

Classification
In the CRF model, the final ice water classification is a maximization problem by assigning a label to each pixel of posterior probability. The global optimization is solved using the Graph cuts (Boykov et al., 2001) algorithm. As shown in 20 water) as background. Each vertex corresponds to a sub-superpixel segment and they are connected by edges with different weight probability, and the ensemble of the edge is named a "cut". Then, the minimum cut problem is to use the least cost of a cut among all cuts separating the vertices, which is equal to the lowest sum of edge weights. The s-t edge weights minimum problem is converted into solving the minimum of the energy function including unary potential f ( ), the statistical distribution potential f ( ) and the fully connected pairwise potential f ( ). It can be expressed as: 5 The inference procedure using the graph cut algorithm is exploited for indicating the label preference. We construct the graph with SVM potential, pairwise potential and the statistical distribution potential, shown in Figure 2. A S/T cut C on a graph with ice and water categories is used to partition the SAR image into two subsets, where S indicate ice label and T indicate water label. For the graph model, the cost of S/T nodes to the pixels is represented as a partial posterior marginal probability, which is derived from the unary potential and the statistical distribution potential in combination. 10 Pixel interaction cost is represented as a joint posterior marginal probability, which is derived from the fully connected pairwise potential. The optimization goal is to find the minimum cost of a cut C = {S, T}, which is a sum of all the costs from partial posterior marginal probabilities and joint posterior marginal probabilities. Then, the minimum cost of label S/T indicates the sub-superpixel to be given either the ice or water label.
In this study, the pre-processing of Sentinel-1 imagery including incidence angle correction and thermal noise removal

Selection of reference incidence angle 20
To choose the best reference angle for calculating the incidence angle corrected 0 , we utilized the MSTA-CRF to calculate the OA (overall accuracy) and STD on the example images in Figure 3. The ice water classification results using difference reference angle ranges from 20° to 40° are shown in Figure 6. We use CV (coefficient of variance) to evaluate the performance of different parameters (Table 4). A lower CV means classification result are more stable. For the reference angle of 21°, 23°, 26°, 29°, 30°, 35°, 36°, the averaged OA are larger than 90%. The reference angle of 26° gets the 25 largest OA but the STD is larger than 10%. As a tradeoff, the CV with the reference angle of 23° is chosen and used for the further sea ice water classification tasks.

Training samples selection for operational ice water classification
For ice-water classification, training samples are needed to train the MSTA-CRF classifier. It is, however, usually 30 difficult to select in-situ measurements as the training samples because they do not cover a significant amount of the satellite scene. For operational application, the method also needs to be stable in time and for different regions without the need for re-training. To deal with the problem of sample selection for ice-water classification, the stability of MSTA-CRF https://doi.org /10.5194/tc-2021-85 Preprint. Discussion started: 30 March 2021 c Author(s) 2021. CC BY 4.0 License. model using different training samples, as well as the performance of ice-water classification using a fixed model should be discussed. As a result, we would like to test whether the choice of different training samples may influence the final ice-water classification result. We also test which of the five statistical distribution are best suited for the MSTA-CRF and will keep only the significant ones.
In MSTA-CRF model, training samples mainly affect the calculation of statistical distribution parameters and the 5 corresponding weight coefficients. As a result, the dynamic of these parameters using different samples is a direct indicator to represent the stability of the adapted model. In this paper, we randomly select 82 models from 2015-2107. The sensitivity of parameter using different training samples is shown in Figure 6. It is clear to see that and distributions are sensitive in some cases in Figure 6(a), but the corresponding weight in Figure 6 Table 5.  Table 5 shows the average accuracy and STD using fixed model and adapted model respectively. The average accuracy using adapted model is slightly higher than for the fixed model since it contains the training samples from the test data, 10 and the STD of the adapted model is lower than for the fixed model. For operational ice-water classification, however, the main task is to obtain the ice type and its distribution in unknown scenes. The performance of fixed model is only slightly lower than the adapted model with 1.5% lower in average accuracy and 1.06% higher in STD. The results indicate that the MSTA-CRF is not sensitive to the change of training samples and thus it can be used as an operational method for icewater classification. 15

Comparison of different statistical distribution models
To evaluate the performance of different statistical distribution models, the reconstructed PDF (probability density function) and CDF (Cumulative Distribution Function) are compared with original SAR data for both ice and OW categories (Figure 7). The histograms for water (OW) and ice (bar plots) in Figure 7(a-f) are obtained by selecting 20 samples of each categories from the training dataset described in Section 3.2. For the PDF fitting curves in Figure 7 Rayleigh, Log-Normal, Alpha-Stable and the mixture model can get good correspondence for the ice histogram. For OW, all the five single distribution based models (a-e) underestimated the OW peak. The mixture model is the only one that can reconstruct the OW peak and shows the best agreement for both classes with an RMSD of about 0.03 (Table 6). From the CDF curves in Figure 7(h)-(i), we can see that the mixture model is most similar to the original SAR (followed by the Log-Normal distribution) data which means that the mixture model has a better performance than single distribution model. 25 The comparison of the Root Mean Square Deviation (RMSD) for the different distributions with the manually classified SAR backscatter (test data) is given in Table 6. The RMSD for the Ice class using Alpha-Stable, Log-Normal and mixture model are similar. For the OW class the RMSD using mixture model is the smallest compared with the other models. Thus we conclude that using a mixture model is beneficial for the ice-water classification as especially the backscatter of the water class can be much better modeled compared to a single distribution model. Since the PDF reconstruction of Gamma 30 and Weibull are not satisfying, we do not use them in the mixture model (see also next section). Thus also the corresponding https://doi.org/10.5194/tc-2021-85 Preprint. Discussion started: 30 March 2021 c Author(s) 2021. CC BY 4.0 License. statistical distribution based CRF algorithm using these two models are not discussed in the following part. 5.

Ice-water classification and validation
The example image shown in Figure 4 acquired on August 25, 2015 was used to evaluate the performance of the proposed MSTA-CRF algorithm. To validate the stability of the training dataset, the training samples were selected randomly from the dataset for 10 times, and the OA and the STD are reported in Table 6. We implemented the general CRF classifier, statistical distribution based CRF as well as the SVM classifier for comparison. The general CRF only contains 5 the unary and pairwise potential (see equation 5 in section 3), the statistical distribution based CRF adds the statistical distribution potential to the general CRF classifier. The SVM classifier transforms the data into a higher dimensional space and implements the "one-against-one" technique by separating the nonlinear data using a RBF(Radial Basis Function) kernel (Zajhvatkina et al., 2017). A features ensemble including the GLCM (gray level co-occurrence matrix) is input into the SVM classifier. Then ice-water classification is calculated pixel by pixel. In the experiment. The GLCM features 10 include mean value, standard deviation, energy, contrast, homogeneity, correlation and entropy. For the GLCM the backscatter values are discretized to 32 gray levels with the window size of 16 × 16 pixel. The separated distance is 8. As shown in Table 7, the MSTA-CRF achieves the highest classification accuracy of 92% among these algorithms.  (Table 6). Moreover, the variation of the OA and STD for the randomly selection training samples is significantly lower than for the other algorithms.

Ice-water classification comparison with ice charts
To validate the MSTA-CRF results we compare them to MET Norway ice charts, which is shown in Figure 9. For our purpose the ice chart classes "open drift ice", "close drift ice", "very close drift ice" and "fast ice" are defined as ice class.
And "open water" , "ice free" and "very open drift ice" are defined as water class. In Figure 9(c), it is clear to see that After the sub-superpixel based ice-water classification, which provides an ice or water class label for every 40x40 m pixel, we calculate the sea ice concentration (SIC) on a 1 km 2 grid. The percentage of ice pixels within a 25x25 cell provides the SIC with 1 km spatial resolution based on the MSTA-CRF model. The classification is based on sub-superpixel (on average 24 pixels), i.e., not every pixel is independent. This leads to a discretization of the SIC for the 1 km grid cells of about 2% to 4% (mind that many of the sub-superpixels are cut in smaller areas by the 1-km gridding). 5 The validation of MSTA-CRF sea ice concentration with ice chart on August 25, 2015 is shown in Figure 10. Similar to the ice-water classification comparison in Figure 9, the ice concentration errors are mainly located in boundary areas, and the overestimated and underestimated areas are mixed with each other Due to the lower spatial resolution of the ice

MSTA-CRF MET_ice_chart
https://doi.org /10.5194/tc-2021-85 Preprint. Discussion started: 30 March 2021 c Author(s) 2021. CC BY 4.0 License. chart, some detailed information are lost. Besides, sea ice drift, which can be pronounced during melting season, may also lead to the difference between classification results and ice chart because their acquisition may be different.
The ice edge difference between the MSTA-CRF SIC result and the MET Norway ice chart are shown in Figure 11.
The MSTA-CRF ice edge is defined as the 15% SIC isocline, and the MET ice edge indicates the boundary line between the "ice" and "water" category, which is defined in Table 2. It is clear to see that these two edges are matching very well 5 in most cases. Only near [81.5N, 5E] and [79.5N, 12.5W] a slight difference can be seen. It also demonstrates that the proposed method is effective for sea ice and water classification during melting seasons. The detailed daily IIEE (Integrated Ice Edge Error) (Goessling, et. al., 2016) analysis will be discussed in Section 5.3.

Monthly Average Classification results 10
We now apply the trained MSTA-CRF classifier to all Sentinel-1 SAR in our Fram Strait target region (Figure 1) for the six years 2015 to 2020 and calculate the accuracy by comparison to the daily MET Norway ice charts (Section 2.2.b).
The monthly combined sea ice and water classification result in Figure 12 shows that the average accuracy for MAST-CRF algorithm is larger than 85%, and the STD is less than 10%, among which the OA in 2018 and 2020 is the highest of over 90%. We can also find that from the months June to September the accuracy in August always is the highest during 15 the six years besides in 2015. Compared with the "OW" error, the "Ice" error is larger. The main misclassified area originates from melting water on fast ice, which occasionally is misclassified as water. Since we defined the very open drift ice as OW category, it may also lead to some misclassification. Besides the MET Norway ice charts, the ASI SIC product from AMSR2 on a 3.125 km grid is also used to validate the performance of MSTA-CRF algorithm. For the comparison, the MSTA-CRF SIC is calculated on a 1 km grid from the ice-water classification (see section 4.4.3). Figure 13 gives two example results of SIC in the same area in 2020. One is acquired on June 1 and the other is on September 29, and these two images are for the same orbits, i.e., show the same area. The MSTA-CRF sea ice concentration contains much more detailed information of sea ice conditions with the spatial 5 resolution of 1 km. Especially for the lower SIC areas of leads and the marginal ice zone (MIZ), the MSTA-CRF SIC can detected more ice fragments. Moreover, MSTA-CRF SIC can also identify fragmented ice areas with many leads and lower SIC (central area in Figure 13e), which is some distance away from the sea ice edge. Overall, the SIC of the two products show very good agreement. This is confirmed in the scatter plots in Figure 13c and 13e, which show that the two SIC results match well near 0% and 100%, and the difference of these two SIC products is less than 10% on average. Most of 10 these areas are in the marginal ice zone or sea ice leads, where the MSTA-CRF SAR SIC has the advantage of a higher spatial resolution and can reproduce a higher SIC variability. There is also a tendency that the difference of classification result is caused by very open drift ice and open drift ice like for the MET Norway ice charts (in Figure 9 and 11). In conclusion, the MSTA-CRF using high resolution Sentinel-1 SAR data can provide much more detailed information of sea ice conditions, while still reproducing similar SIC values as the ASI AMSR2 SIC. 15 Figure 13. Comparison of ASI sea ice concentration and MSTA CRF derived sea ice concentration. The ASI SIC in (a) and (f) are 3.125km resolution, MSTA-CRF SIC in (b) and (e) are 1km resolution and the density of scatter map is in 3.125km by down sampling the MSTA-CRF SIC into 3.125km. (a) Daily sea ice area compared with ASI sea ice concentration (b) Daily IIEE compared with MET Noway ice chart (c) Monthly sea ice area Here we discuss the daily MSTA-CRF SAR SIC time series for the Fram Strait region during summer for the six years 2015 to 2020 and compare it with the ASI AMSR2 SIC (Figure 14). In September 2020, the Arctic has experienced the lowest sea ice extents since the record minimum year of 2012 (https://nsidc.org/arcticseaicenews/charctic-interactivesea-ice-graph/). For the Fram Strait region ([75N 83N], [-15W 15E]; Figure 1), however, the temporal development and 5 seasonal cycle are different. The time series of daily sea ice area in Figure 14a indicate that the sea ice minimum for this region appeared in late August, 2018. During the four months June to September the sea ice area in Fram Strait usually is highest in late June and then decreases until its minimum in mid-August and then increase again from late August to September. Compared with the sea ice area from the ASI SIC product, the MSTA-CRF result is lower in many cases from June to mid July. While during the sea ice decrease period from July to early September, sea ice area calculated from 10 MSTA-CRF is higher than ASI, the main difference is that thin ice may be misidentified as lower concentration using ASI method as well as the lower resolution of the dataset. The two results also shows a significant agreement with the 2 of 0.9976. For the monthly sea ice area in 2018, the sea ice area in August only reaches the half of 2016. that the IIEE compared with MET Norway is can be as large as over 3000km 2 in the beginning of June, and then it decrease until the late August. After that, it stars to increase again. The averaged IIEE only 3300km, which is below 2% of the total sea ice area in the Fram Strait. The monthly average result in Figure 14c shows that sea ice area did not change so much in 2015-2017 and 2019-2020, while it has a significant decrease in 2018, it has decreased 28.75% compared with sea ice area in 2019. The main reason for the sea ice minimal in 2018 may be due to the sea ice export from Fram Strait to lower 5 latitude area which is caused by sea ice drifting.

Summary and Conclusions
We have developed an algorithm for ice-water classification and sea ice concentration (SIC) retrieval using Sentinel-1 Extra Wide swath (EW) mode data acquired over the Fram Strait region during melting seasons (June -September). 10 Especially, during this time of the year sea ice surface melt and wind-roughened surface conditions make ice-water classification challenging. The classification scheme utilizes the SAR backscatter after corrections, which includes incidence angle normalization and thermal noise reduction in a MSTA-CRF (mixture statistical distribution based conditional random fields) framework. The mixture statistical distribution (Alpha-Stable, Log-Normal and Rayleigh) is used to construct an additional probability potential for the classification in addition to the SVM potential and the fully 15 connected pairwise potential. The mixture of statistical distributions can serve as a proxy for turbulent ocean wave conditions or calm water condition and to characterize ice and water. For each sub-superpixel, the proposed MSTA-CRF constructs a fully connected network to eliminate the uncertainty of small fragments of ice caused ice break-up or by the melting water on ice-covered areas during melting seasons. The proposed algorithm has been tested on Sentinel-1 dual polarization SAR scenes for the Fram Strait region. The comparison with MET Norway ice charts shows that the proposed 20 algorithm reaches an overall classification accuracy (OA) of about 90% with STD less than 10%.
To construct a reliable training dataset, 9760 patches from 488 S1 images in 2015-2018 are used for training, we first select twenty samples (ten for sea ice and ten for open water) with the size of 64×64 pixels to construct the training dataset.
The proposed algorithm can distinguish between small open water leads and sea ice and very well preserves the ice edge (averaged integrated ice edge error of below 2% of total ice area) by selecting 1000 patches from the training dataset. 25 Validation of the ice-water classification result was conducted by comparing with ice charts. The results indicate that the main error is open drift ice misclassified as very open drift ice (defined as "water" in this paper). This may have caused some of the discrepancies between our classification results and the MET Norway ice charts.
The acquisition scheme of the Sentinel-1a (launch 04/2014) and b (launch 04/2016) constellation allows very frequent coverage of the Fram Strait region. We have classified the sea ice area from Sentinel-1 almost every day for the six years 30 2015 to 2020 for the four summer months June to September. The sea ice area in Fram Strait shows a clear seasonal cycle in all years with the sea ice area minimum in August about one month earlier than the Arctic-wide sea ice minimum. Out of the four years the year 2018 clearly shows the lowest sea ice area during summer with an especially pronounced decrease from late August to September in 2018. In comparison to ASI AMSR2 microwave radiometer SIC our S1 SIC shows a very similar seasonal cycle and close agreement in SIC values. However, the 1-km SAR SIC provides much more details, winter in a future study the SAR SIC can be compared with SIC from infrared imagers (e.g. Ludwig et al., 2020), which provide a similar spatial resolution but are not available during summer. To extend the MSTA-CRF S1 classification also to areas that are not as well covered by S1 acquisitions as the Fram Strait region, in future, they can be combined with passive microwave SIC datasets in an active-passive microwave fusion framework.

5
Data availability. All processed data can be obtained by contacting the first author. The SAR data used in this paper are downloaded from the European Space Agency website (https://scihub.esa.intdhus/). The Norwegian Ice Service at the Norwegian Meteorological Institute (MET Norway) provided daily sea ice charts for validation (delivered by Nick Hughes).
The daily averaged wind speed, wave height and ice surface temperature are downloaded from ECMWF website.

10
Author contributions. The concept of the study was conceived by YZ and TZ. The SAR processing was done and reported jointly by YZ and TZ. CM, GS, and MH commented the SAR processing and methodology. NH provided the MET Norway ice charts and commented on the experiment part. YZ and GS wrote the conclusions section. All authors commented on the results. All authors contributed to the editing of the text and agreed to the submission.

15
Competing interests. The authors declare that they have no competing interest.