Journal cover Journal topic
The Cryosphere An interactive open-access journal of the European Geosciences Union
Journal topic
The Cryosphere, 13, 1565-1582, 2019
https://doi.org/10.5194/tc-13-1565-2019
The Cryosphere, 13, 1565-1582, 2019
https://doi.org/10.5194/tc-13-1565-2019

Research article 04 Jun 2019

Research article | 04 Jun 2019

# Estimation of turbulent heat flux over leads using satellite thermal images

Estimation of turbulent heat flux over leads using satellite thermal images
Meng Qu1, Xiaoping Pang1,2, Xi Zhao1,2,3, Jinlun Zhang3, Qing Ji1, and Pei Fan1 Meng Qu et al.
• 1Chinese Antarctic Center of Surveying and Mapping, Wuhan University, Wuhan, China
• 2University Corporation for Polar Research, Beijing, China
• 3Applied Physics Laboratory, Polar Science Center, University of Washington, Seattle, WA, USA
Abstract

1 Introduction

Leads are linear structures of the ocean surface within pack ice that are exposed to the atmosphere during an opening event caused by various forces, such as wind and water stresses. In winter, thin ice forms quickly in newly opened leads due to the large temperature difference between the ocean and the atmosphere (Kwok, 2001). The opening of leads breaks the continuity of insulating ice and creates windows for a stronger air–ocean interaction. Newly opened leads are the main source of ice production, brine rejection, and heat transfer from the ocean to the atmosphere (Alam and Curry, 1998). Turbulent heat flux over open water could be 2 orders of magnitude larger than that through mature ice (Maykut, 1978). Although decreasing rapidly with growing ice thickness, ice growth rates can still be an order of magnitude larger for 50 cm thick young ice than for 3 m thick ice (Maykut, 1986). In the central Arctic, open water usually comprises no more than 1 % of the ice pack area during the winter. However, open water, together with thin ice (<1 m) estimated to be 10 % of the whole ice area, contributes to more than 70 % of the upward heat flux (Maykut, 1978; Marcq and Weiss, 2012). A model study shows that an increased lead fraction by 1 % can lead to local air temperature warming up to 3.5 K in winter (Lüpkes et al., 2008).

Leads also allow more surface absorption of radiation due to their lower albedo compared to thick ice. This will accelerate sea ice thinning in summer and delay refreezing in early winter and therefore decrease the mechanical strength of the ice cover and allow even more fracturing, larger drifting speed and deformation, and faster export of sea ice to lower latitudes (Rampal et al., 2009). As the ice pack gets thinner (Kwok and Rothrock, 2009) and more mobile (Spreen et al., 2011), favorable for deformation and opening, networks of more intensive lead with stronger local influence are expected.

Bulk aerodynamic formulae are frequently used in climate models to generalize the turbulent heat flux from open water in Arctic pack ice (Lindsay and Rothrock, 1994; Walter et al., 1995). The bulk formulae attribute heat flux over leads to wind speed, temperature differences between the surface and the atmosphere, and a turbulent transfer coefficient for heat, which is a function of the stability of the near-surface atmosphere and the roughness of the surface. In this approach, Lindsay and Rothrock (1994) estimated sensible heat flux using surface temperature retrieved from the Advanced Very High Resolution Radiometer (AVHRR), while observations suggest that for small leads, down to dozens of meters in width, the discontinuity between leads and pack ice causes the creation of a thermal internal boundary layer (TIBL) in the bottom atmosphere, reducing turbulent heat exchange on the downwind side (Venkatram, 1977; Andreas et al., 1979). Convective plumes formed above leads may further complicate the process within the TIBL (Tetzlaff et al., 2015).

Models were developed for estimation of TIBL thickness and turbulent heat flux over coastal polynyas, leads, and ice edges (Alam and Curry, 1997; Andreas and Cash, 1999; Renfrew and King, 2000; Chechin and Lüpkes, 2017). Chechin and Lüpkes (2017) modeled boundary layer development downwind of the ice edge, potential temperature, and mix-layer height, and wind speed variation was analyzed as well. Renfrew and King (2000) modeled turbulent heat flux over large fetch (5–50 km wide, typical for coastal polynya) during cold-air outbreaks. The dependence of turbulent heat flux on lead width was estimated in several studies (Andreas and Murphy, 1986; Alam and Curry, 1997; Andreas and Cash, 1999). On the basis of the Monin–Obukhov similarity theory and the surface renewal theory, Alam and Curry (1997) estimated turbulent heat flux over leads using an intricate surface roughness model (Bourassa et al., 2001). Sensible heat flux across a single lead is integrated from fetch 0 to fetch X. Andreas and Murphy (1986) calculated transfer coefficient CN10 at 10 m height for turbulent heat in neutral stability, using the nondimensional fetch $-X/L$, where L is the Obukhov length. A maximum CN10 of $\mathrm{1.8}×{\mathrm{10}}^{-\mathrm{3}}$ was found at small fetch, and the value decrease to $\mathrm{1.0}×{\mathrm{10}}^{-\mathrm{3}}$ with increasing $-X/L$. Andreas and Cash (1999) computed lead-average turbulent heat flux using transfer coefficient C as a function of stability parameter $-h/L$, where h is the fetch-dependent height of the TIBL. For small fetch ($-h/L<\mathrm{6}$), turbulent heat is exchanged by mixed free and forced convection, resulting in a large C and higher heat flux.

A power law distribution of lead widths was also reported in various studies (Wadhams, 1981; Wadhams et al., 1985; Lindsay and Rothrock, 1995), indicating that small leads prevail in the Arctic. Impacts of lead width on heat flux were reported by Maslanik and Key (1995) and Marcq and Weiss (2012) using different width distribution models. However, fetch-limited models have not been applied to surface temperature fields retrieved from remote sensing imagery to estimate turbulent heat flux at regional scale, due to the coarse resolution of operational thermal sensors. Fortunately, the launch of Landsat-8 in February 2013 has provided a unique opportunity for the estimation of turbulent heat flux with finer-resolution temperature fields.

In this paper, we derive lead distribution using thermal images from two sensors, Moderate Resolution Imaging Spectroradiometer (MODIS) and Thermal Infrared Sensor (TIRS) aboard Terra and Landsat-8, respectively, at different resolution scales. Then we estimate heat flux over leads with remote sensing temperature fields using both the bulk formulae and a fetch-limited model proposed by Andreas and Cash (1999). With the result, we analyze how the scale property of leads may affect the calculation of heat exchange through leads.

2 Data

Three successive scenes of Level 1 terrain-corrected (L1T) Landsat-8 TIRS images and one corresponding MODIS image acquired on 26 April 2015 were used in this study (Table 1). As shown in Fig. 1, the mosaic image of the three TIRS scenes covers an area of about 98 000 km2 in the marginal ice zone (MIZ) in the east Beaufort Sea, with floes and leads of various lengths and widths spread in the region. We obtained corresponding 10 m wind vector, 2 m air temperature, and dew point temperature from the European Center for Medium-Range Weather Forecasts (ECMWF) ERA-Interim reanalysis dataset. This dataset provides global coverage with a temporal resolution of 3 h; 0.125 grid data are available for download (∼10 km in study area, interpolated from original 0.75 grid). The time difference between reanalysis data and Landsat-8 or MODIS images is within half an hour.

Table 1Satellite images and other data used in this study.

Figure 1Location of study area. Background image is brightness temperature from Moderate Resolution Imaging Spectroradiometer (MODIS) band 31 (∼11µm). Location of Landsat-8 images is marked by a red rectangle.

Willmes and Heinamann (2015) used the MOD29 ice surface temperature (IST) product (Hall and Riggs, 2015) from the National Snow and Ice Data Center (NSIDC) to retrieve leads. The MOD29 product is filtered for cloud contamination using a cloud mask from MOD35. However, inspection of the corresponding MOD29 map of the study area revealed that the pixels within leads marked as clouds are likely open water with ocean fog or plume over the surface (Fett et al., 1997). Apart from that, the study area within the Landsat-8 frame is mostly unobstructed by clouds. To preserve potential lead areas, we applied the NSIDC algorithm (Hall et al., 2001) to thermal images from MODIS L1B to calculate IST instead of using the MOD29. Therefore, no cloud mask procedure was performed in our study.

