Arctic Sea Ice Thickness Estimation from CryoSat-2 Satellite Data Using Machine Learning-Based Lead Detection

Abstract

Satellite altimeters have been used to monitor Arctic sea ice thickness since the early 2000s. In order to estimate sea ice thickness from satellite altimeter data, leads (i.e., cracks between ice floes) should first be identified for the calculation of sea ice freeboard. In this study, we proposed novel approaches for lead detection using two machine learning algorithms: decision trees and random forest. CryoSat-2 satellite data collected in March and April of 2011–2014 over the Arctic region were used to extract waveform parameters that show the characteristics of leads, ice floes and ocean, including stack standard deviation, stack skewness, stack kurtosis, pulse peakiness and backscatter sigma-0. The parameters were used to identify leads in the machine learning models. Results show that the proposed approaches, with overall accuracy >90%, produced much better performance than existing lead detection methods based on simple thresholding approaches. Sea ice thickness estimated based on the machine learning-detected leads was compared to the averaged Airborne Electromagnetic (AEM)-bird data collected over two days during the CryoSat Validation experiment (CryoVex) field campaign in April 2011. This comparison showed that the proposed machine learning methods had better performance (up to r = 0.83 and Root Mean Square Error (RMSE) = 0.29 m) compared to thickness estimation based on existing lead detection methods (RMSE = 0.86–0.93 m). Sea ice thickness based on the machine learning approaches showed a consistent decline from 2011–2013 and rebounded in 2014.

1. Introduction

Sea ice impacts the Earth’s radiation balance because thermal feedback between the Sun and the Earth is highly sensitive to sea ice reflectivity. Thus, Arctic sea ice is considered an important factor in understanding the global climate change process [1]. The reflectivity of sea ice strongly depends on the spatial distribution and extent of the ice [2,3], which have rapidly changed due to global warming over the past two decades [4,5]. Boe et al. [6] used various climate model simulations to predict that the Arctic Ocean would probably be ice-free by the end of the 21st century. Furthermore, several studies have shown that the decline of sea ice is occurring faster than model predictions [7,8]. Thus, there is an increasing need for accurate monitoring of sea ice concentration and thickness to better understand polar and global climate systems and processes.

Sea ice thickness has been measured with various methods. While direct field measurements of sea ice thickness using a submarine upward looking sonar [9,10] or an electromagnetic system (e.g., Airborne Electromagnetic (AEM)-31) [11,12,13] can provide accurate ice thickness information, such techniques can only be applied to local areas within a very limited time frame. Observation of sea ice thickness over vast areas has utilized various space-borne radar and laser altimeter sensors [14,15,16,17,18,19,20]. Laxon et al. [14] retrieved sea ice thickness using European Remote Sensing Satellite-1 (ERS-1) and ERS-2 satellite data based on radar altimetry. Kwok et al. [15] also estimated sea ice thickness using Ice, Cloud and land Elevation Satellite (ICESat) data based on laser altimetry. Unfortunately, the operation of ICESat stopped in 2009 due to the failure of its main instrument. Since the launch of CryoSat-2 in 2010, researchers have developed various methods to use the radar altimetry observations to estimate sea ice thickness from CryoSat-2 data [16,17,18,19].

Sea ice thickness can be estimated from sea ice freeboard based on isostasy [21]. Derivation of sea ice freeboard is an important procedure for estimating the ice thickness by laser or radar altimeter measurements. In particular, identification of leads (i.e., fractures between sea ice floes) is crucial to estimate the freeboard. The height of leads extracted by such measurements enables the calculation of the Local Sea Surface Height (LSSH), and then, freeboard can be estimated using LSSH, actual sea surface height and the surface elevation of the ice extracted by altimetric measurements. Kwok et al. [15] detected leads through direct comparison between the surface elevation profiles extracted by ICESat data and near-coincident Synthetic Aperture Radar (SAR) images. Zwally et al. [22] assumed that the lowest 2% values of the surface elevation profiles from ICESat would correspond to leads. In addition to these relatively simple methods, Farrell et al. [23] proposed a threshold-based method to distinguish leads from ice floes using various parameters extracted from ICESat level 1b data, such as gain, reflectivity, radiance and waveform characteristics. In the case of CryoSat-2, Pulse Peakiness (PP) and Stack Standard Deviation (SSD) parameters are frequently used for lead detection. For example, Ricker et al. [19] used various waveform parameters, such as PP, SSD, stack kurtosis and sea ice concentration, to distinguish leads from ice floes. Although these lead detection methods have been developed in several studies, the determination of ice thickness from CryoSat-2 still suffers from a lack of precise lead discrimination [24]. Simple thresholding methods might not perfectly distinguish leads from ice floes because parameters, such as PP, SSD, stack skewness, stack kurtosis and backscatter sigma-0 (Section 2.1) typically contain aliasing between leads and ice floes, which can result in large errors and uncertainties in sea ice thickness estimates. Therefore, advanced techniques to optimize such thresholds and minimize the associated errors are needed. This study proposes decision trees and random forest machine learning approaches to identify leads and ice floes from CryoSat-2 and Moderate Resolution Imaging Spectroradiometer (MODIS) in order to estimate sea ice thickness.

2. Observational Datasets

