Characterization of ice features in the southwest Greenland ablation zone using multi-modal SAR data

This study investigates ice features of unknown glaciological origin in the ablation zone of southwest Greenland, focusing on the land-terminating Russell Glacier. Using data from an experimental airborne SAR (Synthetic Aperture Radar) campaign of the German Aerospace Center (DLR), the research employs a range of advanced techniques, including SAR polarimetry, interferometry, tomography, and modeling, to characterize these features. The analysis reveals that in low-backscatter areas, surface scattering is dominant with no correlation to topography or surface characteristics. In contrast, surrounding high-backscatter areas are characterized by volume scattering and the presence of subsurface scattering structures. A significant aspect of this study involves comparing the observed ice features with known glaciological phenomena in the ablation zone. However, the combined findings, along with the temporal stability of these features, as seen through annual SAR backscatter analysis, complicate a straightforward glaciological explanation. A first theory involves the presence of a weathering crust causing these low-backscatter features due to residual liquid water content. These findings could improve our understanding of surface and subsurface processes in the ablation zone, contributing to better mass balance assessments.

1 Introduction

The ablation zone of the Greenland Ice Sheet (GrIS) is one of its most dynamic regions, where ice melt drives significant meltwater production, contributing to global sea-level rise and reflecting the impacts of climate change (Box et al., 2012; Fettweis et al., 2013). Situated at the lowest elevations of the GrIS, this zone experiences continuous ice loss that exceeds accumulation. Its bare ice, minimal snow cover, and prominent crevasses make it particularly vulnerable. Understanding these dynamic processes is essential for assessing the overall GrIS mass balance and the broader implications for global climate systems (Mote, 2007; Paterson, 1994). In this context, monitoring variations in ice properties within the ablation zone is crucial. Our observations revealed a puzzling spatial pattern of low-backscatter areas in Synthetic Aperture Radar (SAR) data, suggesting corresponding variations in the ice's dielectric or structural properties. However, no preliminary glaciological explanation matched the SAR observations, rendering their investigation important to potentially enhance our understanding of surface and subsurface properties in the ablation zone.

SAR has become a valuable tool for studying and monitoring glaciers and ice sheets. Currently, monitoring of the GrIS relies heavily on space-borne X-band (8–12 GHz) and C-band (4–8 GHz) systems, such as Sentinel-1 and TanDEM-X, possibly in combination with optical remote sensing data, to create e.g. Digital Elevation Models (DEM) (Krieger et al., 2007), observe flow velocities (Joughin et al., 2017) and to classify glacier zones based on SAR backscatter (Paterson, 1994; Braun et al., 2000). Overall, these higher frequency bands (X- and C-band) primarily capture surface characteristics. In contrast, lower frequency bands such as L-band (1–2 GHz) and P-band (0.3–1 GHz) have the ability to penetrate deeper into the ice and can therefore capture subsurface characteristics or internal structures of glaciers (Dall et al., 2001). However, while we have access to space-borne data at these frequencies from systems like ALOS, their temporal and spatial coverage in the GrIS is limited (Rosenqvist et al., 2007). Only through the analysis of low-frequency data acquired during airborne SAR campaigns have aforementioned low-backscatter areas been found (Banda et al., 2016; Parrella et al., 2021). Moreover, these areas remain almost undetectable in both optical imagery and higher-frequency SAR data, making them fascinating features to investigate in a multi-frequency analysis.

Figure 1 Radar-dark (low-backscatter areas) and radar-bright (high-backscatter areas) features in the ablation zone of southwest Greenland, observed in airborne SAR at P-band with HV polarization (shown here at the Russell Glacier). In the background, a Sentinel-2 optical image provides context for the surrounding area. In the lower left corner, a Greenland map highlights the test site with a red marker.

In this study, we refer to these low-backscatter areas on the glacier as radar-dark features, which are often surrounded by high-backscatter areas, called radar-bright features. Our analysis shows that these radar-dark features exhibit distinct spatial characteristics, typically appearing as round or oval shapes, measuring 100 to 500 m in size. These features often form interconnected patterns that contrast sharply with radar-bright features, creating a clear visual distinction across a section of the Russell Glacier's ablation zone, as observed in airborne SAR data at HV polarization and P-band frequency (Fig. 1). Low-frequency, space-borne SAR data further revealed that these radar-dark and radar-bright features represent a large-scale phenomenon across the ablation zone of southwest Greenland, particularly in low-velocity, land-terminating glaciers. As this study primarily relies on remote sensing data, and lacks ground measurements, the investigation of the unknown glaciological origins of these radar-dark and radar-bright features is based on an analysis of their scattering processes and structures in comparison with known glaciological phenomena.

To analyze the scattering processes of these ice features, we apply advanced SAR techniques: Polarimetric SAR (PolSAR), Interferometric SAR (InSAR), and Tomographic SAR (TomoSAR). PolSAR differentiates scattering mechanisms by distinguishing surface, volume, and dihedral scattering (Cloude and Pottier, 1997). It is used in ice sheet monitoring to detect changes in ice properties, such as surface roughness and volume scattering caused by snow or ice layers (Rott and Davis, 1993; Parrella et al., 2021). InSAR provides critical insights into the vertical structure of ice sheets. Key parameters such as coherence and phase center height are crucial for interpreting the scattering processes of the study area (Rignot et al., 2001). Volumetric decorrelation, a major contributor to total coherence loss, can be linked to the vertical distribution of scatterers and helps to assess the subsurface structure (Bamler and Hartl, 1998; Fischer et al., 2020). TomoSAR provides 3D imaging capabilities (Reigber and Moreira, 2000), enabling detailed exploration of the ice subsurface, effectively mapping internal layers and detecting complex subsurface structures at different heights (Tebaldini et al., 2016; Banda et al., 2016). Moreover, conducting a temporal analysis is essential for determining the origin of these ice features. We can observe annual variations through ALOS backscatter data (Bolon et al., 2007; Ruan et al., 2012), enabling us to monitor the formation and evolution of radar-dark features, providing critical insights into whether these features represent stable structures or transient phenomena. A first polarimetric and tomographic SAR characterization of these radar-dark and radar-bright features was conducted in Pardini et al. (2016) focusing on the scattering properties at one example location without addressing a geophysical explanation. Building on this preliminary study, we perform a more comprehensive investigation.

After determining the scattering processes and temporal changes, we aim to gain a glaciological understanding of radar-dark and radar-bright features observed in our test site, the Russell Glacier’s ablation zone. To achieve this, we compare these ice features with known surface and subsurface processes, related to glaciological phenomena commonly found in the ablation zone. Considered surface phenomena include drained and refrozen supraglacial lakes (Hu et al., 2021; Yang et al., 2021), as well as physical impurities, such as the so-called “dark zone” and cryoconite deposits on the glacier surface (Wientjes and Oerlemans, 2010; Ryan et al., 2018). The subsurface phenomena investigated include seasonal temperature variations and englacial water bodies such as water pockets and channels (Catania et al., 2008; Lampkin and VanderBerg, 2013). We also consider the weathering crust, a near-surface subsurface phenomenon, as observed at Russell Glacier by Cooper et al. (2018). All these phenomena may reduce backscatter through absorption, attenuation, or reflection, and could potentially be related to the radar-dark features. Moreover, by focusing our investigation on data collected during stable, frozen conditions, we aim to minimize the influence of active melt and better isolate the glaciological processes responsible for radar-dark features in SAR imagery.

The paper is organized as follows: Sect. 2 describes the test site and data. Section 3 outlines the methods, focusing on the SAR techniques applied. Section 4 presents the investigation of radar-dark and radar-bright features. Section 5 provides a discussion of the results in relation to glaciological processes, and Sect. 6 concludes with a summary of findings and potential directions for future research.

2 Test Site and Data

This study examines the ablation zone of southwest Greenland, with a focus on the Russell Glacier, also known as K-Transect due to the availability of NASA’s Operation IceBridge data (MacGregor et al., 2021; Studinger et al., 2022) and long-term records from “Programme for Monitoring of the Greenland Ice Sheet” (PROMICE) weather stations (Fausto et al., 2019). Located at approximately 67.1° N and 50.0°  W, the analysed area spans elevations from 200 to 860 m above sea level for about 26 km east-west and 3 km north-south (Fig. 1). The average glacier flow velocity at the eastern edge of the test site is approximately 100 m yr⁻¹ (ENVEO, 2024), which decreases towards the terminus, as typical for land-terminating glaciers in this region (Nagler et al., 2015; van de Wal et al., 2015). Ablation rates vary annually between 2.5 to 5.6 m yr⁻¹ (ENVEO, 2023), although the terminus position has remained stable despite these fluctuations.