The Landsat-8 satellite is in the same near-polar, sun-synchronous, 705 km circular orbit and position as the Landsat-5 satellite decommissioned in 2013. Landsat-8 data are acquired in 185 km swaths and segmented into 185 km ×180 km scenes defined in the second Worldwide Reference System (WRS-2) of path (ground track parallel) and row (latitude parallel) coordinates (Arvidson et al., 2001). The TIRS instrument, a major payload aboard Landsat-8, can observe the ocean surface at a resolution of 100 m by using split-window thermal infrared bands, comparable to MODIS thermal infrared bands, at a resolution of 1000 m. The two narrower thermal infrared channels in the atmospheric window enable application of the widely used split-window algorithm (SWA) in IST retrieval rather than the single-channel method used for TIRS predecessors.

Note that in the L1T product, the TIRS bands at 100 m resolution were resampled to 30 m by cubic convolution and co-registered with the Operational Land Imager (OLI) spectral bands. Apart from the TIRS thermal bands, the top of atmosphere reflectance from the Landsat-8 near-infrared band was used for classification between ice and open water in surface temperature retrieval. A panchromatic band with a resolution of 15 m was used as validation data for lead detection in this study.

3 Method

## 3.1 IST retrieval

Key et al. (1997) developed an SWA for IST retrieval from AVHRR, and the algorithm was then adapted for MODIS thermal images with a different set of coefficients (Hall et al., 2001). The equation is as follows:

$\begin{array}{}\text{(1)}& {T}_{\mathrm{s}}=a+b{T}_{\mathrm{31}}+c\left({T}_{\mathrm{31}}-{T}_{\mathrm{32}}\right)+d\left[\left({T}_{\mathrm{31}}-{T}_{\mathrm{32}}\right)\left(\mathrm{sec}q-\mathrm{1}\right)\right],\end{array}$

where T31 and T32 are brightness temperature from MODIS thermal bands B31 and B32; a, b, c, and d are coefficients developed for specific sensors using a radiance transfer model; q represents the incidence angle; and sec q is the secant of q.

Since there is no special SWA available for sea ice surface temperature retrieval from Landsat-8, a land surface temperature formulation (Du et al., 2015) developed for a wide range of surface types, including ice and snow, was used:

$\begin{array}{ll}{T}_{\mathrm{s}}& ={b}_{\mathrm{0}}+\left({b}_{\mathrm{1}}+{b}_{\mathrm{2}}\frac{\mathrm{1}-\mathit{\epsilon }}{\mathit{\epsilon }}+{b}_{\mathrm{3}}\frac{\mathrm{\Delta }\mathit{\epsilon }}{{\mathit{\epsilon }}^{\mathrm{2}}}\right)\frac{{T}_{i}+{T}_{j}}{\mathrm{2}}\\ \text{(2)}& & +\left({b}_{\mathrm{4}}+{b}_{\mathrm{5}}\frac{\mathrm{1}-\mathit{\epsilon }}{\mathit{\epsilon }}+{b}_{\mathrm{6}}\frac{\mathrm{\Delta }\mathit{\epsilon }}{{\mathit{\epsilon }}^{\mathrm{2}}}\right)\frac{{T}_{i}-{T}_{j}}{\mathrm{2}}+{b}_{\mathrm{7}}\left({T}_{i}-{T}_{j}{\right)}^{\mathrm{2}},\end{array}$

where Ti and Tj are the brightness temperatures measured in channels i (∼11.0µm) and j (∼12.0µm), respectively; ε is the mean emissivity for the two channels (ε=0.5 [εi+εj]); and Δε is the emissivity difference between the channels ($\mathrm{\Delta }\mathit{\epsilon }={\mathit{\epsilon }}_{i}-{\mathit{\epsilon }}_{j}\right)$; bk (k=0, 1,... 7) represents the algorithm coefficients derived from the simulated dataset.

As reported in previous studies (Montanaro et al., 2014a, b, c; Barsi et al., 2014), thermal infrared radiance measured by Landsat-8 TIRS suffers from stray light, which is caused by out-of-field radiance that scatters onto the detectors, adding a nonuniform banding signal across the field of view. The magnitude of this extra signal can be ∼8 % or higher (band 11) and is generally twice as large in band 11 as in band 10. Correction algorithms for this artifact have been developed and applied in the new version of level L1T Landsat-8 data (Montanaro et al., 2015), and the stray light artifact in the current product is reduced by half on average (Gerace and Montanaro, 2017). However, the artifact could be amplified in a surface temperature map when SWA is used, with a temperature error of 0–2 K or more (Gerace and Montanaro, 2017), which would certainly impact lead detection from IST maps. A post-processing procedure utilizing the linear pattern of the stray light artifact is applied to remove this banding noise. First, a median temperature is determined for each image pixel from a long enough along-track-only neighborhood. Then a noise image can be obtained by detrending this median image (Eppler and Full, 1992); thus the surface temperature image from SWA can be improved for lead detection.

In remote sensing images, leads (thin ice and open water) are represented by negative albedo anomalies in the optical range, negative brightness temperature anomalies in the near-infrared (NIR), and positive surface temperature anomalies compared to the surrounding thick ice (Fett et al., 1997). Variance caused by uneven illumination, view angle, and air temperature should also be taken into account.

Willmes and Heinemann (2015) reported the use of surface temperature anomalies to detect leads. The temperature anomaly ΔTs for each IST pixel is defined as a deviation from the median in a square neighborhood; thus temperature variation due to the air temperature field can be removed. This can be expressed in the following equation:

$\begin{array}{}\text{(3)}& \mathrm{\Delta }{T}_{\mathrm{s}}={T}_{\mathrm{s}}-{M}_{{T}_{\mathrm{s}},w},\end{array}$

where ${M}_{{T}_{\mathrm{s}},w}$ represents the median surface temperature in a square neighborhood with a side length of w. This equation was adapted for the Landsat-8 IST map using a median from an along-track-only linear neighborhood to further minimize the stray light artifact. Since median temperature is selected as background temperature, length w should be at least twice as large as the largest lead width within the image area (or along the track) to preserve the lead profile and reduce the temperature gradient caused by air temperature variation across the image.

Generally, surface temperature anomalies for thick ice follow normal distribution with a mean of zero; thus any large deviation from the mean can be assumed as a potential lead and extracted using a proper threshold. Several image-based threshold selection techniques for binary lead segmentation were compared in Willmes and Heinemann (2015), and an iterative threshold selection method (Ridler and Calvard, 1978) was recommended for extracting leads from a temperature anomaly map. Assuming an initial threshold using the mean temperature anomaly (m0) of the whole image, the iterative method seeks a threshold mi which is the mean of averages mA and mB from two clusters divided by mi: leads (A) and pack ice (B). However, any image-based threshold method provides a threshold that can vary significantly due to temperature noise and lead distribution. For consistency in different scales, several threshold methods were compared for lead detection in both MODIS and TIRS temperature maps (see Sect. 5.1), and an iterative threshold method was adopted.

Using width samples crossed by transects, Lindsay and Rothrock (1995) found a mean lead width between 2 and 3 km in the Arctic winter; larger means are found in peripheral seas. We modified the method by using an orthogonal system (vertical, south–north; horizontal, west–east; Fig. 2) to determine lead width for every lead pixel. A minimum lead extent in two orthogonal directions was selected for the pixel; i.e., X= min (x1, x2). Since the orientation of a single lead is unknown, this method tends to overestimate width x due to a mismatch between the preset direction and the orientation of the lead (Key and Peckham, 1991), but the orthogonal system will help contain the error ($X\le \sqrt{\mathrm{2}}x\right)$. Since we assign lead width to every pixel across the lead, length Li for lead width Xi can be calculated as follows:

$\begin{array}{}\text{(4)}& {L}_{i}=\frac{{a}_{\mathrm{0}}^{\mathrm{2}}{N}_{i}}{{X}_{i}}=\frac{{a}_{\mathrm{0}}{N}_{i}}{i},\end{array}$

where a0 is the pixel size; for TIRS, the value is 30 m, for MODIS, 1 km; and Ni is the number of pixels for width Xi=a0i, (i=1, 2, 3...).

Figure 2Detection of lead width using two orthogonal directions. x is the real lead width. Lead extents in orthogonal system in v and h directions are measured as x1 and x2, respectively.

## 3.3 Heat flux model used for lead area

Turbulent heat flux between the Arctic Ocean and the atmosphere, including sensible (Fs) and latent (Fl) heat flux, is mostly dominated by heat flux over open water and thin ice, which constitute leads in pack ice and polynya in coastal areas. Turbulent heat flux over leads can be estimated using bulk aerodynamic formulae or a fetch-limited model developed based on field observations.

### 3.3.1 Bulk aerodynamic formulae

