Spatiotemporal Analysis of Sea Ice Leads in the Arctic Ocean Retrieved from IceBridge Laxon Line Data 2012–2018

Abstract

The ocean and atmosphere exert stresses on sea ice that create elongated cracks and leads which dominate the vertical exchange of energy, especially in cold seasons, despite covering only a small fraction of the surface. Motivated by the need of a spatiotemporal analysis of sea ice lead distribution, a practical workflow was developed to classify the high spatial resolution aerial images DMS (Digital Mapping System) along the Laxon Line in the NASA IceBridge Mission. Four sea ice types (thick ice, thin ice, open water, and shadow) were identified, and relevant sea ice lead parameters were derived for the period of 2012–2018. The spatiotemporal variations of lead fraction along the Laxon Line were verified by ATM (Airborne Topographic Mapper) surface height data and correlated with coarse spatial resolution sea ice motion, air temperature, and wind data through multiple regression models. We found that the freeboard data derived from sea ice leads were compatible with other products. The temperature and ice motion vorticity were the leading factors of the formation of sea ice leads, followed by wind vorticity and kinetic moments of ice motion.

1. Introduction

Arctic sea ice functions as a sensitive indicator of global warming because sea ice responds to even a small increase in temperature [1,2,3]. On the other hand, Arctic sea ice is also an important driver of climate change, and it plays an important role in the Earth’s solar radiation budget. This is due to how sea ice has a significantly higher albedo compared to that of the water surface. Therefore, when the Arctic sea ice starts to melt, the oceans absorb more solar radiation and warm up, accelerating the melting of sea ice in a positive feedback [4].

Among all types of sea ice features, leads have unique characteristics. A lead is an elongated crack in the sea ice developed by the divergence or shear of floating ice floes when moving with currents and winds [5]. Leads vary in width from meters to hundreds of meters depending on their development and the directions of surrounding pressure and tension. Since a lead is physically an open water body, thin ice, or mixed open water and thin ice within (thicker) sea ice floe or between sea ice floes, it allows the direct interaction between the atmosphere and the ocean and is the only (or major) channel in the cold Arctic. Thus, leads play an important role in the local radiation energy budget, ship navigation, and the Arctic sea ice ecosystem [6]. In particular, they dominate the vertical exchange of energy during winter when turbulent heat fluxes over leads can be orders of magnitude larger than that over thick ice. The width of leads and their orientation markedly influence associated vertical sensible and latent heat fluxes and associated cloud formation [7,8]. Recent studies suggest that these fluxes could influence the atmospheric properties tens to hundreds of kilometers downstream [9,10,11]. Even a small fraction of thin ice and open water within the sea ice pack can significantly modify the total energy transfer between the ocean and the atmosphere [12]. Furthermore, leads are elusive and inconsistent features. If sea water temperature drops below around −1.8 °C, the open water within a lead quickly refreezes (in a few hours), and leads will be partly or entirely covered by a thin layer of new ice [13,14,15]. Therefore, leads are an important component of the Arctic surface energy budget, and more quantitative studies are needed to explore and model their impact on the Arctic climate system.

Arctic climate models require a detailed spatial distribution of leads to simulate interactions between the ocean and the atmosphere. Remote sensing techniques can be used to extract sea ice physical features and parameters and calibrate or validate climate models [16]. However, most of the sea ice leads studies focus on low-moderate resolution (~1 km) imagery such as Moderate Resolution Imaging Spectroradiometer (MODIS) or Advanced Very High-Resolution Radiometer (AVHRR) [17,18,19,20], which cannot detect small leads, such as those smaller than 100 m. On the other hand, high spatial resolution (HSR) images such as aerial photos are discrete and heterogeneous in space and time, i.e., images usually cover only a small and discontinuous area with time intervals between images varying from a few seconds to several months [21,22]. Therefore, it is difficult to weave these small pieces into a coherent large-scale picture, which is important for coupled sea ice and climate modeling and verification. Onana et al. used operational IceBridge airborne visible DMS (Digital Mapping System) imagery and laser altimetry measurements to detect sea ice leads and classify open water, thin ice (new ice, grease ice, frazil ice, and nilas), and gray ice [23]. Miao et al. utilized an object-based image classification scheme to classify water, ice/snow, melt ponds, and shadow [24]. However, the workflow used in Miao et al. was based on some independent proprietary software, which is not suitable for batch processing in an operational environment. In contrast, Wright and Polashenski developed an Open Source Sea Ice Processing (OSSP) package for detecting sea ice surface features in high-resolution optical imagery [25,26]. Based on the OSSP package, Wright et al. investigated the behavior of meltwater on first-year and multiyear ice during summer melting seasons [26]. Following this approach, Sha et al. further improved and integrated the OSSP modules into an on-demand service in cloud computing-based infrastructure for operational usage [22].