2.1 Data

2.1.1 SAR Data: The ARCTIC15 Campaign

During the ARCTIC15 airborne campaign led by the German Aerospace Center (DLR) and the Swiss Federal Institute of Technology in Zurich (ETH), fully polarimetric and tomographic SAR acquisitions were performed using DLR’s F-SAR system over the selected test sites at X-, L-, and P-band frequencies (Horn et al., 2008). Table 1 summarizes the acquisition types, look directions, frequencies and resolution for each flight. Figure 2 illustrates the typical SAR side-looking geometry. For another study, C-band data were also acquired and analyzed for glacier facies classification but are not part in this study (Parrella et al., 2021).

Table 1 Overview of flights, and their parameters, conducted during the ARCTIC15 campaign over the Russell glacier.

For each resolution cell within the covered area, after interacting with the ice elements, the transmitted electromagnetic pulses are scattered back to the radar along the range direction, which forms an incidence angle θ with the surface normal. For the considered F-SAR acquisitions θ varies between 25 to 55°, providing a swath width of around 3000 m for a flight height of 3000 m. The azimuth direction is perpendicular to the range-height plane and aligns with the radar’s flight path. As the SAR pulses penetrate into the ice, the incidence angle changes due to refraction from θ to θ_r, which is dependent on the permittivity (ε_r=3.15) and must be accounted for during data analysis.

In the ARCTIC15 TomoSAR acquisitions, multiple tracks, horizontally displaced in a non-uniform way at the same flight height, were flown, as depicted in Fig. 2. The realized displacements provide sensitivity to the distribution of scatterers in the direction perpendicular to the range direction and therefore along the height z, enabling InSAR height estimations and TomoSAR imaging. The vertical wavenumber κ_z describes the interferometric sensitivity to height. It is defined as (Papathanassiou and Cloude, 2001; Reigber and Moreira, 2000; Bamler and Hartl, 1998):

$$\begin{array}{}\text{(1)}& {\mathit{\kappa }}_z={\displaystyle \frac{\mathrm{4}\mathit{\pi }}{\mathit{\lambda }}}{\displaystyle \frac{\mathrm{\Delta }\mathit{\theta }}{\sin \mathit{\theta }}}\end{array}$$

where λ is the radar wavelength, Δθ is the difference in incidence angles between acquisitions, and θ is the incidence angle. For subsurface characterization of ice, κ_zVol accounts for refraction and permittivity effects (Sharma et al., 2013):

$$\begin{array}{}\text{(2)}& {\mathit{\kappa }}_{z\mathrm{Vol}}={\displaystyle \frac{\mathrm{4}\mathit{\pi }\sqrt{{\mathit{\epsilon }}_{\mathrm{r}}}}{\mathit{\lambda }}}{\displaystyle \frac{\mathrm{\Delta }{\mathit{\theta }}_{\mathrm{r}}}{\sin {\mathit{\theta }}_{\mathrm{r}}}}\end{array}$$

where ε_r is the relative permittivity of ice, θ_r is the refracted incidence angle, and Δθ_r is the difference in refracted angles. The Rayleigh limit for the vertical resolution (ρ_z) of TomoSAR reconstructions depends on the maximum κ_zVol across all the tracks (Reigber and Moreira, 2000):

$$\begin{array}{}\text{(3)}& {\mathit{\rho }}_z={\displaystyle \frac{\mathrm{2}\mathit{\pi }}{max{\mathit{\kappa }}_{z\mathrm{Vol}}}}\end{array}$$

Figure 2 SAR side-looking geometry with multiple acquisitions from horizontally displaced tracks at the same height.

Table 2 summarizes the tomographic parameters at the Russell glacier, at L- and P-band, including the number of tracks, displacements, and vertical resolution.

Table 2 Summary of acquisition parameters at the Russell glacier.

Regardless of the frequency used, the minimum track displacement (t. disp. in Fig. 2) was set at 5 m because of limitations in the flight accuracy. As κ_zVol scales with the wavelength, the displacements at P-band were adjusted to provide the same tomographic performance as at L-band (including vertical resolution). Each TomoSAR acquisition was conducted within two hours under stable, frozen conditions with no variability from surface melt or internal structural changes. The slow glacial movement (less than 2.7 m during the entire campaign) at the Russell glacier ensured consistent ice properties across flights.

During one F-SAR flight over the Russel glacier, the P-band radar was operated as a Sounder (i.e. nadir-looking along the vertical direction) collecting data along a transect in full polarimetric mode. The processing of this Sounder data uses coherent summation of successive pulses and fully focused SAR processing. This approach enhances the along-track (azimuth) resolution. Motion compensation, adapted for the nadir-looking geometry, is also applied. In order to avoid excessive across-track (range) clutter and enable the Sounder profile to accurately isolate reflections from shallow subsurface layers we selected the cross-polarisation for our comparison. This methodology is based on an experimental processing chain developed at DLR (Scheiber et al., 2008).

It is worth remarking at this point that both TomoSAR's side-looking and Sounder nadir-looking modes offer distinct sub-surface imaging capabilities. TomoSAR provides vertical backscattering profiles with penetration down to maximum −80 m in P-band at the Russell glacier. In contrast, the Sounder vertical propagation is subject to a lower attenuation than the side-looking slanted one, and penetration in ice increases to depths of several hundred meters.

2.1.2 Additional Datasets

Additional optical and space-borne SAR datasets have been considered to comprehensively analyze the test site at the Russell glacier (Table 3). Optical data includes high-resolution optical imagery from Sentinel-2, providing a way to infer large-scale surface characteristics (Drusch et al., 2012). Additionally, an orthophoto obtained from NASA's Operation IceBridge (MacGregor et al., 2021), an airborne mission delivering detailed, small-scale imagery, allows for detailed examination of surface roughness and topography. Sampled height measurements have been provided by Airborne Topographic Mapper (ATM) lidar instrument from the NASA IceBridge campaign, enabling precise elevation assessments (Studinger et al., 2022). Finally, L-band SAR ALOS-2 data acquired between 2017 and 2024 have been used to detect annual changes of the ice features (JAXA, 2024).

Table 3 Supplementary SAR and optical data utilized at the Russell glacier.

Further data was collected from various sources to support the study. The ENVEO velocity map provides data across Greenland with spatial resolution of 50 m × 50 m (ENVEO, 2024), while PROMICE weather stations at the Russell glacier contributed local measurements of flow velocity, surface melt, surface and subsurface temperatures, and precipitation (Fausto et al., 2019).

3 SAR Methods

This section explains SAR polarimetry for identifying scattering processes, InSAR and TomoSAR to characterize and reconstruct vertical backscattering profiles.

3.1 SAR Polarimetry

Polarimetric SAR data enables qualitative and quantitative characterization of scattering mechanisms (Lee and Pottier, 2017). The systems considered operate in the horizontal (H) – vertical (V) transmit/receive basis, recording single-look complex (SLC) scattering amplitudes $S_{\mathrm{HH}},S_{\mathrm{HV}},S_{\mathrm{VH}},S_{\mathrm{VV}}$, where the subscripts indicate the transmit/receive configuration. In the case of a co-located transmitter and receiver and a reciprocal medium, it becomes S_HV=S_VH. To facilitate the analysis and interpretation of scattering processes, the Pauli scattering vector k_P (Lee and Pottier, 2017) is typically used:

$$\begin{array}{}\text{(4)}& {\mathit{k}}_{\mathrm{P}}={\displaystyle \frac{\mathrm{1}}{\sqrt{\mathrm{2}}}}\left[\begin{array}{c}{S}_{\mathrm{HH}}+S_{\mathrm{VV}}\\ S_{\mathrm{HH}}-S_{\mathrm{VV}}\\ \mathrm{2}{S}_{\mathrm{HV}}\end{array}\right]\text{.}\end{array}$$

S_HH+S_VV corresponds to the contribution of surface scattering (or more in general of odd bounces) to the total, S_HH−S_VV to dihedral scattering (even bounces), and S_HV to volume scattering. The intensities of these components are typically distributed across the channels of a single RGB image to provide a straightforward visual tool for a qualitative separation of different scattering behaviors within a scene (Lee and Pottier, 2017).