Assuming that temperatures in the atmospheric boundary layer are determined by the heat balance over thicker ice and turbulent heat exchange does not vary significantly across the narrow areas of leads, then turbulent heat fluxes are mainly determined by temperature and humidity differences between the surface and atmosphere at reference height r (Maykut, 1978). Sensible heat flux (Fs) and latent heat flux (Fl) are given by the following bulk formulae:

$\begin{array}{}\text{(5)}& {F}_{\mathrm{s}}={\mathit{\rho }}_{\mathrm{a}}{c}_{\mathrm{p}}{C}_{\mathrm{sh}}{u}_{\mathrm{r}}\left({T}_{\mathrm{s}}-{T}_{\mathrm{r}}\right)\end{array}$

$\begin{array}{}\text{(6)}& {F}_{\mathrm{l}}={\mathit{\rho }}_{\mathrm{a}}{L}_{\mathrm{v}}{C}_{\mathrm{le}}{u}_{\mathrm{r}}\left({Q}_{\mathrm{s}}-{Q}_{\mathrm{r}}\right),\end{array}$

where ρa is the air density; cp is the specific heat at constant pressure; Lv is the latent heat of vaporization; ur, Tr, and Qr are wind speed, air temperature, and specific humidity at reference height r=2 m, respectively; Ts is surface temperature; and Qs is specific humidity close to the surface. Assuming that air at the surface of ice or water is always saturated, the specific humidity at the surface can be derived as

$\begin{array}{}\text{(7)}& {Q}_{\mathrm{s}}=\frac{\mathrm{0.622}{e}_{\mathrm{s}\mathrm{0}}}{p-\mathrm{0.378}{e}_{\mathrm{s}\mathrm{0}}},\end{array}$

where p is the air pressure, and es0 represents the saturated vapor pressure at surface temperature Ts:

$\begin{array}{}\text{(8)}& {e}_{\mathrm{s}\mathrm{0}}={e}_{\mathrm{0}}{\mathrm{10}}^{\frac{at}{b+t}},\end{array}$

with e0 representing saturated vapor pressure at 0 C, approximately 6.11 hPa; t is temperature in Celsius; and a and b are coefficients (for water surface, a=7.5, b=237.3 K; for thin ice, a=9.5, b=265.5 K). These equations are also applied for specific humidity at 2 m height using dew point temperature data from ERA-Interim.

Csh and Cle are transfer coefficients for sensible heat and latent heat, calculated using equations from Oberhuber (1988) and Goosse et al. (2001) (see Appendix B). Since the wind speed and air temperature from ERA-Interim are at different heights, a wind profile equation was used to calculate wind speed at 2 m height (Ray et al., 2006):

$\begin{array}{}\text{(9)}& \frac{{u}_{\mathrm{10}}}{{u}_{\mathrm{r}}}=\frac{\mathrm{ln}\mathrm{10}-\mathrm{ln}{Z}_{\mathrm{0}}}{\mathrm{ln}r-\mathrm{ln}{Z}_{\mathrm{0}}},\end{array}$

where u10 and ur are wind speed at 10 and 2 m height, and Z0 is surface roughness length. In our study area, the main direction of wind from the reanalysis dataset is roughly perpendicular to the dominant orientation of leads. Therefore, only the wind magnitude was used in our study.

### 3.3.2 Fetch-limited model

When cold air travels to a warmer surface, a convective atmospheric TIBL forms and thickens with distance downwind of the surface discontinuity or fetch X (Stull, 1988). As the wind blows over water (or thin ice), the near-surface air gets warmer with more vapor, while new ice accumulates at the downwind side of the lead, progressively narrows, and seals the window. Thus, the temperature and humidity differences between the air and the surface decrease. Since sensible and latent heat fluxes are proportional to temperature and humidity differences, respectively, turbulent heat transfer also recedes with increasing lead width or fetch. Another mechanism is described in Esau (2007) for leads 1–10 km wide. Under weak wind conditions (∼2 m s−1), convective overturning prevents cold breezes from penetrating into the lead area, reducing the average turbulent heat flux.

To estimate turbulent heat flux over small leads, fetch-limited models were developed based on a few observations (Andreas and Murphy, 1986; Alam and Curry, 1997; Andreas and Cash, 1999). However, the assumption of universal water surface in leads and the application of sea surface roughness model (Andreas and Murphy, 1986; Alam and Curry, 1997) are not applicable in our case, where open water and thin ice dominate. Since the signal of TIBL is absent in the coarse grid of 2 m air temperature from the ERA reanalysis dataset, the data might not be appropriate to demonstrate the Alam and Curry (1997) model, which relies on the accurate measurement of meteorological parameters, whereas the Andreas and Cash (1999) model is more sensitive to lead width than atmospheric conditions (Marcq and Weiss, 2012). Therefore, only the Andreas and Cash (1999) model was used in our experiment.

Andreas and Cash (1999) gave direct formulations of heat fluxes as a function of lead width X based on data fitting from three sets of measurements. For free convection conditions in large fetch,

$\begin{array}{}\text{(10)}& {F}_{\mathrm{s}\left(X\right)}={C}_{\ast }{\mathit{\rho }}_{\mathrm{a}}{C}_{\mathrm{p}}D\left({T}_{\mathrm{s}}-{T}_{\mathrm{r}}\right)/\mathrm{\Delta }{z}_{\mathrm{T}}\end{array}$

$\begin{array}{}\text{(11)}& {F}_{\mathrm{l}\left(X\right)}={C}_{\ast }{\mathit{\rho }}_{\mathrm{a}}{L}_{\mathrm{v}}{D}_{\mathrm{w}}\left({Q}_{\mathrm{s}}-{Q}_{\mathrm{r}}\right)/\mathrm{\Delta }{z}_{\mathrm{Q}},\end{array}$

where D and Dw are the molecular diffusivities of heat and water vapor in air, respectively, and ΔzT and ΔzQ are length scales for heat and humidity, respectively, which consider the viscosity of air v and buoyancy differences between the surface and reference height r:

$\begin{array}{}\text{(12)}& \mathrm{\Delta }{z}_{\mathrm{T}}={\left(\frac{vD}{\mathrm{\Delta }B}\right)}^{\mathrm{1}/\mathrm{3}}\end{array}$

$\begin{array}{}\text{(13)}& \mathrm{\Delta }{z}_{\mathrm{Q}}={\left(\frac{v{D}_{w}}{\mathrm{\Delta }B}\right)}^{\mathrm{1}/\mathrm{3}}\end{array}$

$\begin{array}{}\text{(14)}& \mathrm{\Delta }B=\frac{\mathrm{g}}{\stackrel{\mathrm{‾}}{T}}\left(\mathrm{\Delta }T+\frac{\mathrm{0.61}\stackrel{\mathrm{‾}}{T}\mathrm{\Delta }Q}{\mathrm{1}+\mathrm{0.61}\stackrel{\mathrm{‾}}{Q}}\right),\end{array}$

where ΔB is the buoyancy difference; g is acceleration due to gravity; ΔT and ΔQ are the difference of temperature and specific humidity between surface and air at reference height r, respectively; and $\stackrel{\mathrm{‾}}{T}$ and $\stackrel{\mathrm{‾}}{Q}$ are the average temperature and specific humidity between them.

The coefficient C is a function of stability, which facilitates the generalization of Eqs. (10) and (11) to the transition between free and forced convection, thus making them applicable to smaller fetch. C is estimated using lead and polynya data:

$\begin{array}{}\text{(15)}& {C}_{\ast }=\frac{\mathrm{0.3}}{\mathrm{0.4}-h/L}+\mathrm{0.15}\end{array}$

$\begin{array}{}\text{(16)}& h=\mathrm{0.82}\mathrm{ln}X+\mathrm{0.02},\end{array}$

where h is the depth of the TIBL in meters as a function of lead width X, and L is the Obukhov length given in Eq. (17); L is a length scale of stability and is negative for unstable stratification, while its magnitude rises with instability.

$\begin{array}{}\text{(17)}& {L}^{-\mathrm{1}}=\mathrm{8.0}\ast \left(\frac{\mathrm{0.65}}{r}+\mathrm{0.079}-\mathrm{0.0043}r\right)\ast R{i}_{\mathrm{b}},\end{array}$

where Rib is the bulk Richardson number:

$\begin{array}{}\text{(18)}& R{i}_{\mathrm{b}}=-\frac{rg}{\stackrel{\mathrm{‾}}{T}}\frac{{T}_{\mathrm{s}}-{T}_{\mathrm{r}}}{{u}_{\mathrm{r}}^{\mathrm{2}}},\end{array}$