Following the previous studies, this paper focuses on the spatiotemporal analysis of sea ice lead distribution through NASA’s Operation IceBridge images, which used a systematic sampling scheme to collect high spatial resolution DMS aerial photos along critical flight lines in the Arctic. A practical workflow was developed to classify the DMS images along the Laxon Line into four classes, i.e., thick ice, thin ice, water, and shadow, and to extract sea ice lead and thin ice during the missions 2012–2018. Finally, the spatiotemporal variations of lead fraction along the Laxon Line were verified by ATM surface height data (freeboard), and correlated with sea ice motion, air temperature, and wind data. The paper is organized as follows: Section 2 provides a background description of DMS imagery, the Laxon Line collection, and auxiliary sea ice data. Section 3 describes the methodology and workflow. Section 4 presents and discusses the spatiotemporal variations of leads. The summary and conclusions are provided in Section 5.

2. Dataset

2.1. IceBridge DMS Images and Study Area

This study uses IceBridge DMS images to detect Arctic sea ice leads along the Laxon Line one day over the course of 7 years in 2012–2018, since these are the longest continuous yearly data available in this Arctic region. The DMS images were collected during the IceBridge sea ice flights using an airborne digital camera. DMS has a high spatial resolution 0.1–2.5 m [27], depending on the aircraft flight height. It has three natural color (red, green, and blue) bands, and each image has a field of view of approximately 400 m by 600 m. The IceBridge campaigns had been designed to survey the Arctic region in March and April since 2009 to partially fill the temporal gap between the ICESat (2003–2009) and ICESat-2 (2018–present) missions.

