Antarctic snow-covered sea ice topography derivation from TanDEM-X using polarimetric SAR interferometry

Single-pass interferometric synthetic aperture radar (InSAR) enables the possibility for sea ice topographic retrieval despite the inherent dynamics of sea ice. InSAR digital elevation models (DEMs) are measuring the radar scattering center height. The height bias induced by the penetration of electromagnetic waves into snow and ice leads to inaccuracies of the InSAR DEM, especially for thick and deformed sea ice with snow cover. In this study, an elevation difference between the satellite-measured InSAR DEM and the airborne-measured optical DEM is observed from a coordinated campaign over the western Weddell Sea in Antarctica. The objective is to correct the penetration bias and generate a precise sea ice topographic map from the single-pass InSAR data. With the potential of retrieving sea ice geophysical information by the polarimetric-interferometry (Pol-InSAR) technique, a two-layer-plus-volume model is proposed to represent the sea ice vertical structure and its scattering mechanisms. Furthermore, a simplified version of the model is derived, to allow its inversion with limited a priori knowledge, which is then applied to a topographic retrieval scheme. The experiments are performed across four polarizations: HH, VV, Pauli 1 (HH + VV), and Pauli 2 (HH − VV). The model-retrieved performance is validated with the optically derived DEM of the sea ice topography, showing an excellent performance with root-mean-square error as low as 0.26 m in Pauli-1 (HH + VV) polarization.

1 Introduction

Sea ice topography is defined as the elevation of the ice volume including the snow cover above the sea level. The sea ice topographic height on spatial scales of meters is dominated by ice ridges, shear zones, and hummocks, due to the forces from ocean winds and currents, together with the blocking effects from the coast and islands (Rampal et al., 2009). Timco and Burden (1997) estimated the ratio of the keel depth (i.e., depth of ice below the seawater) to sail height (i.e., height of ice above the seawater) for both first-year and multi-year sea ice ridges in the Beaufort Sea and highlighted their differences in ridge height and shape. Haas et al. (1999) presented the pressure ridge frequencies to be 3–30 ridges per kilometer over Bellingshausen, Amundsen, and Weddell seas in Antarctica. Tin and Jeffries (2003) indicated that first-year ridges in the Antarctic are flatter and less massive than those in the Arctic. Sea ice ridging height is a crucial parameter to evaluate total ice mass in both polar regions (Hibler et al., 1974; Melling and Riedel, 1995; Lytle et al., 1998; Tin et al., 2003). In the Antarctic, the mean height of the ridges in the Weddell Sea was found to be ∼1.1 m, which is similar to the ridging statistics from the Ross Sea (Lytle and Ackley, 1991), whereas it is considerably less than in the Arctic (Lytle and Ackley, 1991; Dierking, 1995).

The snow layer and ice properties in the Arctic and Antarctic are significantly different due to the diverse growing conditions in the two polar regions (Gloersen, 1992; Walsh, 2009; Sturm and Massom, 2009; Webster et al., 2018). In the Antarctic, the snow depth is reported to be thicker than in the Arctic (Jeffries et al., 1997; Massom et al., 2001; Willatt et al., 2009). When thick enough, the snow will overburden the ice floe and be flooded by seawater, resulting in higher salinity of the snow layer in the Antarctic. Besides, compared to the Arctic, snow on the Antarctic sea ice comprises more heterogeneous layers resulting from highly variable temperature (Massom et al., 2001). The layer heterogeneity in types, density, salinity, and wetness determines the electromagnetic characteristics of the snow. As for ice properties in general, Antarctic sea ice is reported to be thinner (Worby et al., 2008; Kurtz and Markus, 2012; Lindsay and Schweiger, 2015), younger (Webster et al., 2018), and more saline than in the Arctic at comparable age and thickness (Gow et al., 1982). Quantitatively, the mean salinities of the Antarctic first-year ice and multi-year ice profiles are 4.6 ‰ and 3.5 ‰, respectively, whereas the average values are 3 ‰ for the first-year ice and 2−2.5‰ for the multi-year ice in the Arctic (Cox and Weeks, 1974). These variable properties of sea ice, its ridges, and snow cover, at both small and large spatial scales, highlight the challenge and necessity for accurate sea ice topographic information with large spatial coverage and high resolution.

Characterizing sea ice topography is valuable for various geophysical parameters over polar oceans. For instance, the atmospheric drag coefficient over sea ice is an important topography-dependent parameter to understand the interaction at the ice–atmosphere boundary (Garbrecht et al., 2002; Castellani et al., 2014). Interpretation of sea ice topography is also essential in estimating sea ice thickness. Focusing on first-year sea ice in the Alaska region, Tucker and Govoni (1981) observed a square-root relation between the ridge height and thickness, which is further validated by additional in situ observations in Tucker et al. (1984). Petty et al. (2016) presented a detailed characterization of Arctic sea ice topography across both first-year and multi-year sea ice and analyzed the topographic differences between the two ice regimes. A square-root relation function between sea ice topographic height and thickness was established for ice thickness retrieval (Petty et al., 2016). The results demonstrated a maximum ±2 m difference between the measured and predicted ice thickness. Note that the measured thickness ranges from 0 to 8 m with an initial uncertainty of 0.8 m (Petty et al., 2016). In Antarctica, Toyota et al. (2011) demonstrated that the mean ice thickness of snow-covered ice is highly correlated with the sea ice topography. Nowadays, precise characterization of sea ice topography is a topic of active research.

The sea ice topography can be measured by various instruments, such as laser altimeters (Dierking, 1995; Schutz et al., 2005; Abdalati et al., 2010; Farrell et al., 2011, 2020) and stereo cameras using photogrammetric techniques (Dotson and Arvesen., 2014; Divine et al., 2016; Nghiem et al., 2018; Li et al., 2019). However, the major limitation of above measurements is the small spatial coverage. Synthetic aperture radar (SAR) achieves a good balance between wide spatial coverage and high resolution and becomes an invaluable asset for monitoring polar regions thanks to its ability to provide all-weather, day/night imagery. Interferometric SAR (InSAR) offers an opportunity to estimate surface height from two or more image pairs (Rodriguez and Martin, 1992). However, due to the inherent dynamics of sea ice, it is impossible to derive height over sea ice from a single SAR sensor by repeat-pass interferometry, because of its temporal decorrelation. Only single-pass interferometry offers the possibility to characterize the sea ice topography (Dierking et al., 2017). TanDEM-X is a single-pass SAR interferometer developed by the German Aerospace Center (DLR) (Krieger et al., 2007) and is providing high-resolution coregistered data on a global scale. Dierking et al. (2017) demonstrated the theoretical potential of generating sea ice height from single-pass InSAR data and discussed the factors that may impede the accuracy of the retrieval. From TanDEM-X InSAR acquisitions, the derivation of topography over snow-free multi-year sea ice was demonstrated and verified with laser and photogrammetric measurements (Yitayew et al., 2018). Until now, the InSAR technique has become one of the most promising tools for sea ice height estimation.

However, a digital elevation model (DEM) derived with InSAR is affected by the penetration of microwave signals into dry, frozen snow and ice. In fact, an InSAR DEM is actually a measurement of the radar scattering center height, which can be below the surface due to the microwave penetration. This height bias leads to inaccuracies of InSAR DEMs, especially for multi-year sea ice with snow cover. The microwave penetration into snow and ice is described by the electromagnetic penetration depth δ_p. It is determined by the signal extinction coefficient σ in units of decibels per unit length (dB m⁻¹), which indicates the decrease in the signal strength inside the medium. The total electromagnetic loss in a medium consists of both scattering and absorption losses. Scattering loss results from particles of different relative permittivity embedded in a host medium. The absorption loss depends on the imaginary part of the relative permittivity ϵ^′′ (Hallikainen and Winebrenner, 1992). Larger δ_p values are found in multi-year ice due to the smaller ϵ^′′ attribution to reduced brine compared to first-year ice. For sensors operating at X-band, experimental penetration depth for sea ice ranges from about 0.05 to 1 m, depending on the sea ice type, salinity, and temperature (Hallikainen and Winebrenner, 1992). Snow on top causes a greater range of δ_p due to the high sensitivity of ϵ^′′ to water content. Dry and fine-grained snow can have δ_p values up to hundreds of wavelengths (Cloude, 2010).

In this study, we observed an elevation discrepancy between the InSAR DEM and the photogrammetric DEM which were acquired in a coordinated campaign (Nghiem et al., 2018) conducted with DLR's TanDEM-X satellite and the NASA IceBridge aircraft over the sea ice in the western Weddell Sea, Antarctica. The elevation difference reveals the necessity to consider the penetration depth δ_p when retrieving sea ice topography from the InSAR imagery. The objective of this study is to compensate for the penetration depth and thereby obtain a more accurate sea ice topographic map with wide spatial coverage. Note that the studied area is snow-covered sea ice; therefore, the term “sea ice topographic height”, throughout the paper, refers to the sea ice height including snow depth above local sea level.

The estimation of the penetration depth of InSAR signals can be inferred from the interferometric volume decorrelation, which is one of the key components of the interferometric coherence. The volume decorrelation is caused by backscatter contributions from different depths and can be derived from the integral of an assumed vertical scattering distribution function. The investigation of vertical distribution functions for various scattering processes, known as the polarimetric-interferometry SAR (Pol-InSAR) technique (Papathanassiou and Cloude, 2001), is widely applied in retrieving geophysical parameters from natural volumes, such as forests (Kugler et al., 2015), agriculture (Joerg et al., 2018), ice sheets (Fischer et al., 2018), and glaciers (Sharma et al., 2012). To the best of our knowledge, few studies have assessed the potential of retrieving geophysical information by means of the Pol-InSAR technique for sea ice. Dierking et al. (2017) estimated the penetration depth into sea ice volume under the assumption of a uniform lossy volume with an exponential vertical function. However, for snow-covered sea ice, the scattering effects from the snow volume and sub-layers, such as the snow–ice interface, also need to be considered. To achieve an effective estimation of penetration depth, factors including the physical structures, the electromagnetic properties, and the scattering mechanisms within the sea ice volume need to be understood and properly modeled.