where ur is wind speed obtained from Eq. (9). Apart from lead width, meteorological parameters are also important in the model. As shown in Fig. 3, for the narrowest lead from TIRS (X=30 m), turbulent heat flux, especially sensible heat, rises quickly with larger temperature difference and stronger wind. Most importantly, assuming a constant temperature difference and steady crossing wind, heat flux decreases with increasing fetch and becomes saturated for lead width great than 1 km, as depicted in Fig. 4. As the fetch dependence of heat flux over lead is negligible for lead widths greater than 1 km, this model was applied to TIRS data only.

Figure 3Turbulent heat flux rises with increasing temperature difference ΔT and intense wind at lead width of 30 m. Solid and dashed lines represent sensible and latent heat, respectively. Wind speed is illustrated by line color. Clearly, sensible heat flux is basically proportional to ΔT.

Figure 4Turbulent heat flux for each width at wind speed of 5 m s−1. Temperature difference between air and lead surface is marked by line color.

4 Results

## 4.1 Ice surface temperature

IST maps retrieved from MODIS and TIRS using Eqs. (1) and (2) are shown in Fig. 5. The temperature signature of small leads in the northern part of the image area is largely reduced in the MODIS IST map, due to its coarse resolution and heterogeneous pixels, compared to that from TIRS. In addition, the banding effect of stray light is very obvious in the TIRS IST map. This artifact was detected and removed by using a median from the along-track linear neighborhood and detrending the median image (Fig. 6). The corrected TIRS IST map is shown in Fig. 5 for comparison.

Figure 5Ice surface temperature (IST) maps from MODIS and Landsat-8 Thermal Infrared Sensor (TIRS) using split-window algorithms: (a) IST map from MODIS, (b) IST map from Landsat-8 TIRS, and (c) corrected IST map from TIRS.

Although the median and artifact images show a little bias around large leads, the corrected TIRS IST map is very smooth and more suitable for lead detection and heat flux calculation. Scatter plots of IST from MODIS and TIRS before and after correction are shown in Fig. 7. The correlation of IST from two sensors estimated by interpolating MODIS IST to the TIRS scale (30 m) is quite good, with a Pearson coefficient of approximately 0.9 (0.902 and 0.896 before and after correction for stray light, respectively). The primary coefficient of linear regression improved from 1.025 to 1.004 before and after correction, indicating that the corrected TIRS IST is in better agreement with MODIS. However, the root mean square error (RMSE) from regressions increased from 1.216 to 1.233 K. It also reveals that for the 250–270 K temperature range, IST retrieved from TIRS is generally 0.61–0.70 K higher than that from MODIS.

Figure 6Local median and noise image from TIRS IST: (a) along-track median temperature map and (b) noise image by detrending of median temperature map.

Figure 7Correlation between IST from MODIS and Landsat-8 TIRS before and after correction for stray light. Black lines are reference for y=x, and red lines are linear regression lines with a fitting equation. Number density of scattered points is marked by color. (a) Scatter plot of IST from MODIS and Landsat-8 TIRS using split-window algorithm. (b) Scatter plot of IST from MODIS and corrected IST from Landsat-8 TIRS.

## 4.2 Sea ice lead retrieval

Regional temperature anomaly maps calculated from IST maps are shown in Fig. 8. The mean surface temperature anomaly is 0.116 K with a standard deviation (SD) of 1.180 K for MODIS and 0.283 K with a SD of 1.619 K for TIRS.

Figure 8Local temperature anomalies from (a) MODIS and (b) Landsat-8 TIRS.

Binary lead maps were generated using iterative thresholds (Fig. 9). Large floes and small leads dominate the northern part of the images, where temperature is lower, while two very large leads can be observed in the southern portion. The maps illustrate that the lead binary retrieved from MODIS captures major lead structures, but small leads are missed in most cases compared to leads detected from TIRS.

Figure 9Binary lead maps from (a) MODIS and (b) Landsat-8 TIRS.

Lead area estimated from MODIS is 8074.0 km2, which accounts for 8.25 % of the frame area, and from TIRS, 7376.2 km2 and 7.53 %. Validation with Landsat-8 panchromatic images shows that large leads tend to be amplified by blurred mixed pixels along boundaries, while small leads are neglected due to the coarse resolution of MODIS.

Lead width was calculated for every lead pixel in the binary maps from MODIS and TIRS and divided into three categories (Table 2): small leads (width ≤1 km), medium leads (1 km < width  5 km), and large leads (width > 5 km). Although the 1 km resolution is the finest for MODIS thermal images, the 1 km wide lead category should provide a reasonable guess of potential small leads or subpixel leads at MODIS scale (Lindsay and Rothrock, 1995).

Table 2Retrieved leads from MODIS and TIRS and turbulent heat flux estimated using bulk formulae.

The width distribution of leads from MODIS and small leads from TIRS is plotted in Fig. 10 relating to the lengths of leads. Similar to the concept of number density, the length of each lead width can be fitted with a power law distribution, and the exponents from power law fitting are 2.241 and 2.346 for leads from MODIS and TIRS, respectively. The power law distribution indicates that narrow leads are prevalent, while a larger exponent implies that smaller leads are more dominant at TIRS scale.

Figure 10Width distribution of leads from MODIS and TIRS in a log–log plot. Data points from MODIS and TIRS are plotted as orange and blue dots, respectively. Power law fitting is applied. Fitting equations and R squares are shown for comparison.

The total length of leads with various widths is 10 150.3 km from TIRS, including 8502.2 km (83.76 %) from small leads with width no more than 1 km, compared to a total length of 2746.4 km from MODIS, where the narrow leads (1 km wide) only account for 1050.0 km (38.23 %). Total length of leads is underestimated by 72.9 % in MODIS data compared to TIRS data. As for the area of leads, small leads (width ≤1 km) account for 34.54 % of total lead area from TIRS and only 13.00 % of lead area from MODIS (Table 2).

## 4.3 Heat flux over leads

IST, described in Sect. 4.1, and lead width from TIRS (Sect. 4.2) were used in bulk formulae and the fetch-limited model along with ERA-Interim reanalysis data to estimate turbulent heat flux through leads. For consistency, the estimated heat flux is positive from the ocean to the atmosphere.

### 4.3.1 Bulk formulae

Turbulent heat flux over lead area is obtained by summing up sensible and latent heat flux from Eqs. (5) and (6) using IST and lead maps retrieved from MODIS or TIRS (Fig. 11). Table 2 reveals that total heat flux over leads calculated using TIRS IST is 8.40×1011 W over a total area of 7376.2 km2. This is 56.70 % larger than that from MODIS data (5.36×1011 W). About 23 % of the difference can be explained by IST bias between MODIS and TIRS, but most of the difference comes from small leads. Small leads account for 2.16×1011 W (25.75 %) of total heat flux in TIRS data, almost 7 times the heat flux from the narrow lead category in MODIS (3.10×1010 W, 5.79 %).

### 4.3.2 The Andreas and Cash (1999) model

Figure 11Spatial distribution of heat flux derived from MODIS and Landsat-8 using bulk formulae and a fetch-limited model. (a) Turbulent heat flux from MODIS using bulk formulae, (b) turbulent heat flux from Landsat-8 TIRS using bulk formulae, and (c) turbulent heat flux from Landsat-8 TIRS using a fetch-limited model.

Table 3Estimated turbulent heat flux (W) for Landsat-8 TIRS using bulk formulae, the Andreas and Cash (1999) model, and the modified Andreas and Cash model using Obukhov length from Eqs. (B8) and (B13).

Inspection of input data revealed that the 2 m air temperature from ERA-Interim has almost the same mean value around 262 K as the surface temperature from Landsat-8. The temperature difference between air and surface, ΔT, spreads from 1.58 to 12.38 K, with a mean of 4.88 K, along with an average wind speed of about 7 m s−1 at 2 m height over leads. The distributions of air temperature and surface temperature of leads are plotted in Fig. 12. The temperature difference over narrow leads is too small to obtain a robust estimation of turbulent heat flux.

5 Discussion

## 5.1 Threshold method