For distributed scatterers, such as ice volumes, individual SAR resolution cells contain randomly distributed scattering contributions. Their behavior is better described using second-order statistics, with the (3×3) coherency matrix T derived from k_P with $\mathbf{T}=\langle {\mathit{k}}_{\mathrm{P}}{\mathit{k}}_{\mathrm{P}}^{\mathrm{H}}\rangle \text{,}$ where 〈…〉 denotes spatial averaging over range-azimuth cells (i.e. multi-looking), and (.)^H represents the Hermitian operator. A way to perform a quantitative polarimetric data analysis separates three orthogonal elementary scattering mechanisms relying on the eigen-decomposition $\mathbf{T}=\mathbf{U}\mathbf{\Lambda }{\mathbf{U}}^{-\mathrm{1}}$, in which Λ contains the eigenvalues ${\mathit{\lambda }}_{\mathrm{1}}\ge {\mathit{\lambda }}_{\mathrm{2}}\ge {\mathit{\lambda }}_{\mathrm{3}}\ge \mathrm{0}$ on the diagonal and zeros elsewhere, and U=[e₁ e₂ e₃] contains three 3−dimensional eigenvectors on the columns. The eigenvalues are used to calculate the entropy $\mathrm{0}\le H\le \mathrm{1}$, which measures the scattering randomness or complexity of the scattering process (Cloude and Pottier, 1997):

$$\begin{array}{}\text{(5)}& H=-\sum _{i=\mathrm{1}}^{\mathrm{3}}{p}_i{\log }_{\mathrm{3}}{p}_i\text{,}{p}_i={\displaystyle \frac{{\mathit{\lambda }}_i}{{\sum }_{i=\mathrm{1}}^{\mathrm{3}}{\mathit{\lambda }}_i}}\text{.}\end{array}$$

Low H values indicate deterministic, single-scattering processes, whereas high values suggest multiple-scattering behavior. The polarimetric alpha angle (Cloude and Papathanassiou, 1998) can be calculated from the first element e_i1 of each eigenvector as ${\mathit{\alpha }}_i=\arccos \left(\left|e_{i\mathrm{1}}\right|\right)$. Then, a statistical interpretation is obtained by calculating the mean alpha angle as $\mathit{\alpha }={\sum }_{i=\mathrm{1}}^{\mathrm{3}}{p}_i{\mathit{\alpha }}_i$. Low values of α between 0 and 30° typically indicate surface scattering, intermediate values ($\mathrm{30}\mathit{°}\le \mathit{\alpha }\le \mathrm{60}\mathit{°}$) volume scattering, and high values ($\mathrm{60}\mathit{°}\le \mathit{\alpha }\le \mathrm{90}\mathit{°}$) dihedral scattering (Lee and Pottier, 2017).

3.2 SAR Interferometry

InSAR configurations use a limited number (typically two) of SLC SAR images S_i(ω) acquired under slightly different incidence angles (Fig. 2) to characterize the vertical backscattering profile within the same resolution cell. The three-dimensional (complex-valued) unitary vector ω represents the polarization state, defined as S_i(ω)=ω^Hk_Pi. For an interferometric pair ($i=\mathrm{1},\mathrm{2}$), the complex InSAR coherence ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{obs}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })$ is the normalized cross-correlation of both images (Bamler and Hartl, 1998):

$$\begin{array}{}\text{(6)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{obs}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })={\displaystyle \frac{\langle S_{\mathrm{1}}\left(\mathit{\omega }\right)S_{\mathrm{2}}^{*}\left(\mathit{\omega }\right)\rangle }{\sqrt{\langle S_{\mathrm{1}}\left(\mathit{\omega }\right)S_{\mathrm{1}}^{*}\left(\mathit{\omega }\right)\rangle \langle S_{\mathrm{2}}\left(\mathit{\omega }\right)S_{\mathrm{2}}^{*}\left(\mathit{\omega }\right)\rangle }}}\text{.}\end{array}$$

The measured coherence can be decomposed into multiple factors accounting for various sources of decorrelation, including temporal, range spectral, and systematic effects (${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Tmp}}$, ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Sys}}$, ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Rg}}$) (Zebker and Villasenor, 1992). These factors depend on either κ_zVol or ω:

$$\begin{array}{}\text{(7)}& \begin{array}{rl}{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{obs}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })& ={\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Tmp}}\left(\mathit{\omega }\right)\cdot {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Rg}}\left({\mathit{\kappa }}_{z\mathrm{Vol}}\right)\cdot {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Sys}}\left(\mathit{\omega }\right)\\ & \cdot {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Vol}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })\text{.}\end{array}\end{array}$$

Finally, the volumetric decorrelation ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Vol}}({\mathit{\kappa }}_{z\text{Vol}},\mathit{\omega })$ is related to the vertical backscattering profile σ_V(z,ω), where z indicates the height, by means of a Fourier relationship as (Bamler and Hartl, 1998):

$$\begin{array}{}\text{(8)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Vol}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })=e^{i{\mathit{\kappa }}_z{z}_{\mathrm{0}}}{\displaystyle \frac{{\int }_{-\mathrm{\infty }}^{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{V}}(z,\mathit{\omega })e^{i{\mathit{\kappa }}_{z\mathrm{Vol}}z}\mathrm{d}z}{{\int }_{-\mathrm{\infty }}^{\mathrm{0}}{\mathit{\sigma }}_{\mathrm{V}}(z,\mathit{\omega })\mathrm{d}z}}\text{,}\end{array}$$

in which z₀ indicates the topographic height of the glacier surface. σ_V(z,ω) is an electromagnetic quantity that depends on the vertical distribution (density) of surface and subsurface scatterers and their dielectric properties, and on the radar frequency, polarization and incidence angle (Bamler and Hartl, 1998). Moreover, the interferometric phase center height,

$$\begin{array}{}\text{(9)}& h_{\mathrm{pc}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })={\displaystyle \frac{\arg \left({\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{obs}}\right({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega }\left)\right)}{{\mathit{\kappa }}_{z\mathrm{Vol}}}}\text{,}\end{array}$$

where arg (.) is the argument of a complex number, approximates the centroid height of the vertical backscattering profile, especially for small κ_zVol (Dall, 2007). In Sect. 4.1, phase center heights are computed relative to an external DEM, estimating the glacier surface topography, i.e. after compensating the related interferometric phase term from ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{obs}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })$. Moreover, estimating ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Vol}}({\mathit{\kappa }}_{z\mathrm{Vol}},\mathit{\omega })$ is crucial for characterizing sub-surface scattering. For the airborne interferometric acquisitions (Table 2), non-volumetric decorrelation contributions are negligible due to the short acquisition time, large bandwidth, low noise, and high interferometric processing accuracy. Therefore, the complex InSAR coherence is assumed to be solely influenced by volumetric effects (Eq. 8), without any additional compensation for other sources of decorrelation.

3.3 SAR Tomography

TomoSAR techniques reconstruct the vertical backscattering profile by combining SAR images from multiple incidence angles (Reigber and Moreira, 2000). The TomoSAR data vector $\mathit{y}\left(\mathit{\omega }\right)=\left[S_{\mathrm{1}}\right(\mathit{\omega }),S_{\mathrm{2}}(\mathit{\omega }),\mathrm{\dots },S_K(\mathit{\omega }){]}^T$ is constructed, where K is the number of SAR images (Fig. 2), e.g., $K=\mathrm{7},\mathrm{8}$ for the F-SAR data sets at L- and P-band (Table 2). The TomoSAR data vector is then used to compute the (K×K) TomoSAR covariance matrix R(ω) as:

$$\begin{array}{}\text{(10)}& \mathbf{R}\left(\mathit{\omega }\right)=\langle \mathit{y}\left(\mathit{\omega }\right){\mathit{y}}^{\mathrm{H}}\left(\mathit{\omega }\right)\rangle .\end{array}$$

Assuming negligible non-volumetric decorrelation, each element of R(ω) is related to σ_V(z,ω) through a Fourier relationship with the vertical wavenumber (Eq. 2). In absence of scattering models, σ_V(z,ω) is reconstructed as (Lombardini and Reigber, 2003; Stoica and Moses, 2005):

$$\begin{array}{}\text{(11)}& F(z,\mathit{\omega })={\mathit{h}}^{\mathrm{H}}(z,\mathit{\omega })\mathbf{R}\left(\mathit{\omega }\right)\mathit{h}(z,\mathit{\omega }),\end{array}$$

