Estimation of turbulent heat flux over leads using satellite thermal images

Sea ice leads are an important feature in pack ice in the Arctic. Even covered by thin ice, leads can still serve as prime windows for heat exchange between the atmosphere and the ocean, especially in the winter. Lead geometry and distribution in the Arctic have been studied using optical and microwave remote sensing data, but turbulent heat flux over leads has only been measured on-site during a few special expeditions. In this study, we derive turbulent heat flux through leads at different scales using a combination of surface temperature and lead distribution from remote sensing images and meteorological parameters from a reanalysis dataset. First, ice surface temperature (IST) was calculated from Landsat-8 Thermal Infrared Sensor (TIRS) and Moderate Resolution Imaging Spectroradiometer (MODIS) thermal images using a split-window algorithm; then, lead pixels were segmented from colder ice. Heat flux over leads was estimated using two empirical models: bulk aerodynamic formulae and a fetch-limited model with lead width from Landsat-8. Results show that even though the lead area from MODIS is a little larger, the length of leads is underestimated by 72.9 % in MODIS data compared to TIRS data due to the inability to resolve small leads. Heat flux estimated from Landsat-8 TIRS data using bulk formulae is 56.70 % larger than that from MODIS data. When the fetch-limited model was applied, turbulent heat flux calculated from TIRS data is 32.34 % higher than that from bulk formulae. In both cases, small leads accounted for more than a quarter of total heat flux over leads, mainly 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 air–ocean temperature differences and stronger winds.

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.

Since the late 1970s, remote sensing images obtained by satellite sensors, including optical, thermal, and microwave, have been used to detect sea ice leads in the Arctic (Fetterer and Holyer, 1989; Fily and Rothrock, 1990; Fett et al., 1997). Lindsay and Rothrock (1995) promoted the concept of potential open water for lead detection, which requires both temperature and albedo differences between ice surface pixels and open water tie points. Based on different emissivities of thin ice at two microwave frequencies available for the Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E), Röhrs and Kaleschke (2012) developed a retrieval algorithm to estimate Arctic lead concentration, similar to sea ice concentration. The algorithm could provide subpixel information on lead distribution, but the resolution is still too coarse to detect small leads prevailing in pack ice. Willmes and Heinemann (2015) mapped pan-Arctic lead distribution at 1 km resolution using the local temperature anomaly ΔT_s to identify leads from surrounding thick ice. Other remote sensing data, including altimetry, high-resolution optical, and synthetic-aperture radar (SAR) images, were also used to identify leads in limited areas due to constraints of cloud contamination and data acquisition restrictions (Key et al., 1993; Miles and Barry, 1998; Kwok, 2001; Weiss and Marsan, 2004; Wernecke and Kaleschke, 2015; Murashkin et al., 2018).

Regardless of spectral characteristics used for lead detection, the scale dependence of lead statistics was explored in a few studies (Key et al., 1994; Weiss and Marsan, 2004; Marsan et al., 2004). Key et al. (1994) studied the effects of the sensor's field of view (FOV) using degraded optical images from the Landsat Multi-Spectral Scanner (MSS). They suggested that the mean lead width expands, and the lead fraction drops as the pixel size builds up in gradually degraded images. Assuming higher heat flux over narrow leads than wider leads, estimated turbulent heat flux was reduced by 45 % as the FOV was degraded from 80 to 640 m, mainly due to reduced lead fraction.

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 C_N10 at 10 m height for turbulent heat in neutral stability, using the nondimensional fetch $-X/L$, where L is the Obukhov length. A maximum C_N10 of $\mathrm{1.8}\times {\mathrm{10}}^{-\mathrm{3}}$ was found at small fetch, and the value decrease to $\mathrm{1.0}\times {\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 km² 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 1 Satellite images and other data used in this study.

Figure 1 Location 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(\sec q-\mathrm{1}\right)\right],\end{array}$$