Sea ice can be modeled as a multi-layer structure behaving as a mixture of surface, volume, and surface–volume interaction scattering in microwave remote sensing (Nghiem et al., 1995a, b; Albert et al., 2012). When the electromagnetic waves penetrate the volume, the inhomogeneous materials inside the volume (i.e., a mixture of constituents such as brine, ice, and air bubbles) excite the occurrence of volume scattering (Nghiem et al., 1995a). Besides, the surface conditions such as rough interfaces, hummocks, and snow cover can increase the surface scattering at the rough air–snow interface, snow–ice interface, and ice–water interface (Nghiem et al., 1990). In the Arctic first-year thin ice, snow capillary force gives rise to brine wicking, and, consequently, a layer of high-salinity slush ice appears at the snow–ice interface (Reimnitz and Kempema, 1987; Drinkwater and Crocker, 1988; Nghiem et al., 1995a). In the Antarctic, ice-surface flooding widely occurs, resulting from the generally thicker snow layer loading on the thinner ice floes, often followed by freezing of the slush layer at the snow–ice interface (Massom et al., 2001; Jeffries et al., 2001; Maksym and Jeffries, 2000). Even without flooding, the upward wicking of brine from the ice surface can also form a saline layer at the bottom of the snowpack (Massom et al., 2001; Toyota et al., 2011; Webster et al., 2018). The slush layer at the snow–ice interface can induce significant surface scattering and thus has been included in sea ice scattering modeling (Nghiem et al., 1995a, b; Maksym and Jeffries, 2000). Moreover, the surface–volume interaction components (Albert et al., 2012) further complicate the overall scattering mechanisms.

Huang and Hajnsek (2021) investigated the X-band SAR polarimetric behavior for several types of ice over the western Weddell Sea, including new ice, thin ice, thick ice, and deformed ice with ridges. For the area covered by the thick and deformed ice, an empirically inverse relation between the elevation difference (i.e., penetration bias) and the co-polarimetric coherence was observed, indicating that SAR polarimetry carries significant topographic information (Huang and Hajnsek, 2021). Based on Huang and Hajnsek (2021), this study offers a further understanding of the InSAR penetration bias by investigating the polarimetric behavior and exploiting the interferometric volume decorrelation. A novel model is proposed to characterize the scattering processes, and an inversion scheme is developed for height retrieval. Therefore, compared to the previous work, this study is a crucial step forward towards developing an advanced algorithm for sea ice topographic estimation.

In this study, inspired by multi-layer sea ice models utilized for electromagnetic simulation (Nghiem et al., 1990, 1995a, b; Albert et al., 2012), a two-layer-plus-volume model is proposed to relate interferometric coherence to extinction coefficients, layer depths, and layer-to-volume scattering ratios. The model sensitivity to the variation of several parameters is analyzed, and the model accuracy is assessed with various baseline configurations. With the goal to develop and invert the model for sea ice topographic retrieval, the proposed theoretical model is further simplified by reducing the required amount of input parameters. An inversion scheme for topographic retrieval using both the theoretical and simplified model is established. The sea ice topographic height is retrieved in different polarizations over around 50 km × 18 km in the western Weddell Sea. The model-retrieved DEM is validated against a photogrammetric DEM, proving the effectiveness of the proposed model and the inversion scheme.

The organization of the paper is as follows. Section 2 introduces the basic concepts of Antarctic sea ice conditions and the Pol-InSAR technique. Section 3 introduces the data sets and the preprocessing procedures. A two-layer theoretical model and a simplified model are proposed in Sect. 4. A model-inversion scheme is developed in Sect. 5 to achieve the sea ice topographic retrieval. In Sect. 6, the proposed methodology is applied to the study area, and the experimental results are analyzed. More discussions about the model are given in Sect. 7, and the conclusion is drawn in Sect. 8.

2 Basic concepts

2.1 Sea ice conditions

In the Antarctic, the presence of a saline layer at the snow–ice interface due to the flooding or capillary suction of brine from the ice surface has been recognized as a widespread and critical phenomenon (Massom et al., 2001). For thinner ice, flooding may occur when the weight of the snow pushes the ice surface below the water level, yielding a negative freeboard. In this case, as shown in Fig. 1a, seawater infiltrates into the snowpack, floods the ice surface, and creates a high-saline slush layer which may refreeze into snow ice (Lange et al., 1990; Jeffries et al., 1997; Maksym and Jeffries, 2000). Note that snow ice forms from slush which by definition (and ignoring brine wicking or other second-order processes) locates below the water level, and the freeboard of snow ice is generally considered to be zero. The freeboard of snow ice could become positive if there is subsequent ice growth at the bottom of the ice sheet, without new snow accumulation; however, that is a secondary effect as snow thickness will generally tend to increase in the meanwhile. The thickness of snow ice was observed to be ∼42 %–70 % of the total snow accumulation (i.e., the thickness of snow ice plus snow depth) (Jeffries et al., 2001).

Figure 1 Schematic of (a) thinner ice floes flooded by seawater and (b) thicker ice floes without flooding.

For thicker and deformed ice with ridges, less flooding occurs due to the increased buoyancy of the ice mass contained in the ridges (Jeffries et al., 1998). However, even in the absence of flooding, a thin slush layer can also occur due to the capillary suction of brine from the ice surface (Massom et al., 2001; Webster et al., 2018). Besides, the deformed ice in the ridging and rafting area is often poorly consolidated, and thus seawater may reach the snow layer and form a thin slush layer (Maksym and Jeffries, 2000). The sea ice structure for thicker ice without flooding is sketched in Fig. 1b, including snow on top, the ice volume, and a thin and high-saline layer in between.

The condition of flooding can be quantified by a simple hydrostatic balance (Lange et al., 1990):

$$\begin{array}{}\text{(1)}& \begin{array}{rl}& {\mathit{\rho }}_{\mathrm{w}}d={\mathit{\rho }}_{\mathrm{i}}d+{\mathit{\rho }}_{\mathrm{i}}f+{\mathit{\rho }}_{\mathrm{s}}s\\ & t=d+f,\end{array}\end{array}$$

where ρ_w is the seawater density, ρ_i is the ice density, and ρ_s is the snow density; d and f are the thickness of ice below and above the sea level, respectively; t is the thickness of the total ice thickness; and s is the snow depth. In case of a flooding, f should be zero (i.e., d=t), and Eq. (1) becomes

$$\begin{array}{}\text{(2)}& s/t=\left({\mathit{\rho }}_{\mathrm{w}}-{\mathit{\rho }}_{\mathrm{i}}\right)/{\mathit{\rho }}_{\mathrm{s}}\approx \mathrm{0.12}/{\mathit{\rho }}_{\mathrm{s}}\end{array}$$

by assuming ρ_w=1.03 Mg m⁻³ and ρ_i=0.91 Mg m⁻³ (Lange et al., 1990). For snow density ρ_s being 0.3 Mg m⁻³ (Lange et al., 1990), the ratio between snow depth and ice thickness $s/t$ is estimated to be 0.4. The snow depth on Antarctic sea ice during September and November was shown to be below 0.8 m for 99 % of the samples in Webster et al. (2018). This range of snow depth will lead to flooding for ice thickness <2 m.

The relation between ice thickness H_i and surface height h_sur (i.e., ice height above sea surface including snow depth) has been discussed over different regions (Petty et al., 2016; Toyota et al., 2011; Ozsoy-Cicek et al., 2013). Ozsoy-Cicek et al. (2013) showed a linear relation $H_{\mathrm{i}}=\mathrm{2.24}{h}_{\mathrm{sur}}+\mathrm{0.228}$ fitted from large-scale, survey-averaged data over the western Weddell Sea, which is the same region as this study. According to this linear relation, H_i=2 m corresponds to a surface height of ∼0.8 m.

This paper focuses on thicker (>2 m) and deformed ice (Fig. 1b), which is the main ice typology in the studied area. In the following sections, the model and experiments are conducted only for the samples above ∼0.8 m surface height. We assume that samples exceeding this threshold are thicker and deformed ice without flooding. The potential to extend the proposed model to thinner ice scenarios (e.g., Fig. 1a) is discussed in Sect. 7.3.

2.2 Pol-InSAR technique

The complex interferometric coherence ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$ is a measurement of signal correlation between two acquisitions. For single-pass systems, ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$ can be decomposed into a product of terms (Cloude, 2010):

$$\begin{array}{}\text{(3)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}=e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\mathit{\gamma }}_{\mathrm{s}}{\mathit{\gamma }}_{\mathrm{SNR}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}},\end{array}$$

where ϕ₀ is the topographic phase. γ_s is the baseline or surface decorrelation, which depends on the nature of the surface scattering; it can always be removed by employing range spectral filtering and thus is set equal to 1 in this study. γ_SNR denotes decorrelation due to additive noise in the signals. ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$ refers to the complex volume decorrelation. The topographic phase can be converted to topographic height h_topo by

$$\begin{array}{}\text{(4)}& h_{\mathrm{topo}}={\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}}{{\mathit{\kappa }}_z}}=h_{\mathrm{a}}{\displaystyle \frac{{\mathit{\varphi }}_{\mathrm{0}}}{\mathrm{2}\mathit{\pi }}},\end{array}$$

where h_a is the height of ambiguity and κ_z is the vertical wavenumber in free space. Note that h_a corresponds to an interferometric phase change of 2π, and is inversely proportional to the perpendicular baseline between the two acquisitions (Dall, 2007):

$$\begin{array}{}\text{(5)}& h_{\mathrm{a}}={\displaystyle \frac{\mathit{\lambda }H\tan \mathit{\theta }}{\mathrm{2}{b}_{⟂}}},\end{array}$$

where λ is the wavelength, H is orbit height, and b_⟂ is the effective perpendicular baseline of the TanDEM-X bi-static mode.

The magnitude of the ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$ can be corrected for γ_SNR and γ_s by rewriting Eq. (3) as

$$\begin{array}{}\text{(6)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}={\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}}{{\mathit{\gamma }}_{\mathrm{s}}{\mathit{\gamma }}_{\mathrm{SNR}}}}=e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}},\end{array}$$

where ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}$ is the noise-removed interferometric coherence. γ_SNR can be estimated as a function of signal-to-noise ratio (SNR) (Cloude, 2010):

$$\begin{array}{}\text{(7)}& {\mathit{\gamma }}_{\mathrm{SNR}}={\displaystyle \frac{\mathrm{SNR}}{\mathrm{1}+\mathrm{SNR}}}={\displaystyle \frac{S\left(\mathrm{dB}\right)-N\left(\mathrm{dB}\right)}{\mathrm{1}+S\left(\mathrm{dB}\right)-N\left(\mathrm{dB}\right)}},\end{array}$$

with S being the backscattering signal and N being the noise floor (i.e., the noise equivalent sigma zero – NESZ) (Eineder et al., 2008).

In the case of pure surface scattering, the interferometric coherence can be approximated to ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}\approx e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\mathit{\gamma }}_{\mathrm{SNR}}$, assuming that volume scattering can be neglected (i.e., ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}\approx \mathrm{1}$). γ_SNR only contributes to the magnitude of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$; therefore, the InSAR scattering phase center, denoted as $\mathrm{\angle }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$, purely contains the information of topographic phase ϕ₀. In this case, $\mathrm{\angle }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$ can be directly converted to topographic height.

However, in the case of snow-covered, thick and deformed sea ice, when the microwaves penetrate into the snow and ice volume, the inhomogeneous materials inside the volume can excite volume scattering (Nghiem et al., 1995a). Then, the volume decorrelation in Eq. (3) is not 1. Both the topographic phase ϕ₀ and the complex ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$ contribute to the InSAR scattering phase center $\mathrm{\angle }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$. In this case, in order to obtain an accurate topographic phase ϕ₀, ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$ has to be properly modeled and estimated. The main contribution of this paper is the development of a novel two-layer-plus-volume model (Sect. 4) for ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$, which is applied for an improved sea ice topographic retrieval.

The volume decorrelation ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$ depends on the vertical distribution of backscattering σ_v(z) (Cloude, 2010):

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}={\displaystyle \frac{\underset{\mathrm{0}}{\overset{D}{\int }}{\mathit{\sigma }}_{\mathrm{v}}\left(z\right)e^{i{\mathit{\kappa }}_{z\mathit{\_}\mathrm{vol}}z}\mathrm{d}z}{\underset{\mathrm{0}}{\overset{D}{\int }}{\mathit{\sigma }}_{\mathrm{v}}\left(z\right)\mathrm{d}z}},\end{array}$$

where the surface is located at z=0, D is the thickness of volume, and κ_{z_vol} is the vertical wavenumber in the volume (Sharma et al., 2012; Dall, 2007):

$$\begin{array}{}\text{(9)}& {\mathit{\kappa }}_{z\mathit{\_}\mathrm{vol}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{h_{\mathrm{a}\mathit{\_}\mathrm{vol}}}}={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{h_{\mathrm{a}}\frac{\sqrt{{\mathit{ϵ}}^{\prime }-{\sin }^{\mathrm{2}}\mathit{\theta }}}{{\mathit{ϵ}}^{\prime }\cos \mathit{\theta }}}}={\displaystyle \frac{{\mathit{\kappa }}_z}{\frac{\sqrt{{\mathit{ϵ}}^{\prime }-{\sin }^{\mathrm{2}}\mathit{\theta }}}{{\mathit{ϵ}}^{\prime }\cos \mathit{\theta }}}},\end{array}$$

where h_{a_vol} is the height of ambiguity in the volume, θ is the incidence angle (in air), and ϵ^′ is the dielectric constant of the volume and is assumed to be 2.8 (Dierking et al., 2017) throughout this study. The phase of the volume decorrelation $\mathrm{\angle }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$ can be translated to height h_volume as

$$\begin{array}{}\text{(10)}& h_{\mathrm{volume}}={\displaystyle \frac{\mathrm{\angle }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}}{{\mathit{\kappa }}_{z\mathit{\_}\mathrm{vol}}}}.\end{array}$$

In Eq. (8), ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$ can be estimated by choosing an appropriate vertical structural function σ_v(z) and a suitable InSAR baseline configuration using Eq. (9) and then be substituted into Eq. (3) to obtain the topography of snow-covered sea ice.

3 Data sets and preprocessing

This section introduces the campaign, the study area, and the data sets. The InSAR processing and its performance are also described.

3.1 Campaign and study area

A coordinated campaign of NASA's Operation IceBridge (OIB) airborne mission and the DLR's TanDEM-X satellite mission was successfully conducted on 29 October 2017, named the OIB/TanDEM-X Coordinated Science Campaign (OTASC) (Nghiem et al., 2018). The OTASC data have been successfully used in investigating the topography of icebergs (Dammann et al., 2019) and sea ice (Huang and Hajnsek, 2021).

As presented in Fig. 2a, the study area is located in the western Weddell Sea, near the east coast of the Antarctic Peninsula. The TanDEM-X SAR intensity image of the study area is shown in Fig. 2b, where the optical images and the transect of photogrammetric measurements of the sea ice topography are superimposed. From the optical images of the airborne digital camera, it is visible that the study region consists of snow-covered, thick and deformed ice with ridges.

Figure 2 (a) Geolocation of the study area (the grey rectangle). (b) Composite of the SAR intensity image in HH polarization and the airborne digital mapping system (DMS) digital camera images. The green dashed line indicates the 50 km transect of the DMS DEM data. Sub-image A and sub-image B are close-ups of area A and area B with small-scale DMS images (acquired at around 17:50 UTC) superimposed on the large-scale DMS images (acquired at around 22:05 UTC). The red and yellow dots denote the selected reference pairs from the large- and small-scale DMS image, respectively. The green arrow denotes the shift vector, which is used for data coregistration.

3.2 TanDEM-X

The German TanDEM-X mission is a single-pass SAR interferometer operating at X-band at a wavelength of 3 cm. With nearly no temporal gap, TanDEM-X collects two images of the same footprint seen from slightly different viewing angles to generate the topography of the Earth's surface (Krieger et al., 2007). The studied InSAR images were acquired at 23:41 UTC on 29 October 2017 in bi-static mode. The InSAR pair is a dual-polarization (HH and VV) StripMap product. The incidence angle of the scene center is 34.8^∘, and the pixel spacing is 0.9 m × 2.7 m in range and azimuth. The effective perpendicular baseline b_⟂ is 175.7 m, and the along-track baseline b_al is 201.9 m. The height of ambiguity h_a is 32.5 m.

3.3 DMS

The digital mapping system (DMS) is one of the OIB instruments acquiring different data from a digital camera system. This study uses the airborne DMS digital camera images (Dominguez, 2018) and the airborne DMS DEMs (Dotson and Arvesen., 2014).

The DMS digital camera captures natural color and panchromatic imagery, hereby named DMS images for conciseness. The DMS images are geolocated and orthorectified, with a high spatial resolution varying from 0.015 to 2.5 m depending on the flight altitude (Dominguez, 2018). Two types of photography, large- and small-scale DMS images, are obtained during the airborne overflights. The large-scale DMS images over the study area were acquired from 22:01 to 22:07 UTC on 29 October 2017, each with a spatial coverage of around 5.8 km by 8.8 km, shown in Fig. 2b. The small-scale DMS images were captured from 17:45 to 17:52 UTC on 29 October 2017, each with about 400 m by 400 m spatial coverage. The transect of small-scale DMS images is shown in the green dashed line in Fig. 2b and enlarged in Fig. 2b sub-image A and B, where the details of sea ice structure become visible.

The DMS DEM is generated from the small-scale DMS image by a photogrammetric technique and is calibrated with lidar measurements (Dotson and Arvesen., 2014). For the snow-covered sea ice, the DMS DEM measures sea ice height including snow depth. The data are acquired along a 50 km transect with a swath width of 400 m (Dotson and Arvesen., 2014) and 40×40 cm spatial resolution. Note that the temporal gap between the DMS DEM and TanDEM-X SAR acquisitions is about 6 h.

In this study, the DMS DEMs are further processed following four steps: reprojection, mosaicing, geocoding, and calibration. First, the DMS DEMs are reprojected from Antarctic Polar Stereographic to WGS84 spatial reference. The second step is the mosaicing of adjacent files into the 50 km transect with a swath width of 400 m. Third, using the GAMMA software, the merged DMS DEM is geocoded into the SAR coordinate system and resampled into the same resolution (i.e., $\sim \mathrm{10}\times \mathrm{10}$ m in range and azimuth) as the multi-look TanDEM-X image. Finally, as the DMS DEMs are given in meters above the WGS-84 ellipsoid, the sea ice topographic height in this paper is calibrated to the local sea level by selecting the water-surface reference from DMS images. In total, we label nine points as water-surface reference according to the DMS images. Since the interferometric coherence magnitude over water is very low, it can also be used to classify open water (Dierking et al., 2017). All the nine points have an interferometric coherence magnitude below 0.3, which is the threshold of the open-water area mask in Huang and Hajnsek (2021). The average height of the open-water points is subtracted from the DMS DEMs to obtain the sea ice topographic height relative to the local sea level.

3.4 Data coregistration

Due to the inherent dynamics of sea ice and the temporal gap between the DMS DEM and TanDEM-X acquisitions, data coregistration is employed to cancel the shift and thereby ensure a valid pixel-by-pixel comparison.

The large- and small-scale DMS images, although acquired at different times, both clearly reveal the shape and size of ice floes; therefore, they are used to match identical sea ice features, referred to as the “reference pair” in the following. The coregistration is performed by tracking the movement of the selected reference pairs. Specifically, we manually label several pairs of distinguishable features (i.e., the ice floes of a particular shape, or leads) on both the large-scale DMS image and the small-scale DMS image, which is acquired about 4 h 15 min later. By extracting the two geolocations of the reference pair, a shift-velocity vector can be derived. As the temporal difference between the DMS DEM and the SAR acquisition is 5 h 49–56 min, the shift distance can be estimated, assuming constant sea ice motion, based on the shift-velocity vector and utilized for the coregistration.