Freeboard height and ice thickness for March and April in 2011–2014 were calculated from CryoSat-2 data based on machine learning-based lead detection approaches. MODIS and sea ice type data were used as ancillary data when estimating sea ice thickness. The estimated ice thickness was validated using CryoSat Validation experiment (CryoVex) field campaign data (i.e., airborne electromagnetics data) over northwestern Greenland acquired in April 2011.

2.1. CryoSat-2

CryoSat-2 was launched in April 2010 and carries the space-borne Synthetic Aperture Interferometric Radar Altimeter (SIRAL) developed by the European Space Agency (ESA) [25]. SIRAL has a center frequency of 13.575 GHz (K_u-band) and a bandwidth of 320 MHz. It has three operation modes: Low Resolution Mode (LRM), Synthetic Aperture Radar (SAR) and SAR Interferometry (SIN). ESA explains that data collected in SAR and SIN modes are optimized to estimate sea ice thickness because the sensor in the operation modes can measure sea ice characteristics with high spatial resolution comparable to the size of leads [26]. In the CryoSat-2 waveform data, the power of the received microwave signal is recorded in 128 range bins in SAR mode and 512 range bins in SIN mode. The interval of each range bin is almost 1.563 ns (~0.234 m). Detailed specifications of CryoSat-2 are presented in