where T₃₁ and T₃₂ 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}{\displaystyle }{T}_{\mathrm{s}}& {\displaystyle }=b_{\mathrm{0}}+\left(b_{\mathrm{1}}+b_{\mathrm{2}}{\displaystyle \frac{\mathrm{1}-\mathit{\epsilon }}{\mathit{\epsilon }}}+b_{\mathrm{3}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\epsilon }}{{\mathit{\epsilon }}^{\mathrm{2}}}}\right){\displaystyle \frac{T_i+T_j}{\mathrm{2}}}\\ \text{(2)}& {\displaystyle }& {\displaystyle }+\left(b_{\mathrm{4}}+b_{\mathrm{5}}{\displaystyle \frac{\mathrm{1}-\mathit{\epsilon }}{\mathit{\epsilon }}}+b_{\mathrm{6}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\epsilon }}{{\mathit{\epsilon }}^{\mathrm{2}}}}\right){\displaystyle \frac{T_i-T_j}{\mathrm{2}}}+b_{\mathrm{7}}(T_i-T_j{)}^{\mathrm{2}},\end{array}$$

where T_i and T_j 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)$; b_k (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.

3.2 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 ΔT_s 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 (m₀) of the whole image, the iterative method seeks a threshold m_i which is the mean of averages m_A and m_B from two clusters divided by m_i: 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 (x₁, x₂). 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)$. Since we assign lead width to every pixel across the lead, length L_i for lead width X_i can be calculated as follows:

$$\begin{array}{}\text{(4)}& L_i={\displaystyle \frac{a_{\mathrm{0}}^{\mathrm{2}}{N}_i}{X_i}}={\displaystyle \frac{a_{\mathrm{0}}{N}_i}{i}},\end{array}$$

where a₀ is the pixel size; for TIRS, the value is 30 m, for MODIS, 1 km; and N_i is the number of pixels for width X_i=a₀i, (i=1, 2, 3...).

Figure 2 Detection 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 x₁ and x₂, respectively.

3.3 Heat flux model used for lead area

Turbulent heat flux between the Arctic Ocean and the atmosphere, including sensible (F_s) and latent (F_l) 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 (F_s) and latent heat flux (F_l) 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; c_p is the specific heat at constant pressure; L_v is the latent heat of vaporization; u_r, T_r, and Q_r are wind speed, air temperature, and specific humidity at reference height r=2 m, respectively; T_s is surface temperature; and Q_s 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}}={\displaystyle \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 e_s0 represents the saturated vapor pressure at surface temperature T_s:

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

with e₀ 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.

C_sh and C_le 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)}& {\displaystyle \frac{u_{\mathrm{10}}}{u_{\mathrm{r}}}}={\displaystyle \frac{\ln \mathrm{10}-\ln Z_{\mathrm{0}}}{\ln r-\ln Z_{\mathrm{0}}}},\end{array}$$

where u₁₀ and u_r are wind speed at 10 and 2 m height, and Z₀ 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⁻¹), 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(T_{\mathrm{s}}-T_{\mathrm{r}})/\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}}(Q_{\mathrm{s}}-Q_{\mathrm{r}})/\mathrm{\Delta }{z}_{\mathrm{Q}},\end{array}$$

where D and D_w are the molecular diffusivities of heat and water vapor in air, respectively, and Δz_T and Δz_Q 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({\displaystyle \frac{vD}{\mathrm{\Delta }B}}\right)}^{\mathrm{1}/\mathrm{3}}\end{array}$$

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

