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

: 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 workﬂow 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 identiﬁed, 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 veriﬁed 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.


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.

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 Remote Sens. 2021, 13, 4177 3 of 18 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 overall distance of the Laxon Line is around 3398 km, and the distance for the overlapped track through the years is around 2437 km. 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). 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 overall distance of the Laxon Line is around 3398 km, and the distance for the overlapped track through the years is around 2437 km.   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) [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. 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].

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 Na- 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.

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). from SILDAMS (Sea Ice Lead Detection Algorithm) [23,32].

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) [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.

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. 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 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 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.

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.

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 ).
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.
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]: where and 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 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]: shearing de f ormation = ∂F y ∂x 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). 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]. 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)

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. 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. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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  The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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 Remote Sens. 2021, 13, x FOR PEER REVIEW 10 o The classification accuracies were evaluated at the pixel-level, and all calculated curacies are summarized in Table 6. The overall accuracy across the 10 test samples lected by year and illumination conditions was 90.9 ± 3.5%, where the latter number is standard deviation, and the Kappa coefficient was 0.85 ± 0.05. Since sea ice leads w defined as a combination of thin ice and open water, classification accuracy was de The classification accuracies were evaluated at the pixel-level, and all calculated accuracies are summarized in Table 6. 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. Table 6. Pixel-level classification accuracy for each production group. All values except Kappa coefficient are in percentages.

Testing Group
Overall Accuracy Kappa Coef.  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.

UA_
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.

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

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, 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). 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.
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.    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 Figure 7. Scatter plot between ATM freeboard derived by our lead detection and NSIDC freeboard product for every 400 m (2% random selection of the total points).

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 2 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. 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 2 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).
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 dynamicthermodynamic 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. 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.

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

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 highspatial-resolution sea ice fractions were statistically modeled using large scale dynamicthermodynamic 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.