The operational definition of a lead is a fracture or passageway through ice that is navigable by surface vessels (Canadian Ice Service, 2005; World Meteorological Organization, 2014). However, within any optical, thermal, or microwave image, the radiometric signature of a narrow lead with open water may be identical to that of a wider lead with thin ice. In most studies involving the utility of remote sensing data, any linear features of open water or thin ice within pack ice are regarded as leads for convenience (Fetterer and Holyer, 1989; Fily and Rothrock, 1990; Lindsay and Rothrock, 1995). Due to the confusion in the definition of leads in remote sensing images and the need to extract lead signatures from the background, threshold segmentation has been frequently used (Eppler and Full, 1992; Lindsay and Rothrock, 1995; Weiss and Marsan, 2004; Marcq and Weiss, 2012). Willmes and Heinemann (2015) presented several threshold selection techniques for binary lead segmentation. However, thresholds given by image-based methods can vary significantly depending on noise level (caused by air temperature variance) and lead distribution.

In our study, a set of thresholds was tested for extracting leads from temperature anomaly maps, areal fractions of leads from fixed thresholds, SD thresholds, and an iterative threshold are shown in Table 4. The obtained lead fractions are a composite of thresholds and contrast in surface temperature of leads compared to the surrounding ice, i.e., temperature anomaly ΔTS. Since the anomaly maps from the two sensors have different means and standard deviations, mainly due to different definitions of neighborhood in calculating ΔTS, the results from a fixed threshold might be biased. The iterative thresholds from both anomaly maps are a little larger than their first SD thresholds. The difference in lead fractions from the two sensors mainly resulted from mixed pixels at MODIS scale, but the threshold should also be considered. When high thresholds (second and third SD) are applied, the lead fraction extracted from MODIS drops quickly below that from TIRS (as shown in Table 4), and this is consistent with results from Key et al. (1994). While larger thresholds lead to the underestimation of lead distribution, lower thresholds allow more pixels to be detected as leads, also giving rise to false leads caused by air temperature variance.

Validation with Landsat-8 panchromatic images shows that the iterative threshold detects most lead structures (89.5 %) and exhibited better resistance against air temperature noise. Thus, iterative thresholds were selected for lead extraction in this study.

Figure 12Distribution of 2 m air temperature over leads and surface temperature of all leads, small leads with width <1 km, and larger leads with width >5 km.

Lead geometry and distribution in the Arctic have been studied using optical and microwave remote sensing data (Fily and Rothrock, 1990; Lindsay and Rothrock, 1995; Tschudi et al., 2002). A simple one-parameter exponential model was used for the number density distribution of lead width (Key and Peckham, 1991; Key et al., 1994; Maslanik and Key, 1995):

$\begin{array}{}\text{(19)}& {f}_{\left(X\right)}=\frac{\mathrm{1}}{\mathit{\lambda }}{e}^{\frac{-X}{\mathit{\lambda }}},\end{array}$

where λ is the mean lead width. However, a mean lead width can be oversimplified in diverse circumstances. Lindsay and Rothrock (1995) reported the power law distribution of lead width in AVHRR imagery:

$\begin{array}{}\text{(20)}& {N}_{\mathrm{T}\left(X\right)}=a{X}^{-b},\end{array}$

where NT(w) is the number density of leads of width w per kilometer of width increment, and the exponent b indicates the relative frequency of large and small leads, while the coefficient a is directly related to the lead concentration and the range of widths over which the power law is thought to apply. The annual mean of exponent b was found to be 1.60 using AVHRR images (Lindsay and Rothrock, 1995). Larger values of b were reported using data with better resolution: 2 and 2.29 for submarine sonar observation in the Fram Strait (Wadhams, 1981) and the Davis Strait (Wadhams et al., 1985) when a 100 m bin width was used and 2.1–2.6 for 10 m SPOT images in orthographic directions using different thresholds (Marcq and Weiss, 2012). Note that most of these studies used only width samples crossed by limited linear transects.

In our study, although lead width follows the power law distribution at both scales, the fitted exponents vary from 2.241 to 2.346 at resolution from 1 km to 30 m. Since the 30 m L1T images were resampled from the original 100 m TIRS data, the actual distribution of leads less than 100 m wide is debatable. In comparison with Landsat-8 TIRS and panchromatic images, we find that the lead map generated from the MODIS IST data neglects very small leads but overestimates the width of other leads approximately 1 km wide. Overall, the 1 km wide lead category at MODIS scale should provide a reasonable guess of potential small or subpixel leads. The small leads retrieved using TIRS provide a valuable reference for the capacity of MODIS to detect narrow leads.

## 5.3 Comparison of the models

In both the Andreas and Murphy (1986) and Andreas and Cash (1999) models, for reference height r<10 m, the ratio between roughness lengths for momentum and heat, Z0ZT, is assumed to be e2 to calculate Obukhov length L using the Richardson number Rib (see Eq. 17). The calculated Obukhov length L has absolute values about 67 % higher than those using Eqs. (B8) and (B13) from the bulk formulae (Oberhuber, 1988; Goosse et al., 2001). If Eq. (B8) and (B13) were used to solve Obukhov length L and coefficient C in the Andreas and Cash (1999) model, estimated turbulent heat flux would be smaller (Table 3) but still 15.53 % larger than that from the bulk formulae, with an even larger part of the difference from the small lead category (42.48 %, compared to 32.95 % in Sect. 4.3.2).

Our results suggest that the contribution of heat flux from small leads mainly results from their large length, or number density, and vast area instead of efficiency. Though small leads are more efficient for heat exchange between the ocean and the atmosphere, thin ice growing in newly opened leads can quickly cover the exposed ocean surface, thus reducing heat exchange. Moreover, due to the mixture of subpixel leads and thick ice, the surface temperature of some pixels in small leads is much lower than the freezing point.

Nonetheless, our results show that the fetch-limited model could be used to estimate turbulent heat flux on a regional scale with surface temperature fields from remote sensing. However, the fetch-limited model proposed by Andreas and Cash (1999) was based mainly on a few observations over open leads and polynya, while most lead pixels detected using temperature anomalies in our study are likely covered by thin ice (surface temperature <270 K; Fig. 12). Thus, near-surface air temperature with finer resolution is needed for validating the turbulent heat flux estimated using the fetch-limited model.

## 5.4 Heat flux for larger temperature differences

For comparison, a test using preset meteorological conditions was performed using the TIRS lead binary. Assuming the surface temperature in leads is right at the freezing point, with a wind speed of 7 m s−1 at 2 m height and a temperature difference of 5 and 10 K, turbulent heat fluxes from both models were calculated (Fig. 13) and are summarized in Table 5. Note that lead width in Fig. 13 is on a logarithmic scale.

Figure 13Contribution of heat flux from each lead width using bulk formulae and a fetch-limited model. Turbulent heat flux retrieved using a fetch-limited model and bulk formulae are plotted using solid and dashed lines, respectively. Heat flux calculated using satellite surface temperature, air temperature, and wind speed from reanalysis datasets is drawn in orange; simulated heat flux at ΔT=5 and 10 K is in blue and green, respectively.

Table 5Turbulent heat flux (W) for higher temperature differences using Landsat-8 TIRS data and the Andreas and Cash (1999) model.

Clearly, turbulent heat flux estimated using the Andreas and Cash (1999) model is always higher than that using the bulk formulae. For both models, estimated turbulent heat flux with ΔT of 5 or 10 K peaks at lead width of ∼270 m, a smaller width than the 360 m using ΔT obtained from TIRS IST and air temperature from ERA-Interim.

The distribution of turbulent heat flux estimated using bulk formulae with ΔT of 5 and 10 K depends on the areal fraction from each lead category. The contribution from leads with widths greater than 1 km converges to the lower end with fluctuation. As expected, the estimated total heat flux of 1.68×1012 W at ΔT=10 K is about twice as large as that at ΔT=5 K (8.46×1011 W).

When the Andreas and Cash (1999) model was applied, small leads were found to have a larger contribution at higher ΔT, 3.27×1011 W (35.86 %) and 6.66×1011 W (36.57 %) at ΔT=5 and 10 K, respectively, compared to the areal fraction of 34.54 %. More contributions from small leads can be expected at larger temperature differences and stronger wind in winter.

6 Conclusions

Although the same local temperature anomaly and threshold methods were applied, leads retrieved at MODIS and Landsat-8 TIRS resolution scales presented very different geometries and distributions. Within the studied area, the total length of leads is 10 150.3 km from TIRS, including 8502.2 km (83.76 %) from small leads with width less than 1 km. This is in contrast to the total length of 2746.4 km from MODIS, where the narrow leads (1 km wide) only account for 1050.0 km (38.23 %). The total length of leads is underestimated by 72.9 % in the MODIS data. For the area of leads, small leads (width ≤1 km) account for 34.54 % of the total lead area from TIRS but only 13.00 % of the total lead area from MODIS. Although the lead width follows the power law distribution at both scales, the fitted exponents vary from 2.241 to 2.346.