DMS images are collected, processed, and maintained by the Airborne Sensor Facility located at the NASA AMES Research Center. We downloaded the Level 1B geolocated and orthorectified images for the Arctic Laxon Line in spring from 2012 to 2018 from the NASA National Snow and Ice Data Center Distributed Active Archive Center (NSIDC DAAC) (https://nsidc.org/data/iodms1b) (accessed on 6 August 2021). The Laxon Line starts from the Thule Air Base, Greenland to Fairbanks, AK, USA, transiting across the Arctic Ocean (Figure 1). It passes through both multiyear ice (MYI) regions in the north of the Canadian Archipelago and the first-year ice (FYI) regions in northern Alaska. Thus, sea ice data along this line provides useful insights on the transition of sea ice conditions over the Central Arctic in the spring. Furthermore, the IceBridge mission collected data along this track repeatedly every year from 2012 to 2018, which is appropriate for spatiotemporal analysis of sea ice leads. The overall DMS image collection along the Laxon Line is 106,674 aerial photos (1.54 TB) with an overlap of 60–90% along the track. The photo distribution from 2012–2018 is summarized in

Table 1. The DMS images selected for lead detection along the Laxon Line from 2012 to 2018.

Name	Date	Image #	# Image with Sea Ice Leads	Selected/Original Image Size (GB)	Lighting Condition
Flight 12-426-04	14 March 2012	16,544	1066	14.8/260	Cloudy
Flight 13-426-05	21 March 2013	18,480	993	13.8/290	Normal
Flight 14-426-14	14 March 2014	14,322	492	5.2/150	Cloudy
Flight 15-439-08	26 March 2015	20,038	816	9.3/250	Normal
Flight 16-043-08	20 April 2016	15,205	1069	18.4/270	Normal
Flight 17-426-05	10 March 2017	10,939	659	8.67/93	Cloudy
Flight 18-426-38	6 April 2018	11,146	1040	22.2/240	Normal

. The overall distance of the Laxon Line is around 3398 km, and the distance for the overlapped track through the years is around 2437 km.

2.2. Auxiliary Sea Ice Data

2.2.1. AMSR Data

AMSR (Advanced Microwave Scanning Radiometer) is a passive microwave satellite sensor developed by the Japan Aerospace Exploration Agency. Due to its low spatial resolution, the AMSR data can only be used to examine sea ice concentrations at the regional scale. We collected the AMSR-E/AMSR-2 Unified Level 3 daily brightness temperature and sea ice concentration data which has a spatial resolution of 25 km through NSIDC (

Table 2. Auxiliary sea ice datasets.

Product Name	Type	Source	Spatial Resolution	Category
AMSR-E/AMSR2 Unified L3 Daily Brightness Temperatures & Sea Ice Concentration	Passive microwave	NSIDC	25 km	Sea Ice
IceBridge Airborne Topographic Mapper (ATM)	Laser altimeter	NSIDC	~1 m footprint (resampled to 2 m grid)	Sea Ice
Global sea ice type	Sea ice type	EUMETSAT OSI SAF	10 km	Sea Ice
Polar Pathfinder Daily EASE-Grid Sea Ice Motion Vectors	Sea ice motion	NSIDC	25 km	Dynamic
ERA5 (air temperature and wind velocity)	Climate reanalysis	European Centre for Medium-Range Weather Forecasts (ECMWF)	0.25°	Dynamic and thermodynamic

) [28]. The data contain vertically polarized and horizontally polarized brightness temperatures at four frequency channels: 18.7, 23.8, 36.5, and 89.0 GHz. The Arctic sea ice concentration (SIC) was calculated by the NASA Team 2 (NT2) algorithm, which provides <2% of error compared with the high-resolution optical data [29,30,31]. The collected AMSR data coincides with the days of the IceBridge mission from 2012 to 2018, so that the SIC can be compared with that retrieved from the DMS images. Furthermore, the passive microwave data can be used to calculate thin ice concentration (TIC). Röhrs and Kaleschke used brightness temperatures at the vertically polarized 18.7 and 89.0 GHz to identify water and thin ice (i.e., new ice, nilas, and pancake ice) from thick ice, and the sea ice leads and TIC showed a good agreement with the MODIS, Envisat ASAR, and CryoSat-2 data [14]. In this study, we calculated TIC following the Röhrs and Kaleschke’s algorithm. The coarser spatial resolution of 25 km of TIC were compared with our lead and thin ice fractions retrieved from the DMS images.

2.2.2. ATM Surface Height Data (DMS Level)

Our DMS-based lead detection results can be used to cross-validate sea ice freeboard products derived from IceBridge Airborne Topographic Mapper (ATM) Level 1B data [23]. The ATM is an airborne conically scanning laser altimeter with a wavelength of 532 nm. A laser pulse is emitted from the ATM at a rate of 5 kHz, and it has ~1 m of footprint at a typical 500 m altitude above the surface. Since ATM covers exactly the same location and time with the DMS images with a smaller cross-track width (~400 m), DMS images are usually used as good reference for extracting the ATM-based freeboard data [32,33]. In this study, the ATM data are resampled in a 2 m grid and projected to the same projection system as DMS (NSIDC sea ice polar stereographic North) to match the geographical location. After retrieving thin ice and leads through DMS images, we geographically linked the leads with the ATM data to extract freeboard variations along the Laxon Line, and compare with freeboard data derived from SILDAMS (Sea Ice Lead Detection Algorithm) [23,32].

2.3. Oceanic and Atmospheric Geophysical Parameters

NSIDC provides sea ice motion data (nsidc.org/data/NSIDC-0116) with a spatial resolution of 25 km on the Equal-Area Scalable Earth grid [34]. This sea ice motion vector is derived from multiple data sources, including AVHRR, AMSR-E, SMMR, SSMI, SSMI/S satellite sensors, International Arctic Buoy Program (IABP) buoys, and the National Center for Environmental Prediction/National Center for Atmospheric Research (NCEP/NCAR) reanalysis.

We also acquired a global sea ice type product provided by the European Organization for the Exploitation of Meteorological Satellites (EUMETSAT) Ocean and Sea Ice Satellite Application Facility (OSI SAF, www.osi-saf.org) (6 August 2021) [35]. This product assigns different sea ice types, such as multiyear ice (MYI), first-year ice (FYI), and open water, from various satellite data. This is a daily product and has 10 km of spatial resolution.

Other data we used included air temperature (2 m above sea level) and wind velocity data coincident with the DMS images acquired from the European Centre for Medium-Range Weather Forecasts ERA-5 reanalysis. The ERA-5 product has 0.25° spatial resolution and consists of hourly variables, and we integrated this hourly data into daily products and resampled them to 25 km resolution to match the ice motion data. This ERA-5 product was downloaded from the Climate Data Store (cds.climate.copernicus.eu) of the Copernicus Climate Change Service.

In this study, the high spatial resolution lead fractions derived from DMS along the Laxon Line were linearly regressed with the coarse spatial resolution sea ice motion, air temperature, and wind velocity products to identify potential significant drivers.

3. Methods

3.1. Batch Classification Processing Workflow

Since the IceBridge DMS images are highly overlapped along the track (60–90%), we selected one image from every three consecutive images along the Laxon Line to reduce the computation burden. All images in continental land masses and poor-quality images due to overwhelming cloud coverage and low lighting conditions were manually removed, finally generating a collection of sea ice lead images (Figure 2).

The object-based classification scheme was designed based on the color and texture of sea ice features on DMS images. Four sea ice classes were defined: (1) thick ice is usually thick ice or snow-covered ice with a high albedo; (2) thin ice is usually fresh and newly formed ice, which has a smooth surface with a low albedo, since solar radiation is partially absorbed by the water beneath it; (3) open water is dark and smooth due to its strong absorbance of solar radiation; and (4) shadow is within a thick-ice area and is a relative dark feature projecting on the ice surface by surrounding ridges or snow dunes. DMS images collected in different years have different lighting conditions, which affects the image quality (Table 1). Furthermore, even in the same year, the quality of images was quite distinctive due to the local cloud coverage and lighting conditions, as shown in Figure 3. For example, three subgroups were identified in 2012 DMS images: normal images contained regular sea ice scenes with an appropriate exposure and contrast, and all sea ice classes were recognizable by color and texture; gray images were partially cloudy images with a poor lighting condition, so they were relatively dark, and shadows were difficult to detect; and poor images were under extremely poor lighting conditions, and the boundaries between water, thick ice, and thin ice were blurred due to low contrast.

Therefore, training samples were selected using a divide-and-conquer strategy based on image quality. All DMS images taken in 2013, 2015, 2016, and 2018 were under good lighting conditions, and training samples were selected for all four sea ice features. However, the images taken for the other three years were processed in different ways. The training samples for all images taken in 2012, 2014, and 2017 were only selected for thin ice, open water, and thick ice, without considering shadow due to low lighting conditions. Furthermore, the 2012 images were manually classified into three subgroups, i.e., normal, medium, and poor. The 2014 images were manually classified into two subgroups, i.e., normal and medium, and all poor images were abandoned due to serious vignetting, caused by light hitting the lens aperture at a large angle, and significantly reduced brightness values on the four corners of the image. The 2017 images were all classified into the medium subgroup only. In summary, the independent training samples were collected for each subgroup and year for supervised classification.

The OSSP package uses an object-based classification scheme. For each image, the watershed segmentation method is used to convert pixels into objects. Therefore, training samples are labelled at the object level. Only distinctive and typical sea ice objects are selected across the whole scene, and each sea ice class has around 120–250 objects. The attributes of objects such as color values, band ratios, textures, and shape indexes are calculated and served as supervised classification features. Based on these training datasets, the OSSP package uses the random forest classification method to label all unknown objects in DMS images [24,25].

To evaluate the accuracy of classification results, the independent test object samples were also collected.

Table 3. The DMS images selected for lead detection in the Laxon Line from 2012 to 2018.

Testing Group	# Training Image	# Training Object	# Test Image	# Test object
DMS2012_normal	6	50	5	114
DMS2012_medium	7	90	5	94
DMS2012_ poor	7	65	5	124
DMS2013	13	196	7	221
DMS2014_normal	8	106	6	178
DMS2014_medium	6	66	6	119
DMS2015	11	150	9	254
DMS2016	8	144	12	444
DMS2017	12	140	6	150
DMS2018	13	135	9	319

lists the selected image and object numbers for the training and testing process of each classification group. Finally, the confusion matrix was generated at the pixel level and was used for calculating the overall accuracy, user’s accuracy, producer’s accuracy, and Kappa coefficient.

3.2. Sea Ice Leads Parameters Definitions

Based on the classified result in each surface type, we derived the sea ice leads by combining thin ice and open water. Then, the sea ice lead fraction, open water fraction, thin ice fraction, and sea ice concentration were calculated on a per-scene basis. The sea ice lead fraction for each DMS image can be calculated using the following equations:

Sea Ice Lead Fraction (SILF):

SILF = (ThinIce + OpenWater)/(ThickIce + ThinIce + OpenWater + Shadow) ∗ 100,

(1)

where ThinIce, OpenWater, ThickIce, and Shadow are pixel numbers of classified thin ice area, open water, thick ice, and shadow for a DMS image, respectively.

3.3. Spatiotemporal Analysis with Auxiliary Sea Ice Data

The auxiliary sea ice datasets can be used to assess the DMS-based lead detection results to deepen the understanding of the formation mechanism of leads. In this research, first, our lead detection result was used to determine local sea reference height and calculate the sea ice freeboard. This retrieved freeboard was compared with the existing NSIDC freeboard data at the scale of 400 m [36]. Furthermore, the coincident AMSR thin ice concentration (TIC) data, and the geophysical atmosphere and ocean data, such as air temperature, wind velocity, and sea ice motion, were compared with the lead fraction results.

Based on our DMS lead detection algorithm, sea ice freeboards were retrieved from the ATM lidar data using the same method as in [32]. Specifically, we removed variations in the instantaneous sea surface height by subtracting geoid and ocean tide height. Then, we calculated the freeboard by subtracting locally determined leads surface height (z_shh) from the corrected height (H_corr).

Freeboard = H_corr − z_shh,

(2)

where z_shh is determined from the sets of individual lead elevation estimates through ordinary kriging. We calculated the mean freeboard for each DMS image (around 400 m by 600 m) and resampled the value to 400 m resolution. On the other hand, Kurtz et al. used an automated lead detection algorithm through the minimal signal transform [23,32] and then retrieved the freeboard at the resolution of 400 m. Therefore, the two products can be compared and cross-verified at this scale.

TIC could be calculated from the AMSR as described in Röhrs and Kaleschke [14] with a rather coarse spatial resolution of 25 km. This AMSR-based TIC represents the existence of open water and thin ice on sea ice leads. This TIC is conceptually equivalent to our SILF. Since the AMSR and DMS have different resolutions and geographical coverage, they cannot be compared directly. Therefore, we resampled and averaged the DMS-based ice lead fractions for every 25 km grid cell to match the spatial resolution of AMSR data, as shown in Figure 4. Then, the mean of sea ice lead fractions within the range of each 25 km block was calculated.

Furthermore, the 25 km resampled lead fractions were also correlated with other 25 km resolution sea ice and atmospheric data including NSIDC sea ice motion, ERA5 air temperature, and wind velocity. Since kinetic moments of sea ice movement can play an important role in formations of leads, four kinetic moments or tensions were calculated from the NSIDC sea ice motion data by the following equations [37]:

$$divergence=\frac{\partial F_x}{\partial x}+\frac{\partial F_y}{\partial y}$$

(3)

$$vorticity=\frac{\partial F_y}{\partial x}-\frac{\partial F_x}{\partial y}$$

(4)

$$shearing deformation=\frac{\partial F_y}{\partial x}+\frac{\partial F_x}{\partial y}$$

(5)

$$stretching deformation=\frac{\partial F_x}{\partial x}-\frac{\partial F_y}{\partial y}$$

(6)

where $F_x$ and $F_y$ refer to the velocity of sea ice along the x and y axes, respectively. Divergence is a measure of parcel area change without the change of orientation or shape, and vorticity is a measure of orientation change without area or shape change. Shearing and stretching deformation are measures of shape change produced by differential motions parallel and normal to the boundary, respectively [37].

Finally, based on the assumption that these atmosphere and sea ice variables for a series of the previous days would contribute to the formation of sea ice leads, the average of these dynamic and thermodynamic variables up to 30 successive days before the DMS acquisition day were calculated (

Table 4. Variables for the multiple linear regression models.

Department	Factors	Description
Sea Ice Leads	mean_leads	Mean lead fraction for 25 km segment
Temperature	tmpXX	Averaged air temperature for last XX days (e.g., tmp03 means average temperature of last 1, 2, 3 days)
Wind	U10_XX	Averaged u-component of wind velocity for last XX days
	V10_XX	Averaged v-component of wind velocity for last XX days
	wind_XX	Averaged wind velocity for last XX days (e.g., wind_10 means wind velocity for last 10 days)
Sea Ice Motion	u_ice_XX	Averaged u-component of ice velocity for last XX days (e.g., u_ice_10 means u-velocity for last 10 days)
	v_ice_XX	Averaged v-component of ice velocity for last XX days (e.g., v_ice_10 means v-velocity for last 10 days)
	vel_ice_XX	Averaged ice velocity for last XX days (e.g., v_ice_10 means ice velocity for last 10 days)
	divXX	Averaged divergence of sea ice motion for last XX days (e.g., div10 means divergence for last 10 days)
	vorXX	Averaged vorticity of sea ice motion for last XX days (e.g., vor10 means vorticity for last 10 days)
	shrXX	Averaged shearing deformation of sea ice motion for last XX days (e.g., shr10 means shearing deformation for last 10 days)
	stcXX	Averaged stretching deformation of sea ice motion for last XX days (e.g., stc10 means stretching deformation for last 10 days)

). By comparing these variables and the lead fractions, we hoped to identify the potential contribution of these explanatory variables to lead formation.

Multiple linear regression (MLR) was used for modelling the mean lead fractions in terms of large-scale sea ice dynamic–thermodynamic variables, including the NSIDC sea ice motion data with four kinetic moments, ERA-5 air temperature, and wind velocity data. The forward and backward stepwise regression methods were used to identify the most important explanatory variables. This strategy refers to the process of building a regression model by adding or removing explanatory variables in a stepwise manner until the predicted variable does not change significantly [38].

4. Result and Discussion

4.1. Classification Result

A total of 106,674 DSM images along the Laxon Line from 2012–2018 were processed, and a total of 6135 images with sea ice leads were visually selected (Table 1). All images were classified through the OSSP package integrated in the ArcCI online service [22].

Six classified images in 2012 are shown in

Table 5. Comparison of original 2012 DMS images and classified results for three subgroups. Two samples were selected for each subgroup.

	Sample Result 1		Sample Result 2
	Raw Image	Classified Result	Raw Image	Classified Result
Normal
Medium
Poor

. The first row shows the classification results for the subgroup of normal images, the second row for the medium images, and the third row for the poor images. All six images were selected to show a variety of sea ice features under different lighting conditions. The classified results illustrate four sea ice classes: open water, shadow, thin ice, and thick ice.

The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in

Table 6. Pixel-level classification accuracy for each production group. All values except Kappa coefficient are in percentages.

Testing Group	Overall Accuracy	Kappa Coef.	UA_Thick **	UA_Thin	UA_Shadow	UA_Water	PA_Thick **	PA_Thin	PA_Shadow	PA_Water
DMS2012_normal	88.9	0.83	88.0	91.7	83.8	nan *	98.4	94.2	63.8	nan
DMS2012_medium	93.6	0.85	97.3	85.0	nan	95.5	93.8	93.1	nan	97.5
DMS2012_poor	93.8	0.86	95.0	96.0	nan	61.9	98.9	81.2	nan	94.9
DMS2013	96.4	0.95	92.2	100.0	99.4	95.5	99.7	96.5	88.3	99.9
DMS2014_normal	88.0	0.82	74.7	86.2	93.9	98.0	97.1	81.3	99.7	89.0
DMS2014_medium	93.7	0.89	91.7	96.3	nan	97.1	100.0	75.7	nan	97.1
DMS2015	86.4	0.78	86.6	83.5	98.6	93.4	99.8	80.9	82.2	57.9
DMS2016	87.9	0.83	82.1	89.3	95.0	95.7	99.4	68.8	89.7	90.2
DMS2017	86.7	0.75	87.4	82.8	nan	99.4	97.6	76.5	nan	60.7
DMS2018	93.5	0.88	91.9	96.5	95.2	97.9	98.5	79.1	89.4	98.4
Average Accuracy	90.9	0.84	88.7	90.7	94.3	92.7	98.3	82.7	85.5	87.3

*,** User’s accuracy and producer’s accuracy for each classified ice type represented as UA_XX, and PA_XX, and XX could be thick ice, thin ice, shadow, or open water.

. The overall accuracy across the 10 test samples selected by year and illumination conditions was 90.9 ± 3.5%, where the latter number is the standard deviation, and the Kappa coefficient was 0.85 ± 0.05. Since sea ice leads were defined as a combination of thin ice and open water, classification accuracy was determined by these two classes. The user’s accuracy for thin ice and water were 90.7 ± 5.9% and 92.7 ± 11.0%, respectively. The low accuracy of 61.9% for open water in the 2012 poor subgroup was due to the confusion between water and thin ice under extremely poor lighting conditions.

4.2. Overall Integrated Statistical Analysis and Trend of Sea Ice Leads and Freeboard

4.2.1. Sea Ice Leads Fraction, Area, and Frequency

Figure 5a shows the averaged lead fraction for every 25 km along the Laxon Line. Relatively large lead fractions (>15%) were only observed near the Beaufort Sea area (track distance > 1200 km) in 2013, 2014, and 2016, where they were generally located in the FYI region or transition region between FYI and MYI. However, the smaller lead fraction region in the central Arctic (track distance < 1200 km) was primarily covered by MYI and thick ice. Although these observations of one day per year for seven years cannot represent the overall continuous spatiotemporal variations of lead fraction, this general spatial pattern agrees with that of previous lead studies [5,18,19,39]. Figure 5b portrays the averaged area of individual leads for the 25 km track segment, and Figure 5c portrays the ratio of the number of lead-included images to the total number of images for the 25 km segment. The lead fraction (Figure 5a) was determined by the individual lead area (Figure 5b) and the frequency of leads (Figure 5c). For example, although large leads were observed in 2013 for 0–500 km (Figure 5b), lead frequency for this part was low (Figure 5c) due to the small number of large leads. As a result, the averaged lead fraction for this segment was not high because of the low lead frequency. In addition, the lead frequency in 2018 for 1000–2500 km was relatively high, but the averaged lead fraction was not so high due to the large number of small leads.

4.2.2. Retrieval of Freeboard

Based on the DMS lead detection result, we calculated the 400 m mean sea ice freeboard from the ATM surface height data (Figure 6). The MYI area (near central Arctic Ocean) at track distance <1200 km showed a higher freeboard (i.e., thicker ice) compared to that of the FYI area (near the Beaufort Sea with a track distance beyond 1200 km). As shown in

Table 7. ATM sea ice freeboard retrieved from the DMS lead detection.

Year	FYI	MYI	Total
2013	0.263	0.519	0.409
2014	0.277	0.339	0.320
2015	0.275	0.470	0.407
2016	0.335	0.398	0.354
2017	0.211	0.467	0.366
2018	0.320	0.505	0.414

, the FYI area always showed a lower freeboard than the MYI area. In addition, the freeboard retrieved from our lead detection shows a good correlation with the ATM freeboard product provided by NSIDC [32]—correlation coefficient (R) was 0.832, and root mean square difference (RMSD) was 0.105 m (

Table 8. R and RMSE between our freeboard estimation and NSIDC freeboard estimation.

Year	R	RMSD (m)
2013	0.928	0.089
2014	0.907	0.063
2015	0.755	0.140
2016	0.784	0.114
2017	0.742	0.119
2018	0.869	0.082
Total	0.832	0.105

). It is also noted that 2015, 2016, and 2017 showed relatively lower R and higher root mean square error (RMSE) than the other years (Table 8 and Figure 7), which might be due to the lower classification accuracy of these years (Table 6). Some misclassified leads can make substantial differences in estimation of sea surface height, eventually leading to the differences between our freeboard estimation and the NSIDC freeboard products. Nevertheless, the freeboard differences between MYI and FYI and the cross-validation with the NSIDC freeboard product showed that our lead detection result was reasonable and compatible with other lead detection products.

4.3. Sea Ice Lead Fraction Modelling with Auxiliary Sea Ice Product

In general, March and April have the lowest lead fraction and lead frequency in a year because of the highly packed sea ice conditions [5,23]. Since the OIB missions were conducted during these months of packed sea ice, the widths of individual observed leads were usually less than 1 km. Indeed, as shown in Figure 5b, most leads had less than 0.1 km² of area, which accounts for a tiny portion of the entire 25 × 25 km grid cells. Hence, it is reasonable that the DMS-based lead detection and AMSR-based TIC were not highly correlated (R~0.21, Figure 8), because narrow leads are hardly detected by the coarse resolution satellite data [14,40]. For example, we found that most of AMSR-based TIC along the track was zero and AMSR-based SIC was 100% even though the DMS images clearly showed leads in that area.

Figure 9 shows the lead fractions and related dynamic and thermodynamic variables at the scale of 25 km on the same days that DMS images were taken from 2012 to 2018. In general, the lead fractions did not show significant correlation with any single auxiliary variable or kinetic property from sea ice motion data. This is reasonable because (1) these ancillary data have 25 km spatial resolution, which is much coarser than the spatial resolution of the DMS image; (2) the DMS images have only ~500 m of width, representing only a small portion along the Laxon Line; and (3) the formation of sea ice leads results from the accumulative and complex effects of multiple dynamic and thermodynamic variables, rather than just one variable.

Although the DMS images have different spatial scale with the ancillary datasets, we attempted to explore the potential relationship the DMS-based lead fractions and sea ice dynamic and thermodynamic variables from the ancillary datasets. Assuming that (1) these variables are the results of the large-scale atmosphere and ocean circulation and (2) the combination of these variables somehow affects the formation of leads, we normalized all explanatory variables and constructed a series of multiple-variables linear regression models, as shown in Equation (7).

$$\mathrm{SILF}={\displaystyle \sum }_{k=0}^n{a}_k{x}_k$$

(7)

where $x_k$ is one of the normalized dynamics-thermodynamic variables, and $a_k$ are corresponding coefficients.

The lead fraction variable is the mean of all DMS image-based lead fractions within a 25 km block. On the other hand, all dynamic-thermodynamic variables, including four kinetic moments from the NSIDC sea ice motion data, ERA5 air temperature, and wind velocity data, were averaged by 1, 2, 5, 10, 20, and 30 days prior to the date when the DMS image was taken, considering the accumulative effects of these explanatory variables.

After exploring all possible multiple linear regression models, we found that dynamic-thermodynamic variables integrated by 10 days showed the highest correlation coefficient. Therefore, these explanatory variables were used to reconstruct the linear regression models using the forward and backward stepwise regression approach. The coefficients of all normalized explanatory variables for all models are illustrated in

Table 9. Selected variables and coefficients in 14 stepwise linear regressions.

Year	Approach	R2	Tmp10	U10_10	V10 _10	Wind_10	U_Ice_10	V_Ice_10	Vel_Ice_10	Div10	Vor10	Shr10	Stc10	Constant
2012	Forward	0.26	/	/	/	/	−0.39	−0.38	0.16	−0.10	−0.08	/	/	0.41
	Backward	0.26	0.10	/	/	/	−0.34	−0.19	/	−0.12	/	/	/	0.31
2013	Forward	0.48	−1.19	/	/	0.35	−6.46	−2.78	9.51	/	−0.01	−0.14	/	0.60
	Backward	0.48	−1.18	/	/	0.35	−6.44	−2.75	9.45	/	/	−0.15	/	0.08
2014	Forward	0.87	4.61	−5.60	−0.97	1.09	1.24	15.31	−12.98	/	0.89	−0.55	/	−2.08
	Backward	0.87	4.64	−5.37	/	/	1.16	13.34	−11.25	−0.16	0.87	−0.59	/	−1.94
2015	Forward	0.34	/	/	−0.53	/	−1.35	/	1.19	0.15	0.14	0.28	−0.33	0.40
	Backward	0.34	/	/	−0.53	/	−1.35	/	1.19	0.15	0.14	0.28	−0.33	0.40
2016	Forward	0.29	/	−0.79	/	/	/	/	0.29	0.30	−0.39	0.57	0.15	0.21
	Backward	0.34	0.67	−4.62	−0.53	4.09	/	/	/	/	−0.36	0.46	/	0.22
2017	Forward	0.66	−1.17	−6.54	−3.08	6.77	2.98	−0.09	−2.01	−0.19	/	/	/	1.50
	Backward	0.66	−1.15	−6.57	−3.11	6.86	3.02	/	−2.09	−0.19	/	/	/	1.45
2018	Forward	0.30	0.34	−1.40	−1.40	1.83	/	/	/	/	−0.03	/	−0.31	0.45
	Backward	0.30	0.34	−1.31	−1.33	1.72	/	/	/	/	/	/	−0.32	0.42

. There were 11 thermodynamic-dynamic variables, including one thermodynamic variable (temperature), six dynamic variables (velocity of wind and ice motion), and four kinetic moments caused by ice motion.

The forward and backward stepwise regression models for each year identified different sets of explanatory variables. Both 2012 models identified ice motion velocity and divergence as the significant explanatory variables. The 2013 models mainly identified the ice motion velocity and temperature variables. Other than ice motion velocity and temperature, the 2014 models included wind velocity at u-direction, and the correlation coefficient was significantly higher than that of other models. The 2015 models emphasized the functions of wind and ice motion velocity. The 2016 forward model identified more kinetic moments, but the backward model emphasized wind velocity, which represents the possible correlation among these variables. Finally, the 2017 and 2018 models showed significant influence of wind velocity and temperature.

Except for that of 2014, all other models had only moderate correlation, and R² ranged from 0.26 to 0.66. This was because (1) the sea ice fractions were derived from high spatial resolution DMS images, and the dynamic-thermodynamic variables had a much coarser resolution of 25 km; (2) the atmospheric and oceanic dynamics that contribute to lead formation can occur in a much smaller scale (<25 km scale), which cannot be captured by coarse resolution products; and (3) the uncertainty of the DMS-based lead detection (accuracy of 90%) can be carried and exaggerated in the data fusion and resampling process.

Based on all stepwise regression results, the relative explanatory variable importance could be ranked based on their frequencies in a total of 14 regression models (Table 9), as summarized in Figure 10. It showed that temperature and ice motion vorticity were the leading factors of the formation of sea ice leads, followed by wind vorticity and kinetic moments or tensions of ice motion. However, since this result is only based on stepwise regression of several available variables, it cannot clearly explain the detailed mechanism of lead formation that is a complex combination of multiple ocean and atmospheric parameters. In addition, it is noted that the spatial resolution of the variables can be too coarse to represent the formation of leads in the DMS image scale. Therefore, more comprehensive studies are needed to clearly understand small-scale lead formations in the future.

5. Conclusions

This research demonstrates a scientific case study for sea ice lead detection during 2012–2018 along the IceBridge Laxon Line. To address the lack of standard image processing workflow for sea ice parameter extraction from massive and long-term HSR imagery, a practical object-based image classification workflow was implemented based on the OSSP package to extract multiscale multitype sea ice features and to calculate sea ice lead fractions and freeboard parameters. These sea ice products could be directly used to validate other coarse resolution remote sensing images/products. Furthermore, the high-spatial-resolution sea ice fractions were statistically modeled using large scale dynamic-thermodynamic models.

We found that thick ice, thin ice, water, and shadow could be successfully classified using an object-based classification algorithm or the OSSP package with reasonable overall accuracies of 86.4–96.4%. The sea ice lead fractions along the Laxon Line could be calculated for each DMS image accordingly. The temporal and spatial distribution of leads were verified by ATM surface height data and an independent freeboard product. Finally, the lead fractions were aggregated and modelled with 25 km resolution dynamic and thermodynamic variables including sea ice motion, air temperature, and wind data. All stepwise linear regression models had medium to high correlation coefficients. It seems that temperature and ice motion vorticity were the leading factors of the formation of sea ice leads, and each year could have different dominant factors. The results could provide insightful understanding of the mechanism of sea ice leads, which is useful for climate modelling.

In the future, novel image classification algorithms such as deep learning could be used to improve the traditional machine learning methods. The methods can be extended to other sea ice regions and data types. The results and parameters derived from this study can help the sea ice community to better understand the mechanisms driving sea ice variability so that they can be better represented in climate models.