Simultaneous estimation of wintertime sea ice thickness and snow depth from space-borne freeboard measurements

A method of simultaneously estimating snow depth and sea ice thickness using satellite-based freeboard measurements over the Arctic Ocean during winter was proposed. The ratio of snow depth to ice thickness (referred to as α) was defined and used in constraining the conversion from the freeboard to ice thickness in satellite altimetry without prior knowledge of snow depth. Then α was empirically determined using the ratio of temperature difference of the snow layer to the difference of the ice layer to allow the determination of α from satellite-derived snow surface temperature and snow–ice interface temperature. The proposed method was evaluated against NASA's Operation IceBridge measurements, and results indicated that the algorithm adequately retrieves snow depth and ice thickness simultaneously; retrieved ice thickness was found to be better than the methods relying on the use of snow depth climatology as input in terms of mean bias. The application of the proposed method to CryoSat-2 radar freeboard measurements yields similar results. In conclusion, the developed α-based method has the capacity to derive ice thickness and snow depth without relying on the snow depth information as input for the buoyancy equation or the radar penetration correction for converting freeboard to ice thickness.

1 Introduction

Satellite altimeters have been used to estimate sea ice thickness for nearly 2 decades (Laxon et al., 2003, 2013; Kwok et al., 2009). The altimeters do not measure sea ice thickness directly but measure the sea ice freeboard which is then converted to sea ice thickness with assumptions, for example, regarding the snow depth, snow–ice densities, and radar penetration (Ricker et al., 2014). We hereafter refer to this procedure as “freeboard to thickness conversion”.

Generally, there are two types of satellite altimeters measuring different sea ice freeboards. (1) Lidar altimeters such as NASA's ICESat (Zwally et al., 2002) and ICESat-2 (Markus et al., 2017) missions measure the total freeboard (F_t), which is the height from the sea surface in leads to the snow surface. (2) Radar altimeters such as ESA's CryoSat-2 (CS2) (Wingham et al., 2006) measure the radar freeboard (F_r), which is the difference in the radar ranging between the sea surface and the radar scattering horizon. By applying two correction terms regarding the wave propagation speed change in the snow layer (F_c) and displacement of the scattering horizon from the ice surface (F_p), the radar freeboard is converted to the ice freeboard (F_i), which is the height from the sea surface to the snow–ice interface (F_i). Several studies indicate that the radar scattering horizon is at or above the snow–ice interface depending on ice type and snow–ice conditions (Nandan et al., 2017; Armitage and Ridout, 2015; Willatt et al., 2011; Tonboe et al., 2010). However, the radar scattering horizon is often treated as the snow–ice interface (Kurtz et al., 2014; Kwok and Cunningham, 2015; Hendricks et al., 2016; Guerreiro et al., 2017, Tilling et al., 2018). The three different freeboards are indicated in Fig. 1.

Figure 1 Schematic diagram of a typical snow–ice system during the winter. Snow depth (h_s), ice thickness (H_i), total freeboard (F_t), radar freeboard (F_r), and ice freeboard (F_i) are indicated. Correction terms regarding the wave propagation speed change in the snow layer (F_c) and the displacement of the scattering horizon from the ice surface (F_p) are indicated by blue arrows. The red line denotes a typical temperature profile with air–snow interface temperature (T_as), snow–ice interface temperature (T_si), and ice–water interface temperature (T_iw).

For both lidar and radar altimeters, snow depth (h_s) is required as an input to constrain the freeboard to thickness conversion; thus, the conversion results are highly dependent on snow depth (Ricker et al., 2014; Zygmuntowska et al., 2014; Kern et al., 2015). The buoyancy equation used in the freeboard to thickness conversion describes the balance between buoyancy and the weight of snow and ice. For a given freeboard, snow–ice densities, and assumptions on radar penetration of the snow layer, sea ice thickness (H_i) is a function of h_s. According to Zygmuntowska et al. (2014), up to 70 % of uncertainty in the freeboard to thickness conversion stems from the poorly constrained snow depth. However, mapping the Arctic-scale snow depth distribution is challenging. The most commonly used snow depth information necessary for the freeboard to thickness conversion is the modified version of the snow depth climatology by Warren et al. (1999) (hereafter W99). W99 is based on in situ measurements from Soviet drifting stations (1954–1991) mostly of multiyear ice (MYI). Kurtz and Farrell (2011) compared W99 with Operation IceBridge (OIB) snow depth measurements in 2009 and claimed that W99 was still valid in the MYI region and significantly different from OIB snow depth on first-year ice (FYI). Based on that study, modified W99 (hereafter MW99) was developed, which halves W99 snow depth in regions covered by FYI. MW99 is often used in ESA's CryoSat-2 (CS2) ice thickness products available from the Center for Polar Observation and Modeling Data Portal (CPOM-UCL; Laxon et al., 2013), the Alfred Wegener Institute (AWI; Ricker et al., 2014), and the National Snow and Ice Data Center (NSIDC; Kurtz and Harbeck, 2017).

However, the use of MW99 for the freeboard to thickness conversion understandably yields a substantial error, considering that W99 is climatology and not actual snow depth. This is because the actual snow depth distribution is subject to the year-to-year variation of the snow–ice system; thus, the climatology based on the 37-year measurements of snow depth would deviate significantly from the actual distribution (Webster et al., 2014). Accordingly, such deviation causes errors in the estimation of ice thickness. Thus, additional snow observations covering both MYI and FYI at the Arctic basin scale would be ideal as a replacement of MW99.

There have been various approaches aimed at obtaining the snow depth distribution over the Arctic scale using satellite observations. Markus and Cavalieri (1998) developed an algorithm based on the brightness temperatures (TBs) of the Special Sensor Microwave/Imager (SSM/I) based on the negative correlation of the snow depth with the spectral gradient ratio between 18 and 37 GHz of vertically polarized TBs on the Antarctic FYI. Comiso et al. (2003) have updated the coefficients of this algorithm for the Advanced Microwave Scanning Radiometer for the Earth Observing System (AMSR-E). However, snow depth retrieval using this algorithm is relatively less accurate when the MYI fraction within the grid cell is significant (Brucker and Markus, 2013). Recently, Rostosky et al. (2018) suggested a new method: using the lower frequency pair of 7 and 19 GHz to overcome this limitation. Nonetheless, estimating the basin-scale snow depth distribution seems to be a difficult task.

There are other approaches involving the use of the lower frequency measurements at L-band. Using Soil Moisture Ocean Salinity (SMOS) measurements, Maaß et al. (2013) found that 1.4 GHz TB depends on the snow depth through the insulation effect of the snow layer, and they determined snow depth by matching TBs simulated with a radiative transfer model (RTM) with SMOS-measured TBs. Zhou et al. (2018) simultaneously estimated the sea ice thickness and snow depth by adding additional laser altimeter freeboard information, improving the Maaß et al. (2013) approach. However, both of these RTM-based approaches require a priori information on ice properties (e.g., temperature and salinity profiles).

Other satellite remote sensing approaches include the snow depth retrieval using dual-frequency altimetry (Guerreiro et al., 2016; Lawrence et al., 2018; Kwok and Markus, 2018), multilinear regression (Kilic et al., 2019), and a neural network approach (Braakmann-Folgmann and Donlon, 2019). In spite of promising results, the dual frequency altimetry method is available only for regions where two altimeters overlap with each other, reducing a great deal of spatial coverage. On the other hand, the regression/neural network methods based on AMSR-2 TBs are prone to the overfitting problem, limiting their applications to other microwave sensors.

Here, let us switch our point of view to solving the buoyancy equation instead of retrieving snow depth directly. Remember that there are two unknowns (snow depth and ice thickness) in the buoyancy equation for given snow–ice densities, freeboard, and assumptions on radar penetration of the snow layer. The attempt so far has been to add one constraint (snow depth information) to the buoyancy equation for solving ice thickness. However, if a particular relationship between two unknowns is available, it can be used to constrain the equation yielding both ice thickness and snow depth simultaneously.