where h^H(z,ω) is a (polarization- and height-dependent) coefficient vector preserving backscattering at height z while attenuating other contributions. The analysis in Sect. 4.3.1 follows the Capon method (Lombardini and Reigber, 2003; Stoica and Moses, 2005):

$$\begin{array}{}\text{(12)}& \mathit{h}(z,\mathit{\omega })={\displaystyle \frac{{\mathbf{R}}^{-\mathrm{1}}\left(\mathit{\omega }\right)\mathit{a}\left(z\right)}{{\mathit{a}}^{\mathrm{H}}\left(z\right){\mathbf{R}}^{-\mathrm{1}}\left(\mathit{\omega }\right)\mathit{a}\left(z\right)}},\end{array}$$

where the K-dimensional steering vector $\mathit{a}\left(z\right)={\left[\begin{array}{c}\mathrm{1},e^{-j{\mathit{\kappa }}_{z\mathrm{Vol},\mathrm{2}}z},\mathrm{\dots },e^{-j{\mathit{\kappa }}_{z\mathrm{Vol},K}z}\end{array}\right]}^T$ is a function of the interferometric phase generated by a scattering contribution at height z across all the acquisitions with respect to a reference one.

Here, κ_zVol,k is the vertical wavenumber for the kth track with respect to a reference one. The Capon formulation adapts for each height to the data through R(ω), hence it leads to reconstructions of σ_V(z,ω) with vertical resolution better than the limit ρ_z in Eq. (3) imposed by the track geometry. The inherent loss of radiometric accuracy (Lombardini and Reigber, 2003; Cazcarra-Bes et al., 2020) induced by the increase of resolution is deemed of secondary importance for the analysis in this paper.

The TomoSAR data vector can also be processed in order to isolate scattering contributions within a height interval and reconstruct their vertical backscattering profile as detailed above. This filtering is here performed using a (K×K)-dimensional (complex-valued) matrix G (Lombardini and Pardini, 2009; Joerg et al., 2017):

$$\begin{array}{}\text{(13)}& {\mathit{y}}_f\left(\mathit{\omega }\right)=\mathbf{G}\mathit{y}\left(\mathit{\omega }\right).\end{array}$$

The derivation of the matrix filter G requires first of all the definition of a height interval in which backscattering contributions are left undistorted, and a height interval in which backscattering contributions are cancelled. Then, the elements of G are calculated from a least-squares optimization involving a set of steering vectors calculated for sampled heights in the undistorted and attenuated intervals. The complete derivation and formulation is reported in Joerg et al. (2017).

4 Characterization of Ice Features

This section investigates the radar-dark and radar-bright features in our SAR data, examining their scattering processes, surface and subsurface characteristics, and temporal stability to understand their origins and potential glaciological formation processes.

4.1 Analysis of Polarimetric and Interferometric Data

We begin by examining scattering mechanisms of radar-dark and radar-bright features at the testsite, using Sentinel-2 optical imagery and multi-frequency Pauli representations (X-, L-, and P-band).

Figure 3 (a) Optical data (Sentinel-2, 4 May 2017) in geo-coordinates and polarimetric data in range-azimuth acquired in (b) X-, (c) L- band with a yellow rectangle indicating the area for surface characterization discussed in Sect. 4.2, and a white cross marks the location of the KAN L weather station, as detailed in Sect. 4.3.2, (d) P-band at the Russell glacier and shown in the Pauli basis (HH+VV (surface scattering) in blue, HH-VV (dihedral scattering) in red, HV (volume scattering) in green). For an additional comparison, the appendix (Fig. B1) contains geocoded backscatter images in HV polarization from different frequency bands (X-, L, and P-band), presented alongside a Sentinel-2 optical image captured during the summer.

The Pauli (Eq. 4) representations, in range-azimuth coordinates, offer qualitative insights into scattering mechanisms, visualizing surface scattering (HH+VV) in blue, dihedral scattering (HH-VV) in red, and volume scattering (HV) in green (Sect. 3.1). By combining optical and SAR data, we integrate visual and scattering information to enhance our understanding of the site's characteristics. The test site was processed in two parts for X- and L-band and color-scaled separately, then recombined for analysis.

Figure 3 shows the extent of the test site, with the glacier terminus at the western end, surrounded by rocky terrain. The minimal snow cover on the glacier reveals the overall topography and surface roughness, including small and large crevasses in the optical image in Fig. 3a. In the X-band Pauli image in Fig. 3b, these crevasses are represented as green lines. Their appearance may shift to red, indicating dihedral scattering, depending on their orientation relative to the radar's flight direction. Bright green areas usually signify volume scattering, associated with high surface roughness, highly crevassed areas, or potentially shallow penetration into the ice. L-band data (Fig. 3c) penetrates deeper into the ice (Table 5), providing more subsurface information while reducing surface detail. Here, radar-dark features dominated by surface scattering sharply contrast with radar-bright features characterized by volume scattering. In P-band data (Fig. 3d), the scattering patterns are similar to L-band, with enhanced subsurface contrast, though at a lower resolution (Table 1). By comparing HV backscatter in Fig. 1 with the Pauli representation in Fig. 3c, radar-dark features are linked to low backscatter and surface scattering, while radar-bright features correspond to high backscatter and volume scattering. However, not all volume scattering areas correspond strictly to radar-bright features, as highly crevassed areas may also appear as stronger volume scattering than surrounding regions.

Since Pauli representations are qualitative, we further analyze scattering entropy H and mean alpha angle α to quantify scattering mechanisms for radar-dark and radar-bright features across X-, L-, and P-band (Eq. 5), where both features exhibit distinct characteristics with dominant scattering mechanisms (Table 4). Radar-dark features have low entropy values in L- and P-band (0.0 to 0.6) and low mean alpha angles (0 to 30°), indicating surface-dominated scattering with minimal volume contribution. X-band, however, shows a broader range for H and α, likely due to different interactions with surface roughness. In contrast, radar-bright features have higher entropy (up to 1.0) and mean alpha angles (up to 60°) across all frequency bands, indicating multiple scattering mechanisms, including volume scattering, that are typical of more heterogeneous or layered subsurface structures (Parrella et al., 2021).

It is important to note that these parameters are typically influenced by variations in the incidence angle. In our case, both H and α show a slight systematic trend: lower values occur at steeper incidence angles (near range), gradually increasing toward the far range. However, as found in the analysis, entropy and mean alpha angle remain relatively stable overall. This suggests that, although a minor trend exists, the polarimetric characteristics used to differentiate radar-dark and radar-bright features are far less affected by incidence angle than by the spatial variation induced by the two radar features. This is in accordance with the H and α analysis in Parrella et al. (2021), which indicates that the incidence angle dependency is stronger in the percolation zone and less in the ablation zone. We therefore do not perform a dedicated incidence-angle analysis, and instead emphasize that this stability reinforces the robustness of our interpretation of the dominant scattering mechanisms.

Table 4 Entropy H and mean alpha angle α statistics for radar-dark and radar-bright features across frequency bands. The corresponding maps are provided for reference in the Appendix (Fig. B2).

After using PolSAR to identify scattering mechanisms, we shift our focus to InSAR, which complements this analysis by providing phase center heights (Eq. 9) from two acquisitions with a 10 m horizontal displacement between their tracks, allowing us to gain vertical information that PolSAR cannot capture. A change over different horizontal displacements and therefore a change in κ_zVol values (Eq. 2) will be shown in Sect. 4.4 for InSAR coherences.

Table 5 Phase center heights (PCH) for radar-dark and radar-bright features across frequency bands at an incidence angle between 30 to 40°.

Table 5 demonstrates that radar-dark features consistently exhibit shallow signal penetration depths, with negative phase center heights across all frequency bands and polarizations, typically ranging between −6 to 0 m. This limited depth aligns with PolSAR analysis, which indicates dominant surface scattering with minimal contributions from subsurface structures or volume scattering. In contrast, in L- and P-band, radar-bright features reach significantly lower phase center heights, extending from −50 to −10 m. This trend can also be confirmed across all incidence angles within the testsite. However, the penetration depth is far more sensitive to changes in the incidence angle than in the polarimetric analysis. Overall, a higher incidence angle generally leads to less radar penetration depth. This is because the signal travels a longer path through the ice, which increases attenuation from both absorption and volume scattering. Despite this, the relative difference in penetration between radar-bright and radar-dark features remains consistent across all incidence angles.