When bulk aerodynamic formulae are applied to the reanalysis dataset, the heat flux estimated using TIRS data is 8.40×1011 W, 56.70 % larger than that from MODIS data (5.36×1011 W). About 23 % of the difference can be explained by IST bias between MODIS and TIRS, but most of the difference comes from small leads. Small leads account for 2.16×1011 W (25.75 %) of the total heat flux over all leads in the TIRS data, almost 7 times the heat flux from the narrow lead category in MODIS (3.10×1010 W, 5.79 %).

The turbulent heat flux over leads estimated from the TIRS data by the Andreas and Cash (1999) model is 1.11×1012 W, 32.34 % higher than that from bulk formulae (8.40×1011 W). In both cases, small leads account for about a quarter of the total heat flux in both models, due to the large area, though the heat flux estimated using the fetch-limited model is 41.39 % larger. A greater contribution from small leads can be expected with larger temperature differences and stronger wind conditions. A near-surface air temperature with finer resolution is still needed for the validation of turbulent heat flux estimated using the fetch-limited model before extensive application.

Appendix A: Validation using Landsat-8 panchromatic images

Top of atmosphere reflectance from Landsat-8 panchromatic images was corrected for solar zenith angle and mosaicked for validation. Using Jenks' natural breaks classification method (Jenks, 1963), panchromatic pixels were classified into three surface categories: open water and thin ice, refrozen leads, and pack ice. In terms of turbulent heat flux, only pixels in the open water and thin ice category were regarded as leads. As can be seen in Table A1, the producer's accuracy of lead detection using the iterative threshold is 89.5 %, with an omission error of 10.5 % and a commission error of 16.1 %.

Table A1Leads and pack ice pixels detected by Landsat-8 TIRS and panchromatic images at 15 m resolution.

Appendix B

Equations used for turbulent heat flux estimation using bulk formulae (Large and Pond, 1981, 1982; Oberhuber, 1988; Goosse et al., 2001; Marcq and Weiss, 2012) are as follows:

$\begin{array}{}\text{(B1)}& {c}_{\mathrm{sh}}=\mathrm{0.0327}\frac{k}{\mathrm{ln}\left(r/{z}_{\mathrm{0}}\right)}{\mathrm{\Phi }}_{\mathrm{sh}}={c}_{\mathrm{shN}}{\mathrm{\Phi }}_{\mathrm{sh}}\end{array}$

$\begin{array}{}\text{(B2)}& {c}_{\mathrm{le}}=\mathrm{0.0346}\frac{k}{\mathrm{ln}\left(r/{z}_{\mathrm{0}}\right)}{\mathrm{\Phi }}_{\mathrm{le}}={c}_{\mathrm{leN}}{\mathrm{\Phi }}_{\mathrm{le}}\end{array}$

$\begin{array}{}\text{(B3)}& {\mathrm{\Phi }}_{\mathrm{sh}}=\frac{\sqrt{{c}_{\mathrm{M}}/{c}_{\mathrm{MN}}}}{\mathrm{1}-{c}_{\mathrm{shN}}{k}^{-\mathrm{1}}{C}_{\mathrm{MN}}^{-\mathrm{1}/\mathrm{2}}{\mathrm{\Psi }}_{\mathrm{H}}}\end{array}$

$\begin{array}{}\text{(B4)}& & {\mathrm{\Phi }}_{\mathrm{le}}=\frac{\sqrt{{c}_{\mathrm{M}}/{c}_{\mathrm{MN}}}}{\mathrm{1}-{c}_{\mathrm{leN}}{k}^{-\mathrm{1}}{C}_{\mathrm{MN}}^{-\mathrm{1}/\mathrm{2}}{\mathrm{\Psi }}_{\mathrm{L}}}\text{(B5)}& & \sqrt{\frac{{c}_{\mathrm{M}}}{{c}_{\mathrm{MN}}}}=\frac{\mathrm{1}}{\left(\mathrm{1}-\sqrt{{c}_{\mathrm{MN}}}{k}^{-\mathrm{1}}{\mathrm{\Psi }}_{\mathrm{M}}\right)}\text{(B6)}& & {c}_{\mathrm{MN}}=\frac{{k}^{\mathrm{2}}}{{\left(\mathrm{ln}\left(\frac{r}{{z}_{\mathrm{0}}}\right)\right)}^{\mathrm{2}}}\text{(B7)}& & {u}_{\ast }^{\mathrm{2}}={c}_{\mathrm{M}}{u}_{\mathrm{r}}^{\mathrm{2}}\text{(B8)}& & {T}_{\mathrm{0}}={T}_{\mathrm{r}}\left(\mathrm{1}+\mathrm{2.2}×{\mathrm{10}}^{-\mathrm{3}}{T}_{\mathrm{r}}{Q}_{\mathrm{r}}\right).\end{array}$

Surface roughness lengths for momentum are given as

$\begin{array}{}\text{(B9)}& {z}_{\mathrm{0}}=\mathrm{0.032}\frac{{u}_{\ast }^{\mathrm{2}}}{g}.\end{array}$

For unstable conditions,

$\begin{array}{ll}\text{(B10)}& & {\mathrm{\Psi }}_{\mathrm{H}}\left(A\right)={\mathrm{\Psi }}_{\mathrm{L}}\left(A\right)=\mathrm{2}\mathrm{ln}\left(\frac{\mathrm{1}+{A}^{\mathrm{2}}}{\mathrm{2}}\right)& {\mathrm{\Psi }}_{\mathrm{M}}\left(A\right)=\mathrm{2}\mathrm{ln}\left(\frac{\mathrm{1}+A}{\mathrm{2}}\right)+\mathrm{ln}\left(\frac{\mathrm{1}+{A}^{\mathrm{2}}}{\mathrm{2}}\right)\\ \text{(B11)}& & -\mathrm{2}\mathrm{arctan}A+\frac{\mathit{\pi }}{\mathrm{2}}\text{(B12)}& & A={\left(\mathrm{1}-\mathrm{16}\left(r/L\right)\right)}^{\mathrm{1}/\mathrm{4}}\text{(B13)}& & r/L=\frac{\mathrm{100}r}{{T}_{\mathrm{0}}{u}_{\mathrm{r}}^{\mathrm{2}}}\left(\left({T}_{\mathrm{s}}-{T}_{\mathrm{r}}\right)+\mathrm{2.2}×{\mathrm{10}}^{-\mathrm{3}}{T}_{\mathrm{0}}^{\mathrm{2}}\left({Q}_{\mathrm{s}}-{Q}_{\mathrm{r}}\right)\right).\end{array}$
Appendix C

Constants used in IST calculation from Landsat-8 TIRS (Du et al., 2015) are as follows.

1. ASTER emissivity library (Skoković et al., 2014):
${\mathit{\epsilon }}_{\mathrm{water},\mathrm{10}}=\mathrm{0.991}$; ${\mathit{\epsilon }}_{\mathrm{water},\mathrm{11}}=\mathrm{0.986}$; ${\mathit{\epsilon }}_{\mathrm{snow}/\mathrm{ice},\mathrm{10}}=\mathrm{0.986}$; ${\mathit{\epsilon }}_{\mathrm{snow}/\mathrm{ice},\mathrm{11}}=\mathrm{0.959}$
${\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}_{\mathrm{water}}=\mathrm{0.9885}$; Δεwater=0.005
${\stackrel{\mathrm{‾}}{\mathit{\epsilon }}}_{\mathrm{snow}/\mathrm{ice}}=\mathrm{0.9725}$; $\mathrm{\Delta }{\mathit{\epsilon }}_{\mathrm{snow}/\mathrm{ice}}=\mathrm{0.027}$.

2. NIR reflectance threshold for classification between water and ice/snow: 0.1.

3. Water vapor content from MOD05: <2.5 g cm−2.

4. bi: b0∼7: [−2.78009, 1.01408, 0.15833, −3.4991, 4.04487, 3.55414, −8.88394, 0.09152].

5. RMSE: 0.34 K.

Constants used in turbulent heat flux estimation are as follows.

• Air pressure: p=101 kPa.

• Air density: ρa=1.3 kg m−3.

• Kinematic viscosity of air: $v=\mathrm{1.31}×{\mathrm{10}}^{-\mathrm{5}}$ m2 s−1.