To identify such a relationship, this study examines the vertical thermal structure within the snow–ice layers observed by drifting buoys. The vertical thermal structure of a snow–ice system in winter is rather simple; the temperature profile of the snow–ice system can be assumed to be piecewise linear, as illustrated in Fig. 1. Therefore, the temperatures at three interfaces can represent the thermal state of the snow–ice system fairly well: they are (1) air–snow interface temperature (T_as), (2) snow–ice interface temperature (T_si), and (3) ice–water interface temperature (T_iw). T_iw is assumed to be nearly constant at the freezing temperature of seawater (Maaß et al., 2013), implying that the two other interface temperatures (T_as and T_si) are sufficient to describe the thermal structure of the system.

Based on this thermal structure, there is a constraint relating the snow depth and ice thickness. In identifying this constraint, conductive heat flux is assumed to be continuous through the snow–ice interface (Maykut and Untersteiner, 1971), implying that conductive heat fluxes within the snow and ice layers are the same under the steady state assumed in the given thermal structure. As the conductive heat flux is proportional to the bulk temperature difference of the layer divided by its thickness, it is possible to deduce the relationship between snow depth and ice thickness from the given thermal structure.

Once the relationship is obtained, then it is possible to apply it at the Arctic Ocean basin scale because the thermal structure can be resolved from satellites, as shown in the recently available basin-scale and long-term satellite-derived interface temperatures (Dybkjær et al., 2020; Lee et al., 2018). In determining the snow depth along with the ice thickness, instead of using the snow depth as an input to solve for the ice thickness, we intend to (1) examine the relationship between the vertical thermal structure of a snow–ice system (T_as and T_si) and the thicknesses of the snow and ice layer (h_s and H_i) using buoy measurements and (2) retrieve the sea ice thickness and the snow depth simultaneously by applying their relationship to the freeboard to thickness conversion as a constraint, thus replacing the snow depth information. The result may reduce uncertainty in the freeboard to ice thickness conversion by replacing the currently used snow depth climatology.

2 Method

Here we provide the theoretical background of how the snow–ice thickness ratio (α=h_s ∕ H_i) can be related to T_as and T_si. Then, after empirically determining the relationship of α to T_as and T_si from buoy-measured temperature profiles, α obtained from satellite observed T_as and T_si is then used to constrain the conversion from freeboard to ice thickness over the Arctic Ocean during winter.

2.1 Theoretical background

We intend to find a relationship between snow depth and ice thickness in terms of the vertical thermal structure of the snow–ice system. Because the temperature gradients within the snow and ice layers are linked to both temperature and thickness, we focus on the temperature gradient. Owing to the physical reasoning that the conductive heat flux is continuous across the snow–ice interface (Maykut and Untersteiner, 1971), the following relationship is valid at the snow–ice interface:

$$\begin{array}{}\text{(1)}& {\left.k_{\text{snow}}{\displaystyle \frac{\partial T_{\text{snow}}}{\partial z}}\right|}_{z=\mathrm{0}}={\left.k_{\text{ice}}{\displaystyle \frac{\partial T_{\text{ice}}}{\partial z}}\right|}_{z=\mathrm{0}}.\end{array}$$

In Eq. (1), the subscripts snow and ice denote their respective layers, while T, k, and z denote temperature, thermal conductivity, and depth, respectively. The snow–ice interface is defined as z=0. Assuming a piecewise linear temperature profile within the snow–ice layer, Eq. (1) can be rewritten as follows:

$$\begin{array}{}\text{(2)}& k_{\text{snow}}{\displaystyle \frac{T_{\text{as}}-T_{\text{si}}}{h_{\text{s}}}}=k_{\text{ice}}{\displaystyle \frac{T_{\text{si}}-T_{\text{iw}}}{H_{\text{i}}}},\end{array}$$

where subscripts as, si, and iw denote the air–snow, snow–ice, and ice–water interface, respectively, and H_i and h_s denote the sea ice thickness and snow depth as in Fig. 1. Introducing a variable α, which is the snow–ice thickness ratio, Eq. (2) becomes the following:

$$\begin{array}{}\text{(3)}& \mathit{\alpha }={\displaystyle \frac{h_{\text{s}}}{H_{\text{i}}}}={\displaystyle \frac{k_{\text{snow}}}{k_{\text{ice}}}}{\displaystyle \frac{\mathrm{\Delta }{T}_{\text{snow}}}{\mathrm{\Delta }{T}_{\text{ice}}}}.\end{array}$$

Here, ΔT denotes the temperature difference between the top and bottom of each of the snow and ice layers (i.e., $\mathrm{\Delta }{T}_{\text{snow}}=T_{\text{as}}-T_{\text{si}}$, $\mathrm{\Delta }{T}_{\text{ice}}=T_{\text{si}}-T_{\text{iw}}$). As explained in detail in Sect. 2.3, α can be used to constrain the freeboard to thickness conversion. Thus, once α is known, both snow depth and ice thickness can be simultaneously estimated from altimeter-measured freeboard instead of using snow depth data for ice thickness retrieval.

2.2 Empirical determination of “α-prediction equation” from buoy measurements

To obtain α, the conductivity ratio (k_snow ∕ k_ice) should be known even if the temperature difference ratio (ΔT_snow ∕ ΔT_ice) is given. In this study, instead of using the conventional conductivity ratio found in the literature, it is empirically determined using buoy-measured α and ΔT_snow ∕ ΔT_ice. Thus, the interface should be defined and determined from buoy-measured temperature profiles, which show a piecewise linear temperature profile as shown in Fig. 1.

The buoy-measured temperature profiles at the vertical resolution of 10 cm are used in this study (Sect. 3.1). Although the instrument initially sets the zero-depth reference position to be approximately at the snow–ice interface, the reference position can deviate from the initial location if the ice deforms or if the snow refreezes after the temporary melt into snow ice. In addition, the interfaces (air–snow, snow–ice, and ice–water) may be located in between measurement levels in a 10 cm spacing. Therefore, an interface searching algorithm is developed to determine three interfaces (y_as, y_si, y_iw) and their respective temperatures (T_as, T_si, T_iw) by extrapolating each piecewise linear temperature profile iteratively.

Figure 2 The flow chart of the interface searching algorithm: y_i and T_i denote the position and temperature of a data point in the temperature profile, y_as,y_si, and y_iw denote the position of the interfaces, and T_layer denotes a set of temperature data points.

The interface searching algorithm iterates three processes to find the location and temperature of each interface: it (1) divides the temperature profile into four layers using the most recently available locations of the three interfaces, (2) finds a linear regression line of the temperature profile at each layer, and (3) updates the location and temperature of each interface by finding an intersection between two adjacent regression lines. The algorithm fails if the temperature profile is far from linear or if the thickness of a certain layer is too thin to have fewer than two data points. More detailed procedures for determining the interface are provided in Fig. 2 as a flow chart. The outputs are T_as, T_si, T_iw, H_i ($H_{\text{i}}=y_{\text{as}}-y_{\text{si}}$), and h_s ($h_{\text{s}}=y_{\text{si}}-y_{\text{iw}}$). Examples of the interface searching results for 15 d averaged temperature profiles are shown in Fig. 3. The algorithm works adequately for both Cold Regions Research and Engineering Laboratory Ice Mass Balance (CRREL-IMB; Fig. 3a–c) and Surface Heat Energy Budget of the Arctic (SHEBA; Fig. 3d–f) buoy data.

Figure 3 Examples of interface searching results with an averaging period of 15 d: (a) 2012G period: 2, (b) 2013F period: 8, (c) 2014G period: 1, (d) Q2 period: 6, (e) R4 period: 6, and (f) SEA period: 10. The period number indicates the sequential 15 d period from 1 November (e.g., period: 2 denotes a time-averaging period of 16 to 30 November). Blue dots are buoy-measured temperature profiles, and red lines are regression lines. Dashed black lines indicate the intersections between adjacent regression lines.

Since T_as, T_si, T_iw, H_i, and h_s can be obtained from the previous interface determination with buoy data, the calculation of ΔT_snow ∕ ΔT_ice and α is straightforward. An empirical relationship can then be obtained by relating α to ΔT_snow ∕ ΔT_ice by running a regression model, and details are given in Sect. 4. However, for the time being, we assume that the regression equation (referred to as an “α-prediction equation”, which will be discussed in Sect. 4) is used to predict α from ΔT_snow/ΔT_ice.

2.3 Simultaneous estimation of ice thickness and snow depth from satellite-based freeboard using α

In this section, we describe how α can be used to constrain the freeboard to thickness conversion. Based on the assumed hydrostatic balance, ice thickness can be obtained from satellite-borne total freeboard or ice freeboard as follows:

$$\begin{array}{}\text{(4)}& {\displaystyle }{H}_{\text{i}}={\displaystyle \frac{{\mathit{\rho }}_w}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}}}{F}_{\text{t}}-{\displaystyle \frac{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{s}}}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}}}{h}_{\text{s}},\text{(5)}& {\displaystyle }{H}_{\text{i}}={\displaystyle \frac{{\mathit{\rho }}_w}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}}}{F}_{\text{i}}+{\displaystyle \frac{{\mathit{\rho }}_{\text{s}}}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}}}{h}_{\text{s}}.\end{array}$$

Here, ρ_w, ρ_i, and ρ_s denote the bulk densities of the water, ice, and snow layers, respectively. Ice freeboard is obtained from radar freeboard by applying two correction terms regarding the change in the wave propagation speed in snow layer (F_c) and the displacement of the scattering horizon from the ice surface (F_p) (Kwok and Cunningham, 2015):

$$\begin{array}{}\text{(6)}& F_{\text{i}}=F_{\text{r}}+\left(F_{\text{c}}-F_{\text{p}}\right).\end{array}$$

The correction terms are expressed in the following equations (Armitage and Ridout, 2015; Kwok and Markus, 2018):

$$\begin{array}{}\text{(7)}& {\displaystyle }{F}_{\text{c}}=\left({\mathit{\eta }}_{\text{s}}-\mathrm{1}\right)f{h}_{\text{s}},\text{(8)}& {\displaystyle }{F}_{\text{p}}=\left(\mathrm{1}-f\right)h_{\text{s}}.\end{array}$$

Here, η_s denotes the refractive index of the snow layer and f denotes the radar penetration factor (Armitage and Ridout, 2015), which is the depth of the radar scattering horizon relative to the snow depth (e.g., f=1 if the radar scattering horizon is at the snow–ice interface and f=0 if the radar scattering horizon is at the air–snow interface). A combination of Eqs. (6)–(8) yields the following relationship:

$$\begin{array}{}\text{(9)}& F_{\text{i}}=F_{\text{r}}+\left(f{\mathit{\eta }}_{\text{s}}-\mathrm{1}\right)h_{\text{s}}.\end{array}$$

Ice freeboard in Eq. (5) can be substituted by radar freeboard and snow depth using Eq. (9), i.e.,

$$\begin{array}{}\text{(10)}& H_{\text{i}}={\displaystyle \frac{{\mathit{\rho }}_w}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}}}{F}_{\text{r}}+{\displaystyle \frac{\left(f{\mathit{\eta }}_{\text{s}}-\mathrm{1}\right){\mathit{\rho }}_w+{\mathit{\rho }}_{\text{s}}}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}}}{h}_{\text{s}}.\end{array}$$

According to Eq. (10), the ice thickness can be estimated from the radar freeboard and the snow depth. Note that Eq. (10) becomes equivalent to the equation for the total freeboard (Eq. 4 if f=0, i.e., if there is no radar penetration into the snow layer). With the use of α, defined in Eq. (3), Eqs. (4) and (10) become the following:

$$\begin{array}{}\text{(11)}& {\displaystyle }{H}_{\text{i}}={\displaystyle \frac{{\mathit{\rho }}_w}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}+\mathit{\alpha }\left({\mathit{\rho }}_w-{\mathit{\rho }}_{\text{s}}\right)}}{F}_{\text{t}},\text{(12)}& {\displaystyle }{H}_{\text{i}}={\displaystyle \frac{{\mathit{\rho }}_w}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}-\mathit{\alpha }\left\{\left(f{\mathit{\eta }}_{\text{s}}-\mathrm{1}\right){\mathit{\rho }}_w+{\mathit{\rho }}_{\text{s}}\right\}}}{F}_{\text{r}}.\end{array}$$

From Eqs. (3), (11), and (12), it is evident that the snow depth and ice thickness can be simultaneously estimated from the freeboards once α, ρ, f, and η_s are known.

In order to obtain α from satellite measurements of T_as and T_si, we need to calculate the temperature difference ratio (ΔT_snow ∕ ΔT_ice). For the calculation, T_iw is set to be −1.5 ^∘C. The freezing temperature of seawater is often assumed to be −1.8 ^∘C; however, the value of −1.5 ^∘C has been chosen based on the buoy observations. A sensitivity test indicated that the influence of a 0.3 ^∘C difference in the freezing temperature on α was negligible (e.g., approximately 1.2 % difference for typical interface temperatures of $T_{\text{as}}=-\mathrm{30}$ ^∘C and $T_{\text{si}}=-\mathrm{20}$ ^∘C). The α values are calculated only at the pixel whose monthly sea ice concentration (SIC) is greater than 95 % and rejected if T_as is warmer than T_si. The densities are prescribed with those used for OIB data processing: ρ_s, ρ_i, and ρ_w are 0.320, 0.915, and 1.024 g cm⁻³, respectively (Kurtz et al., 2013). Although ρ_s varies seasonally (Warren et al., 1999) and ρ_i is greater for MYI than FYI (Alexandrov et al., 2010), we use the same densities as those of OIB data because we intend to compare outputs against OIB data. In solving Eq. (12), cases showing negative ice thickness ($\mathit{\alpha }\ge {\mathit{\alpha }}_{\text{crit}}=\mathrm{0.291}$ for the given densities and radar penetration factor) are rejected. Radar penetration factor f is set to be 0.84 for CS2 (Armitage and Ridout, 2015), and η_s is parameterized as a function of the snow density, i.e., ${\mathit{\eta }}_{\text{s}}=(\mathrm{1}+\mathrm{0.51}{\mathit{\rho }}_{\text{s}}{)}^{\mathrm{1.5}}$ (Ulaby et al., 1986).

Before the Arctic-basin-scale retrieval, ice thickness is estimated from OIB total freeboard measurements using Eq. (11) and from OIB-derived radar freeboards (Sect. 3.3) using Eq. (12) and using satellite-derived α as a constraint. At the same time, the corresponding snow depth is derived by multiplying the obtained sea ice thickness and α (Eq. 3). Sea ice thicknesses are also calculated from Eqs. (4) and (10) using MW99 as the snow depth to examine how simultaneous retrievals compare with ice thickness estimation using MW99. To differentiate various outputs, obtained snow depth and ice thickness are expressed with nomenclature such as (constraint, freeboard source). For example, the snow depth estimated from satellite-derived α and OIB total freeboard is referred to as h_s (α^sat, $F_{\text{t}}^{\text{OIB}}$), and sea ice thickness from the MW99 and OIB radar freeboard is referred to as H_i ($h_{\text{s}}^{\text{MW99}}$, $F_{\text{r}}^{\text{OIB}}$). Finally, ice thickness and snow depth are estimated from CS2 radar freeboard (Sect. 3.4) over the Arctic Ocean.

3 Data

Here we provide detailed information on the datasets used for the development of the retrieval algorithm, evaluation, and application at the Arctic Ocean basin scale.

3.1 CRREL and SHEBA buoy data

To determine the empirical relationship between α and ΔT_snow ∕ ΔT_ice using Eq. (3), we need information regarding h_s, H_i, T_as, T_si, and T_iw (as depicted in Fig. 1). These are sourced from temperature profiles observed by buoys deployed over the Arctic as part of the Surface Heat Energy Budget of the Arctic (SHEBA) campaign (Perovich et al., 2007) and the Cold Regions Research and Engineering Laboratory Ice Mass Balance (CRREL-IMB) buoy program (Perovich et al., 2019). Those buoy observations are stored for further analysis if there are no missing records over the entire period ranging from November to March of the following year. Detailed information regarding ice type and initial snow–ice thickness at deployment locations is given in Table 1.

Table 1 Information on the measurement sites of buoys whose observations were used in this study.

^∗ The initial snow depth and ice thickness at the SHEBA sites are average values of all thickness gauge measurements in the corresponding site because there was one thermistor string but several thickness gauges at each measurement site.

Time averages of temperature profiles are used as input for the interface searching algorithm (described in Sect. 2.2) to meet the required near-equilibrium states (e.g., linear temperature profile). However, because of the possibility that the results are dependent on the averaging period, we examine the results using various averaging periods from 1 to 30 d.

3.2 Satellite-derived skin and interface temperatures