Overall, this significant difference suggests that an inhibiting factor limits further penetration into the ice for radar-dark features. It may be due to increased attenuation caused by changes in permittivity, likely from water content or impurities in the ice, a lack of scatterers deeper in the ice, or smooth surfaces reflecting the radar signal, limiting penetration through strong surface reflectance (Dall et al., 2001), which will be explored further in Sect. 4.2 and 4.3.

4.2 Surface Characterization

4.2.1 Analysis of Topography and Surface Roughness

After reviewing PolSAR and InSAR scattering processes, we analyzed whether radar-dark features correlate with surface processes by examining glacier topography and roughness.

Surface roughness, derived from a high-resolution orthophoto (Fig. 4a) and ATM lidar data from NASA IceBridge (Studinger et al., 2022), highlights surface variability in a central sample area (yellow rectangle in Fig. 3c). HV L-band backscatter (Fig.  4b), processed with a 10 m × 10 m multi-looking window, was used to locate these features. The temporal proximity of L-band and IceBridge data (Tables 1 and 3) ensures minimal surface changes, allowing for better comparability. Additionally, the reference DEM, acquired during the same F-SAR campaign, was analyzed to assess surface topography across the sample area (Fig. 4c). Given the overall slope in the sample area, the DEM was de-trended using a 50 m × 50 m window smoothing filter, and the difference from the original DEM was computed to emphasize small-scale topographic variations (Fig. 4d). The orthophoto (Fig. 4a) reveals a snow-covered glacier surface, where the typical winter snow height ranges from 0.5 to 1.5 m (Fausto et al., 2019). It is characterized by grid-like surface roughness patterns, which may result from intersecting strain fields or meltwater channels (Hambrey and Lawson, 2000; Yang et al., 2021), or alternatively, may represent wind-formed features. The ATM lidar tracks overlaid on the HV L-band backscatter image (Fig. 4b) show elevation differences along the flight lines (500 to 540 m elevation) similar to the reference DEM in Fig. 4c, illustrating slope changes and overall surface variability within the sample area. While the lidar tracks suggest subtle elevation gradients across the area, the de-trended DEM (Fig. 4d) emphasizes small-scale topographic variations, revealing a pattern broadly resembling the surface roughness observed in the orthophoto. However, despite these surface features, no direct spatial correlation is evident between the grid-like roughness patterns and radar-dark or radar-bright features, presented in the HV L-band data (Fig. 4b), suggesting that the observed radar signature is not primarily controlled by surface roughness or topography. Also at a larger scale over the complete test site, no correlation was found between small-scale topography and radar-dark features, see Fig. B3.

Figure 4 Surface characterization of the sample area: (a) orthophoto by NASA IceBridge, (b) L-band, HV backscatter (grey-scale) with overlaid point-wise elevation of ATM lidar data (color-scale) from NASA IceBridge, (c) reference DEM derived from DLR’s F-SAR data, (d) de-trended DEM based on (c). The DEM and de-trended DEM of the full test site are included in the Appendix (Fig. B3) for contextual reference.

4.2.2 Analysis of Snow Albedo Dynamics

Building upon the previous analysis, examining snow albedo dynamics provides deeper insight into surface composition. Higher albedo values indicate cleaner, brighter surfaces, such as fresh snow, which reflect more sunlight, whereas lower values correspond to darker, impurity-rich surfaces that absorb more heat (Navari et al., 2021). During the campaign acquisition period, albedo values measured at the KAN L Automatic Weather Station ranged from 0.7 to 0.9, indicating the presence of a snow cover (Fausto et al., 2019), consistent with the optical orthophoto shown in Fig. 4a.

Tedstone and Cook (2020) investigated snow albedo in relation to surface types and utilized Sentinel-2 data from summer 2017 to generate a comprehensive albedo map of the Russell Glacier, exhibiting an overall albedo range of 0.1 to 0.7. Radar-dark features emerge with their distinct shapes (Fig. 1) as areas of slightly lower albedo compared to the surrounding radar-bright features within the sample area, indicated by a red rectangle (Fig. 5). However, closer to the terminus, distinguishing albedo variations unrelated to topography and surface roughness becomes increasingly challenging.

Figure 5 Albedo map modified after Tedstone and Cook (2020), derived from Sentinel-2 imagery acquired on 21 July 2017. Low-albedo areas indicate cloud coverage, while the shape of radar-dark features is distinctly visible within the red rectangle as regions with lower albedo.

Furthermore, Ryan et al. (2018) classified Russell Glacier surface types into distinct categories during summer 2015, identifying two key types pertinent to this study: clean ice and ice containing uniformly distributed impurities. Although their point-based measurements do not encompass the entire surface, their findings reinforce the hypothesis that radar-bright features may correspond to clean ice, while radar-dark features are likely associated with impurity-rich areas.

Impurity-rich areas have also been associated with the formation and persistence of a weathering crust, as demonstrated by Tedstone et al. (2020), which plays a critical role in enhancing meltwater retention and altering ice surface properties, contributing to the weathering crust's development. The presence of residual liquid water within this crust, even during winter, may lead to near-surface attenuation of the radar signal, thereby explaining the observed SAR characteristics. This relationship suggests that radar-dark features may indicate areas of crust development, while radar-bright features likely represent cleaner ice, free of residual liquid water, which permits deeper radar penetration (Sect. 4.1). If this hypothesis holds, then SAR analysis significantly enhances the visibility of areas with weathering crust.

4.3 Subsurface Characterization

4.3.1 Analysis of Tomographic and Sounder Data

For the tomographic analysis, multiple frequency bands were investigated to gain insights into subsurface scattering at the Russell Glacier. X-band data primarily shows the surface and its roughness for both radar-bright and radar-dark features, with negligible penetration. In contrast, L- and P-band tomograms confirm significant differences in penetration depth (Sect. 4.1) and subsurface scattering scenarios between the two radar features. Therefore, we further analyze P-band since it provides the deepest penetration into the ice, which allows for the most detailed analysis of the glacier's subsurface.

The RGB Pauli representation (Fig. 6a) clearly illustrates the transitions between radar-bright and radar-dark features along the white transect line, where Sounder and tomographic profiles are generated. Figure 6b shows HV Sounder data in P-band. Polarization variability in TomoSAR is illustrated by comparing HH and HV, as shown in Fig. 6c–d using Capon tomograms (Eqs. 11 and 12), with VV omitted due to its similarity with HH. In Fig. 6e, another HV Capon tomogram is shown, generated using a matrix filter (Eq. 13), which removes scattering within a height interval of −4 to 2 m relative to the reference DEM, enhancing subsurface scattering visibility (Sect. 3.3). All Capon tomograms in Fig. 6c–e are normalized for each tomographic profile along azimuth, ranging from minimum to maximum values (0 to 1), offering a qualitative representation of the vertical backscattering distribution across the transect line.

In the Capon tomograms, a clear separation between the two ice features can be observed, especially in HV (Fig. 6d). Here, radar-bright features exhibit a well-defined, continuous subsurface scattering layer at heights between −15 to −50 m, which remains consistent and prominent across the transect line, indicating a robust subsurface structure. Tomographic analyses at various locations confirmed the presence of this consistent subsurface scattering layer across the entire test site. The HH tomogram (Fig. 6c) reveals a similar subsurface scattering layer, though with more pronounced surface scattering that reduces the visibility of the subsurface layer. For radar-dark features, the HV tomogram shows primarily surface scattering with minimal subsurface scattering. The HH tomogram further emphasizes the scattering at the surface, further diminishing the visibility of any possible subsurface scattering structures. Due to the poor visibility of these potential subsurface scattering structures for radar-dark features in both HH and HV polarizations, we applied the matrix filter to enhance subsurface scattering in Fig. 6e. This removal of surface scattering primarily affects the radar-dark features, as surface scattering components in HV for radar-bright features are already minimal. Here, the weak subsurface scattering, which was already present in Fig. 6d, but so weak that it is suppressed during normalization, is now revealed for radar-dark features. This finding indicates the existence of a subsurface layer, a characteristic previously associated primarily with radar-bright features, although it remains less pronounced.

Figure 6 (a) RGB Pauli representation with the white transect line for Sounder and tomographic profiles, (b) Sounder data (nadir-looking), (c) Capon tomogram (side-looking) in HH, (d) Capon tomogram (side-looking) in HV, and (e) Capon tomogram (side-looking) in HV with surface scattering components removed. (f) Polarimetric parameters (entropy and mean alpha angle) along the white transect line. All plots are shown in P-band. The height for (c)–(e) is set with respect to the surface elevation from the reference DEM. All heights are scaled assuming a permittivity of 3.15.