• Molecular diffusivities of heat in the air: $D=\mathrm{1.86}×{\mathrm{10}}^{-\mathrm{5}}$ m2 s−1.

• Molecular diffusivities of water vapor in the air: ${D}_{\mathrm{w}}=\mathrm{2.14}×{\mathrm{10}}^{-\mathrm{5}}$ m2 s−1.

• Specific heat at constant pressure: cp=1004 J kg−1 K−1.

• Latent heat of vaporization or sublimation: ${L}_{\mathrm{water}}=\mathrm{2.51}×{\mathrm{10}}^{\mathrm{6}}$ J kg−1, ${L}_{\mathrm{ice}}=\mathrm{2.86}×{\mathrm{10}}^{\mathrm{6}}$ J kg−1.

• Reference height: r=24 m. Von Karman constant: k=0.4.

• Gravitational constant: g=9.8 m s−2.

• Salinity of seawater in the Beaufort Sea: Sw=27.947 (‰).

• Freezing point of seawater:

• ${T}_{\mathrm{S}\mathrm{0}}=\mathrm{273.15}-\mathrm{0.0137}-\mathrm{0.05199}$ Sw−0.00007225 ${S}_{\mathrm{w}}^{\mathrm{2}}=\mathrm{271.68}$ K.

Data availability
Data availability.