For applying the buoy-based α-prediction equation in retrieving the snow–ice thicknesses over the Arctic Ocean, satellite-derived T_as and T_si data are necessary. In this study, T_as is obtained from the Arctic and Antarctic ice Surface Temperatures from thermal Infrared satellite sensors version 2 (AASTI-v2) data (Dybkjær et al., 2020), and the monthly mean for the 1982–2015 period is obtained from daily products. AASTI-v2 T_as is derived from the CM SAF cLouds, Albedo and surface RAdiation dataset from AVHRR (Advanced Very High Resolution Radiometer) data edition 2 (CLARA-A2) dataset (Karlsson et al., 2017), which is based on the algorithm described in Dybkjær et al. (2018). Information on the validation of this product is found in Dybkjær and Eastwood (2016). It is available in a 0.25^∘ grid format; however, because other satellite datasets such as SIC are available in a 25 km Polar Stereographic SSM/I grid, AASTI-v2 data are re-gridded in the same 25 km grid format. This reformatted AASTI-v2 dataset is called “satellite skin temperature”.

T_si is obtained from snow–ice interface temperature (SIIT) produced by Lee et al. (2018) over 30 years (1988–2017) of winters (December to February) using SSM/I and the Special Sensor Microwave Imager/Sounder (SSMIS) homogenized TBs (Berg et al., 2018). The daily data are in the 25 km grid format. Lee et al. (2018) reported that the satellite-derived T_si is consistent with snow–ice interface temperatures observed by CRREL-IMB buoys, with a correlation coefficient, bias, and RMSE of 0.95, 0.15, and 1.48 K, respectively. In this study, we also produced T_si for March using the same algorithm as Lee et al. (2018) for evaluating results against OIB data, which are mostly collected during spring. Monthly composites are constructed by averaging daily data for grid cells where the data frequency is over 20 d. This product is called “satellite interface temperature”.

3.3 OIB data

In this study, OIB snow depth ($h_{\text{s}}^{\text{OIB}}$) and total freeboard ($F_{\text{t}}^{\text{OIB}}$) are used as a reference in the evaluation of snow depth and ice thickness retrieved from the developed algorithm. NASA's OIB is an aircraft mission, and it measures snow depth and total freeboard over the Arctic using snow radar, Digital Mapping System (DMS), and Airborne Topographic Mapper (ATM) (Kurtz et al., 2013). OIB ice thickness is derived from measured snow depth and total freeboard for the given snow and ice densities using Eq. (4). In this study, the OIB radar freeboard ($F_{\text{r}}^{\text{OIB}}$) is derived from $F_{\text{t}}^{\text{OIB}}$ and $h_{\text{s}}^{\text{OIB}}$ using the combined relationship of $F_{\text{i}}=F_{\text{t}}-h_{\text{s}}$ and Eq. (9) as follows:

$$\begin{array}{}\text{(13)}& F_{\text{r}}^{\text{OIB}}=F_{\text{t}}^{\text{OIB}}-h_{\text{s}}^{\text{OIB}}-\left(f{\mathit{\eta }}_{\text{s}}-\mathrm{1}\right)h_{\text{s}}^{\text{OIB}}.\end{array}$$

Because the main objective of using OIB data is to evaluate the relative performance of the simultaneous retrieval method when the method is applied to CS2 data, the radar penetration factor (f) for OIB data processing is also set to be 0.84. In the data processing chain, $h_{\text{s}}^{\text{OIB}}$ is removed if it is smaller than the given uncertainty level of the dataset (∼5.7 cm) or it is larger than the total freeboard $F_{\text{t}}^{\text{OIB}}$.

The 5 years of OIB data during the 2011–2015 period are utilized in this study. The level 4 dataset (Kurtz et al., 2015) during the 2011–2013 period and Quick Look dataset during the 2014–2015 period are obtained from the NSIDC website (see the Data Availability section). Because we use the November–March period for the buoy analysis, only March OIB data are considered for the evaluation. The OIB data are also reformatted into the 25 km grid format by averaging pixel-level OIB observations on the 25 km grid.

3.4 CS2 data

For examining the Arctic Ocean basin distribution of ice thickness and snow depth, CS2 freeboard measurement summary data are used (Kurtz and Harbeck, 2017). They are monthly mean composites of CS2 ice freeboard data in the 25 km Polar Stereographic SSM/I grid format covering the entire Arctic and available from September 2010. Detailed descriptions of the retracker algorithm used in this dataset are found in the study by Kurtz et al. (2014). The dataset also includes MW99 ($h_{\text{s}}^{\text{MW99}}$) and W99 snow density climatology used for producing the ice freeboard.

The CS2 ice freeboard data ($F_{\text{i}}^{\text{CS2}}$) distributed by NSIDC (Kurtz and Harbeck, 2017) assumed that the radar scattering horizon is at the snow–ice interface and applied a wave propagation speed correction. However, the correction was made using $h_{\text{s}}^{\text{MW99}}$ and W99 snow density climatology with an erroneous form of $h_{\text{c}}=(\mathrm{1}-{\mathit{\eta }}_{\text{s}}^{-\mathrm{1}}$) h_s instead of the proper form of $h_{\text{c}}=({\mathit{\eta }}_{\text{s}}-\mathrm{1}$) h_s (Mallett et al., 2020). Thus, at this point, it is straightforward to derive the CS2 radar freeboard by removing the correction term, as in the following equation:

$$\begin{array}{}\text{(14)}& F_{\text{r}}^{\text{CS2}}=F_{\text{i}}^{\text{CS2}}-\left(\mathrm{1}-{\mathit{\eta }}_{\text{s}}^{-\mathrm{1}}\right)h_{\text{s}}^{\text{MW99}}.\end{array}$$

Here, η_s was parameterized as a function of the snow density, i.e., ${\mathit{\eta }}_{\text{s}}=(\mathrm{1}+\mathrm{1.7}{\mathit{\rho }}_{\text{s}}+\mathrm{0.7}{\mathit{\rho }}_{\text{s}}^{\mathrm{2}}{)}^{\mathrm{0.5}}$ (Tiuri et al., 1984), and ρ_s is taken from the W99 climatology based on Kurtz and Harbeck (2017). Then CS2 ice thickness is reproduced from $F_{\text{r}}^{\text{CS2}}$ and $h_{\text{s}}^{\text{MW99}}$ by using Eq. (10) with the constant densities and the radar penetration factor described in Sect. 2.3. Those $h_{\text{s}}^{\text{MW99}}$ and $H_{\text{i}}(h_{\text{s}}^{\text{MW99}}$, $F_{\text{r}}^{\text{CS2}}$) values are used for comparison with results from our simultaneous method.

3.5 Sea ice concentration

Calculation of α is done for those pixels whose monthly SIC is greater than 95 % (as described in Sect. 2.3). To determine pixels that meet this SIC criterion, “bootstrap sea ice concentrations from Nimbus-7 SMMR [Scanning Multichannel Microwave Radiometer] and DMSP [Defense Meteorological Satellite Program] SSM/I-SSMIS version 3” produced by Comiso (2017) are used. This SIC dataset is provided in the 25 km Polar Stereographic SSM/I grid format.

Figure 4 (a) Scatterplots of the temperature difference ratio of the snow and ice layer (ΔT_snow ∕ ΔT_ice) and the snow–ice thickness ratio (α). Color denotes the collected year of buoy data. The red lines are the regression lines (defined in Eq. 15). (b) The scatter plot of observed and regressed α.

4 Results

4.1 The empirical relationship between α and ΔT_snow ∕ ΔT_ice

We examine variables (i.e., T_as, T_si, T_iw, H_i, and h_s) obtained from buoy observations by applying the interface searching algorithm. In the scatter plot of weekly averaged ΔT_snow ∕ ΔT_ice vs. α (Fig. 4a), it appears that α linearly increases with ΔT_snow ∕ ΔT_ice when the ratio is smaller than 1.8, but the linear slope becomes smaller when ΔT_snow ∕ ΔT_ice is larger than 1.8. This pattern of the slopes is found to be nearly invariant from year to year, as is observed in the different colors appearing in the entire range of ΔT_snow ∕ ΔT_ice in Fig. 4a. We also found that this slope pattern is of a consistent nature even for different datasets; two different datasets (red points for SHEBA and other points for CRREL) covering various ranges of ΔT_snow ∕ ΔT_ice show similar distributions along the two different slopes. Thus, the slope pattern is not due to different data sources or different data periods. Further analysis of the two slopes is found in Appendix A.