To further understand these subsurface scattering structures for both ice features, the Sounder data (Fig. 6b) serves as a valuable comparison. For radar-bright features, the Sounder HV data aligns well with the findings in Fig. 6d, revealing a consistent strong radar return starting at heights between −15 to −50 m, closely matching the main subsurface scattering layer observed in the tomograms. The downwards extent of this layer in the Sounder data could be due to off-nadir returns of the same, rather compact, layer as in the tomograms, or due to an actual vertical extent of the strong radar returns. For radar-dark features, the Sounder data reveals much deeper locations of the strong radar return, down to −100 m, which occurs often deeper than what is detectable in the matrix-filtered tomogram in Fig. 6e. This suggests that radar-dark features may indeed contain deeper subsurface structures that are weak but still observable, though the evidence is less conclusive compared to the radar-bright features. The same subsurface features were also observed in data from the Center for Remote Sensing of Ice Sheets (CRESIS) airborne accumulation Sounder (Paden et al., 2014), a radar system specifically designed to measure fine-scale layering. This system provides a finer vertical resolution than the F-SAR Sounder within the upper several hundred meters of the ice sheet (Arnold et al., 2019). The consistency between the two datasets serves as a valuable cross-validation, confirming the reliability of our F-SAR Sounder data.

Overall, the tomographic and Sounder analysis reveals distinct subsurface characteristics between radar-bright and radar-dark features and confirm the PolSAR (Fig. 6f) and InSAR understanding of dominant surface scattering at the radar-dark features and subsurface volume scattering at the radar-bright features. The presence of weak but observable subsurface scattering in the radar-dark features indicates that englacial water bodies, such as water pockets or channels are unlikely to be the cause of these features. If significant liquid water were present, the radar signal would be fully absorbed, preventing any subsurface scattering response. However, a small amount of liquid water within the ice could still explain an attenuation of the signal, causing the low-backscatter areas observed in Fig. 1 for the radar-dark features. In contrast, the well-pronounced subsurface scattering layer within the radar-bright features suggests these areas are frozen, with no water inclusions or remaining liquid water content, at least down to the depth of the main subsurface scattering layer.

4.3.2 Analysis of Temperature Profiles

Subsurface temperature profiles recorded by the KAN L Automatic Weather Station (Fausto et al., 2019) at the Russell glacier offer valuable insights into seasonal and vertical temperature variations within the radar-bright feature area.

Figure 7 Monthly subsurface temperature profiles for the year 2015 by the KAN L Automatic Weather Station (Fausto et al., 2019). Its location within the test site is indicated in Fig. 3.

Nine temperature measurements were taken at depths ranging from −10 to 0 m, relative to the ice surface, as illustrated in Fig. 7. The exact depths may vary due to accumulation and ablation processes (Fausto et al., 2019). Figure 7 highlights notable temperature variations at all available measurement depths during 2015. Below, subsurface temperatures are expected to be more consistent, suggesting a stabilization of thermal conditions at deeper levels. However, the lack of deeper measurements between −50 to −15 m presents a limitation for comparison, as this range corresponds to the main subsurface scattering layer identified in SAR tomography and Sounder data (Sect. 4.3.1). Despite this constraint, it is reasonable to assume that the ice in the upper layers undergoes significant seasonal temperature fluctuations, whereas the conditions below remain comparatively stable. This transition zone may represent a thermal boundary that is detectable in both tomography and Sounder data, manifesting as the main subsurface scattering layer.

4.4 InSAR Modelling of Scattering Structures

Building on the surface and subsurface characterization of radar-dark and radar-bright features, this modeling approach aims to transfer the current understanding of the scattering structures into a quantitative InSAR structure model at example locations. For this purpose, established InSAR structure model components are combined to simulate coherence magnitudes and comparing them with those measured from different track combinations in both L- and P-band. This exploits the sensitivity of InSAR coherences to volume decorrelation due to the vertical distribution of scattering (Eq. 8), as explained in Sect. 3.2. A model-based analysis of the entire test site goes beyond the scope of this paper, therefore, an exemplary location, defined as a single, multi-looked pixel, was selected for both radar-dark and radar-bright features. Detailed formulas and modeling parameters used in this analysis are provided in Appendix A. For the modeling of the radar-bright features, a two-component model is applied, one for the surface and one for the subsurface, due to the clear separation of scatterers in the tomographic analysis (Sect. 4.3.1). Hereby, the surface component is modeled using a rectangular function with a height of 2 m around the ice surface, which is based on the reference DEM, to account for the surface roughness on the glacier (Sect. 4.2). The subsurface component is modeled using a Uniform Volume (UV) model, which assumes a uniform distribution of scatterers within the volume (Zebker and Hoen, 2000; Fischer et al., 2020), starting at a height of −15 m for L- and P-band. The coherence magnitude of this two-component model is simulated across different κ_zVol values, representing various track displacements, and incorporates a surface-to-volume ratio m(ω) that changes with polarization. Figure 8a–b present the results of modeling radar-bright features in L-band, while Fig. 8c–d show the corresponding results for P-band. In both cases, the simulated coherence is represented by a line, and the measured coherences by discrete points. The general drop of coherence with increasing track displacement and κ_zVol shows the expected volume decorrelation, whereas the undulating behavior is typical for two scattering components with a certain vertical distance (Cloude and Papathanassiou, 2003). The two-component model generally provides a good fit across both frequency bands and polarizations.

Figure 8 Modeling of coherence magnitudes at the representative location of radar-bright features in different polarizations (HH and HV) for L-band (a–b) and P-band (c–d). The points represent the measured coherence, while the lines represent the simulated coherence.

Modeling the radar-dark features presents additional challenges due to their less distinct subsurface scattering structures. Previous analyses suggest that radar-dark features are primarily dominated by surface scattering, with minimal subsurface contributions. These subsurface contributions could be enhanced in Fig. 6e, which reveals a weak subsurface scattering layer at similar heights to those observed for radar-bright features. However, the Sounder data in Fig. 6b suggests a lower height for the subsurface scattering layer than in the radar-bright features. In this case, a three-component model can match the surface and subsurface scattering processes observed in both frequency bands and polarizations. The surface component is modeled identically as for the radar-bright features, supporting that surface properties are independent of the formation of radar-dark and radar-bright features. The main subsurface scattering layer, detected in Fig. 6e, is also modeled using a UV model with adjusted parameters (see Appendix A). Additionally, to fit the measured coherences, we added a weak subsurface layer at −5 m, even though this is generally not indicated by the tomograms.

Figure 9 Modeling of coherence magnitudes at the representative location of dark features in different polarizations (HH and HV) for L-band (a–b) and P-band (c–d). The points represent the measured coherence, while the lines represent the simulated coherence.

Overall, Fig. 9 demonstrates a good match between the modeled and measured data for both frequency bands and polarizations, providing insights into the scattering processes of radar-dark features. The success of this approach lies in its ability to fit both radar-dark and radar-bright features using the same structural model across all tested polarizations and frequencies, with the primary difference being the surface-to-volume scattering ratio. In the L-band HH polarization (Fig. 9a), coherence values remain consistently high, gradually decreasing from 1.0 to 0.8 as κ_zVol increases, as expected for dominant surface scattering, although two subsurface components are required to match this specific decrease. In contrast, the P-band HH polarization (Fig. 9c) exhibits a steeper decline in coherence, indicating a greater volume decorrelation due to deeper penetration into the ice – a trend consistent with physical expectations. For HV polarization in both frequency bands (Fig. 9c–d), the coherence shows significant variability over changing κ_zVol values, suggesting a complex subsurface structure. We purposely focused the simulation on capturing the higher-order oscillations of the coherences in the HV polarization to understand the required complexity in the vertical scattering structure. Moreover, the difference in scattering behavior between polarizations is more pronounced for radar-dark features, where HH polarization shows predominantly surface scattering, while HV polarization clearly exhibits strong volume decorrelation. The surface-to-volume ratio values from Tables A1 and A2, which were used to generate the modeled lines in Figs. 8 and 9, behave as expected. Specifically, higher values were estimated for radar-dark features, for HH polarization compared to HV, and for L-band (shorter wavelength) relative to P-band.