In the studied image, the focus is the 50 km transect (green dashed line in Fig. 2b) of the DMS flight track. The transect is divided into 50 segments, and each segment contains 11×100 pixels in range and azimuth corresponding to about 110×1000 m area. For each segment, several reference pairs are selected and labeled (marked in red and yellow dots on the large- and small-scale DMS images, respectively, in Fig. 2b sub-image A and B). The shift-velocity vectors are calculated and annotated by the green arrows. Then, data coregistration is conducted by multiplying the derived shift-velocity vector with the temporal gap for each segment respectively. The coregistered results of all segments are confirmed by the visualization of the DMS images and the SAR images. Among the 50 segments, 12 segments which still contain residual mis-coregistration induced by the sea ice non-linear movement or rotation are excluded and will not be used in the following experiments. A total of 76 % of the segments from the whole SAR scene are accepted as correctly coregistered segments in this study.

3.5 InSAR processing

The TanDEM-X InSAR pair is already coregistered and common spectral band filtered in range and azimuth (Duque et al., 2012). The remaining InSAR processing includes interferogram generation, flat earth removal, interferogram filtering, low-coherence area mask, phase unwrapping, and phase-to-height conversion; see details in Huang and Hajnsek (2021). All steps are carried out with the GAMMA software.

In single-pass interferometry, two simultaneous observations, denoted as s₁ and s₂, are made. The complex interferogram γ and interferometric phase ϕ_γ can be calculated as (Cloude, 2010)

$$\begin{array}{}\text{(11)}& {\displaystyle }& {\displaystyle }\mathit{\gamma }=s_{\mathrm{1}}{s}_{\mathrm{2}}^{*},\text{(12)}& {\displaystyle }& {\displaystyle }{\mathit{\varphi }}_{\mathit{\gamma }}=\arg \left\{s_{\mathrm{1}}{s}_{\mathrm{2}}^{*}\right\},\end{array}$$

where symbol (^*) denotes the complex conjugate.

The interferometric coherence is estimated by (Cloude, 2010)

$$\begin{array}{}\text{(13)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}={\displaystyle \frac{<s_{\mathrm{1}}{s}_{\mathrm{2}}^{*}>}{\sqrt{<s_{\mathrm{1}}{s}_{\mathrm{1}}^{*}><s_{\mathrm{2}}{s}_{\mathrm{2}}^{*}>}}},\end{array}$$

where the symbol $<.>$ denotes an ensemble average. Here, a 4×12 window in azimuth and slant range, corresponding to about 10 m × 10 m spatial size, is applied to estimate ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}$ for four polarizations: HH, VV, Pauli 1 (HH + VV), and Pauli 2 (HH − VV). Areas with $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}\right|$ less than 0.3 are masked out and will not be considered in the following processing. For conciseness, only the interferometric coherence in HH polarization is shown in Fig. 3a. The $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}\right|$ histograms for the four polarizations are plotted in Fig. 4, from which the interferometric decorrelation varying among different polarizations can be observed. $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}\right|$ values for HH and VV polarizations show small differences, mainly lying in a range of 0.6–0.8. The Pauli-1 polarization has the highest $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}\right|$ of 0.7–0.8 with a narrow distribution, whereas the Pauli-2 polarization shows the lowest values with a wider spread of the coherence, which is mainly due to the lower SNR. The observed interferometric decorrelation indicates the necessity to consider the volume scattering contributing to the InSAR decorrelation.

Figure 3 (a) Magnitude of the interferometric coherence $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}\right|$ and (b) InSAR DEM h_InSAR for HH polarization.

Figure 4 Magnitude of the interferometric coherence $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{InSAR}}\right|$.

The InSAR DEM h_InSAR is derived for the four polarizations. Again, only the HH polarization is shown in Fig. 3b for conciseness. The comparison between h_InSAR for the four polarizations and the DMS DEM h_DMS along the flight track (the green dashed line in Fig. 2) is shown in Fig. 5, with a maximum elevation difference around 2 m. The differences of h_InSAR across the four polarizations are illustrated in Fig. 6, where the height differences mostly lie in the range of −0.5 to 0.5 m. The InSAR-derived heights from Pauli-1 and HH polarizations reveal the most similar values, while Pauli-2 and HH polarizations show the largest height difference, which is in accordance with the wider spread distribution of InSAR coherence in the Pauli-2 channel (see Fig. 4). The root-mean-square errors (RMSEs) between h_InSAR and h_DMS are averaged to be ∼1.10 m for the four polarizations, indicating that the penetration of electromagnetic waves into snow and ice should be properly considered and corrected for sea ice topographic retrieval, at least for the deformed thick ice with snow cover in this study.

Figure 5 The InSAR-derived height profiles (h_InSAR) and the DMS DEM (h_DMS). Each profile represents the height along a 1×5000-pixel section at the center of the coregistered segment. The mis-coregistered and h_DMS-below-0.8 m samples are excluded from the plots.

Figure 6 Height difference between the InSAR-derived height (h_InSAR) in HH polarization and the other three polarizations.

It can be summarized that both the interferometric decorrelation and the elevation difference between the InSAR DEM and the DMS DEM highlight the necessity for developing an appropriate method aimed at an accurate sea ice topography retrieval.

4 Model development

This section proposes a two-layer-plus-volume model to describe the interferometric coherence of sea ice and snow on top. Simulations are performed to analyze the model sensitivity and accuracy by varying parameter sets and baseline configurations. In a separate step, the model is further simplified for the practical purpose of deriving sea ice topography with limited a priori parameter knowledge.

4.1 Composite coherence model for sea ice

The sea ice volume has been modeled as a multi-layer structure in microwave remote sensing (Nghiem et al., 1990; Albert et al., 2012). We propose a two-layer-plus-volume model considering snow cover, ice volume, and seawater, illustrated in Fig. 7, behaving as a mix of surface and volume scattering under radar illumination. The uppermost surface (i.e., snow–air interface) is located at z₀.

Surface scattering is considered to originate mainly from two interfaces, named the top layer and the bottom layer, respectively. The top layer located at z₁ is the snow–ice interface, which can induce significant surface scattering due to a slush layer with high permittivity (Hallikainen and Winebrenner, 1992; Maksym and Jeffries, 2000). This slush layer is widespread on the Antarctic sea ice and increases the radar backscattering as well as limits the signal penetration compared to a smooth and dry snow–ice interface. As long as the slush layer has a small vertical extent, it is irrelevant for the Pol-InSAR scattering structure model, regardless of whether the top layer represents the snow–ice interface, the snow–slush interface, or both. The position of the bottom layer (z₂) could be somewhere inside the ice volume or at the ice–water interface. The vertical distributions of the top and bottom surface can be modeled as two Dirac delta functions at the specific layer position. An additional parameter, the layer-to-volume scattering ratio, accounts for the (relative) scattering from these interfaces, depending, e.g., on roughness and dielectric contrast (Fischer et al., 2018). Hence, the surface scattering component in the context of interferometry is modeled as (Fischer et al., 2018)

$$\begin{array}{}\text{(14)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{Layer}}=e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\displaystyle \frac{m_{\mathrm{1}}{e}^{i{\mathit{\kappa }}_{z\mathit{\_}\mathrm{vol}}{z}_{\mathrm{1}}}+m_{\mathrm{2}}{e}^{i{\mathit{\kappa }}_{z\mathit{\_}\mathrm{vol}}{z}_{\mathrm{2}}}}{m_{\mathrm{1}}+m_{\mathrm{2}}}},\end{array}$$

where m₁ and m₂ are the layer-to-volume ratio of the top and bottom layer, respectively.

The volume scattering is attributed to the constituents in the snow (from z₀ to z₁) and ice (from z₁ to z₂) volumes (Hallikainen and Winebrenner, 1992). Both volumes are assumed to be uniform volumes, which means that the scattering coefficient per unit volume and the extinction coefficient σ have no spatial variation. In this case, the vertical structure function σ_v(z) becomes exponential. The ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}$ for a uniform volume model can be formulated as (Papathanassiou and Cloude, 2001)

$$\begin{array}{}\text{(15)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}(\mathit{\sigma },D)={\displaystyle \frac{\underset{\mathrm{0}}{\overset{D}{\int }}{e}^{\frac{\mathrm{2}\mathit{\sigma }z}{\cos {\mathit{\theta }}_r}}{e}^{i{\mathit{\kappa }}_{z\mathit{\_}\mathrm{vol}}z}\mathrm{d}z}{\underset{\mathrm{0}}{\overset{D}{\int }}{e}^{\frac{\mathrm{2}\mathit{\sigma }z}{\cos {\mathit{\theta }}_r}}\mathrm{d}z}},\end{array}$$

where θ_r is the incidence angle in the volume, σ is the constant extinction coefficient, and D is the volume thickness. The corresponding volume coherences can be derived according to Eq. (15), denoted as ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}})$ and ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{2}},z_{\mathrm{12}})$ for the snow and ice volumes, respectively, where $z_{\mathrm{01}}=z_{\mathrm{0}}-z_{\mathrm{1}}$ is the thickness of the snow volume and $z_{\mathrm{12}}=z_{\mathrm{1}}-z_{\mathrm{2}}$ is the thickness of the ice volume.

Figure 7 Schematic of the proposed two-layer-plus-volume model for sea ice.

For the overall two-layer-plus-volume model, the interferometric coherence can be given as a combination of volume and surface effects which are described by Eqs. (14) and (15), respectively. As represented in Fig. 7, if we set z₀ as the origin of the coordinate and thus to be 0, with z₁ and z₂ being the position of two layers with negative values, the composite interferometric coherence is postulated to be

$$\begin{array}{}\text{(16)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}=\\ & e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\displaystyle \frac{\mathit{\alpha }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}\left({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}}\right)+e^{i{\mathit{\varphi }}_{\mathrm{1}}}(\mathrm{1}-\mathit{\alpha }){\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}\left({\mathit{\sigma }}_{\mathrm{2}},z_{\mathrm{12}}\right)+m_{\mathrm{1}}{e}^{i{\mathit{\varphi }}_{\mathrm{1}}}+m_{\mathrm{2}}{e}^{i{\mathit{\varphi }}_{\mathrm{2}}}}{\mathrm{1}+m_{\mathrm{1}}+m_{\mathrm{2}}}}\\ & =e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}\left({\mathit{\sigma }}_{\mathrm{1}},{\mathit{\sigma }}_{\mathrm{2}},\mathit{\alpha },m_{\mathrm{1}},m_{\mathrm{2}},z_{\mathrm{1}},z_{\mathrm{2}}\right),\end{array}\end{array}$$