Table 1. The specifications of CryoSat-2 (Synthetic Aperture Interferometric Radar Altimeter (SIRAL). LRM, Low Resolution Mode; SIN, SAR Interferometry.

CryoSat-2
Center frequency	13.575 GHz
Bandwidth	320 MHz
Pulse Repetition Frequency (PRF)	1.97 kHz (LRM)/18.181 kHz (SAR and SIN)
Pulse duration	44.8 ms
Samples in echo	128 (LRM and SAR)/512 (SIN)
Antenna footprint	0.29 km
Range bin sample	0.4684 (LRM)/0.2342 m (SAR and SIN)

.

CryoSat-2 (Baseline B) and CryoSat-2 Level 1B (L1B) waveform data were used to estimate surface elevation (ftp://science-pds.cryosat.esa.int). Five parameters (i.e., SSD, stack skewness, stack kurtosis, PP and backscatter sigma-0) were used to distinguish leads from ice floes and ocean, as they can represent surface characteristics, such as surface roughness and dielectric property. SSD was available from L1B data; stack skewness and stack kurtosis were available from Level 2I (L2I) data, which are provided by ESA. SSD is the variation of the stacked power distribution with an incidence angle [26]. Stack skewness and stack kurtosis measure the asymmetry and peakedness of the range stacked power distribution, respectively [25]. PP is commonly used to identify leads and ice floes [16,18,19,27]. The equations to calculate SSD, stack skewness and kurtosis, as well as PP are summarized in

Table 2. The equations to retrieve Stack Standard Deviation (SSD), stack skewness and kurtosis, as well as Pulse Peakiness (PP) [26,28].

Parameter	Equation
SSD	$\frac{1}{2}\frac{{\sum }_{i=1}^{N}S{P}^2\left(i\right){\sum }_{i=1}^{N}S{P}^2\left(i\right)}{{\sum }_{i=1}^{N}S{P}^4\left(i\right)}$
Stack skewness	$\frac{\frac{1}{N}{\sum }_{i=1}^N{\left(SP\left(i\right)-\mu \right)}^3}{{[\frac{1}{N-1}{\sum }_{i=1}^N{\left(SP\left(i\right)-\mu \right)}^2]}^{2/3}}$, $\mathsf{\mu }=\frac{1}{N}\sum _{i=1}^{N}SP\left(i\right)$
Stack kurtosis	$\frac{\frac{1}{N}{\sum }_{i=1}^N{\left(SP\left(i\right)-\mu \right)}^4}{{[\frac{1}{N-1}{\sum }_{i=1}^N{\left(SP\left(i\right)-\mu \right)}^2]}^3}-3$, $\mathsf{\mu }=\frac{1}{N}\sum _{i=1}^{N}SP\left(i\right)$
Pulse peakiness	$\frac{k \text{x }{P}_{max}}{{\sum }_{i=1}^n{p}_i}$ n = 128 (SAR) and 512 (SIN)

where SP stands for integrated stacked power that is not obtainable in the L1b data. The integrated stacked power is the summation of each single look echo power. $p_{max}$ is the maximum power of the waveform from L1b data, and $p_i$ is the power of i-th range bin. k is a multiplying factor based on the assumption that the waveform is almost centered in the range bins. A k value of 1 was used in this study following [29].

.

The radar backscatter sigma-0 (i.e., backscatter coefficient) from Level 2 (L2) data, documenting the observed surface, is a function of dielectric properties, the radar frequency, incidence angle, the target surface roughness, geometric shape and volume scattering [25]. The SAR L1b waveforms can be converted into watts using power scaling parameters that are available in the L1b product. The radar equation is solved using transmit power, range and instrument gain and bias correction to retrieve backscatter sigma-0. A bias correction value is then applied to remove any residual bias [30]. Since these parameters are sensitive to change in surface condition, they can be used to discriminate leads from ice floes.

2.2. MODIS

MODIS onboard the Terra and Aqua satellites, which were launched in 1999 and 2002, respectively, has 36 spectral bands from 0.4–14.4$\mathsf{\mu }$m and plays a vital role in observing the Earth’s environments, such as the land, lower atmosphere and oceans. MODIS images are an ideal way to separate leads and ice floes because of the albedo difference between the two. MOD02QKM, one of the MODIS L1B products, is a calibrated and geolocated dataset with two bands (0.645$\mathsf{\mu }$m and 0.858$\mathsf{\mu }$m) at a 250-m ground sample distance. In this study, training data of leads, ice and ocean for the machine learning models were extracted from MOD02QKM images through visual interpretation based on reflectance differences.

2.3. Sea Ice Type

The European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) Ocean & Sea Ice Satellite Application Facility (OSI SAF) provides sea ice type data (http://osisaf.met.no/p/ice/) with 10-km resolution. The sea ice type includes First-Year Ice (FYI) and Multi-Year Ice (MYI) based on differences in ice surface roughness. The sea ice type was used as an input variable to calculate sea ice thickness from ice freeboard.

2.4. Airborne Electromagnetics Data

The CryoVex airborne and field campaign was conducted to validate the measurements of CryoSat-2. As a part of the campaign, sea ice thickness was measured with an Airborne Electromagnetic (AEM)-bird sensor onboard the Alfred-Wegener Institutes (AWI) Polar-5 aircraft. AEM uses electric conductivity differences between sea water and ice to measure sea ice thickness with an accuracy of $\pm$0.1 m over level ice [13,31]. From 14–17 April 2011, AEM measured four tracks of ice thickness around the Lincoln Sea. Considering the length of the tracks and sample size, two of the data tracks were used to validate the sea ice thickness estimated from CryoSat-2 in this study.

3. Sea Ice Thickness Estimation and Machine Learning Algorithms for Lead Detection

3.1. Sea Ice Thickness Estimation

Estimation of the snow-covered Arctic sea ice thickness from CryoSat-2 measurements is based on the assumption of hydrostatic equilibrium [18] (Figure 1). If the sea ice freeboard ($h_{fb}$) is accurately determined from altimeter measurements, the freeboard can be directly converted into sea ice thickness by Equation (1).

$$h_{si}=\frac{{\rho }_{sw }}{{\rho }_{sw}-{\rho }_{si}}{h}_{fb}+\frac{{\rho }_{s }}{{\rho }_{sw}-{\rho }_{si}}{h}_s$$

(1)

where ${\rho }_{sw}$, ${\rho }_{si}$ and ${\rho }_s$ are the density of sea water, sea ice and snow, respectively, and $h_s$ is the snow depth. Although the density parameters and snow depth should be observed concurrently with the altimeter measurements to best estimate the ice thickness, this is challenging due to the extreme weather conditions of the Arctic Ocean. Thus, studies have used typical values based on field measurements or numerical simulation. For example, Giles et al. [32] and Wadhams [10] used the density of sea water, sea ice and snow as 1023.8 ± 3, 915.1 ± 5 and 319.5 ± 3 kg/m³, respectively, from field observations. In this study, 916.7 kg/m³ and 882 kg/m³ were used as the density of FYI and MYI, respectively, according to Alexandrov et al. [33]. Snow depth simulated by the Warren 99 (hereafter W99) climatology model of Warren et al. [34] has been widely applied. However, the original W99 data only capture the seasonal variability of snow depth. Kurtz and Farrell [35] thus applied a modification to the snow depth data to reflect the significant decline in MYI over the last few years. Kurtz and Farrell [35] suggested reducing the snow depth over FYI by 50%. In this study, FYI and MYI were discriminated by the ice type products derived by EUMETSAT OSI SAF, and we used the typical density and snow depth values derived by Kurtz and Farrell [35].

As mentioned above, it is important to determine the sea ice freeboard from CryoSat-2 data in order to successfully estimate the ice thickness. Figure 2 shows the procedure to determine the sea ice freeboard from CryoSat-2 L1B data. Initially, the surface height of the sea ice (i.e., the distance between the sea ice surface and the WGS84 ellipsoid) is estimated by Equation (2).

$${\eta }_{sea ice}=H_{sat}-R_{win}-R_{err}-\mathsf{\Delta }R$$

(2)

where $H_{sat}$ is the height of the satellite platform mass above the WGS84 ellipsoid. $R_{win}$ is the window delay field, which means the distance between the mid-point of the range bin (i.e., 64th range bin in SAR mode and 256th range bin in SIN mode) of the waveform data and the satellite platform. $R_{err}$ is a range correction term associated with the phase range due to geophysical properties, such as atmospheric effects. These variables are given in CryoSat-2 L1B data; detailed descriptions and processing methods of the variables are well explained in [25]. $\mathsf{\Delta }R$ is another correction term derived by various retracking methods [36,37,38,39]. The aim of these terms is to determine the range offset between the mid-point of the range bin and a realistic range point of the leading edge of sea ice. The retracking method used in this study is the Threshold First Maximum Retracker Algorithm (TFMRA) [16,18,19].

Davis [37] introduced the threshold retracking concept, which is useful for measuring the surface elevations of ice sheets or sea ice from radar altimeter data [14,19]. The retracking method determines the range point of the leading edge between the threshold level and the range point of the first maximum power peak. The threshold level (${\rho }_l$) is determined by:

$${\rho }_l=P_n+\alpha \left(P_{max}-P_n\right) here, P_n=\frac{1}{5}{\displaystyle \sum }_{i=\overline{n}}^{\overline{n}+4}{P}_i$$

(3)

where $P_n$ is the thermal noise of the CryoSat-2 system. $\alpha$ is the threshold value, the percentage of the maximum waveform amplitude above the thermal noise. $P_{max}$ is the first maximum power of the waveform. $\overline{n}$ is the range bin of the first unaliased waveform. $P_i$ is the power at the $i$-th range bin of the waveform. Finally, the retracking point ($n_r$) as the leading edge is estimated using the following equation.

$$n_r=\left(\widehat{n}-1\right)+\frac{{\rho }_l-P_{\widehat{n}-1}}{P_{\widehat{n}}-P_{\widehat{n}-1}}$$

(4)

where $\widehat{n}$ is the first range point exceeding the threshold level. It is essential to detect the first peak in the range bin. Here, Rose [16] indicated that the maximum power peak in the range bin may not be the first peak due to time delay effects of complicating factors, such as multiple scattering (i.e., in the surface) and volume scattering. Thus, the true range point (i.e., local maxima in the waveform) is detected by the peak detection algorithm, which identifies the range point using derivatives of the waveform signal [19,40]. Lastly, $\mathsf{\Delta }R$, the retracking correction, is calculated using Equation (5).

$$\mathsf{\Delta }R=C_{2m}\left(n_r- n_{tr}\right)$$

(5)

where $n_r$ is the retracking point and $n_{tr}$ is an on-board retracking point. $C_{2m}$ is a factor to convert from range bins to meters, which is 23.24 cm/bin for CryoSat-2. Figure 3 shows an example of the peak detection algorithm. This figure illustrates that the range point of the first maximum power (the open square) was found prior to the maximum power peak, and the retracking point (the open circle) was determined between the range point of the first maximum power and the threshold level (the dotted line). While various threshold values$\left(\alpha \right)$ have been used in the literature, several studies have found that thresholds of 40% and 50% give the best result for determining the leading edge of the ice floe [16,40]. A threshold of 40% was used in the retracking method in this study as it was frequently used in the literature [16,19].

The next step removes the distance between the actual sea surface and the WGS84 ellipsoid from the surface height (η_{sea_ice}) in order to estimate the sea ice freeboard. In general, the actual Sea Surface Height (SSH) is estimated from the sum of the mean SSH and local Sea Surface Height Anomaly (SSHA). Mean SSH data were obtained from the Technical University of Denmark 10 (DTU10) product [41]; local SSHA data were derived from the proposed lead detection method (Section 3.2) that extracts leads from CryoSat-2 data. The SSHA observations were discontinuous because leads were detected at irregular intervals, thus linear interpolation and low-pass filtering were applied to make spatially-continuous SSHA. The local SSHA was used to remove the surface height from the mean SSH at the leads. Although many studies have tried to develop effective lead detection methods, it is still very difficult to accurately identify leads due to limited reference data, the irregular shape and size of leads and the characteristics leads share with ice or ocean. To overcome these challenges and correctly identify leads, this study proposed a novel lead detection method explained in Section 3.2.

A correction to the sea ice thickness estimates from freeboard should be applied to account for the penetration of microwave radiation on snow and lower propagation speed in the snow pack. First of all, while typical microwave pulses do not penetrate the snow surface when the snow layer is wet during the melting season, it is well known that a K_u-band microwave penetrates the air/snow interface of dry and cold snow during the freezing season [42,43,44]. This complexity makes it difficult to determine the optimum penetration depth to correct the sea ice freeboard. Nevertheless, Laxon et al. [18] believed that microwaves fully penetrate the snow layer. Since the speed of microwave is typically lower in the snow pack [45], it should also be corrected. However, given the uncertainty in these corrections, we did not apply the correction terms in order to enable consistent comparison with Laxon et al. [18] and Ricker et al. [19], who did not apply these corrections.

3.2. Machine Learning Algorithms for Lead Detection

In order to detect leads using machine learning approaches, reference samples were extracted using MODIS data. All 5-min MOD02QKM images above latitude 65°N in March and April 2011–2014 were downloaded. Cloud-free images were selected through visual interpretation. A total of 48 cloud-free March and April images were selected from MOD02QKM between 2011 and 2014 to clearly identify sea ice, leads and ocean based on visual interpretation of the images (Figure 4). However, visual interpretation with MODIS is not always reliable because the leads in the MODIS images could refreeze, and new thin ice is formed. CryoSat-2 paths were geolocated over the MODIS images to extract five parameters (i.e., SSD, stack skewness, stack kurtosis, PP and backscatter sigma-0) for each class (i.e., lead, sea ice and ocean). The time difference between CryoSat-2 paths and MODIS images was set to within 30 min (12 min on average) to minimize sampling errors as sea ice sometimes moves fast. Since there were more leads found in the Arctic in April than March, the number of samples for April was larger than that for March. It should be noted that we could not extract samples all over the Arctic Region because spatiotemporal coincidence between CryoSat-2 and MODIS was limited during the given time period. Lead reference samples were not collected when the size of leads was smaller than 250 m considering the movement velocity of sea ice.

Since the characteristics of the sea ice surface have monthly and annual variations, three schemes were examined to develop machine learning models in this study. The first scheme was Classification of Monthly data (CM), which used the reference samples by month regardless of year and developed the machine learning models for both months (i.e., March and April). The second was Classification of Annual data (CA), which divided the samples by year and developed the models separately for each year (i.e., 2011, 2012, 2013 and 2014). Individual Classifications (IC) used all reference data to develop the machine learning models to consider the tradeoff between transferability and accuracy.

In order to detect leads, we used two rule-based machine learning approaches: decision trees and random forest. Decision trees are one of the most widely-used machine learning algorithms for inductive inference [46,47,48]. To implement decision trees, See 5.0 was used. See 5.0 recursively splits training data into subdivisions based on a set of attributes defined at each node in a tree [49]. An attribute is selected at each node and two branches that descend from that node use a value of the attribute as a threshold. Selecting an attribute (i.e., STD, stack skewness, stack kurtosis, PP or backscatter sigma-0 in this study) at each node is crucial for successful classification. In general, statistical properties, such as information gain or the Gini index, are used to choose an appropriate attribute in decision trees. See 5.0 uses information gain to select which candidate attribute is used at each node. See 5.0 has been widely used for various remote sensing applications, including land cover/land use classification [50,51,52], climate region delineation [53], vegetation species mapping [54,55], ice mapping [56] and change detection [57,58]. Using a See 5.0 decision tree has some advantages. First, it provides a non-parametric classification, and thus, it does not require any assumptions in terms of the distribution of training data. See 5.0 can also handle non-linear relationships between classes and features, even with missing values. In addition, See 5.0 transforms a decision tree into a series of production rulesets, which makes it easier and more straightforward for human interpretation of the results.

Random forest uses an ensemble approach that combines a boosting sampling strategy and Classification And Regression Trees (CART) [59] to improve the weaknesses of a single CART such as overfitting and sensitivity to training data configuration. CART uses a Gini index to measure impurity from training samples, while See 5.0 uses the concept of entropy. The Gini index is defined as shown in Equation (6)

$$\text{Gini index}\left(\mathrm{S}\right)=1-{\sum }_{i=1}^c{p}_i^2$$

(6)

where c is the number of classes and $p_i$ is the proportion of S belonging to class i. The Gini gain is used to identify the most appropriate attribute at each node. Since it is similar to the information gain, it is defined by replacing the entropy with the Gini index in the Equation (6). However, a single CART is often unstable and tends to overfit training data. Bagging can overcome such weaknesses by creating n independent trees and help minimize errors that can be caused from unstable classifiers [60]. Random forest produces numerous independent trees through two bagging-based randomization processes: (1) using a random subset of training data for each tree; and (2) using a random subset of input variables at each node of a tree. Breiman [58] pointed out that it is not necessary to use a separate dataset for model validation, as random forest uses out-of-bag data (i.e., training data that are not used) for internal cross-validation. A majority voting strategy is used to combine the results from multiple classifiers to determine the final class for a given sample. In addition, random forest provides the relative importance of a variable using out-of-bag data when the variable is permuted. Because of these strengths, random forest has proven robust in various remote sensing applications [61,62,63,64,65,66,67,68].

4. Results and Discussion

4.1. Typical Waveform over Leads, Ice Floes and Ocean

Radar signals of each of the three target features have different characteristics because of the impact of several factors, especially surface roughness, on the signals. In particular, flat surfaces produce strong signals, and rough surfaces produce weak signals. The shape of the typical waveform of ice floes is similar to that of ocean (Figure 5a,b). The sea ice waveform has large variation as it contains both diffuse (e.g., from ice floes and ridges) and specular (e.g., from lead and new ice) signals. In particular, since the surface of MYI is rougher than that of FYI, more diffuse reflection occurs on MYI. Leads have a typical specular reflection and symmetric waveform because they are relatively flat and there is little surface wave in leads (Figure 5c).

4.2. Characteristics of Five Parameters Based on CryoSat-2 Waveform

Figure 6 depicts the box plots of the five parameters (i.e., SSD, stack skewness, stack kurtosis, PP and backscatter sigma-0) by feature (i.e., leads, sea ice and ocean) using the reference samples (refer to

Table 3. Reference data used in the machine learning models by scheme and target feature. CM, Classification of Monthly data; CA, Classification of Annual data; IC, Individual Classification.

Scheme	Target Feature (Number of Observations)
	Leads	Sea Ice	Ocean
CM (March)	331	660	724
CM (April)	641	1284	1220
CA (2011)	179	357	357
CA (2012)	458	919	919
CA (2013)	209	419	420
CA (2014)	126	249	248
IC	972	1944	1944

). Among the three target features, ocean showed the narrowest distribution for all parameters, except SSD. This is because the ocean surface is relatively homogeneous. Since the state of the sea ice surface varies significantly, all parameters resulted in a wide distribution in the sea ice plots. Waveform over ocean generally has higher backscattered signal intensity than that over sea ice because of the higher diffuse reflection of ocean. However, the backscatter sigma-0 value of ocean was lower than that of sea ice (Figure 6). This might be because the ocean samples were mostly collected around the Svalbard islands, where strong winds frequently cause high waves, which may reduce the backscattered intensity. Leads showed large variation for all parameters because the size and shape of leads were diverse with different neighboring environments, such as sea ice melting states, combined with the samples being collected in March and April across multiple years, which undoubtedly increased the variation of the parameter values. The range stacked power of single look echoes from leads are similarly high because leads are relatively flat with little waves, which makes SSD of leads low with a narrow distribution in Figure 6. On the other hand, SSDs for sea ice and ocean have a broad distribution due to the surface roughness and waves. The median values of each parameter seem to distinguish lead, sea ice and ocean. However, the distribution of the parameter values of sea ice and leads partly overlapped, possibly due to off-nadir observations of CryoSat-2. This implies that simple thresholding approaches are not suitable to clearly identify leads from sea ice.

4.3. Comparison of Lead Detection Performance

Both See 5.0 and random forest produced similar classification results for the three features.

Table 4. Accuracy assessment results of See 5.0 and random forest by scheme through 10-fold cross-validation. The overall accuracy in percentage averaged for 10 folds is provided.

Scheme	See 5.0	Random Forest
CA (March)	99.50	99.43
CA (April)	93.47	90.40
CM (2011)	92.49	94.80
CM (2012)	94.87	96.96
CM (2013)	95.07	95.48
CM (2014)	94.40	91.60
IC	94.20	94.05

summarizes the overall accuracy by model and scheme through 10-fold cross-validation. All of the cases produced very high overall accuracy (>90%). The most common misclassification for both approaches was between leads and sea ice, possibly due to sampling around the boundaries between them. Since CM and CA resulted in varied accuracy patterns and did not produce significantly higher accuracy than IC, we focused the following discussion on IC. Using IC can reduce temporal variability by including all samples in the subsequent analyses, including sea ice freeboard and thickness estimation.

Table 5. Relative variable importance (i.e., contribution) to lead detection using See 5.0 and random forest by IC.

	SSD	Stack Skewness	Stack Kurtosis	PP	Sigma-0
See 5.0 Usage (%)	0	100	21	40	100
Random forest mean accuracy decrease (%)	19.97	20.44	36.75	20.50	97.72

presents relative variable importance for lead classification by model when using IC. While stack skewness and sigma-0 were used at every node in See 5.0, SSD was not used at all. For the random forest analysis, sigma-0 was identified as the predominant contributing variable, followed by stack kurtosis, PP and stack skewness. Similar to See 5.0, SSD was the least contributing variable to lead detection in random forest. Stack skewness was useful because it was able to distinguish ocean from leads and sea ice with very low error. Interestingly, backscatter sigma-0 was considered a critical parameter for lead detection in both See 5.0 and random forest, but it has not been used in previous studies for lead detection. Backscatter sigma-0 represents not only surface roughness, but also dielectric properties, radar frequency, incidence angle and geometric shape, while the other parameters are mainly sensitive to surface roughness [39].

Table 6. An example of threshold-based rules produced by See 5.0 using IC to classify leads, sea ice and ocean.

	SSD	Skewness	Kurtosis	PP	Sigma-0
Lead		>0.73	>17.53		>27.8
	≤25.6	>0.73	≤17.53		>27.8
Sea ice		≤0.73			≤14.89
		≤0.73		>0.043	14.89$<$ Sigma-0 $<16.48$
		>0.73			≤27.8
	≤25.6	>0.73	≤17.53		27.8 $<$ Sigma-0 $<$ 31.47
Ocean		≤0.73		≤0.043	>14.89
		$\le 0.73$		$>0.043$	$>14.89$
	$>\mathbf{25.6}$	$>\mathbf{0.73}$	$\le \mathbf{17.53}$		$>\mathbf{27.8}$

summarizes threshold-based rules produced using See 5.0 by IC. Previous studies have used SSD, for example, Laxon et al. [18] and Ricker et al. [19] used SSD <4 as one of the conditions to detect leads. However, SSD was not primarily used in the threshold-based rules in this study. It should be noted that the rules were an integration of multiple factors, which implies that the simple thresholding approaches might over- or under-estimate leads resulting in uncertainty in sea ice thickness estimation.

Figure 7 shows two examples of identifying leads using four different lead detection methods. A simple thresholding approach based on PP and SSD (i.e., PP > 0.25 and SSD < 4 for leads and PP < 0.45 and SSD > 4 for ice floes) used in Rose [16] was adopted to identify leads (Figure 7a,e). The simple thresholding method resulted in somewhat over-identification of leads; some leads were mistakenly found on the ice. Laxon et al. [18] also used a similar thresholding approach based on PP [69] and SSD (leads: PP > 18 and SSD < 4; ice floes: PP < 9 and SSD > 4), which also resulted in overestimation of leads on the ice (Figure 7b,f). Although PP and SSD are considered useful parameters for lead detection, simple thresholding based on just two parameters appears insufficient for effectively distinguishing leads from ice. Since surface height on leads is considered as LSSH, if leads are identified on the ice, then LSSH would be overestimated, which would result in an increased bias towards smaller freeboard and thinner sea ice estimates. On the other hand, the two machine learning approaches applied—See 5.0 decision trees and random forest—resulted in improved lead detection and less overestimation of leads compared with the existing approaches (Figure 7c,g,d,h, respectively). Different lead detection approaches were quantitatively assessed and compared (

Table 7. The error matrix based on the See 5.0-based lead detection results for IC.

Reference Classified as	Lead	Sea Ice	Sum	User’s Accuracy (%)
Lead	36	5	41	87.8
Sea ice	6	192	198	96.7
Sum	42	197	239
Producer’s accuracy (%)	85.7	97.5
Overall accuracy (%)	95.4
Kappa coefficient (%)	84

,

Table 8. The error matrix based on the random forest-based lead detection results for IC.

Reference Classified as	Lead	Sea Ice	Sum	User’s Accuracy (%)
Lead	36	2	38	94.7
Sea ice	6	195	201	97.0
Sum	42	197	239
Producer’s accuracy (%)	85.7	98.9
Overall accuracy (%)	96.2
Kappa coefficient (%)	86.4

,

Table 9. The error matrix based on the lead detection results by the approach of Rose [16].

Reference Classified as	Lead	Sea Ice	Sum	User’s Accuracy (%)
Lead	36	28	64	56.2
Sea ice	6	169	175	96.7
Sum	42	197	239
Producer’s accuracy (%)	85.7	85.8
Overall accuracy (%)	85.7
Kappa coefficient (%)	59.3

and

Table 10. The error matrix based on the lead detection results by the approach of Laxon et al. [18].

Reference Classified as	Lead	Sea ice	Sum	User’s Accuracy (%)
Lead	41	45	86	47.7
Sea ice	1	152	152	99.3
Sum	42	197	239
Producer’s accuracy (%)	97.6	77.2
Overall accuracy (%)	80.7
Kappa coefficient (%)	53

) using error matrices. Since the approaches by Rose [16] and Laxon et al. [18] considered ice floes and leads only, i.e., excluding ocean, the accuracy assessment was conducted without ocean samples for consistent comparison between the proposed approaches and the existing literature. The machine learning approaches to lead detection resulted in higher overall accuracy and Kappa coefficients than the approaches used by Rose [16] and Laxon et al. [18]. Both See 5.0 and random forest produced high producer’s accuracy for leads and sea ice. However, the user’s accuracy for leads, as well as overall accuracy and kappa coefficient of random forest were slightly higher than those of See 5.0. Random forest uses an ensemble approach based on numerous independent trees through randomization, which can avoid problems associated with sampling biases. On the other hand, See 5.0 uses only a single tree, but provides more straightforward rules to understand the results at the cost of possible overfitting and sampling biases. Based on the lead detection results in this study, both See 5.0 and random forest can be used to identify leads with minimal difference in the performance. However, in order to analyze the physical meaning among the parameters for lead detection, See 5.0 would be better, as it provides rulesets in simple forms, compared to the ensemble results of random forest. While all four lead detection methods have high producer’s accuracy for both leads and sea ice, the existing approaches (i.e., Rose [16] and Laxon et al. [18]) produced much lower user’s accuracy for leads than the proposed methods, implying the overestimation of leads.

Figure 8 shows the comparison of the SSHA and freeboard from the different lead detection methods. SSHA is a subtraction of LSSH from the mean SSH, representing the relative vertical location of leads. The proposed machine learning-based lead detection methods (Figure 8a,b) detected fewer leads than the existing methods (Figure 8c,d) with few leads above the latitude of 87°N, where leads are rarely found in April. Almost all of the leads detected by the proposed See 5.0 and random forest approaches were also detected by the approaches from Laxon et al. [18] and Rose [16]. SSHA was linearly interpolated and smoothed using a 3 × 3 (pixel) moving average filter. The freeboard, a derivation of surface height from the sum of SSHA and mean SSH, was smoothed by a 30 × 30 (pixel) moving average filter to remove signal noise. While the overall shape of the freeboard lines with latitudes looks similar, the average freeboards by approach—See 5.0, random forest, Rose [16] and Laxon et al. [18]—were 0.095 m, 0.092 m, 0.089 m and 0.090 m, respectively. The average freeboards of Rose [16] and Laxon et al. [18] were relatively underestimated compared to the freeboards by See 5.0 and random forest because of their over-identification of leads on the ice, especially over higher latitudes (>85°N; Figure 8c,d). Laxon et al. [18] found lower SSHA (Figure 8h) between 76°–80°N that appeared to be ocean. However, since the proposed machine learning-based lead detection approaches discriminate ocean from sea ice and leads, ocean was excluded in the SSHA.

4.4. Spatial Distribution of Arctic Sea Ice Freeboard and Thickness

Figure 9 and Figure 10 present the IC ice freeboard and thickness maps, respectively, for March and April from 2011–2014 using a polar stereographic projection with a 25 × 25 km² grid. A typical MYI zone near the Canadian Archipelago and northwestern Greenland has relatively thick freeboard and thickness, except for 2012. The retrieval of freeboard and thickness in these regions appears to show relatively poor LSSHA due to the very limited leads in the regions. On the other hand, sea ice freeboard and thickness around the Kara Sea and Laptev Sea were consistently stable and low for all cases during March and April from 2011–2014. Interestingly, unlike other years, sea ice thickness was relatively high in the central Arctic in 2012, while it was generally low in the typical MYI zone. The annual variability of sea ice thickness was high on the MYI zones, compared to the FYI zones. The amount of sea ice freeboard and thickness apparently diminished from 2011–2013.

Laxon et al. [18] determined Arctic sea ice thickness from February–March 2012 using CryoSat-2 data. Although it was averaged for two months, the overall distribution of sea ice thickness over MYI zones was similar to the results of this study. Sea ice freeboard and thickness maps for March 2013 from this study were slightly different from the results in Ricker et al. [19]. This is possibly because the two studies used different smoothing approaches to waveform data, lead detection methods and gridding approaches to CryoSat-2 track data. Farrell et al. [23] showed two-month averaged sea ice freeboard maps from 2003–2008 using ICESat data. However, sea ice freeboard from ICESat will be different from the sea ice freeboard determined in this study because the height of the sea ice freeboard derived by laser altimetry (i.e., from ICESat) includes snow depth on the sea ice. Farrell et al. [23] observed a slightly thicker sea ice freeboard between 2003 and 2008 (up to 0.75 m) than the 2011–2014 period in the present study. The major differences were found in the Canadian Archipelago, Northern Greenland and the central Arctic, where the freeboard was observed as being high in Farrell et al. [23], while it decreased from 2011–2013 in this study.

4.5. Comparison with AEM-Bird Data

The monthly ice thickness that was averaged by grid using the novel machine learning-based lead detection approaches was compared to the averaged AEM-bird data collected 15 and 17 April 2011 during the CryoVex campaign, as well as the thickness derived by the lead detection methods from Rose [16] and Laxon et al. [18]. This comparison considered the three schemes and two lead detection methods (Figure 11). The sea ice thickness determined using See 5.0 with CA produced the best validation performance on both days with r ~ 0.83 and root mean square error (RMSE) ~0.29 m (Figure 11c).

While Rose [16] and Laxon et al. [18] produced similar performance, the sea ice thickness from See 5.0 generally showed better performance than that derived using the existing methods. However, there are several uncertainty factors for sea ice thickness estimation using CryoSat-2 measurements [19], including: (1) the error of range measurements from CryoSat-2; (2) the uncertainty of detection of leads, resulting in over- or under-estimation of SSHA; (3) the uncertainty of mean scattering horizon in the snow cover; and (4) the uncertainty of snow depth and the density of snow, ice and sea water. Zygmuntowska et al. [70] estimated the uncertainty of Arctic sea ice thickness and volume in terms of sea ice density and snow depth by using the Monte Carlo approach. They revealed that using snow loading (i.e., W99) produced higher uncertainty with respect to the estimation of sea ice thickness than using mean density. In the present study, we used the densities from Alexandrov et al. [33] that do not have year-to-year variability. In order to more accurately estimate sea ice thickness, such changes should be carefully considered, which requires further examination. Ricker et al. [19] analyzed random and systematic uncertainties of Arctic sea ice thickness from CryoSat-2 using the partial derivative of Equation (1) based on the assumption of hydrostatic equilibrium. They showed that random uncertainty affects the estimation of sea ice thickness less than systematic uncertainty caused by the selection of a retracker threshold and the unknown penetration level of the signals on snow. To remove systematic uncertainty caused by the choice of a retracker threshold, Kurtz et al. [17] used a waveform fitting approach to retrieve sea ice freeboard. Any of the above-mentioned factors could result in uncertainty in this study. The lead detection models proposed in this study produced higher accuracy than the existing approaches for lead detection, which implies a possible reduction of the uncertainty caused by the second factor.

In order to examine the influence of snow penetration on the thickness estimation, we conducted a simple sensitivity analysis on snow penetration by testing different penetration ratios with the assumption that radar signals penetrate into the snow depth with the same rate over the entire Arctic region. The results showed that higher accuracy (i.e., lower RMSE) was achieved with increasing penetration depth ratios. Nevertheless, it is difficult to quantify how many centimeters of snow the K_u-band penetrates simply because the snow penetration depth highly depends on the spatiotemporal distribution of snow and whether it is dry or wet. In order to further enhance the sea ice freeboard and thickness produced in this study, snow penetration depth should be considered.

5. Conclusions

In this study, a novel machine learning-based lead detection approach was proposed to quantify Arctic sea ice freeboard and thickness from CryoSat-2 data. The estimated sea ice thickness was validated with AEM-bird data. Accurate lead detection is crucial in estimating LSSH, which is essential to retrieve the freeboard and thickness [16,23]. The results showed that the proposed lead detection approach successfully estimated the sea ice thickness, compared to the existing methods. The overall accuracies by the proposed lead detection methods—decision trees (See 5.0) and random forest—were 95.4% and 96.2%, respectively, which were higher than those produced using the existing methods.

A total of five parameters were used to detect leads, including SSD, stack skewness, PP, stack kurtosis and backscatter sigma-0. Among the parameters, backscatter sigma-0, which prior methods had not considered, played a significant role in determining the threshold-based rules to distinguish leads from ice floes. The lead detection models developed by year or month (i.e., CM and CA) did not produce better performance than the combined model that used all samples for March and April from 2011–2014. This suggests that sea ice thickness in other months such as May or June could be retrieved when additional reference samples from the months were combined with the existing data. That way, a standard lead detection model can be proposed, which can be applied for any year and month. The results also showed that Arctic sea ice freeboard and thickness consistently decreased from 2011–2013, especially in the Canadian Archipelago region, but rebounded in 2014. Future research includes developing a machine learning-based lead detection model that can be applied to any year and month and modeling snow depth penetration using CryoSat-2 data.