In summary, modeling of radar-bright features confirmed the surface and subsurface scattering structures identified in previous analyses. Moreover, we confirmed that the surface characteristics are independent of the presence of radar-bright or radar-dark features. However, for the radar-dark features, the subsurface scattering structure remains inconclusive, with tomographic, Sounder, and modeling analyses each suggesting different heights for potential subsurface scattering layers. Nonetheless, the UV model, along with an additional weak subsurface layer, matched the measured coherences very well as a two-component, subsurface scattering structure for the representative location of the radar-dark features.

4.5 Temporal Stability of Ice Features

Building on the surface and subsurface characterization of ice features, we introduce a temporal component to further investigate their formation processes or origins. Hereby, we focus on the temporal changes of radar-dark features within our test site at the Russell glacier, using time series data from ALOS-2 shown in Table 3 (JAXA, 2024). Figure 10 presents backscatter images in HV, showing the test site under frozen conditions in 2017, 2023, and 2024 (Fig. 10a–c). A backscatter threshold of −10 dB is used to outline the radar-dark features, which are overlaid in Fig. 10d for 2023 and 2024. In Fig. 10e–g, we highlight three sample areas that demonstrate variations in these features from 2017 to 2024, offering insights into potential changes over time. The overall description of the radar-bright and radar-dark features in the HV backscatter images (Fig. 10a–c) is similar to Fig. 1. The radar-dark features are visible as low-backscatter areas on the glacier. The radar-bright features, representing high-backscatter areas, surround the radar-dark features on the glacier. A visual comparison of the three HV backscatter images (Fig. 10a–c) suggest that radar-dark features remain stable in shape, size and position over the observed period and almost no change can be detected. By overlaying the outlines of the radar-dark features from 2023 to 2024, some changes, although subtle, are visible. Comparing the velocity maps from ENVEO (2024) with Fig. 10d confirms that the glacier moves slowly, with an average velocity of about 100 m yr⁻¹ at the eastern side of the test site (van de Wal et al., 2015), slowing down towards the terminus. Hereby, the radar-dark features generally move at the same velocity and in the same direction as the glacier. Analyzing the sample areas (Fig. 10e–g) shows the different changes between 2017 and 2024. In Fig. 10e, a radar-dark feature observed in 2017 (shown in blue) has moved consistent with the glacier’s velocity by 2024 (shown in red), with a slight increase in size. Additionally, two smaller features in 2017 have merged, forming a larger composite feature by 2024. In Fig. 10f, three distinct radar-dark features are visible in 2017. By 2024, the two larger features have nearly disappeared, while the smallest feature remains stable and continues to move with the glacier. In Fig. 10g, radar-dark features have increased in size and moved with the glacier over the observed period, with new radar-dark features also forming in this time period. These results show that radar-dark features move with the glacier’s velocity and direction, showing generally a high temporal stability, but also varying behaviors over time. These behaviors include no change, as well as significant increases or decreases in their size. Given the significant surface glacier melt (van de Wal et al., 2015; ENVEO, 2023), as described in Sect. 2, their consistent shape and persistence raise questions about their origin and formation processes.

Figure 10 ALOS-2 HV backscatter data from (a) 2017, (b) 2023 and (c) 2024, (d) the change and movement of the radar-dark features between 2023 in blue and 2024 in red; (e–f) three sample areas to display the variations of change between 2017 in blue and 2024 in red.

5 Discussion

5.1 Radar-Dark Features

Radar-dark features exhibit dominant surface scattering across multiple frequencies (X-, L-, and P-band) (Sect. 4.1), with no link to current topography or surface roughness, suggesting a subsurface origin (Sect. 4.2). They generally have an oval or round shape, with a mean size of several hundred meters, and often form interconnected patterns (Fig. 1). InSAR observables (Sect. 4.1) and SAR tomography (Sect. 4.3.1) confirm scattering mainly at the surface with weak contributions from the subsurface, which also could be modeled in Sect. 4.4. Moreover, analysis of ALOS-2 data revealed that the radar-dark features move with the glacier flow and remain relatively stable in shape, with only minor changes over time (Sect. 4.5), despite the large surface melt. The combined SAR analyses suggest that deeper penetration in L- and P-band is obstructed, likely due to high signal attenuation just below the surface. The attenuation is probably caused by small amounts of liquid water within the ice. However, significant water inside the glacier, such as englacial water bodies, would fully absorb the signal, eliminating any subsurface scattering contributions.

5.2 Radar-Bright Features

Radar-bright features are characterized by high backscatter in both L- and P-band, with volume scattering as the dominant scattering mechanism (Sect. 4.1). InSAR analysis (Sect. 4.1) and SAR tomography confirm deeper penetration and a prominent subsurface scattering layer between −15 to −50 m, extending across the entire test site. Moreover, Sounder data confirmed this subsurface scattering layer for the radar-bright features (Sect. 4.3.1). Modeling results further suggest that weak surface contributions exist, but they are less significant compared to the dominant subsurface scattering in the radar-bright features (Sect. 4.4), especially in HV polarization.

5.3 Preliminary Glaciological Interpretation

5.3.1 Comparison of Radar-Dark and Radar-Bright Features

Radar-dark and radar-bright features share several key characteristics, including a lack of correlation with current glacier topography or surface roughness (Sect. 4.2.1). Additionally, both ice features are similarly affected by the high surface melt rate of approximately 2.5 to 5.6 m yr⁻¹ (van de Wal et al., 2012; ENVEO, 2023). Investigated glaciological phenomena such as a direct effect of drained or refrozen supraglacial lakes (Hu et al., 2021; Yang et al., 2021), the dark zone – an area located further inside the Greenland ice sheet – and cryoconites, which may exist but are limited to very small scales and do not reach the size of the radar-dark features (Wientjes and Oerlemans, 2010; Ryan et al., 2018), as well as englacial water bodies (Catania et al., 2008; Lampkin and VanderBerg, 2013), have been ruled out as potential causes.

The most plausible explanation remaining is that the radar-dark features are influenced by a weathering crust (Cooper et al., 2018). The weathering crust, a porous layer formed near the glacier surface through meltwater infiltration and refreezing processes, plays a crucial role in shaping the near-surface properties of the ice (Cooper et al., 2018) and causes lower albedo values in summer due to a higher fraction of impurities (Sect. 4.2.2). This process may explain the formation and behavior of radar-dark features, contributing to their distinctive characteristics. During winter, the weathering crust can trap residual liquid water in the ice just below the surface, leading to signal attenuation, which would reduce the subsurface scattering contribution and overall backscatter in the radar-dark features significantly.

While the weathering crust theory aligns with most observations, it does not explain the shape and long-term stability of radar-dark features (Sect. 4.5). A possible explanation for the stable shape could be a connection to past supraglacial lake extent, as the radar-dark features strongly resemble lake outlines. The increased deposit of sediments at the location of past supraglacial lakes (Selmes et al., 2013; Leidman et al., 2023) might enhance the formation process of the weathering crust (Cooper et al., 2018). In this sense, past supraglacial lakes could drive the spatial pattern of the weathering crust. As lakes tend to form in consistent locations, their residuals, shaped by glacier flow, could produce the spatial patterns observed in our data. However, no direct evidence was found to support this theory.

In contrast, radar-bright features could correspond to bare, clean ice without the presence of the weathering crust. In these areas, there is minimal or no meltwater infiltration into the ice subsurface, resulting in mostly surface runoff and solid refreezing in winter, with no residual liquid water content (Cooper and Smith, 2019). This allows for high penetration of the SAR signal into the ice, particularly at lower frequencies, as shown in the interferometric (Sect. 4.1) and tomographic analysis (Sect. 4.3.1) of the radar-bright features. The enhanced penetration results in strong (volume) backscatter from any structural or dielectric heterogeneities in the ice subsurface of the radar-bright features.

5.3.2 Interpretation of the Main Subsurface Scattering Layer