where ϕ₁=κ_{z_vol}z₁, ϕ₂=κ_{z_vol}z₂, and σ₁ and σ₂ are extinction coefficients of snow and ice volume, respectively, in the unit of nepers per meter (Np m⁻¹). Note that $\mathit{\sigma }\left(\mathrm{dB}{\mathrm{m}}^{-\mathrm{1}}\right)=\frac{\mathrm{10}}{\ln \mathrm{10}}\mathit{\sigma }\left(\mathrm{Np}{\mathrm{m}}^{-\mathrm{1}}\right)=\mathrm{4.343}\mathit{\sigma }\left(\mathrm{Np}{\mathrm{m}}^{-\mathrm{1}}\right)$. The volume coherences of snow ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}})$ and ice ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{2}},z_{\mathrm{12}})$ can be obtained according to Eq. (15). Weight parameter α ($\in [\mathrm{0},\mathrm{1}]$) represents the proportion of the snow volume scattering in the combined (snow and ice) volume scattering.

4.2 Analysis of model sensitivity

The prediction of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ from the proposed two-layer-plus-volume model by Eq. (16) requires seven parameters: extinction coefficients σ₁ and σ₂ (dB m⁻¹), layer-to-volume ratios m₁ and m₂, layer positions z₁ and z₂ (m), and weight parameter α. However, it is impossible to estimate all unknowns based on only two observables: the phase and magnitude of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}$. Therefore, the necessary simplification in terms of model parameters should be considered for the model inversion, which is addressed in Sects. 4.4 and 5.

Figure 8 Simulation of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ for the proposed model by varying (a) snow extinction σ₁, (b) ice extinction σ₂, (c) weight parameter α, (d) top-layer position z₁, (e) top-layer layer-to-volume ratio m₁, and (f) bottom-layer layer-to-volume ratio m₂.

The simulation in this section aims at reducing the number of unknowns of the model by selecting the parameters which induce minor variance of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$. The sensitivity of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ to various parameters is presented in Fig. 8, where the radius and angular rotation corresponds to the coherence magnitude and phase, respectively. The phase can be translated to height via Eq. (4). It shows the complex ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ as a function of the ice-volume height ($h_{\mathrm{v}}=z_{\mathrm{1}}-z_{\mathrm{2}}$, ranging from 0 to −5 m) by varying only one parameter and keeping the others constant. The κ_z for the studied image is 0.28 rad m⁻¹ following Eq. (9).

Figure 8a and b show the loci obtained for α=0.5, m₁=0.5, m₂=0.5, and $z_{\mathrm{1}}=-\mathrm{0.15}$ m with increasing σ₁ and σ₂, respectively. The snow extinction coefficient σ₁ depends on the electromagnetic wave's frequency, snow temperature, volumetric water content, snow density, and the shape of the ice particles. At 10 GHz frequency, the snow extinction coefficient was measured to be 1–10 dB m⁻¹ (Haykin et al., 1994). This range of values is considered for σ₁ in Fig. 8a. Sea ice consists of pure ice, brine inclusions, and air bubbles. The properties and geometry of these constituents together with the environmental conditions influence the sea ice extinction coefficient σ₂. Experimental values of sea ice extinction coefficient at 10 GHz are given in (Hallikainen and Winebrenner, 1992), ranging from 10 to 200 dB m⁻¹ and covering different types of sea ice (cf. Fig. 8b). Figure 8c shows the loci obtained for m₁=0.5, m₂=0.5, σ₁=2 dB m⁻¹, and σ₂=20 dB m⁻¹, with α varying from 0 to 1.

As illustrated in Fig. 8a–c, the simulated ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ is marginally sensitive to the variance of σ₁, σ₂, and α, suggesting the possibility to fix them as constant to reduce the model complexity. This marginal sensitivity of both volume contributions can be understood by looking at their individual contributions to Eq. (16). The complex coherence of the snow volume ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}})$ can be calculated by Eq. (15) with thickness z₀₁=15 cm, and its magnitude and phase can be denoted as $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}\right({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}}\left)\right|$ and $\mathrm{\angle }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}})$, respectively. Then, the phase center location of the snow volume alone can be calculated by Eq. (10). Across the range of σ₁ (i.e., 1–10 dB m⁻¹), the snow volume has an individual coherence magnitude (i.e., $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}\right({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}}\left)\right|$) close to unity and phase center height varying from −6 to −7 cm. Therefore, it acts almost as a constant contribution, which is additionally located close to the Dirac delta of the top layer $m_{\mathrm{1}}{e}^{i{\mathit{\varphi }}_{\mathrm{1}}}$ in the complex unit circle. Similarly, the ice volume ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{2}},z_{\mathrm{12}})$ has an individual coherence magnitude of almost unity and a phase center height between −15 and −33 cm for the investigated range of ice extinction coefficients. Therefore, its effects have a limited variability and are also more or less aligned with the top layer. These observations are only valid for the investigated κ_{z_vol} and might differ for baselines larger than usual for TanDEM-X.

Figure 8d shows the loci obtained for α=0.5, m₁=0.5, m₂=0.5, σ₁=2 dB m⁻¹, and σ₂=20 dB m⁻¹, with snow depth z₁ varying from 0 to −0.5 m. The influence of the snow depth on ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ is not negligible. One way to address this is by using a priori knowledge from external sources. With fixed values of α, σ₁, σ₂, and z₁, the loci with different m₁ and m₂ values are illustrated in Fig. 8e and f, respectively. The simulated ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ shows sensitivity to the layer-to-volume ratio of the top layer in Fig. 8e. However, estimations of m₁ from observations are challenging due to the insufficient measurements of the sea ice conditions over the study area. Therefore, m₁ is approximated to be a constant value in the proposed theoretical model. A simplified model, which avoids estimating m₁, will be introduced in Sect. 4.4. For the layer-to-volume ratio of the bottom layer, m₂ induces significant variance to ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$, indicating m₂ as the most deterministic parameter which should be properly estimated to ensure the accuracy of model inversion.

4.3 Assessment of model accuracy

The observed interferometric coherence ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}$ can be biased by a residual non-volumetric decorrelation component γ_res, even after accounting for γ_s and γ_SNR by means of Eq. (6). γ_res can induce further errors when performing the model inversion for height estimation (Kugler et al., 2015). Therefore, based on Eq. (16), this potential error term is considered as

$$\begin{array}{}\text{(17)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}=e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}{\mathit{\gamma }}_{\mathrm{res}}.\end{array}$$

In the inversion, the height uncertainty depends on the magnitude of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}$ (i.e., $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}\right|$) and the InSAR baseline configuration (i.e., κ_{z_vol}). In this subsection, a Monte Carlo simulation is performed to assess the height uncertainty with various $\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}\right|$ and κ_{z_vol} values. The estimation of coherence itself has a variance due to its stochastic nature and the number of looks (i.e., the size of the coherence estimation window). In other words, the estimation accuracy of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}$ depends on the standard deviation of its magnitude and phase, which are defined by the statistical distribution and the number of looks for multi-look SAR data (Kugler et al., 2015; Touzi and Lopes, 1996; Lopes et al., 1992). The statistical distribution of coherence magnitude and phase can be given as follows.

The probability density function (pdf) of coherence magnitude γ is obtained as (Touzi and Lopes, 1996)

$$\begin{array}{}\text{(18)}& \begin{array}{rl}P\left(\mathit{\gamma }\right)=& \mathrm{2}(N-\mathrm{1}){\left(\mathrm{1}-D^{\mathrm{2}}\right)}^N\mathit{\gamma }{\left(\mathrm{1}-{\mathit{\gamma }}^{\mathrm{2}}\right)}^{N-\mathrm{2}}\\ & F\left(N,N;\mathrm{1},D^{\mathrm{2}}{\mathit{\gamma }}^{\mathrm{2}}\right),\end{array}\end{array}$$

where N is the number of looks, F is a hypergeometric function, and D is the expectation value of coherence level.

The pdf of sample coherence phase ϕ follows (Lopes et al., 1992)

$$\begin{array}{}\text{(19)}& \begin{array}{rl}P\left(\mathit{\varphi }\right)& ={\displaystyle \frac{{\left(\mathrm{1}-D^{\mathrm{2}}\right)}^N}{\mathrm{2}\mathit{\pi }}}{[}_{\mathrm{3}}{F}_{\mathrm{2}}\left(\mathrm{1},N,N;\mathrm{0.5},N;D^{\mathrm{2}}{\cos }^{\mathrm{2}}(\mathit{\varphi }-\mathit{\beta })\right)\\ & +k^{\prime }D\cos (\mathit{\varphi }-\mathit{\beta }){\times }_{\mathrm{3}}{F}_{\mathrm{2}}(\mathrm{1.5},N+\mathrm{0.5},N+\mathrm{0.5};\mathrm{1.5},\\ & N+\mathrm{0.5};D^{\mathrm{2}}{\cos }^{\mathrm{2}}(\mathit{\varphi }-\mathit{\beta })\left)\right],\end{array}\end{array}$$

where $k^{\prime }=\mathrm{\Gamma }\left(\mathrm{0.5}\right)\mathrm{\Gamma }(N+\mathrm{0.5})/\mathrm{\Gamma }\left(N\right)$, ₃F₂ is a generalized hypergeometric function, and β is the mean phase difference.