$$\begin{array}{}\text{(14)}& \mathrm{\Delta }B={\displaystyle \frac{\mathrm{g}}{\stackrel{\mathrm{‾}}{T}}}\left(\mathrm{\Delta }T+{\displaystyle \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 }={\displaystyle \frac{\mathrm{0.3}}{\mathrm{0.4}-h/L}}+\mathrm{0.15}\end{array}$$

$$\begin{array}{}\text{(16)}& h=\mathrm{0.82}\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({\displaystyle \frac{\mathrm{0.65}}{r}}+\mathrm{0.079}-\mathrm{0.0043}r\right)\ast R{i}_{\mathrm{b}},\end{array}$$

where Ri_b is the bulk Richardson number:

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

where u_r 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 3 Turbulent 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 4 Turbulent heat flux for each width at wind speed of 5 m s⁻¹. 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 5 Ice 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 6 Local median and noise image from TIRS IST: (a) along-track median temperature map and (b) noise image by detrending of median temperature map.

Figure 7 Correlation 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 8 Local 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 9 Binary lead maps from (a) MODIS and (b) Landsat-8 TIRS.

Lead area estimated from MODIS is 8074.0 km², which accounts for 8.25 % of the frame area, and from TIRS, 7376.2 km² 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 2 Retrieved 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 10 Width 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×10¹¹ W over a total area of 7376.2 km². This is 56.70 % larger than that from MODIS data (5.36×10¹¹ 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×10¹¹ 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×10¹⁰ W, 5.79 %).

4.3.2 The Andreas and Cash (1999) model

As we can see in Fig. 11 and Table 3, total heat flux over leads estimated by the Andreas and Cash (1999) model is 1.11×10¹² W, 32.34 % higher than that from bulk formulae, i.e., 8.40×10¹¹ W, among which 32.95 % of the difference comes from the small lead class. In both cases, small leads account for a quarter or more of total heat flux over all leads in both models, due to the large area, though the heat flux estimated using the fetch-limited model is 3.06×10¹¹ W, 41.39 % larger than the 2.16×10¹¹ W from bulk formulae. For comparison, the estimated heat fluxes from medium and large lead classes also increased by 38.95 % and 28.10 %, respectively, when the Andreas and Cash (1999) model was applied. However, the contribution of turbulent heat flux from large leads is reduced from 34.17 % to 32.68 %, while the contribution from small leads increased from 25.75 % to 27.50 %. Nonetheless, the fact that large leads with widths greater than 5 km account for 27.16 % of total lead area but contribute more than 32 % of total heat flux over leads is somehow contradictory to the fetch-limited theory.

Figure 11 Spatial 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 3 Estimated 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⁻¹ 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 ΔT_S. Since the anomaly maps from the two sensors have different means and standard deviations, mainly due to different definitions of neighborhood in calculating ΔT_S, 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 12 Distribution 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.

Table 4 Threshold candidates for lead detection and corresponding lead fractions.

5.2 Lead width

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)}={\displaystyle \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 N_T(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, Z₀∕Z_T, is assumed to be ∼e² to calculate Obukhov length L using the Richardson number Ri_b (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⁻¹ 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 13 Contribution 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 5 Turbulent 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×10¹² W at ΔT=10 K is about twice as large as that at ΔT=5 K (8.46×10¹¹ W).

When the Andreas and Cash (1999) model was applied, small leads were found to have a larger contribution at higher ΔT, 3.27×10¹¹ W (35.86 %) and 6.66×10¹¹ 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×10¹¹ W, 56.70 % larger than that from MODIS data (5.36×10¹¹ 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×10¹¹ 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×10¹⁰ W, 5.79 %).

The turbulent heat flux over leads estimated from the TIRS data by the Andreas and Cash (1999) model is 1.11×10¹² W, 32.34 % higher than that from bulk formulae (8.40×10¹¹ 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 A1 Leads 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)}& {\displaystyle }{c}_{\mathrm{sh}}=\mathrm{0.0327}{\displaystyle \frac{k}{\ln \left(r/z_{\mathrm{0}}\right)}}{\mathrm{\Phi }}_{\mathrm{sh}}=c_{\mathrm{shN}}{\mathrm{\Phi }}_{\mathrm{sh}}\end{array}$$

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

$$\begin{array}{}\text{(B3)}& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{sh}}={\displaystyle \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)}& {\displaystyle }& {\displaystyle }{\mathrm{\Phi }}_{\mathrm{le}}={\displaystyle \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)}& {\displaystyle }& {\displaystyle }\sqrt{{\displaystyle \frac{c_{\mathrm{M}}}{c_{\mathrm{MN}}}}}={\displaystyle \frac{\mathrm{1}}{\left(\mathrm{1}-\sqrt{c_{\mathrm{MN}}}{k}^{-\mathrm{1}}{\mathrm{\Psi }}_{\mathrm{M}}\right)}}\text{(B6)}& {\displaystyle }& {\displaystyle }{c}_{\mathrm{MN}}={\displaystyle \frac{k^{\mathrm{2}}}{{\left(\ln \left(\frac{r}{z_{\mathrm{0}}}\right)\right)}^{\mathrm{2}}}}\text{(B7)}& {\displaystyle }& {\displaystyle }{u}_{\ast }^{\mathrm{2}}=c_{\mathrm{M}}{u}_{\mathrm{r}}^{\mathrm{2}}\text{(B8)}& {\displaystyle }& {\displaystyle }{T}_{\mathrm{0}}=T_{\mathrm{r}}\left(\mathrm{1}+\mathrm{2.2}\times {\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)}& {\displaystyle }{z}_{\mathrm{0}}=\mathrm{0.032}{\displaystyle \frac{u_{\ast }^{\mathrm{2}}}{g}}.\end{array}$$