The subsurface scattering layer in the radar-bright features is prominently pronounced in both L- and P-band tomograms at heights between −15 to −50 m, particularly in HV polarization (Fig. 6). The presence and height of this subsurface scattering layer are confirmed by comparisons with Sounder data, particularly at the locations of the radar-bright features (Sect. 4.3.1). This consistency across different viewing geometries – from side-looking SAR to nadir-looking Sounder – confirms that the observed scattering is a genuine subsurface feature, not an artifact of surface clutter or the specific SAR acquisition geometry. Furthermore, this layer's widespread and reliable presence is corroborated by additional tomographic analysis across multiple locations and headings within the test site, reinforcing its existence as a widespread and reliable feature. For radar-dark features, the subsurface scattering layer only becomes detectable after techniques are applied to remove dominant surface scattering in the tomographic analysis (Sect. 4.3.1). Modeling confirms the presence of this layer, although it appears weaker and is located deeper than in the radar-bright features (Sect. 4.4). Despite these differences, it is reasonable to assume that the main subsurface layer is present across the entire test site. Therefore, the consistent presence of this subsurface scattering layer needs an exploration of its origin. Englacial channels are ruled out due to their typically localized nature (Catania et al., 2008). Additionally, a hypothesis that this layer originated from previous glacier zones is dismissed, as such layers would likely show slope and height variations from glacier movement and melting, which are not observed (Florentine et al., 2018). Our current hypothesis is that above this subsurface scattering layer, seasonal warming allows for temperature fluctuations, while below lies impermeable ice (Sect. 4.3.2). The related transition in dielectric properties may cause the scattering. Variations in height of this layer may result from factors such as changes in the onset of impermeable ice, or variations in latent heat transmission (Florentine et al., 2018). Overall, this explanation would align with observed data and provides a coherent understanding of the subsurface scattering layer’s characteristics and distribution across the test site, even though the available temperature measurements do not reach the depths of the observed subsurface layer (Sect. 4.3.2). The exact details of this theory require further dedicated investigation with additional ground measurements in future research.

6 Conclusions

This study investigated the surface and subsurface characteristics of radar-bright and radar-dark features in the ablation zone of southwest Greenland, focusing on the Russell glacier. Using multi-frequency SAR data, we demonstrated that both radar-dark and radar-bright features are detectable in L- and P-band SAR images and can be effectively characterized using different advanced SAR techniques such as polarimetry, interferometry and tomography.

Radar-dark features were identified as low-backscatter surface scattering elements, typically convex-shaped and several hundred meters in size, often forming interconnected patterns. No immediate, consistent, and directly observable correlation with surface topography or roughness was found. Their location and size exclude the “dark zone” or cryoconite deposits as determining factors. Temporal analysis using ALOS-2 data showed that radar-dark features move with the glacier flow while remaining remarkably stable in shape and size, despite the large surface melt.

In contrast, radar-bright features were characterized by high backscatter and volume scattering, with the main subsurface scattering layer located at heights of −15 to −50 m in both L- and P-band. This layer's consistency over the complete testsite rules out englacial channels and formations from prior glacier zones as its origin. Glaciologically, this subsurface scattering layer likely forms due to seasonal temperature changes occurring above this layer, with impermeable ice lying beneath.

A preliminary glaciological interpretation suggests that radar-dark features are associated with the presence of a weathering crust, a porous layer that allows meltwater infiltration and refreezing. Through this process, residual liquid water content remains just below the ice surface during frozen conditions, leading to surface scattering and SAR signal attenuation. Further, their shape strongly suggests a connection to past supraglacial lake extent, where the increased deposit of sediments might enhance the formation process of the weathering crust. Therefore, the weathering crust best explains the characteristics of the radar-dark features, while radar-bright features likely represent bare, clean ice without a weathering crust.

While our SAR and Sounder data has provided significant insights, further validation using ground measurements such as ice core drillings, is essential to fully understand the formation processes of these features. Eventually, long-wavelength SAR could become a monitoring tool for weathering crust formation and near-surface meltwater storage in the ablation zone, contributing to the understanding of mass balance and runoff processes.

Appendix A: Scattering Modeling Formulas and Parameters

The coherence magnitude for radar-bright features is modeled using a two-component approach that accounts for both surface and subsurface scattering contributions. The coherence magnitude as a function of polarization w is given by:

$$\begin{array}{}\text{(A1)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Vol}}\left(\mathit{\omega }\right)=e^{i{\mathit{\kappa }}_z{z}_{\mathrm{0}}}{\displaystyle \frac{{\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Subsurface}}+m\left(\mathit{\omega }\right){\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Surface}}}{\mathrm{1}+m\left(\mathit{\omega }\right)}}\end{array}$$

Here, $e^{i{\mathit{\kappa }}_z{z}_{\mathrm{0}}}$ represents the phase shift due to the reference surface height z₀, which typically corresponds to the ice surface as defined by the reference DEM height. The surface-to-volume ratio m(ω) weights the surface and subsurface contributions, with its value depending on polarization, which affects the balance between surface and subsurface scattering. The surface component ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Surface}}$ is modeled using a rectangular function to account for surface roughness at the representative location (Cloude and Papathanassiou, 2003):

$$\begin{array}{}\text{(A2)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Surface}}=e^{i\left(z_{\mathrm{0}}+\frac{h_v{\mathit{\kappa }}_z}{\mathrm{2}}\right)}{\displaystyle \frac{\sin \left(\frac{h_v{\mathit{\kappa }}_z}{\mathrm{2}}\right)}{\frac{h_v{\mathit{\kappa }}_z}{\mathrm{2}}}}\end{array}$$

In this equation, the sinc function represents the Fourier transform of a rectangular structure function, incorporating z₀ for the reference surface height, h_v as the height of the surface scattering layer, and κ_z as the vertical wavenumber.

The subsurface component ${\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Subsurface}}$ is modeled using a Uniform Volume (UV) model, which assumes a uniform distribution of scatterers and captures diffuse scattering in subsurface ice layers (Zebker and Hoen, 2000; Fischer et al., 2020). The subsurface model is expressed as:

$$\begin{array}{}\text{(A3)}& {\stackrel{\mathrm{̃}}{\mathit{\gamma }}}_{\text{Subsurface}}=e^{i{\mathit{\kappa }}_{z\mathrm{Vol}}{z}_{\mathrm{ul}}}{\displaystyle \frac{\mathrm{1}}{\mathrm{1}+\frac{i{d}_{\text{Pen}}\left(\mathit{\omega }\right){\mathit{\kappa }}_{z\mathrm{Vol}}}{\mathrm{2}}}}\end{array}$$

In this formula, z_ul represents the upper limit of the subsurface scattering layer. The parameter d_Pen(ω), which varies with polarization, represents the penetration depth of the radar signal into this volume component.

The optimization process adjusts the parameters m(ω) and d_Pen(ω) with flexibility across different polarizations and frequency bands (Table A1). This approach enables the model to accurately fit the observed coherence magnitudes by effectively capturing the subsurface scattering structures within the ice. The height of the rectangular function for the surface scattering layer is h_v=2 in all cases. The upper limit z_ul for the UV model is set as −15 m for both L- and P-band consistently across different polarizations.

Table A1 Modeling parameters for the radar-bright features, specifying the values for L- and P-band in both HH and HV polarizations.

For radar-dark features, the scattering model consists of three components: a surface component, a weak subsurface scattering layer (not detected in Fig. 6d–e), and a deeper subsurface scattering layer (detected in Fig. 6e). The surface component is modeled using a rectangular function identical to that of radar-bright features, demonstrating the independence of surface roughness between the two features (Sect. 4.2). The weak subsurface layer is represented by a Dirac delta function at −5 m, accounting for additional scattering observed during the modeling approach. The deeper, main subsurface scattering layer is modeled using a UV model, with the upper limit adjusted to capture the deeper subsurface characteristics specific to radar-dark features. In this case, the upper limit for the UV model is set at −48 m in both frequency bands. The modeling parameters for radar-dark features, including the surface-to-volume ratio m(ω), Dirac-to-volume ratio m(ω), and UV penetration depth d_Pen(ω), are detailed in Table A2.

Table A2 Modeling parameters for the radar-dark features, specifying the values for L- and P-band in both HH and HV polarizations.

Appendix B: Additional Figures

B1 Backscatter Maps at Different Frequency Bands

Figure B1 (a) Optical – true color – image of the test site acquired by Sentinel-2 on 21 July 2017 and backscatter maps in HV polarization for different frequency bands: (b) X-band, (c) L-band, and (d) P-band.

B2 Entropy and Mean Alpha Angle Maps at Different Frequency Bands

Figure B2 Polarimetric decomposition results of the F-SAR data across X, L, and P-band: (a, c, e) Entropy (H) [–] at X, L, and P-band, respectively. (b, d, f) Mean alpha angle (α) [°] at X, L, and P-band, respectively. (g) Incidence angle [°] at P-band.

B3 Digital Elevation Model (DEM) and de-trended topography of the study area

Figure B3 (a) Reference DEM derived from DLR’s F-SAR data, (b) de-trended DEM based on (a), (c) Pauli RGB at P-band.