The simulation is a four-step procedure. First, the complex value of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ is calculated for the designed two-layer-plus-volume model with specific parameters (σ₁, σ₂, α, m₁, m₂, z₁, z₂) and a given κ_{z_vol}. The surface phase ϕ₀ is assumed to be 0. Second, ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{sim}}$ is obtained by ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{sim}}={\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}{\mathit{\gamma }}_{\mathrm{res}}$ with γ_res=0.98 (according to Kugler et al., 2015). Next, a set of (N_s=10 000) ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{sim}}_i}$ complex samples is generated via Eq. (18) (for magnitude) and Eq. (19) (for phase) with $D=\left|{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{sim}}\right|$ and $\mathit{\beta }=\mathrm{\angle }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{sim}}$. Finally, for each simulated ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{sim}}_i}$, the volume height h_vi is estimated by the inversion of Eq. (16) with the specific parameters and κ_{z_vol} of the simulation and compared with the input $h_{\mathrm{v}}=z_{\mathrm{1}}-z_{\mathrm{2}}$ in the first step. The bias $B_{{\mathrm{\Delta }}_h}=|\mathbb{E}(h_{\mathrm{v}i})-h_{\mathrm{v}}|$ and the standard deviation ${\mathit{\sigma }}_{{\mathrm{\Delta }}_h}=\sqrt{\mathbb{E}\left[\right(h_{\mathrm{v}i}-h_{\mathrm{v}}{)}^{\mathrm{2}}]}$ are calculated to quantify the estimation accuracy of the model.

Figure 9 (a–c) Relative bias $B_{{\mathit{\delta }}_{\mathrm{h}}}/h_{\mathrm{v}}$ and (d–f) relative standard deviation ${\mathit{\sigma }}_{{\mathit{\delta }}_{\mathrm{h}}}/h_{\mathrm{v}}$ of the obtained ice volume h_v for the proposed model. Fixed parameters: (a–f) σ₁=2 dB m⁻¹, m₁=0.3, $z_{\mathrm{1}}=-\mathrm{0.15}$ m, and α=0.5; variable parameters: (a, d) σ₂=20 dB m⁻¹, m₂=0.5, (b, e) σ₂=20 dB m⁻¹, m₂=2, and (c, f) σ₂=100 dB m⁻¹, m₂=0.5. Blue line indicates the threshold of $B_{{\mathit{\delta }}_{\mathrm{h}}}/h_{\mathrm{v}}=\mathrm{25}$ % or ${\mathit{\sigma }}_{{\mathit{\delta }}_h}/h_{\mathrm{v}}=\mathrm{25}$ %.

By fixing the model parameters (σ₁, σ₂, α, m₁, m₂, z₁) and varying z₂, $B_{{\mathrm{\Delta }}_{\mathrm{h}}}$ and ${\mathit{\sigma }}_{{\mathrm{\Delta }}_{\mathrm{h}}}$ are functions of ice-volume height $h_v=z_{\mathrm{1}}-z_{\mathrm{2}}$ and volume-corrected vertical wavenumber κ_{z_vol}. The simulation procedure is performed for h_v ranging from 0 to 5 m and κ_{z_vol} ranging from 0.02 to 0.4 rad m⁻¹. The number of looks N is set to be the same value as the experimental data. The bias and the standard deviation relative to volume height are shown in Fig. 9. The plots illustrate that for a specific baseline geometry with a given κ_{z_vol}, the model performance is superior for a certain range of volume height, shown in the blue curve, indicating the 25 % threshold, in Fig. 9. For volume heights lower than at the blue curve, the bias and variance are larger than 25 %, leading to a lower precision of model inversion. With different σ₂ values, there are no obvious distinctions of $B_{{\mathrm{\Delta }}_{\mathrm{h}}}/h_{\mathrm{v}}$ between Fig. 9a and c, as well as ${\mathit{\sigma }}_{{\mathrm{\Delta }}_{\mathrm{h}}}/h_{\mathrm{v}}$ between Fig. 9d and f, respectively. It suggests that the model accuracy is marginally sensitive to ice extinction coefficient σ₂. The layer-to-volume ratio of the bottom layer m₂ plays a key role in model accuracy. From Fig. 9b and e, the $B_{{\mathrm{\Delta }}_{\mathrm{h}}}/h_{\mathrm{v}}$ and ${\mathit{\sigma }}_{{\mathrm{\Delta }}_{\mathrm{h}}}/h_{\mathrm{v}}$ are smaller than those from Fig. 9a and d due to the larger m₂, indicating the higher model accuracy in this case. The 25 % error threshold provides a criteria for selecting the best baseline geometry for the application. For example, with a specific parameter set as Fig. 9b and e, the κ_{z_vol} needs to be larger than 0.40 rad m⁻¹ to ensure an effective inversion for ice-volume thickness less than 0.85 m. Since the κ_{z_vol} of the studied SAR image is 0.28 rad m⁻¹, in order to achieve an 25 % error accuracy, the ice volume needs to be thicker than a certain value depending on m₂. This certain value ranges from 1.1 to 2.7 m for m₂ being 0.5–2; see Fig. 9d and e. The above assessment also indicates the potential of applying the proposed model to achieve a more accurate result using a larger baseline configuration. Note that the baseline should not be too large since it results in stronger interferometric decorrelation which contaminates the topographic information.

4.4 Model simplification

The proposed theoretical model ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ given in Eq. (16) contains seven parameters: σ₁, σ₂, α, m₁, m₂, z₁, and z₂, requiring necessary a priori knowledge of the test site. However, such a priori knowledge is scarce due to the sparse ground measurement of the Antarctic sea ice, therefore impeding the practical application of the proposed model. As described in Sect. 4.2, the contributions of the snow volume ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}})$ and the ice volume ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}({\mathit{\sigma }}_{\mathrm{2}},z_{\mathrm{12}})$ to the theoretical model ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ show a limited sensitivity for the TanDEM-X acquisition geometry within the investigated range of extinction coefficients σ₁ and σ₂. Additionally, their individual coherence loci are located close to the Dirac delta of the top layer $m_{\mathrm{1}}{e}^{i{\mathit{\varphi }}_{\mathrm{1}}}$ in the unit circle. Therefore, the theoretical model can be approximated by merging the contributions of the snow volume, the ice volume, and the top layer into one Dirac delta. This simplified model can be given as

$$\begin{array}{}\text{(20)}& \begin{array}{rl}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}=\\ & e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\displaystyle \frac{\mathit{\alpha }{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}\left({\mathit{\sigma }}_{\mathrm{1}},z_{\mathrm{01}}\right)+e^{i{\mathit{\varphi }}_{\mathrm{1}}}(\mathrm{1}-\mathit{\alpha }){\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{v}}\left({\mathit{\sigma }}_{\mathrm{2}},z_{\mathrm{12}}\right)+m_{\mathrm{1}}{e}^{i{\mathit{\varphi }}_{\mathrm{1}}}+m_{\mathrm{2}}{e}^{i{\mathit{\varphi }}_{\mathrm{2}}}}{\mathrm{1}+m_{\mathrm{1}}+m_{\mathrm{2}}}}\\ & \approx e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\displaystyle \frac{\mathrm{1}\cdot e^{i{\mathit{\varphi }}_{\mathrm{1}}}+m\cdot e^{i{\mathit{\varphi }}_{\mathrm{2}}}}{\mathrm{1}+m}}\\ & =e^{i{\mathit{\varphi }}_{\mathrm{0}}}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{S}}\left(m,z_{\mathrm{1}},z_{\mathrm{2}}\right),\end{array}\end{array}$$

where ϕ₁=κ_{z_vol}z₁, ϕ₂=κ_{z_vol}z₂, z₁ and z₂ are the position of the top layer and the bottom layer, respectively, and m is the layer-to-layer ratio.

Compared to the theoretical model in Eq. (16), the simplified model in Eq. (20) only has three parameters, remarkably improving the applicability in practice. Figure 10 illustrates the simulations of ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ and ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{S}}$ according to Eqs. (16) and (20), respectively. For ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$, the parameters are set as σ₁=2 dB m⁻¹, σ₂=20 dB m⁻¹, $z_{\mathrm{1}}=-\mathrm{0.15}$ m, and $m_{\mathrm{1}}=m_{\mathrm{2}}=\mathrm{0.5}$. For ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{S}}$, z₁ is also set to be −0.15 m, and m varies from 0.3 to 0.4. As we can see, for both coherence magnitude and phase, the simplified model can achieve comparable results to the theoretical model by assuming appropriate m values.

Figure 10 Comparison of complex coherence ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}}$ from the theoretical model (${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$, thick black line) and the simplified model (${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{S}}$, colored lines). (a) Magnitude of the modeled coherence. (b) Phase (^∘) of the modeled coherence.

5 Model inversion

In order to apply the simplified model and theoretical model to geophysical parameter retrieval, a methodology is developed for the model inversion. The objective is to estimate the topographic phase ϕ₀ and thus generate the sea ice topographic height with snow depth for the whole SAR image. The model inversion includes three main steps, illustrated in Fig. 11. As explained in Sect. 2.1, in order to select the ice that is deformed and thick without seawater flooding, the samples with height above 0.8 m, which are 83 % of the coregistered data set, are selected.

As shown in step 1 in Fig. 11, the h_DMS is converted to phase ϕ_DMS via Eq. (4) and used as a priori knowledge. For the simplified model, z₁ is set to be −0.18 m according to the snow depth provided by the AMSR-E/AMSR2 Unified Level-3 Daily data set (Meier et al., 2018). The AMSR-E/AMSR2 data set provides snow depth over sea ice as 5 d running averages. For the studied area on the campaign date, the averaged snow depth was measured to be 18 cm. For the theoretical model, as discussed in Sect. 4.2, since the simulated ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{mod}\mathrm{\_}\mathrm{T}}$ shows marginal sensitivity to the variance of snow layer extinction coefficient σ₁, ice layer extinction coefficient σ₂, and weight parameter α, these three parameters are fixed to constants. For the snow-covered sea ice of the studied area, the σ₁ and σ₂ are assumed to be 2 and 20 dB m⁻¹, respectively, referring to experimental values (Cox and Weeks, 1974; Hallikainen and Winebrenner, 1992). The snow depth z₁ is also set to −0.18 m. The layer-to-volume ratio of the top layer m₁ is set to 0.3. With the above specific parameters, the m₂ (also m for the simplified model) values can be derived by the inversion of the proposed model according to Eqs. (16) or (20) for the theoretical and the simplified model, respectively.

Figure 11 Flow chart of the proposed inversion method.