For unstable conditions,

$$\begin{array}{ll}\text{(B10)}& {\displaystyle }& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{H}}\left(A\right)={\mathrm{\Psi }}_{\mathrm{L}}\left(A\right)=\mathrm{2}\ln \left({\displaystyle \frac{\mathrm{1}+A^{\mathrm{2}}}{\mathrm{2}}}\right){\displaystyle }& {\displaystyle }{\mathrm{\Psi }}_{\mathrm{M}}\left(A\right)=\mathrm{2}\ln \left({\displaystyle \frac{\mathrm{1}+A}{\mathrm{2}}}\right)+\ln \left({\displaystyle \frac{\mathrm{1}+A^{\mathrm{2}}}{\mathrm{2}}}\right)\\ \text{(B11)}& {\displaystyle }& {\displaystyle }-\mathrm{2}\arctan A+{\displaystyle \frac{\mathit{\pi }}{\mathrm{2}}}\text{(B12)}& {\displaystyle }& {\displaystyle }A={\left(\mathrm{1}-\mathrm{16}\left(r/L\right)\right)}^{\mathrm{1}/\mathrm{4}}\text{(B13)}& {\displaystyle }& {\displaystyle }r/L={\displaystyle \frac{\mathrm{100}r}{T_{\mathrm{0}}{u}_{\mathrm{r}}^{\mathrm{2}}}}\left(\left(T_{\mathrm{s}}-T_{\mathrm{r}}\right)+\mathrm{2.2}\times {\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⁻².

4. bi: b_0∼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⁻³.

• Kinematic viscosity of air: $v=\mathrm{1.31}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹.

• Molecular diffusivities of heat in the air: $D=\mathrm{1.86}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹.

• Molecular diffusivities of water vapor in the air: $D_{\mathrm{w}}=\mathrm{2.14}\times {\mathrm{10}}^{-\mathrm{5}}$ m² s⁻¹.

• Specific heat at constant pressure: c_p=1004 J kg⁻¹ K⁻¹.

• Latent heat of vaporization or sublimation: $L_{\mathrm{water}}=\mathrm{2.51}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹, $L_{\mathrm{ice}}=\mathrm{2.86}\times {\mathrm{10}}^{\mathrm{6}}$ J kg⁻¹.

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

• Gravitational constant: g=9.8 m s⁻².

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

• Freezing point of seawater:

• $T_{\mathrm{S}\mathrm{0}}=\mathrm{273.15}-\mathrm{0.0137}-\mathrm{0.05199}$ S_w−0.00007225 $S_{\mathrm{w}}^{\mathrm{2}}=\mathrm{271.68}$ K.