Taking such a two-slope pattern with ΔT_snow ∕ ΔT_ice into account, we introduce a piecewise linear function that may express the slope pattern, i.e.,

$$\begin{array}{}\text{(15)}& y=\left\{\begin{array}{c}{a}_{\mathrm{1}}x+b_{\mathrm{1}}x\le x_{\mathrm{0}}\\ a_{\mathrm{2}}x+b_{\mathrm{2}}x>x_{\mathrm{0}}\end{array}\right.,x_{\mathrm{0}}={\displaystyle \frac{b_{\mathrm{1}}-b_{\mathrm{2}}}{a_{\mathrm{2}}-a_{\mathrm{1}}}}.\end{array}$$

In Eq. (15), x and y correspond to ΔT_snow ∕ ΔT_ice and α, respectively, and x₀ is the point where the slope transition takes place. Applying Eq. (15) to data points from buoy-based variables, the regression coefficients (a₁, b₁, a₂, b₂) and transition point (x₀) are determined by minimizing the total variance – the obtained regression line is plotted in Fig. 4a. The α value is predicted using the determined regression equation (hereafter referred to as the α-prediction equation) and compared to the original α values to see how well the regression was performed. The comparison of α with predicted values in Fig. 4b shows that the regression equation is well fitted because of the zero bias and 91.9 % of explained variance.

Figure 5 (a) The regression coefficients (a₁, b₁, a₂, b₂) in Eq. (15) and (b) the error statistics of the regression with averaging periods from 1 to 30 d.

Table 2 Coefficients of the regression equation for averaging periods of 1, 7, 15, and 30 d; a₁, b₁, a₂, b₂, and x₀ are given in Eq. (15).

Although the slope pattern discussed with Eq. (15) and Fig. 4 is based on the weekly averages, the slope pattern seems to be consistent among the data averaging periods except for an averaging period shorter than 5 d. Regressions in the form of Eq. (15) are performed with buoy data averaged with different averaging periods to understand the slope pattern. The regression coefficients and transition points for the chosen averaging periods are examined, and results for four averaging periods are given in Table 2. Detailed information on the coefficients and associated statistics varying with the averaging period is given in Fig. 5. The positions of slope change (x₀) are located at approximately 1.8, delineating a nearly invariant slope pattern regardless of different data averaging periods. Figure 5a shows that coefficients do not vary much with different averaging periods, while coefficients of the first part of the regression line (a₁ and b₁, x≤x₀) vary less than those of the second part (a₂ and b₂, x>x₀). The regression equations show that the explained variance (R²) rises quickly when the averaging period is longer but levels off when data are averaged over a period that is longer than 7 d. The bias appears to be near zero over the various averaging periods. Thus, the regression performance is found to be comparable if data are averaged over a period that is longer than 1 week.

4.2 Evaluation against OIB estimates

According to the regression results, it is possible to estimate α from ΔT_snow ∕ ΔT_ice. Since ΔT_snow ∕ ΔT_ice can be calculated from the satellite skin and interface temperature (as described in Sect. 3.2), the corresponding α can be estimated from satellite measurements. Thus, we are able to simultaneously retrieve sea ice thickness and snow depth from altimeter-based freeboard measurements following Eqs. (11) and (12). We test and evaluate this simultaneous retrieval approach using OIB data. Accordingly, ice thickness and snow depth are simultaneously estimated from OIB freeboard measurements and evaluated against the OIB snow depth ($h_{\text{s}}^{\text{OIB}}$) and ice thickness ($H_{\text{i}}^{\text{OIB}}$).

To calculate α, a data averaging period must be selected. Considering that the monthly composite of satellite freeboard measurements is needed to retrieve snow–ice thickness at the Arctic basin scale, it seems appropriate to use the monthly averaging period to calculate the monthly α distribution. Thus, we use the monthly averaged satellite temperatures and the coefficients for the 30 d averaging period (Table 2) to calculate α.

We simultaneously retrieved H_i and h_s for March of each year during the 2011–2015 period from the reformatted OIB freeboard measurements (Sect. 3.3), together with satellite-derived α (α^sat). As expressed in Eqs. (11) and (12), two different ice thickness retrievals are possible depending on the use of the freeboard type (i.e., total freeboard F_t vs. radar freeboard F_r). Two accordingly associated retrievals of snow depth are available. Retrieved results of ice thickness (H_i) and snow depth (h_s) from the use of OIB total freeboard and radar freeboard are given in the first and second row of Fig. 6, respectively. Corresponding OIB measurements are given at the bottom of Fig. 6. The comparison between any snow–ice retrievals and OIB measurements appear to be consistent with each other for both snow depth and ice thickness in terms of magnitudes and distribution.

Figure 6 Simultaneously retrieved ice thickness and snow depth from OIB total and radar freeboards in March of the 2011–2015 period. Corresponding OIB products are at the bottom.

To compare the results quantitatively, scatterplots comparing retrievals against OIB measurements are made, along with statistics for the snow depth and ice thickness retrievals, in the top four panels of Fig. 7. The two top-left panels are derived from the use of OIB total freeboard ($F_{\text{t}}^{\text{OIB}}$), while the two top-right panels are derived from the OIB radar freeboard ($F_{\text{r}}^{\text{OIB}}$). The comparison is done only for pixels for which all four products (i.e., snow–ice thicknesses from two different freeboards) are available. This indicates that the snow depth from the total freeboard (top left) is fairly consistent with the OIB snow depth with a correlation coefficient of 0.73 and with a near-zero bias. The retrieved ice thickness from the total freeboard (middle left) appears to be consistent with OIB ice thickness with a correlation coefficient of 0.93 and a bias around 8.5 cm. The RMSEs for snow depth and ice thickness are 6.8 cm and 44.3 cm, respectively. Based on the comparison results, Eq. (15) obtained from buoy measurements can be successfully implemented with space-borne total freeboard measurements for the simultaneous retrieval of snow depth and ice thickness.

Figure 7 Scatter plots between OIB products and the simultaneously retrieved snow depth and ice thickness from OIB total and radar freeboards during the March 2011–2015 period. Corresponding ice thicknesses estimated from MW99 are in the third row. The red lines are linear regression lines.

Following Eq. (12), snow depth and ice thickness retrievals are made from the use of radar freeboard measurements, and results are presented in the two top-right panels in Fig. 7. On the one hand, the comparison of obtained ice thickness against the OIB ice thickness indicates that the retrieved ice thickness shows a similar quality as that retrieved from the total freeboard measurements. On the other hand, snow retrievals from the radar freeboard show more scattered features compared to snow retrieval results from the total freeboard. The more scattered features found in the snow depth from the radar freeboard are likely due to the greater sensitivity of the retrieved α and the prescribed densities, as noted in Eq. (12). Note that Eq. (12) has a smaller denominator than that of Eq. (11). Results of the associated sensitivity analysis can be found in Appendix B.

We now examine how the use of MW99 for retrieving sea ice thickness from ICESat and CS2 measurements compares with results from our simultaneous method. To do so, OIB-measured total freeboard and radar freeboard are converted into ice thickness using MW99 as the input to solve Eqs. (4) and (10). In this study, these two ice thickness retrievals with the use of MW99 are referred to as “ICESat-like” thickness and “CS2-like” thickness, respectively, and their comparisons are now observed in the two panels at the bottom of Fig. 7. According to our analysis, ICESat-like thickness tends to underestimate the ice thickness by about 47.9 cm when MW99 is used in comparison to OIB thickness, and CS2-like ice thickness shows an overestimate of about 25.5 cm. Nevertheless, their correlation coefficients and RMSEs are similar to the results obtained from the α method.

The better agreement of H_i from the simultaneous method with $H_{\text{i}}^{\text{OIB}}$ may be due to the fact that the simultaneously estimated h_s is more consistent with $h_{\text{s}}^{\text{OIB}}$ ($h_{\text{s}}^{\text{MW99}}$ is likely larger than $h_{\text{s}}^{\text{OIB}}$, as shown in Fig. S1). Note that all inputs are the same except the snow depth. The negative bias of ICESat-like thickness and positive bias of CS2-like thickness reflect expected responses in different signs to the same snow depth error, as shown in different signs in the last terms of Eqs. (4) and (10) (also note Eq. B2 in Appendix B). Because of this reasoning, if there are decreasing trends not only in ice thickness but also in snow depth, the decreasing trend of ice thickness estimated from the constant snow depth will be diminished in radar while being amplified in lidar. Because of this, the construction of the ice thickness (or volume) trend from the two different satellite altimeters would be problematic if MW99 is used for the freeboard to thickness conversion. For example, it would be hard to compare the sea ice thickness records estimated from ICESat and CS2 observations and to extend the current ice thickness record from CS2 with NASA's recently launched ICESat-2, which carries a lidar altimeter, for the same reason.

4.3 Simultaneous retrieval of ice thickness and snow depth from CS2 measurements

We have demonstrated that the method of simultaneously retrieving the sea ice thickness and snow depth was successfully implemented with OIB measurements. Now we extend the proposed approach to satellite freeboard measurements. Here the method is tested with CS2 freeboard measurements, solving for H_i in Eq. (12), and α is obtained from the collocated satellite skin and interface temperature data.

Figure 8 Geographical distributions of observed CS2 radar freeboard (F_r) and estimated snow–ice thickness ratio (α), ice thickness (H_i), and snow depth (h_s) from December 2013 to March 2014. Gray areas in the second row denote where α retrievals failed because T_as was warmer than T_si.

Monthly means of CS2-estimated freeboard (F_r), retrieved α, ice thickness (H_i), and snow depth (h_s) for December 2013 to March 2014 are given in Fig. 8. The geographical distribution of α indicates that α is largest in January and becomes smaller during the following months. Geographically, there seems to be no particular distribution of α between months, although interestingly the lowest α values are always found over the north of the Canadian Archipelago, and the western part of the Arctic Ocean shows α values that are generally larger than those over the eastern part.

Retrieved ice thickness from the CS2 freeboard (F_r) using obtained α is presented in the third row of Fig. 8. As expected and as noted in Eq. (12), H_i shows a similar geographical distribution to radar freeboard (the first row). The thickest area is located north of the Canadian Archipelago where the ice appears thicker than 4 m. On the other hand, most of the FYI thickness appears to range from 1.0 to 2.0 m. The snow depth h_s is obtained by multiplying α by H_i (in 2nd and 3rd rows), following Eq. (3), and results are shown in the bottom row. The obtained snow distribution indicates that thicker (thinner) snow areas are generally coincident with thicker MYI (thinner FYI) areas. Such a similarity should be consistent with the notion that MYI should accumulate more precipitation than FYI because of its longer existence.

To assess the accuracy of CS2 retrievals, reference snow depth and ice thickness collocated with CS2 freeboard in space and time are necessary. However, different from simultaneous retrievals from OIB freeboards in Sect. 4.2, the evaluation with the required matching data may not be possible from the monthly composite of CS2 data used in this study. Here, instead of using monthly collocated match-up data, an indirect way is used to examine the accuracy of CS2 retrievals. We do so by examining whether the relationship between the simultaneous method and the MW99 method, based on retrievals from the OIB freeboard, can be reproduced by CS2-based retrievals. If similar results are obtained, respective accuracies can be deduced against those noted from the evaluation against OIB measurements.

Figure 9 Comparison of simultaneously retrieved snow depth and ice thickness to those from the MW99 method. (a) Snow depth from OIB radar freeboard, (b) snow depth from CS2 radar freeboard, (c) ice thickness from OIB radar freeboard, and (d) ice thickness from CS2 radar freeboard.

The relationships which can be obtained from the analysis in Sect. 4.2 – i.e., h_s (α^sat, $F_{\text{r}}^{\text{OIB}}$) vs. $h_{\text{s}}^{\text{MW99}}$ and H_i(α^sat, $F_{\text{r}}^{\text{OIB}}$) vs. H_i ($h_{\text{s}}^{\text{MW99}}$, $F_{\text{r}}^{\text{OIB}}$) – are compared with the relationships found in the current results in Fig. 8 – i.e., h_s (α^sat, $F_{\text{r}}^{\text{CS2}}$) vs. $h_{\text{s}}^{\text{MW99}}$ and H_i(α^sat, $F_{\text{r}}^{\text{CS2}}$) vs. H_i ($h_{\text{s}}^{\text{MW99}}$, $F_{\text{r}}^{\text{CS2}}$); the results are presented in Fig. 9. Observably, the relationships from CS2 freeboard data (Fig. 9b, d) are very similar to the relationship obtained from the comparison results from OIB measurements (Fig. 9a, c). This similarity of the slope strongly indicates that the CS2-based sea ice thickness from the current α method has similar accuracy to that found in the evaluation against OIB measurements (Sect. 4.2). Further uncertainty estimates for CS2-derived products can be found in Appendix C.

5 Conclusions and discussion

A new approach towards simultaneously estimating snow depth and ice thickness from space-borne freeboard measurements was proposed and tested using OIB data and CS2 freeboard measurements. In developing the algorithm, the vertical temperature slopes were assumed to be linear within the snow and ice layers so that a continuous heat flux could be maintained in both layers. This assumption allowed for the description of the snow–ice vertical thermal structure with snow skin temperature, snow–ice interface temperature, the water temperature at the ice–water interface, snow depth, and ice thickness. Based on the continuous heat transfer assumption, the snow–ice thickness ratio (α=h_s ∕ H_i) was introduced and could then be embedded into the freeboard to ice thickness conversion equations. Thus, information on both ice thickness and snow depth can be derived once α is known in the case of the availability of a freeboard without relying on the snow depth information as an input for the conversion from freeboard to ice thickness. From the drifting buoy measurements of the temperature profile, snow depth, and ice thickness over the Arctic Ocean, we demonstrated that α can be reliably determined using the ratio of the vertical difference of the snow layer temperature to the vertical difference of ice layer temperature (ΔT_snow ∕ ΔT_ice). An empirical regression equation was obtained for predicting α from three interface temperatures.

Before applying the α-prediction equation to simultaneously retrieve the ice thickness and snow depth from satellite-borne freeboard measurements, the algorithm was evaluated using OIB measurements in conjunction with satellite-derived snow skin temperature and snow–ice interface temperature. The evaluation of results demonstrated that our proposed algorithm adequately retrieved both parameters simultaneously. As a matter of fact, the ice thickness results were more accurate than they were from the current retrieval methods relying on the input of snow depth (this time MW99 snow climatology) in terms of mean bias. It should be noted that in this case, snow depth is a retrieval product instead of being input for the freeboard to ice thickness conversion adopted by CS2 or ICESat retrieval. The application was finally made for the retrieval of the snow depth and ice thickness from CS2 radar freeboard measurements from December 2013 to March 2014 using α as a constraint. Results showed that the quality of the obtained ice thickness was similar to that obtained from evaluation results against OIB measurements. Retrieved snow depth distributions were also found to be consistent with expectations.

In the retrieval process, we may be concerned about the applicability of the algorithm developed with buoy observations representing the point measurements to the larger spatial and temporal scales of satellite measurements. This concern may be relevant upon observing the range of α values. The α value in the satellite's monthly and 25 km × 25 km spatial scales was found to be generally smaller than 0.2. The smaller range of α compared to that shown in the buoy analysis results is likely due to the scale differences, indicating that extreme α values often shown in buoy measurements (due to very thick snow and/or very thin ice) may never be observed in satellite measurements. However, the range may not be a problem because the relationship (Eq. 3) expresses the thermal equilibrium condition described by the temperature at three interfaces, the ratio of snow and ice thickness, and the ratio of thermal conductivity between snow and ice. Considering that the algorithm is based on the equilibrium conditions, results should be valid regardless of spatial and temporal scales if the prerequisite equilibrium conditions are met. Apparently, buoy observations contain so many different cases that equilibrium conditions are met with different thermal and physical conditions of the snow–ice system. The sound evaluation results and the consistency between OIB and CS2 ice thickness retrieval results, which are subject to different scales, all suggest that the point-measured α-prediction equation can apply to satellite measurements.

Overall, the developed α-based method yields ice thickness and snow depth without relying on a priori “uncertain” snow depth information (MW99), which results in uncertainty in the ice thickness retrieval. The proposed method applies to both lidar and radar altimeter data, although lidar-based altimeter data tend to offer relatively more suitable snow depth information with a smaller RMSE. We expect to continuously monitor the Arctic-scale snow depth and ice thickness by applying the proposed α method to total freeboard observations from the recently launched ICESat-2 and using temperature observations from the upcoming MetOp-SG meteorological imager (MetImage), the microwave imager (MWI), and the proposed Copernicus Imaging Microwave Radiometer (CIMR).

Appendix A: Physical interpretation of the piecewise linearity between α and ΔT_snow ∕ ΔT_ice

The relationship found between α and ΔT_snow ∕ ΔT_ice showed a piecewise linearity which is almost invariant to the data averaging period. Because the slope change is attributable neither to different data sources nor to different data periods, it is likely caused by the physical properties of the snow and ice, as shown in Fig. A1. If the slope change is caused by the snow–ice condition, there will be a significant difference in snow–ice properties between the two parts showing different slopes. Here we examine the possibility of different physical properties causing the difference in slopes. Through this comparison using buoy data, we may identify important properties that might be responsible for the piecewise linearity.

Figure A1 Distribution of physical variables on scatterplots of the temperature difference ratio of snow and ice layer (ΔT_snow ∕ ΔT_ice) and the snow–ice thickness ratio (α). Color denotes the value of physical variables: (a) ice thickness (H_i), (b) snow depth (h_s), (c) air–snow interface temperature (T_as), (d) snow–ice interface temperature (T_si), (e) temperature difference within the snow layer ($\left|\mathrm{\Delta }{T}_{\text{snow}}\right|$), and (f) temperature difference within the ice layer ($\left|\mathrm{\Delta }{T}_{\text{ice}}\right|$).

Figure A2 Histogram of estimated (a, b) k_ice and (c, d) k_snow. The top and bottom rows denote the first and the second parts, respectively. The size of the bins is 0.05 W K⁻¹ m⁻¹.

First, the averages of basic properties available from buoy measurements are compared. They include ice thickness, snow depth, snow–ice interface temperature, ice temperature – $T_{\text{ice}}=(T_{\text{as}}+T_{\text{si}})/\mathrm{2}$ – and so on. The comparison revealed that snow–ice system within the first part (x≤x₀) is found to consist of relatively thicker ice (mean value: 1.84 m), thinner snow (0.29 m), and colder ice (−9.13 ^∘C), while the second part (x >x₀) is found to consist of relatively thinner ice (1.10 m), thicker snow (0.46 m), and warmer ice (−5.00 ^∘C). In general, a thicker snow or ice layer exhibits a greater temperature difference from the top to the bottom of the layer. There is no significant difference between the air–snow interface temperature (T_as) in the two slope parts.

The thermal conductivities, k_snow and k_ice, are also compared because what connects α and ΔT_snow ∕ ΔT_ice is the ratio of thermal conductivities. Before showing the results, we describe how to calculate k_snow and k_ice. First, the thermal conductivity ratio is calculated from buoy-measured variables (i.e., T_as, T_si, T_iw, h_s, and H_i) using Eq. (3). Because the underlying physics in k_snow are significantly more complex, k_ice is estimated first, and then k_snow is obtained by multiplying the calculated k_ice by k_snow ∕ k_ice. To calculate k_ice, the parameterization of Maykut and Untersteiner (1971), which describes k_ice as a function of salinity and temperature, is used:

$$\begin{array}{}\text{(A1)}& k_{\text{ice}}=\mathrm{2.03}+\mathrm{0.117}{\displaystyle \frac{S_{\text{ice}}}{T_{\text{ice}}}}.\end{array}$$

Here, S_ice and T_ice are the salinity (in parts per thousand, ppt) and temperature (in Celsius) of sea ice, respectively. For the calculation, S_ice is estimated according to the empirical relationship between sea ice thickness and mean salinity from Cox and Weeks (1974) as follows:

$$\begin{array}{}\text{(A2)}& S_{\text{ice}}=\left\{\begin{array}{c}\mathrm{14.24}-\mathrm{19.39}{H}_{\text{i}},H_{\text{i}}\le \mathrm{0.4}\text{m}\\ \mathrm{7.88}-\mathrm{1.59}{H}_{\text{i}},H_{\text{i}}>\mathrm{0.4}\text{m}\end{array}\right..\end{array}$$

Although Trodahl et al. (2001) reported that k_ice depends on depth and temperature, here we do not estimate accurate thermal conductivities but attempt to examine the physical consequences of the total ice layer.

The calculated thermal conductivities are presented in Fig. A2. The calculated k_ice ranges from 1.8 to 2.0 W K⁻¹ m⁻¹ (two left panels in Fig. A2). These values are consistent with the in situ measurements by Pringle et al. (2006). The mean values of k_ice of the first part (1.96 W K⁻¹ m⁻¹) and the second part (1.88 W K⁻¹ m⁻¹) show almost no difference. The calculated k_snow ranges from 0.2 to 1.05 W K⁻¹ m⁻¹ (two right panels in Fig. A2). This range is consistent with reported values in Sturm et al. (1997). The first part shows the greater spread in the distribution of k_snow compared to the second part. The mean k_snow values are 0.44 and 0.27 for the first part and second part, respectively.

As a significant difference is observed in k_snow, we would like to find a possible reason for this difference. To do so, we should first review the factors determining k_snow; they are density, temperature, and crystal structure (Sturm et al., 1997). Snow is a mixture of ice particles and air, and air has lower thermal conductivity than ice. Thus, snow with a relatively lower density including a greater portion of air should have relatively lower thermal conductivity. Besides, the thermal conductivity of ice particles depends on the temperature, and the path of heat transfer depends on the crystal structure which describes how the particles are connected. The heat transfer occurs not only by conduction but also by water vapor latent heat transportation and convection through the pore spaces (Sturm et al., 2002), which are hard to quantify explicitly. These two factors are closely related to the temperature gradient (or difference) imposed within the snow layer.

Based on this knowledge, we can infer the condition of the snow layer of the two parts. The relatively higher and varying k_snow of the first part would be related to the compaction process resulting in high density and metamorphic diversity, which changes the crystal structure. According to Sturm et al. (2002), the value of k_snow of a hard wind slab is up to 0.5 W m⁻¹ K⁻¹, while that of k_snow of depth hoar is below 0.1 W m⁻¹ K⁻¹. On the other hand, the lower and nearly constant k_snow of the second part implies that the snow layer of the second part would consist of fresh and dry snow having relatively lower density and a relatively lower likelihood of experiencing particular metamorphism.

In summary, it is concluded that the physical properties of snow and ice can account for the piecewise linearity based on the differences in the physical properties between the first and second parts. Especially the thermal conductivity of the snow, k_snow, seems to play an important role. Nevertheless, further analysis is required to fully understand this phenomenon.

Appendix B: Sensitivity test for the proposed method

Here we present results of a sensitivity test for showing how the snow depth and ice thickness retrieval results are dependent on the uncertainties in α. To do so, the uncertainty in the snow depth (Δh_s) due to the α error (i.e., Δα) and associated ice thickness error (ΔH_i) are estimated. From this sensitivity test, we expect to understand why the simultaneous method for the radar freeboard shows more scattered features than those from the lidar total freeboard.

First, Δh_s is defined by the difference of retrieved h_s with error (α+Δα) and without error (α):

$$\begin{array}{}\text{(B1)}& \mathrm{\Delta }{h}_{\text{s}}=\left\{\begin{array}{l}{h}_{\text{s}}\left(\mathit{\alpha }+\mathrm{\Delta }\mathit{\alpha },F_{\text{t}}\right)-h_{\text{s}}\left(\mathit{\alpha },F_{\text{t}}\right)\left(\mathrm{using}{F}_{\text{t}}\right)\\ h_{\text{s}}\left(\mathit{\alpha }+\mathrm{\Delta }\mathit{\alpha },F_{\text{r}}\right)-h_{\text{s}}\left(\mathit{\alpha },F_{\text{r}}\right)\left(\mathrm{using}{F}_{\text{r}}\right)\end{array}\right..\end{array}$$

Then Δh_s can be converted to the error in the ice thickness (ΔH_i) using the following equation derived from Eq. (10):

$$\begin{array}{}\text{(B2)}& \mathrm{\Delta }{H}_{\text{i}}={\displaystyle \frac{\left(f{\mathit{\eta }}_{\text{s}}-\mathrm{1}\right){\mathit{\rho }}_w+{\mathit{\rho }}_{\text{s}}}{{\mathit{\rho }}_w-{\mathit{\rho }}_{\text{i}}}}\mathrm{\Delta }{h}_{\text{s}}=\left\{\begin{array}{c}-\mathrm{6.46}\mathrm{\Delta }{h}_{\text{s}}\left(\mathrm{using}{F}_{\text{t}}\right)\\ \mathrm{3.44}\mathrm{\Delta }{h}_{\text{s}}\left(\mathrm{using}{F}_{\text{r}}\right)\end{array}\right..\end{array}$$

Because H_i and h_s are the combination of freeboard and α, as in Eqs. (3), (11), and (12), we only examine the uncertainty with some typical sea ice types. Here, physical states for thicker ice (type A), moderate ice (type B), and thinner ice (type C) are chosen, which are summarized in Table B1. Typical values for those three types are shown in the scatterplots of OIB-based α^OIB vs. $F_{\text{t}}^{\text{OIB}}$ and of satellite-based α^sat vs. $F_{\text{r}}^{\text{CS2}}$ – Fig. B1. It is shown that the majority of data points are located around type B, followed by type A. There seems to be a very small portion of total samples showing values around type C.

With $\mathrm{\Delta }\mathit{\alpha }=\pm \mathrm{0.05}$, which is the root mean square difference (RMSD) value between α^OIB and α^sat, Δh_s and ΔH_i are estimated for three ice types. Table B2 summarizes the results and shows that $\left|\mathrm{\Delta }{h}_{\text{s}}\right|$ is within 8 cm and that it tends to decrease as the ice becomes thinner when the current method is applied to the total freeboard. On the other hand, the use of radar freeboard shows that $\left|\mathrm{\Delta }{h}_{\text{s}}\right|$ tends to be more sensitive for the same Δα. Especially the sensitivity of type C is the greatest. This is because the denominator of Eq. (12) becomes smaller when α approaches α_crit, resulting in an unstable solution. For the ice thickness, $\left|\mathrm{\Delta }{H}_{\text{i}}\right|$ is smaller when the total freeboard is used since ΔH_i is proportional to Δh_s. However, the gap between the results from the two freeboards has narrowed because H_i from the total freeboard is more sensitive than the radar freeboard to Δh_s, according to Eq. (B2). The sensitivity characteristics shown here are consistent with the analysis results given in Sect. 4.2. Because there is a much smaller number of data points belonging to type C, at least in the data used for this study, the overall sensitivity would likely be in between types B and A.

It is also of importance to ask what degree of retrievals was yielded successfully. In this study, cases showing T_as>T_si or retrieved α≥α_crit are considered to be failures. Statistics on success/fail ratio of α retrieval for December–March of the 2011–2015 period are provided in Table B3. Overall, the success ratio was over 82 % in December−February, while it was reduced to ∼74 % in March. Most of the failures appear to be associated with cases showing the temperature inversion (i.e., T_as>T_si), whose areas are shaded with gray in the α distributions of Fig. 8. Those failure areas are generally found around the marginal ice zones and in the east of Greenland. On the other hand, there was a near-zero failure (0.02 % of total pixels) for retrieved α≥α_crit. This near-zero failure implies that almost all calculated α values meet the satisfactory condition after the removal of cases showing the temperature inversion. It may be concluded that the calculated α appears to be physically reasonable (i.e., α<α_crit) as long as presumed thermodynamic conditions are met.

Table B1 The physical state of typical cases of points A, B, and C.

Table B2 Errors of snow depth (Δh_s) and ice thickness (ΔH_i) for snow depth to ice thickness ratio error (Δα) of ±0.05.

^∗ Retrieval fail occurs if $\mathit{\alpha }+\mathrm{\Delta }\mathit{\alpha }>{\mathit{\alpha }}_{\text{crit}}$ (α_crit=0.291 for ρ_s=320 kg m⁻³, ρ_I=915 kg m⁻³, ρ_w=1024 kg m⁻³, and f=0.84).

Figure B1 Locations of physical states for typical types (A, B, C) on the freeboard-thickness ratio space. Blue dots are from (a) OIB data and (b) retrieved thickness ratio and CS2 radar freeboard.

Table B3 Statistics of success/fail ratios of α retrieval for winter 2011–2015.

α_crit=0.291 for ρ_s=320 kg m⁻³, ρ_i=915 kg m⁻³, ρ_w=1024 kg m⁻³, and f=0.84.

Appendix C: Uncertainty estimation for CS2 retrievals

Although the sensitivity test regarding uncertainty of satellite-derived α has been conducted in Appendix B, the uncertainty of CS2 freeboard measurements and prescribed parameters should be considered as well for the satellite snow depth and ice thickness estimates. To do so, a simple propagation analysis of errors is performed, regarding the uncertainty of satellite products (α^sat and $F_{\text{r}}^{\text{CS2}}$) and prescribed parameters (ρ_i, ρ_s, and f). Uncertainty due to the variability of ρ_w is neglected (Kurtz and Harbeck, 2017; Hendricks et al., 2016; Ricker et al., 2014). Here we assume that α^sat and $F_{\text{r}}^{\text{CS2}}$ are not correlated and have no systematic bias. Such an assumption may not be true in the real world. However, it allows us to estimate the retrieval uncertainty from satellite-derived products with a certain limit. Uncertainty of ice thickness can be estimated by the following Gaussian error propagation equation:

$$\begin{array}{}\text{(C1)}& {\mathit{ϵ}}_{y,\text{total}}^{\mathrm{2}}=\sum _x{\mathit{ϵ}}_y{\left(x\right)}^{\mathrm{2}}.\end{array}$$

Here, ε_y,total denotes the total uncertainty of retrieved variable y (h_s or H_i) and ε_y(x) denotes the uncertainty of y related to input variable x (α,F_r, ρ_i, ρ_s, or f). The uncertainties on the right-hand side are obtained by the following equation:

$$\begin{array}{}\text{(C2)}& {\mathit{ϵ}}_y\left(x\right)={\displaystyle \frac{\partial y}{\partial x}}{\mathit{\sigma }}_x={\displaystyle \frac{y\left(x+\mathit{\delta }\right)-y\left(x\right)}{\mathit{\delta }}}{\mathit{\sigma }}_x.\end{array}$$

Here, σ_x denotes the uncertainty of x, and δ is set to be 10⁻⁶ for the numerical calculation of the partial derivative using Eqs. (3) and (12); σ_α is estimated to be an RMSD value between α^OIB and α^sat, σ_Fr is given by Kurtz and Harbeck (2017), and σ_f is adopted from Armitage and Ridout (2015). Uncertainties of snow–ice densities are from the relevant literature (Alexandrov et al., 2020; Hendricks et al., 2016; Kern and Spreen, 2015; Ricker et al., 2014; Warren et al., 1999). Those values are summarized in Table C1.

Using Eqs. (C1) and (C2), uncertainties of snow depth and ice thickness retrievals can be estimated. Ice thickness uncertainty estimates are presented in Fig. C1. Total uncertainty of ice thickness estimates ranges from 0.8 to 2.0 m. Generally, F_r-related uncertainty in the third row is greater than α-related uncertainty in the second row. Snow depth uncertainty estimates are presented in Fig. C2. Total uncertainty of snow depths ranges from 0.04 to 0.4 m. In the case of the snow depth, α-related uncertainty is greater than F_r-related uncertainty. Both uncertainties of ice thickness and snow depth are greater for MYI regions than FYI regions. It is thought that the improvement of accuracy in satellite-derived temperatures can reduce the snow depth uncertainty, while the improvement of freeboard accuracy can reduce the ice thickness uncertainty. Uncertainties induced from densities and radar penetration factors are found to be relatively smaller than uncertainties related to α and F_r (shown in Figs. S2 and S3).

Table C1 Values and uncertainties of input variables for uncertainty estimation.

Figure C1 Geographical distributions of sea ice thickness uncertainty: (first row) total uncertainty, (second row) α-related uncertainty, and (third row) F_r-related uncertainty for the period from December 2013 to March 2014.

Figure C2 Geographical distributions of snow depth uncertainty: (first row) total uncertainty, (second row) α-related uncertainty, and (third row) F_r-related uncertainty for the period from December 2013 to March 2014.