Landsat-8 L1T images are available at the U.S. Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center (http://glovis.usgs.gov/, last access: 25 May 2019; Zanter, 2019). MODIS L1B images are available at NASA Level-1 and Atmosphere Archive & Distribution System Distributed Active Archive Center (LAADS DAAC) (https://ladsweb.modaps.eosdis.nasa.gov/, last access: 25 May 2019; MCST, 2015). ERA-Interim reanalysis datasets are available at the European Centre for Medium-Range Weather Forecasts (https://www.ecmwf.int/, last access: 25 May 2019; Berrisford et al., 2011).

Author contributions
Author contributions.

XP and XZ designed the experiments, and MQ carried them out. JZ provided valuable instructions on data acquisition and editing of the paper. QJ and PF helped to develop the model code. MQ prepared the paper with contributions from all co-authors.

Competing interests
Competing interests.

The authors declare that they have no conflict of interest.

Acknowledgements
Acknowledgements.

The authors acknowledge the NASA Goddard Space Flight Center, the U.S. Geological Survey (USGS), and the European Center for Medium-Range Weather Forecasts (ECMWF) for providing the images and datasets used in this study. We thank the anonymous reviewers and handling editor Christian Haas for their valuable comments, which helped improve our paper.

Financial support
Financial support.

This research has been supported by the National Natural Science Foundation of China (grant nos. 41876223, 41576188, and 41606215) and the National Key Research and Development Program of China (grant nos. 2016YFC1402704 and 2018YFC1407100). Jinlun Zhang was supported by the NOAA Climate Program Office (grant no. NA15OAR4310170).

Review statement
Review statement.

This paper was edited by Christian Haas and reviewed by two anonymous referees.

References

Alam, A. and Curry, J. A.: Determination of surface turbulent fluxes over leads in Arctic sea ice, J. Geophys. Res.-Oceans, 102, 3331–3343, 1997.

Alam, A. and Curry, J. A.: Evolution of new ice and turbulent fluxes over freezing winter leads, J. Geophys. Res.-Oceans, 103, 15783–15802, 1998.

Andreas, E. L. and Cash, B. A.: Convective heat transfer over wintertime leads and polynyas, J. Geophys. Res.-Oceans, 104, 25721–25734, 1999.

Andreas, E. L. and Murphy, B.: Bulk transfer coefficients for heat and momentum over leads and polynyas, J. Phys. Oceanogr., 16, 1875–1883, 1986.

Andreas, E. L., Paulson, C. A., William, R. M., Lindsay, R. W., and Businger, J. A.: The turbulent heat flux from Arctic leads, Bound.-Lay. Meteorol., 17, 57–91, 1979.

Arvidson, T., Gasch, J., and Goward, S. N.: Landsat 7's long-term acquisition plan – An innovative approach to building a global imagery archive, Remote Sens. Environ., 78, 13–26, 2001.

Barsi, J. A., Schott, J. R., Hook, S. J., Raqueno, N. G., Markham, B. L., and Radocinski, R. G.: Landsat-8 thermal infrared sensor (TIRS) vicarious radiometric calibration, Remote Sens., 6, 11607–11626, 2014.

Berrisford, P., Dee, D., Poli, P., Brugge, R., Fielding, K., Fuentes, M., Kallberg, P., Kobayashi, S., Uppala, S., and Simmons, A.: The ERA-Interim archive, version 2.0, ERA report series, 1, ECMWF, Shinfield Park, Reading, UK, 23 pp., available at: https://www.ecmwf.int/ (last access: 25 May 2019), 2011.

Bourassa, M., Vincent, D., and Wood, W.: A sea state parameterization with nonarbitrary wave age applicable to low and moderate wind speeds, J. Phys. Oceanogr., 31, 2840–2851, 2001

Canadian Ice Service: MANICE: Manual of standard procedures for observing and reporting ice conditions, Ice Services Branch, Atmospheric Environment Service, Ice Centre Environment Canada, Ottawa, Ontario, Canada, 2005.

Chechin, D. G. and Lüpkes, C.: Boundary-layer development and low-level baroclinicity during high-latitude cold-air outbreaks: A simple model, Bound.-Lay. Meteorol., 162, 91–116, 2017.

Du, C., Ren, H., Qin, Q., Meng, J., and Zhao, S.: A practical split-window algorithm for estimating land surface temperature from Landsat 8 data, Remote Sens., 7, 647–665, 2015.

Eppler, D. T. and Full, W. E.: Polynomial trend surface analysis applied to AVHRR images to improve definition of arctic leads, Remote Sens. Environ., 40, 197–218, 1992.

Esau, I. N.: Amplification of turbulent exchange over wide Arctic leads: Large-eddy simulation study, J. Geophys. Res.-Atmos., 112, D08109, https://doi.org/10.1029/2006JD007225, 2007.

Fett, R. W., Englebretson, R. E., and Burk, S. D.: Techniques for analyzing lead condition in visible, infrared and microwave satellite imagery, J. Geophys. Res.-Atmos., 102, 13657–13671, 1997.

Fetterer, F. M. and Holyer, R. J.: A Hough Transform Technique for Extracting Lead Features from Sea Ice Imagery, in: Geoscience and Remote Sensing Symposium, 1989, IGARSS'89, 12th Canadian Symposium on Remote Sensing, 1989 International, Vol. 2, 1125–1128, IEEE, 1989.

Fily, M. and Rothrock, D. A.: Opening and closing of sea ice leads: Digital measurements from synthetic aperture radar, J. Geophys. Res.-Oceans, 95, 789–796, 1990.

Gerace, A. and Montanaro, M.: Derivation and validation of the stray light correction algorithm for the thermal infrared sensor onboard Landsat 8, Remote Sens. Environ., 191, 246–257, 2017.

Goosse, H., Campin, J. M., Deleersnijder, E., Fichefet, T., Mathieu, P. P., Maqueda, M. M., and Tartinville, B.: Description of the CLIO model version 3.0., Institut d'Astronomie et de Géophysique Georges Lemaitre, Catholic University of Louvain, Belgium, 2001.

Hall, D. K. and Riggs, G.: MODIS/Terra Sea Ice Extent 5-Min L2 Swath 1 km, Version 6, Boulder, Colorado USA, NASA National Snow and Ice Data Center Distributed Active Archive Center, https://doi.org/10.5067/MODIS/MOD29.006, 2015.

Hall, D. K., Riggs, G. A., Salomonson, V. V., Barton, J. S., Casey, K., Chien, J. Y. L., DiGirolamo, N. E., Klein, A. G., Powell, H. W., and Tait, A. B.: Algorithm theoretical basis document (ATBD) for the MODIS snow and sea ice-mapping algorithms, NASA GSFC, 45, 2001.

Jenks, G. F.: Generalization in statistical mapping, Ann. Assoc. Am. Geogr., 53, 15–26, 1963.

Key, J. and Peckham, S.: Probable errors in width distributions of sea ice leads measured along a transect, J. Geophys. Res.-Oceans, 96, 18417–18423, 1991.

Key, J., Stone, R., Maslanik, J., and Ellefsen, E.: The detectability of sea-ice leads in satellite data as a function of atmospheric conditions and measurement scale, Ann. Glaciol., 17, 227–232, 1993.

Key, J., Maslanik, J. A., and Ellefsen, E.: The effects of sensor field-of-view on the geometrical characteristics of sea ice leads and implications for large-area heat flux estimates, Remote Sens. Environ., 48, 347–357, 1994.

Key, J. R., Collins, J. B., Fowler, C., and Stone, R. S.: High-latitude surface temperature estimates from thermal satellite data, Remote Sens. Environ., 61, 302–309, 1997.

Kwok, R.: Deformation of the Arctic ocean sea ice cover between November 1996 and April 1997: a qualitative survey, IUTAM symposium on scaling laws in ice mechanics and ice dynamics, Springer, Dordrecht, 315–322, 2001.

Kwok, R. and Rothrock, D. A.: Decline in Arctic sea ice thickness from submarine and ICESat records: 1958–2008, Geophys. Res. Lett., 36, L15501, https://doi.org/10.1029/2009GL039035, 2009.

Large, W. G. and Pond, S.: Open ocean momentum flux measurements in moderate to strong winds, J. Phys. Oceanogr., 11, 324–336, 1981.

Large, W. G. and Pond, S.: Sensible and latent heat flux measurements over the ocean, J. Phys. Oceanogr., 12, 464–482, 1982.

Lindsay, R. W. and Rothrock, D. A.: Arctic sea ice surface temperature from AVHRR, J. Climate, 7, 174–183, 1994.

Lindsay, R. W. and Rothrock, D. A.: Arctic sea ice leads from advanced very high resolution radiometer images, J. Geophys. Res.-Oceans, 100, 4533–4544, 1995.

Lüpkes, C., Vihma, T., Birnbaum, G., and Wacker, U.: Influence of leads in sea ice on the temperature of the atmospheric boundary layer during polar night, Geophys. Res. Lett., 35, L03805, https://doi.org/10.1029/2007GL032461, 2008.

Marcq, S. and Weiss, J.: Influence of sea ice lead-width distribution on turbulent heat transfer between the ocean and the atmosphere, The Cryosphere, 6, 143–156, https://doi.org/10.5194/tc-6-143-2012, 2012.

Marsan, D., Stern, H., Lindsay, R., and Weiss, J.: Scale dependence and localization of the deformation of Arctic sea ice, Phys. Rev. Lett., 93, 178501, 2004.

Maslanik, J. A. and Key, J.: On treatments of fetch and stability sensitivity in large-area estimates of sensible heat flux over sea ice, J. Geophys. Res.-Oceans, 100, 4573–4584, 1995.

Maykut, G. A.: Energy exchange over young sea ice in the central Arctic, J. Geophys. Res.-Oceans, 83, 3646–3658, 1978.

Maykut, G. A.: The surface heat and mass balance, The geophysics of sea ice, Springer, Boston, MA, 395–463, 1986.

Miles, M. W. and Barry, R. G.: A 5-year satellite climatology of winter sea ice leads in the western Arctic, J. Geophys. Res.-Oceans, 103, 21723–21734, 1998.

MCST (MODIS Characterization Support Team): MODIS 1km Calibrated Radiances Product. NASA MODIS Adaptive Processing System, Goddard Space Flight Center, USA, available at: https://ladsweb.modaps.eosdis.nasa.gov/ (last access: 25 May 2019), 2015.

Montanaro, M., Gerace, A., Lunsford, A., and Reuter, D.: Stray light artifacts in imagery from the Landsat 8 Thermal Infrared Sensor, Remote Sens., 6, 10435–10456, 2014a.

Montanaro, M., Levy, R., and Markham, B.: On-orbit radiometric performance of the Landsat 8 Thermal Infrared Sensor, Remote Sens., 6, 11753–11769, 2014b.

Montanaro, M., Lunsford, A., Tesfaye, Z., Wenny, B., and Reuter, D.: Radiometric calibration methodology of the Landsat 8 thermal infrared sensor, Remote Sens., 6, 8803–8821, 2014c.

Montanaro, M., Gerace, A., and Rohrbach, S.: Toward an operational stray light correction for the Landsat 8 Thermal Infrared Sensor, Appl. Opt., 54, 3963–3978, 2015.

Murashkin, D., Spreen, G., Huntemann, M., and Dierking, W.: Method for detection of leads from Sentinel-1 SAR images, Ann. Glaciol., 59, 124–136, 2018.

Oberhuber, J. M.: An atlas based on the “Coads” data set: The budgets of heat, buoyancy and turbulent kinetic energy at the surface of global ocean, Rep. 15. Max-Planck-Inst. for Meteorol., Hamburg, Germany, 199 pp., 1988.

Rampal, P., Weiss, J., and Marsan, D.: Positive trend in the mean speed and deformation rate of Arctic sea ice, 1979–2007, J. Geophys. Res.-Oceans, 114, C05013, https://doi.org/10.1029/2008JC005066, 2009.

Ray, M. L., Rogers, A. L., and McGowan, J. G.: Analysis of wind shear models and trends in different terrains, University of Massachusetts, Department of Mechanical and Industrial Engineering, Renewable Energy Research Laboratory, 2006.

Renfrew, I. A. and King, J. C.: A simple model of the convective internal boundary layer and its application to surface heat flux estimates within polynyas, Bound.-Lay. Meteorol., 94, 335–356, 2000.

Ridler, T. W. and Calvard, S.: Picture thresholding using an iterative selection method, IEEE Trans. Syst. Man. Cybern., 8, 630–632, 1978.

Röhrs, J. and Kaleschke, L.: An algorithm to detect sea ice leads by using AMSR-E passive microwave imagery, The Cryosphere, 6, 343–352, https://doi.org/10.5194/tc-6-343-2012, 2012.

Skoković, D., Sobrino, J. A., Jimenez-Munoz, J. C., Soria, G., Juşien, Y., Mattar, C., and Cristóbal, J.: Calibration and validation of land surface temperature for Landsat-8 TIRS sensor, in: LPVE (Land Product Validation and Evolution, ESA/ESRIN, 2014.

Spreen, G., Kwok, R., and Menemenlis, D.: Trends in Arctic sea ice drift and role of wind forcing: 1992–2009, Geophys. Res. Lett., 38, L19501, https://doi.org/10.1029/2011GL048970, 2011.

Stull, R. B.: An introduction to boundary layer meteorology (Vol. 13), Springer Science and Business Media, 1988.

Tetzlaff, A., Lüpkes, C., and Hartmann, J.: Aircraft-based observations of atmospheric boundary-layer modification over Arctic leads, Q. J. Roy. Meteorol. Soc., 141, 2839–2856, 2015.

Tschudi, M. A., Curry, J. A., and Maslanik, J. A.: Characterization of springtime leads in the Beaufort/Chukchi Seas from airborne and satellite observations during FIRE/SHEBA, J. Geophys. Res.-Oceans, 107, 8034, https://doi.org/10.1029/2000JC000541, 2002.

Venkatram, A.: A model of internal boundary-layer development, Bound.-Lay. Meteorol., 11, 419–437, 1977.

Wadhams, P.: Sea-ice topography of the Arctic Ocean in the region 70 W to 25 E, Phil. Trans. Roy. Soc. London A, 302, 45–85, 1981.

Wadhams, P., McLaren, A. S., and Weintraub, R.: Ice thickness distribution in Davis Strait in February from submarine sonar profiles, J. Geophys. Res.-Ocean, 90, 1069–1077, 1985.

Walter, B. A., Overland, J. E., and Turet, P.: A comparison of satellite-derived and aircraft-measured regional surface sensible heat fluxes over the Beaufort Sea, J. Geophys. Res.-Oceans, 100, 4585–4591, 1995.

Weiss, J. and Marsan, D.: Scale properties of sea ice deformation and fracturing, Compt. Rend. Phys., 5, 735–751, 2004.

Wernecke, A. and Kaleschke, L.: Lead detection in Arctic sea ice from CryoSat-2: quality assessment, lead area fraction and width distribution, The Cryosphere, 9, 1955–1968, https://doi.org/10.5194/tc-9-1955-2015, 2015.

Willmes, S. and Heinemann, G.: Pan-Arctic lead detection from MODIS thermal infrared imagery, Ann. Glaciol., 56, 29–37, 2015.

World Meteorological Organization.: WMO – No. 259: Sea-Ice Nomenclature, volumes I, II and III, Oostende, Belgium, 2014.

Zanter, K.: Landsat 8 data users handbook, LSDS-I574, Version 4.0, U.S. Geological Survey, Sioux Falls, South Dakota, USA, 115 pp., available at: http://glovis.usgs.gov/, last access: 25 May 2019.