Figure 12 The relation between the bottom-layer layer-to-volume ratio m₂ (respectively m) and the co-polarization (coPol) coherence γ_coPol for the theoretical model (a–d) and the simplified model (e–h). (a, d) HH polarization, (b, f) VV polarization, (c, g) Pauli 1 (HH + VV polarization), and (d, h) Pauli 2 (HH − VV polarization). Note the different y-axis scaling for m₂ and m.

Next, since m₂ or m is the most deterministic parameter in the respective models, the aim is to estimate m₂ or m from the SAR observations (step 2 in Fig. 11). In addition to interferometry, which provides topographic information, polarimetry reveals information on the scattering processes and is a useful tool to characterize sea ice properties. Among several polarimetric signatures, the co-polarization (coPol) coherence γ_coPol is a measurement of the degree of electromagnetic wave depolarization between HH and VV polarizations caused by both the rough surface scattering and the volume scattering (Kasilingam et al., 2001). γ_coPol was demonstrated to be a crucial signature in sea ice characterization (Kim et al., 2011; Wakabayashi et al., 2004; Huang and Hajnsek, 2021). γ_coPol can be calculated as (Lee and Pottier, 2009)

$$\begin{array}{}\text{(21)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\mathrm{coPol}}={\mathit{\gamma }}_{\mathrm{coPol}}\cdot e^{i{\mathit{\varphi }}_{\mathrm{coPol}}}={\displaystyle \frac{<s_{\mathrm{VV}}{s}_{\mathrm{HH}}^{*}>}{\sqrt{<s_{\mathrm{VV}}{s}_{\mathrm{VV}}^{*}><s_{\mathrm{HH}}{s}_{\mathrm{HH}}^{*}>}}},\end{array}$$

where s_HH and s_VV are single-look complex images in HH and VV polarization, respectively. The symbol $<.>$ denotes an ensemble average. A 4×12 window in azimuth and slant range is applied to estimate γ_coPol. It is found that m₂ (also m) is inversely related to γ_coPol for the four polarizations, shown in Fig. 12, enlightening us to derive an empirical function between the parameter m₂ (also m) and γ_coPol. As shown in step 2 in Fig. 11, the linear functions for the different polarizations are derived by least-squares fitting and are detailed in Fig. 12. Note the almost identical correlation coefficients of the fitted linear function for the theoretical model (first row) and the simplified model (second row) in Fig. 12. This underlines that the theoretical model can be approximated by the simplified model. The parameters of the linear functions for m₂ and m are of course different. Then, the fitted functions are applied to estimate ${\widehat{m}}_{\mathrm{2}}$ (also $\widehat{m}$), which will be used as an input to perform the model inversion for the whole image, including the area without DMS measurements.

Finally, as shown in the step 3 in Fig. 11, for the TanDEM-X data without a priori knowledge of DMS measurement, the γ_coPol together with the derived linear function is utilized to estimate ${\widehat{m}}_{\mathrm{2}}$ (also $\widehat{m}$) for each pixel. With the specific parameter setting, the estimated ${\widehat{m}}_{\mathrm{2}}$ (also $\widehat{m}$), and the ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{{\mathrm{InSAR}}^{\prime }}$ from InSAR pairs, the topographic phase ϕ_{0_est} can be retrieved by solving Eqs. (16) or (20) and then converted to height h_mod in meters via Eq. (4). The area overlaid by the DMS flight track is used to verify the model-inversion result quantitatively and visually, which will be illustrated in the next section.

6 Experimental results

In this section, the proposed two-layer-plus-volume model and its simplified version are inverted to estimate sea ice topography following the developed scheme. Note that the retrieved sea ice topography refers to the sea ice height including the snow depth above the local sea level. Both visual and quantitative analyses are given to evaluate the retrieval performance.

6.1 Retrieval performance of the simplified model

The sea ice topographic retrievals based on the simplified model (Fig. 11) are performed for the four polarizations (HH, VV, Pauli 1, and Pauli 2), respectively. Because of the marginal visual distinction among HH, VV, and Pauli-1 polarizations, only the Pauli-1 polarization result is presented for conciseness. The quantitative evaluation is given for the four polarizations.

The model-retrieved sea ice topography in Pauli-1 polarization is shown in Fig. 13. The strip between the grey lines is the area covered by the DMS DEM, which is superimposed on the model-retrieved result (h_{mod_S}) with the same color map. In general, the retrieved height varies from 0.8 to 3 m across the whole image, showing a good agreement with the height range obtained by DMS measurements. Three areas are selected and enlarged for detailed analyses. Each area contains 40×500 pixels (corresponding to 400×5000 m area) along the range and azimuth direction, respectively. Figure 13b is the close-up of area 1, where several sea ice areas are higher than 2.5 m. The model-retrieved height (outside the grey lines) shows good continuity with the DMS measurements (between the grey lines), indicating the effectiveness of the proposed method. Area 2 is mainly covered with ice lower than 2 m. The sea ice topographic retrieval of area 2 is shown in Fig. 13c, where the consistency between the model-retrieved Pol-InSAR DEM and the DMS DEM is again visually verified. Area 3 (Fig. 13d), including sea ice in the range of 2–2.5 m, shows the preservation of continuous sea ice features as well. Besides, the relative retrieval bias ϵ, which can be calculated as $\mathit{ϵ}=|h_{\mathrm{mod}\mathrm{\_}\mathrm{S}}-h_{\mathrm{DMS}}|/h_{\mathrm{DMS}}$, is used to quantify the retrieval accuracy. In Fig. 13e–g, ϵ values over area 1–3 are below 25 % for most parts, whereas only a few parts, often near the masked-out regions (transparent pixels), present higher ϵ values. Note that the masked-out regions refer to water and thinner ice areas with height less than 0.8 m. The averaged ϵ values over area 1–3 are about 19 %, 14 %, and 15 %, respectively, and the averaged ϵ value is 18 % for the whole image, achieving the theoretical 25 % error accuracy derived in Sect. 4.3. Therefore, from the comparison with photogrammetric measurements, the model-retrieved Pol-InSAR DEM demonstrates a good visual and quantitative agreement with the DMS DEM.

Figure 13 Sea ice topographic retrieval with the simplified model (h_{mod_S}). The transect from the DMS DEM is plotted between grey lines. Note that the heights below 0.8 m are set to be transparent. (a) The whole studied SAR image. (b–d) Close-up of areas 1–3. (e–g) Relative retrieval bias $\mathit{ϵ}=|h_{\mathrm{mod}\mathrm{\_}\mathrm{S}}-h_{\mathrm{DMS}}|/h_{\mathrm{DMS}}$ of area 1–3.

Figure 14 (a) Sea ice height profiles from DMS measurement (blue) and model inversion (red). Each profile represents the height along a 1×5000-pixel section at the center of the coregistered segment. (b) Height difference between the DMS measurement and the simplified model-derived height (blue), theoretical model-derived height (red), or original InSAR-derived height (green). The mis-coregistered and h_DMS-below-0.8 m samples are excluded from the plots. The black line plots 0 m as a reference.

The comparison of the height profiles in Pauli-1 polarization along the DMS strip is shown in Fig. 14a. The simplified-model-retrieved DEM h_{mod_S} and the DMS DEM h_DMS along the transect are plotted in the red and blue line, respectively. As observed from Fig. 14a, the model-retrieved height has, in general, a good capture of the topographic variation. The height difference (h_{diff_S}) between the h_DMS and the simplified model-retrieved height is shown in Fig. 14b (blue). Compared with the larger height difference (h_{diff_ori}) between the h_DMS and the initial InSAR DEM in Fig. 14b (green), the surface-elevation bias is properly compensated for by using the simplified model. In Pauli-1 polarization, the RMSE of the model-retrieved height is ∼0.26 m relative to the DMS DEM. Compared to the original InSAR-derived height with RMSE of ∼1.10 m, Fig. 14 reveals the pronounced improvement of applying the proposed model to estimate sea ice topographic height, considering that the RMSE in the DMS DEM is already 0.2 m to start with (Dotson and Arvesen., 2014).

Table 1 summarizes the performances between the retrieved height from the simplified model and h_DMS for the four polarizations. The RMSE values across the HH, VV, and Pauli-1 polarizations are similar, ranging from 0.2637 to 0.2764 m. The larger RMSE value in the Pauli-2 channel likely results from the lower SNR values. The polarization-independent performances in HH, VV, and Pauli 1 further reveal the inherent property of the studied sea ice to be a polarization-independent volume among co-polarization channels in X-band radar frequencies.

Table 1 The RMSE of model-retrieved height for the four polarizations.

6.2 Retrieval performance of the theoretical model

For the theoretical model which requires necessary a priori knowledge (e.g., snow density, ice salinity, and air temperature) to determine the input parameters, the inversion performance depends on the study area and the sea ice structural characteristics. In this subsection, the inversion of the theoretical model is achieved by fixing parameters which have marginal effects on the model predictions and by estimating the layer-to-volume ratio of the bottom layer from a polarimetric signature: The coPol coherence. With the specific parameter set assumed in Sect. 5, the theoretical model is inverted according to the proposed method (Fig. 11) to retrieve the sea ice topography. The height difference (h_{diff_T}) between the h_DMS and the theoretical model-retrieved height is shown in Fig. 14b (red), visualizing the similar performance of the theoretical and the simplified model. The results in Table 1 also show that the retrieval accuracy in terms of RMSE of the theoretical model is almost identical to the simplified model. It demonstrates that the theoretical model can adequately correct the penetration bias of InSAR signals and achieve an effective sea ice height retrieval from dual-polarization single-pass interferometric data, fulfilling the primary goal of this study.

The comparable performance of the theoretical and simplified model convinces the effectiveness of employing the simplified model to achieve an accurate sea ice topographic retrieval. In the cases when the ground measurements are sparse, the simplified model requires only one parameter (i.e., snow depth over sea ice), significantly reducing the model complexity and improving the applicability in practice.

Figure 15 Yellow area: sea ice height profiles from the simplified model in Pauli-1 polarization with z₁ from −0.05 to −0.75 m. Blue dash line: sea ice height profiles from DMS measurements. Each profile represents the height along a 1×5000-pixel section at the center of the coregistered segment. The mis-coregistered and h_DMS-below-0.8 m samples are excluded from the plot.

7 Discussion

7.1 Model complexity and observation space

The proposed two-layer-plus-volume model includes seven parameters with the assumption of their independence of polarization. In the case of a polarization-dependent volume, two more parameters are introduced for each volume and thus further complicate the inversion. In order to achieve the model inversion, one method is to develop a simplified model as presented in Sect. 6.1, which accurately approximates the behavior of the theoretical model and requires only the snow depth as an input parameter. Another method to achieve the theoretical model inversion is to increase the observation space to full polarization and/or multi-baseline configurations. Acquisitions of full-polarization data improve the inversion capability over single-polarization or dual-polarization configurations. For instance, dual-baseline quad-polarization data provide 12 independent observables and thereby offer an opportunity to theoretically invert a model with a maximum of 12 parameters. In Sect. 4.3, we have illustrated the theoretical performance of the proposed model with various baseline configurations and obtained a certain range which can ensure high inversion accuracy. It reveals the potential to establish an inversion scheme by combining observations from a range of different κ_z (i.e., the vertical wavenumber in free space), where larger values of the effective perpendicular baseline b_⟂, corresponding to larger κ_z values, are expected to improve height retrieval accuracy. With quad-polarization and multi-baseline data acquired over sea ice in the future, developing a refined inversion scheme for more diverse scattering scenarios and thinner sea ice heights will be promising.

7.2 Influence of snow depth on experimental result

In Sect. 4.2, we demonstrated that the influence of snow depth on the simulated coherences is not negligible and stated that external data of snow measurements should be used in the model. In this study, snow depth is assumed to be invariant across the scene due to the limited spatial resolution of available snow measurements. Therefore, a constant value of $z_{\mathrm{1}}=-\mathrm{0.18}$ is used in the retrieval. Actually, snow on sea ice undergoes temporal- and spatial-variant processes and is strongly coupled with atmospheric, oceanic, and ice conditions. Thus, a single value is not representative of the actual spatial snow depth distribution. In order to assess the impact of the snow depth on the experimental results, we perform the whole inversion scheme with various inputs of z₁. During September and November, the snow depth on Antarctic sea ice is reported to be maximum ∼1 m and mainly 0–0.8 m (Webster et al., 2018). Therefore, z₁ values ranging from −0.05 to −0.75 m are selected. For each pixel, we retrieve heights using this range of z₁ values, shown as the yellow area in Fig. 15. Δh_{mod_S} is defined as the difference between the maximum and the minimum retrieved height of every pixel. The distribution of Δh_{mod_S} along the transect is presented in Fig. 16, where Δh_{mod_S} has a range of 0.07–1.09 m with an average of 0.31 m, indicating the fluctuation of model-retrieved height by using different snow depths. This analysis with various snow depth assumptions can help to constrain possible model-retrieved topographies, and Δh_{mod_S} can be a quantitative indicator for the uncertainty of the retrieved height in the absence of high-resolution snow depth data.

Figure 16 The distribution of Δh_{mod_S} along the transect.

7.3 Model extension to other ice conditions

The proposed model was proven to be effective in a specific area covered by thick and deformed ice with snow cover in the western Weddell Sea. The extension of the proposed model to other ice types under different environmental conditions needs further research and suitable data.

In order to apply the model over younger and thinner sea ice, the first challenge is the severe misregistration between SAR images and reference measurements due to the stronger dynamics of thinner sea ice. Reduced SAR backscattering intensity corresponding to thinner and smoother sea ice further complicates the data coregistration. Besides, the achievable height sensitivity for thin ice is also a major limitation of InSAR/Pol-InSAR-derived sea ice DEMs with current SAR systems. In this study, the proposed method can achieve sea ice topographic retrieval with an RMSE of 0.26 m for thick and deformed ice; however, this accuracy is insufficient for thinner ice whose height above sea level is only tens of centimeters or even less. Moreover, an additional volume, i.e., snow ice formed by flooding, should be considered when extending the proposed model to a thinner ice area. The presence of snow ice is a challenging retrieval scenario and not only for Pol-InSAR applications. The zero freeboard of snow ice, or the negative freeboard of ice in case of a slush cover, can cause significant difficulties when converting radar or laser altimeter measurements of snow or ice freeboard to ice thickness because the isostatic equations fail in such a scenario. Past studies showed that snow ice contributes an average of 8 % of the total volume in the Weddell Sea (Lange et al., 1990). A greater amount of snow ice, which accounts for 12 %–36 % of the total mass, was reported in the Ross, Amundsen, and Bellingshausen seas (Jeffries et al., 2001). Although the snow ice has a higher salinity than the ice below, there could be still some penetration into the ice volume below. A three-layer model, which includes snow, snow ice, and ice layers, could be feasible to correct the InSAR phase center and retrieve surface height for thinner ice in the Antarctic (as illustrated in Fig. 1a). However, the three-layer model involving more parameters than the proposed theoretical model can only be inverted by increasing the observation space to full polarization and multi-baseline configurations. In the case of single-baseline or dual-polarization configurations where observables are limited, an intelligent model with fewer parameters is worthy of investigation in the future.

In order to assess the transferability of the proposed model to the Arctic regions, further validations, including coregistered SAR images and topographic reference (e.g., optical/lidar measurements), are needed considering the significant difference of ice and snow properties between the Arctic and Antarctic. Ancillary measurements (e.g., snow depth, temperature, and salinities of ice and snow) at a wide range of ice conditions in both polar regions are crucial to understand the properties of various typologies of ice and therefore are valuable for extending the model to general applicability. Part of the ancillary data (i.e., snow depth and ice freeboard height) would be available in OTASC Level-4 products in the future, offering us an opportunity to interpret sea ice electromagnetic properties.

Given the sea ice scenario in this study, we assume that the snow condition over the thick and deformed ice is dry, and the X-band microwaves penetrate both the snow and ice layer. For wet-snow-covered sea ice, the penetration capability of X-band is limited (Hallikainen and Winebrenner, 1992), and therefore the proposed approach in this study cannot be applied. In the case of wet snow, as well as medium wet snow, the snow–air interface influenced by the surface roughness needs to be taken into account as it changes the X-band SAR backscatter (Nandan et al., 2016; Dufour-Beauséjour et al., 2020).

In this study, the InSAR pair in StripMap mode covers 19 km × 50 km in the SAR ground range and azimuth direction, respectively, providing a unique three-dimensional (3D) topographic map rather than a narrow transect from lidar or photogrammetric measurements. TanDEM-X has a regular revisit cycle of 11 d over the Arctic and a larger revisit time due to the particular satellite position configuration required over Antarctica. The current SAR satellites, such as X-band TanDEM-X and COSMO-SkyMed, C-band Sentinel-1 and RADARSAT Constellation, as well as the future X-band LOTUSat, L/S-band NISAR, and L-band ROSE-L, will together achieve a long-term sea ice topographic monitoring in both polar regions. Synergistic use of different SAR satellites offers more extensive spatial coverage and shorter revisit times than a single platform. In the future, the joint use of multi-frequency SAR imagery could develop a better understanding of sea ice properties and processes (Dierking and Davidson, 2021), which would be indispensable for retrieving sea ice topography at a more comprehensive range of ice conditions.

8 Conclusions

In this study, the potential to retrieve sea ice topography with the Pol-InSAR technique was validated with single-pass interferometric SAR data and airborne photogrammetric measurements over the thick (>2 m) and deformed sea ice with snow cover in the western Weddell Sea. The DMS DEM reveals that the sea ice topography along the flight track varies from 0 to 2.68 m, with the average height being 1.27 m. The average elevation difference between the conventional InSAR DEM and the DMS DEM is ∼1 m in the four investigated polarizations (HH, VV, Pauli 1, and Pauli 2), suggesting the demand for a valid method to obtain sea ice topography and to correct for the penetration of the microwave signals into the sea ice. By exploiting the interferometric coherence, a two-layer-plus-volume model was proposed to characterize the sea ice vertical scattering structure, and an inversion scheme was developed for height retrieval. The model's theoretical accuracy was assessed for various vertical wavenumber values to ensure 25 % error accuracy at the employed baseline configuration with κ_z=0.28 rad m⁻¹. The assessment of the model's theoretical accuracy showed the potential to apply the model to multi-baseline configurations, giving the ability to adjust a sensor to the particular type of sea ice. For instance, a configuration with κ_z=0.40 rad m⁻¹ ensures an effective inversion for ice-volume thickness less than 0.85 m.

The proposed theoretical model requires seven input parameters depending on the environmental conditions over the test site, which are unavailable in many practical applications. In order to reduce the model complexity and improve the model applicability, a simplified model requiring only the input of snow depth was proposed based on the analyses of the model sensitivity to different parameterizations. For the theoretical and simplified model, the layer-to-volume ratio of the bottom layer and the layer-to-layer ratio were observed to be inversely correlated to an essential polarimetric signature: the coPol coherence. This relationship was exploited in the inversion scheme by estimating those parameters from the coPol coherence with a fitted linear function. Note that the proposed models and inversion scheme in this study were developed over a specific area containing thick and deformed ice; therefore, they cannot be directly applied to sea ice areas covering various ice types. With more coregistered data acquired from SAR and reference DEMs in the future, extension of the model to regions covered by different sea ice types is worthy of further investigation.

The effectiveness of both the theoretical and simplified models and the proposed inversion scheme were verified with the DMS measurements for the sea ice height above 0.8 m, corresponding to the thick (>2 m) and deformed ice without flooding. The model-retrieved sea ice topography achieved a RMSE as low as 0.26 m, which is significantly better than the RMSE of 1.10 m of the conventional InSAR DEM. This indicates the capability to correct for the microwave signal penetration and to generate a precise wide-swath topographic map from dual-polarization single-pass InSAR data. The polarization-independent volume property of sea ice in the co-polarization channels in X-band radar frequency, which was concluded from the similar retrieval performance across HH, VV, and Pauli-1 polarizations, gave insights to develop superior models for height retrieval in the future. Future work will include an investigation of the possibility of sea ice topographic retrieval for various types of sea ice, such as thin ice and newly formed ice.