Quantifying the impact of X-band InSAR penetration bias on elevation change and mass balance estimation

Abstract

Interferometric synthetic aperture radar (InSAR) data suffer from an elevation bias due to signal penetration into the firn and ice surface, rendering the height information unusable for elevation and mass-change detection. This study estimates the penetration bias in X-band InSAR data to quantify its impact on elevation and mass-change detection and to demonstrate the applicability of TanDEM-X digital elevation models (DEMs) for cryosphere research. To achieve this, a multiple linear regression model is applied to a time series of four TanDEM-X DEMs acquired between 2010 and 2018 over the Sverdrup Glacier basin (SGB), Devon Ice Cap, Canada. The resulting penetration corrected TanDEM-X DEMs agreed to within ±14 cm of spatially and temporally coincident precise in situ kinematic dGPS data (±10 cm RMSE). Additionally, multi-year estimations of mass change for the SGB derived from differencing TanDEM-X DEMs over multi-year periods between 2010 and 2018, showed good agreement with mean deviation of 338 ± 166 mm w.e. with independent measurements of mass change derived from annual in situ surface mass balance over the same time periods. The results show that the penetration bias can vary significantly, leading to random under- and overestimations in the detection of elevation and mass changes.

1. Introduction

Ongoing global warming is causing accelerating melting of ice sheets, ice caps and glaciers worldwide (IPCC, 2022). The dramatic changes in the cryosphere have enormous impacts on humans and the environment that will continue to intensify in the future. These impacts are not limited to the cryosphere but threaten many millions of people in the world's coastal regions due to the associated rise in sea level (Wong and others, 2014). There is a persistent need for reliable and consistent information on the current state and the evolution of the ice sheets, ice caps and glaciers worldwide. In this context, satellite remote sensing enables not only a cost-efficient data acquisition but also regular collection of information on outlying and inaccessible areas. With respect to height and mass changes of ice sheets, ice caps, and glaciers and several remote sensing techniques have been widely used to date: (1) stereo-photogrammetry (Berthier and Brun, 2019; Howat and others, 2019; Dehecq and others, 2020), (2) radar altimetry (Nilsson and others, 2016; Sørensen and others, 2018; Jakob and others, 2021), (3) gravimetry (Tapley and others, 2019; Sasgen and others, 2020; Velicogna and others, 2020), and (4) laser altimetry (Pritchard and others, 2009; Sørensen and others, 2011; Smith and others, 2023).

These techniques have advantages but they have their limitations. Some have low-spatial resolution, others are based on point measurements (i.e. no areal coverage), and they may depend on illumination and weather. This poses a considerable limitation to monitoring the condition and evolution of ice sheets, ice caps and glaciers. In this regard, the single-pass InSAR could complement the existing methods to address the need for area-wide elevation information with high spatial resolution and thus to provide a more comprehensive view of the spatial and temporal dynamics of ice masses worldwide.

Since the launch of the TanDEM-X (TerraSAR-X-Add-on for Digital Elevation Measurements) mission (Krieger and others, 2007) in 2010, single-pass InSAR offers synchronous acquisition of SAR images from different spatial positions and enables an area-wide as well as spatially high-resolution reconstruction of the topography of the ice sheets without accuracy limitations due to temporal decorrelations or atmospheric disturbances (Krieger and others, 2013). The area-wide, high-precision terrain information of the TanDEM-X mission is consistently available worldwide with a spatial resolution of 0.4″ and allows the derivation of height and thus mass changes of the ice sheets on the basis of multitemporal acquisitions (Berthier and others, 2023). Especially for monitoring the narrow, fast-flowing outlet glaciers of the ablation zones with the highest rates of elevation change (Vaughan and others, 2013; Cooper and Smith, 2019), single-pass InSAR with its area-wide high-resolution and high-accuracy elevation information provides important advantages.

A key challenge to exploiting single-pass InSAR for monitoring ice sheets is the penetration of the radar signals into the snow and ice masses and the resulting elevation bias (Dall, 2007). The so-called penetration bias is defined as the difference between the height measured by the interferometer via the position of the scattering center and the real surface height (Weber Hoen and Zebker, 2000). Depending on the physical properties of the snow, firn and ice layers, such as stratigraphy, surface roughness or water content, the radar signals are scattered at different depths below the ice sheet surface up to several meters (Abdullahi and others, 2019; Fischer and others, 2019; Rott and others, 2021). The scattering behavior and thus also the penetration bias can vary strongly locally and they are subject to inter- and intrannual variations (Rignot and others, 1993; Liu and others, 2006). In addition, the backscattering behavior of the radar signals and thus the penetration bias is determined by the properties of the interferometer, such as frequency, polarization or acquisition geometry (Dall, 2007). The complex interaction of these factors and the resulting continuously varying penetration bias considerably complicates the analysis of the X-band InSAR data with respect to changes in glacier and ice sheet topography (Dehecq and others, 2016; Leinss and Bernhard, 2021; Wessel and others, 2021).

To address this problem, previous studies on the derivation of elevation and mass change using single-pass InSAR are often based on a spatially and temporally limited data basis with the assumption of negligible penetration of the radar signals (Abdel Jaber and others, 2019; Farías-Barahona and others, 2020) and/or consider the penetration bias as an error contribution in the uncertainty estimation of the mass balance calculations (Braun and others, 2019; Seehaus and others, 2019; Krieger and others, 2020). However, there are approaches to estimate the penetration bias in single-pass InSAR data that aim either to model the scattering behavior of the radar signals (Weber Hoen and Zebker, 2000; Dall, 2007; Fischer and others, 2020; Rott and others, 2021) or to determine the penetration bias more generally with respect to zones of different backscattering (Rignot and others, 2001; Rizzoli and others, 2017). In most cases, inter- and intraannual variations are not considered. Moreover, the proposed approaches require prior knowledge about the physical properties of the snow and ice surface, such as information on firn and ice density or the presence and number of glacier zones. This imposes a huge limitation on large-scale applicability due to the lack of availability of such additional information. However, an accurate estimation of the penetration bias is crucial to determining accurate elevation and therewith mass change estimates for bi-static interferometric data.

In this context, the current study aims to quantify the impact of X-Band InSAR penetration bias on elevation and mass-change detection using TanDEM-X Digital Elevation Models (DEM). The objectives of this study are: (i) the prediction of penetration bias based on the approach of Abdullahi and others (2019) and (ii) the detection of basin-wide elevation and mass change based on uncorrected as well as corrected TanDEM-X DEMs in order to quantify the impact of X-band InSAR penetration bias on geodetic mass balance estimates.

2. Study area and data

2.1. Study area

The study area is the Sverdrup Glacier basin (SGB), which occupies 765 km² (Randolph Glacier Inventory, Version 6.0, 2017) of the northwest sector of the Devon Ice Cap (DIC) (Fig. 1). The DIC covers ~14 000 km² of eastern Devon Island, which is located in the southeast region of the Queen Elizabeth Islands, Nunavut, Canada. The DIC has a dome-like shape with ice divides running towards east, north, and south from the summit at 1921 m a.s.l. (Boon and others, 2010). Excluding the stagnant southwest arm (Burgess and Sharp, 2004), ~90% of the ice cap area is drained by 15 major dynamic basins, 14 of which flow towards tidewater terminating glaciers. Iceberg calving accounts for ~0.4 Gt of mass loss annually from 13 tidewater glaciers (Van Wychen and others, 2017), with 30% discharged from the Belcher Glacier alone (Burgess and Sharp, 2004). Climatically, DIC is located in the polar desert climatic region, where annual precipitation is often <200 mm w.e. with minimal interannual variations, and surface air temperatures exceed the freezing point only during two to three months in summer (June–August) (Gardner and Sharp, 2007). The intensity and duration of summer melt represents the primary control on total mass balance of glaciers and ice caps in this region (Koerner, 2005; Millan and others, 2017).

Fig. 1. Map of the study area and data coverage. (a) Overview of the DIC with its basins, the outline of the study area and locations of the mass balance poles and automatic weather stations along the Devon NW transect. (b) TanDEM-X data coverage over the study area. (c) Reference data coverage over the study area. (d) Location of Devon Island within the Queen Elizabeth Islands in the Canadian Arctic Archipelago.

The SGB is situated between ~900 and 1800 m a.s.l. on the ice cap proper, below which the basin-ice converges towards the head of the Sverdrup Glacier flowing ~25 km mountain valley to sea-level. Systematic measurements of surface mass balance were initiated within the SGB in 1960 as part of a fully integrated research program conducted by the Arctic Institute of North America (AINA) on Devon Island from 1959–1961 (Apollonio and others, 1961). Since the AINA program ended in 1962, annual measurements of surface mass balance along the original transects have been maintained by the Government of Canada (Koerner, 2005; Burgess, 2017). Consistent with the other long-term glacier monitoring sites in the Canadian high Arctic, results from the surface mass balance surveys on the DIC indicate a trend of increasingly negative surface mass balance beginning in the early 1990s, becoming significantly more negative in the early 2000s (Burgess and Danielson, 2022). Since 2005, all monitored ice caps and glaciers in the Canadian high Arctic have been losing mass due to melt-runoff at rates four times higher than in the previous four decades (Sharp and others, 2011).

2.2. Data

2.2.1. TanDEM‒X

A time series of four TanDEM-X datasets (Fig. 1b, Table 1) acquired between 2010 and 2018 was used in order to derive changes in elevation and mass of the SGB in the northwest sector of the DIC. To maximize signal coherence and consistency between image datasets used, all radar data were acquired in the frozen seasons of November, December, and May, when near-surface hourly averaged temperatures did not exceed 0°C. The crossing equator orbit direction was ascending for the acquisitions from 2010 and 2012, while the other two datasets were acquired with a descending orbit. The distance between satellites varied from one acquisition date to the next with a perpendicular baseline ranging from 81 to 187 meters. Incidence angles between acquisitions however, remained similar (Table 1).

Table 1. TanDEM-X datasets

The TanDEM-X datasets used in this study provide nearly complete spatial coverage across the SGB for the years 2010, 2013 and 2018. Based on the data available, spatial coverage during the frozen season in 2012 was limited to 78% of the SGB. Figure 2 shows the hypsometric coverage of the SGB by the TanDEM-X datasets. It can be assumed that all datasets are representative for the entire SGB.

Fig. 2. Hypsometric coverage of the TanDEM-X datasets compared to the hypsometric coverage of the SGB.

Based on bistatic co-registered single-look complex SAR data (CoSSCs), the Operational Integrated TanDEM-X Processor (ITP) of the German Aerospace Center (DLR) was used to generate pre-calibrated, geocoded single-scene digital elevation models (DEMs) for each acquisition date (Fritz and others, 2011; Rossi and others, 2012). In addition to the DEMs, interferometric coherence and backscatter intensity layers were employed to estimate the penetration bias and to quantify its impact on elevation and mass change. All TanDEM-X datasets were resampled bilinearly to a pixel spacing of 12 m to match the independent pixel spacings of the interferometric processing filters and to reduce the random height error by averaging.

2.2.2. Reference data

For validation purposes of penetration bias estimation as well as mass balance calculation, several sources of point-based and gridded elevation information were used in this study (Fig. 1c, Table 2): Airborne laser measurements from NASA's Operation IceBridge mission (Studinger and others, 2010), differential Global Positioning System (dGPS) data collected during field campaigns, and DEMs based on optical stereo imagery from the U.S. National Science Foundation, National Geospatial-Intelligence Agency ArcticDEM initiative (Polar and Geospatial Center, 2022; Porter and others, 2022). All datasets were rasterized and resampled to 12 m aligned to the same grid as the TanDEM-X data projected to Universal Transverse Mercator (UTM) Zone 17N in ellipsoidal heights (World Geodetic System (WGS84)). An overview of the reference datasets used in this study is given in Table 2.

Table 2. Reference datasets

IceBridge Data. Airborne laser altimeter measurements of the NASA Operation IceBridge Mission were acquired by the NASA Airborne Topographic Mapper (ATM) scanning lidar system. In this study, the L1B product containing geolocated surface elevations processed with aircraft dGPS, corrected for biases including pitch and roll was used. The spatial resolution varies with flying altitude and the lidar scanner configuration. Operating at a typical height of 500 m above the ground, the ATM illuminates a swath ~140 m wide with a footprint size of 1–3 m, an along-track separation of 2 m with a laser pulse rate of 5 kHz, the ATM scanning lidar system achieves a vertical accuracy <10 cm (Martin and others, 2012; Studinger, 2013).

dGPS Data. Reference elevation measurements of the Devon Ice Cap surface were collected via kinematic dGPS snowmobile surveys over, (1) a 45 km by 6 km grid (16–21 May 2011) and (2) five small scale transects (100 × 100 m, 5 lines at 20 m interval) along the Devon NW transect in April 27–28, 2018. The dGPS survey data were collected at 1 s observation rates by two dual frequency (L1/L2) Ashteck Z-Xtreme^TM recievers. The dGPS data were post-processed to an x, y, z accuracy of ±10 cm using the Canadian Spatial Reference System – Precise Point Positioning service, and corrected for antenna height offsets. Data gaps due to GPS signal loss occurred along 8 km, or ~2% of the grid surveyed due to antenna wobble over rough ice cap surfaces.

ArcticDEM. The ArcticDEM consists of stereoscopic DEMs extracted from optical satellite data at centimeter spatial resolution (0.32–0.5 m) from WorldView-1, WorldView-2 and WorldView-3, as well as GeoEye-1 satellite sensors. Stereo image data were acquired at 2-m spatial resolution during the summer season between 2007 and 2021. The individual DEMs, called strip files were processed using the SETSM algorithm with a vertical accuracy of about ±2 m (Noh and Howat, 2015). For validation purposes, three ArcticDEM strips acquired on 30 April 2013, 5 April 2014 and 8 April 2014, were vertically and horizontally co-registered to the corresponding TanDEM-X DEMs of the same frozen season over stable terrain according to Nuth and Kääb (2011). Comparison between values of the ArcticDEM and TanDEM datasets yielded mean deviations of −0.02 ± 0.68, 0.01 ± 0.96 and 0.01 ± 0.74 m, for the three time periods respectively.

Glaciological Mass Balance Estimations. As part of the long-term glacier monitoring activities within the Government of Canada, the National Glaciology Project collects glaciological measurements across the northwest sector of the DIC (Devon NW, Fig. 1a). The Devon NW transect comprises 44 stakes along a ~60 km transect that extend from the ice cap summit region (~1800 m a.s.l.) to the terminus of the Sverdrup Glacier, with a separate branch extending to the west margin (~1200 m a.s.l.). The glaciological surface mass balance surveys are conducted each spring (April–May) according to the Stratigraphic System (Cogley and others, 2011) whereby net annual basin-wide surface mass balance is derived from the density corrected difference in stake height above the last summer surface, over two successive years. For example, net balance for the mass balance year extending from late summer 2020 to late summer 2021 is calculated from pole measurements obtained during the 2021 and 2022 spring field campaigns. Annual winter balance is measured at each pole as the depth and density of the winter snowpack which overlies the last summer surface. The area-weighted sum of net annual, summer, and winter balance calculated at 100 m elevation band intervals provide basin-wide values of annual and seasonal mass balance for each measurement year. Surface mass balance stake measurements are augmented with meteorological data collected from five automatic weather stations along the Devon NW transect. Hourly soundings of snow/ice surface height change (±1 cm), and air temperature provide data critical for interpreting the timing and magnitude of the mass balance processes (i.e. melt and accumulation) across the ice cap surface (Burgess, 2017).

3. Methods

Figure 3 depicts the datasets used and workflow of the current study. During preprocessing, the elevation models were vertically calibrated, interferometric coherence was adjusted and backscatter intensity was radiometrically calibrated. Subsequently, the multiple linear regression model from Abdullahi and others (2019) based on interferometric coherence and backscatter intensity was used to estimate the X-band penetration bias in order to quantify its impact on elevation and mass-change detection. The correctness of the estimated penetration bias and thus, correction of the TanDEM-X elevation models was evaluated using point-based and gridded reference height datasets. Moreover, mass changes derived from the uncorrected and the corrected TanDEM-X DEMs were compared with the independent estimations of mass balance for the SGB.

Fig. 3. Workflow of this study.

3.1. Pre–processing of TanDEM–X data

3.1.1. Absolute height calibration of TanDEM–X DEMs

After DEM generation within the ITP, a block adjustment procedure was applied to the pre-calibrated, geocoded single-scene DEMs through the DEM Mosaicking Calibration Processor (MCP) in order to correct for systematic absolute height offsets due to residual baseline and orbit inaccuracies (Gruber and others, 2012). For this purpose, a least-squares block adjustment procedure based on ground control points over stable terrain from ICESat (Ice, Cloud and land Elevation Satellite) GLA14 products and tie points for each single-scene DEM was used.

Comparison with independent validation ICESat measurements over off-glacier terrain after operational calibration showed high height accuracy for the DEMs from 2010 and 2012 (mean offset <0.4 m, Std dev. <1.4 m). Due to the lack of a sufficient number of ice-free ICESat ground control points in 2013 and 2018, an adapted calibration procedure similar to Wessel and others (2016), which solely relies on tie points over rocky terrain with the well-calibrated adjacent data takes, was applied and also yielded high height accuracies (mean offset <0.3 m, Std dev. <1.4 m).

3.1.2. Adjustment of interferometric coherence

The interferometric coherence is defined as the normalized cross-correlation coefficient between the interferometric image pair and describes the degree of correlation between the two radar images (Moreira and others, 2013). When the radar wave penetrates into a firn and ice volume, the signals are backscattered from different depths and the measured InSAR height corresponds to the height of the mean scattering center below the surface. The deeper the signals penetrate into the firn and ice volume, the more the two radar images disagree and the higher the loss of coherence. This interaction is commonly referred to as volume decorrelation (Rizzoli and others, 2022). However, the correlation between the vertical position of the scattering center (i.e. the magnitude of the penetration bias) and the volume decorrelation depends not only on the characteristics of the firn and ice volume (e.g., grain size, density, liquid water content, and stratigraphy) but also on the acquisition parameters such as frequency and acquisition geometry (e.g., effective baseline of the InSAR system and incidence angle) (Weber Hoen and Zebker, 2000; Martone and others, 2016; Rizzoli and others, 2022). For the TanDEM–X datasets acquired under acquisition geometries given in Table 1, the interferometric coherence layers were adjusted according to Abdullahi and others (2019) in order to ensure that the estimation model retains its validity without further calibration to the dataset of the current study.

3.1.3. Radiometric calibration

Following Ashcraft and Long (2005), the backscatter coefficient σ ⁰ is sensitive to snow grain size, wetness, and subsurface features and thus correlates with the penetration bias, which depends on the snow and ice properties. Therefore, the backscatter intensity was radiometrically calibrated according to

(1)$$\sigma ^0 = 10\cdot \log _{10}( {k_s\cdot {\vert {AMP} \vert }^2-NEBN} ) \cdot \sin \theta _{loc}$$

where AMP is the backscatter intensity in digital numbers (pixel intensity value); k _s is the calibration and processor scaling factor for the SAR signals annotated in the supplied metadata (Wessel, 2018), NEBN is the Noise Equivalent Beta Nought, which represents the influence of different noise contributions to the signal and is annotated in the supplied metadata in the form of polynomials over range with azimuth time tags describing the noise power, and θ _loc is the local incidence angle (Fritz and Eineder, 2008; Airbus, 2014).

3.2. Estimation of TanDEM–X penetration bias

A multiple linear regression model based on interferometric coherence and backscatter intensity according to Abdullahi and others (2019) was applied to the TanDEM-X datasets used in this study (i.e. the adjusted interferometric coherence and backscatter coefficient) to estimate the penetration bias for each acquisition date. The multiple linear regression model was calibrated and validated using a dataset acquired over the northern Greenland ice sheet in 2012. The calibration and validation dataset included observed penetration bias values (i.e. the difference between TanDEM-X DEM and IceBridge laser-altimeter measurements) as the dependent variable and corresponding interferometric coherence and mean backscatter coefficient values as independent variables at each IceBridge footprint. In total, the dataset comprises 87 097 data samples based on TanDEM-X acquisitions from 2 April and 10 April 2012 and IceBridge measurements collected from 30 March to 16 May 2012. The dataset covers a 36 km wide and 405 km long swath from the coast in the north at about 65 m a.s.l. southward in the direction of the center of the ice sheet up to about 2265 m a.s.l.. The data encompass all glacier facies (according to Benson (1962): dry snow zone, percolation zone, wet snow zone, superimposed ice zone, ice zone) and cover a range of penetration bias values from 0 to −8.5 m. For model calibration, 75% of the data were randomly selected, while the remaining 25% were used to assess the model accuracy (Fig. 4).

Fig. 4. Accuracy of the multiple regression model for penetration bias estimation fitted based on a dataset comprising about 65 000 observed penetration biases (i.e. difference between TanDEM-X DEM and IceBridge laser measurements height) samples over the northern Greenland ice sheet. (a) Comparison of observed and estimated penetration bias per IceBridge footprint; and (b) the corresponding distribution of residuals. (modified from Abdullahi and others (2019)).

Figure 4 illustrates the model accuracy in terms of the comparison of observed against estimated penetration bias values. The model achieved estimations with a standard error of 0.69 m and explains almost 70% of the variance of TanDEM-X penetration bias according to the coefficient of determination (R ²) (Fig. 4a). Observations and estimations are highly correlated with a correlation coefficient of 83%. Both independent variables (i.e. interferometric coherence and backscatter coefficient) are highly significant (p < 0.001) and the residuals possess a median value of −0.01 m, with lower and upper quartiles of −0.44 and 0.41 m, respectively (Fig. 4b).

Since the calibration data acquired over the northern Greenland ice sheet cover all glacier facies occurring on DIC (Koerner, 2005) and include the entire elevation range of the SGB (i.e. sea level to 1800 m a.s.l.), the model derived by Abdullahi and others (2019) was applied to the TanDEM-X datasets used in this study without any further adjustment.

3.3. Elevation change and geodetic mass balance

The calibrated DEMs and penetration biases estimated were bilinearly resampled from 6 m to 12 m resolution, and snapped to the same grid as the TanDEM-X DEMs, in order to match the independent pixel spacing of the interferometric processing filters as mentioned above. Subsequently, DEM differencing was performed to derive both the uncorrected and corrected elevation changes. All possible combinations for the time series (2010–2012, 2010–2013, 2010–2018, 2012–2013, 2012–2018, and 2013–2018) (Table 1) were considered.

Since the absolute height calibration within the MCP achieves high height accuracy (Wessel and others, 2018) it can be assumed that the relative vertical height bias between the DEMs has been sufficiently reduced such that further co-registration of the DEMs prior to differentiation is not required. Table 3 summarizes the mean errors and the std dev of the differences between the TanDEM-X DEMs over flat (slope < 20°) ice-free stable terrain, indicating no significant vertical offsets.

Table 3. Mean error (ME) and Std dev. (SD) of the difference between the TanDEM-X DEMs over flat (slope <20°) ice-free stable terrain

Before converting height changes into mass changes, the data gaps were filled with the hypsometric mean values of the elevation changes per 100 m elevation band. As described in Section 2.2.1., it can be assumed that the TanDEM-X data coverage is representative of the entire SGB.

Subsequently, the uncorrected and corrected volume changes were calculated on the basis of the elevation changes according to

(2)$$\Delta V = r^2\mathop \sum \limits_{i = 1}^N \Delta h_i$$

where ΔV is the total volume change over the basin in m³, r is the cell size in m, N is the number of pixels Δh is elevation change at pixel i in m and converted to mass changes assuming a mean density of 850 ± 60 kg m⁻³ according to Huss (2013)

(3)$$B_{geo} = \Delta V\cdot \rho $$

where B _geo is mass change in kg and ρ is the density in kg m⁻³.

3.4. Uncertainty assessment

3.4.1. Geodetic mass balance

The variance in elevation difference between the TanDEM-X DEMs over ice-free stable terrain is used to assess the uncertainty of the derived elevation changes. Layover, shadow, and regions with low coherence values <0.3 were excluded include only valid values in the uncertainty assessment. The Std dev. of the elevation difference σ ^Δh (SD, Table 3) over the ice-free and stable terrain regions ~ the SGB are used as an indicator for the uncertainty of elevation difference measures. This uncertainty is assumed to be representative not only for the glacier-off regions but also for the on-glacier regions. Since the elevation and mass balance changes are integrated over the entire SGB, for uncertainty estimation of the averaged elevation difference spatial correlation must be taken into account (Rolstad and others, 2009). The uncertainty of the spatial average of the elevation change rate σ ^(Δh/Δt) was assessed according to

(4)$$\sigma ^{{{\Delta h} \over {\Delta t}}} = \displaystyle{{\sigma ^{\Delta h}} \over {( \Delta t/\sqrt N ) }}$$

where Δt is the time period between the two DEMs and N is the number of uncorrelated observations included in the mean elevation difference (Abdel Jaber and others, 2019).

To determine the number of uncorrelated observations, multiple ranges of spatial correlation in the elevation differences were determined by means of empirical variograms obtained on the ice-free and stable terrain and fitted with a double spherical variogram according to Hugonnet and others (2022).

Following the rule of error propagation, the uncertainty of mass change rates σ ^(ΔM/Δt) is calculated according to

(5)$$\sigma ^{{{\Delta M} \over {\Delta t}}} = \left\vert {\mu^{{{\Delta M} \over {\Delta t}}}} \right\vert \cdot \sqrt {{\left({\displaystyle{{\sigma^{( \Delta h/\Delta t) }} \over {\mu^{( \Delta h/\Delta t) }}}} \right)}^2 + {\left({\displaystyle{{\delta A} \over A}} \right)}^2 + {\left({\displaystyle{{\delta \rho } \over \rho }} \right)}^2} $$

where μ ^(ΔM/Δt) is the mass change rate in Gt a⁻¹, δA the uncertainty related to the basin area in m², and δρ the uncertainty related to the volume-to-mass conversion in kg m⁻³.

3.4.2. Glaciological mass balance

To assess the uncertainty of the glaciological mass balance, the most recent estimate of the uncertainty (random and systematic) of ±250 mm w.e. a⁻¹ (Adams, 1966; Cogley and others, 1996) was applied, consistent with other studies on Devon by Koerner (1970) and on global data sets by Zemp and others (2009). Furthermore, an uncertainty of ±68 km² of the basin area was taken into account for the conversion of the glaciological specific mass balance into Gt according to Burgess and Sharp (2004).

4. Results

4.1. Snow surface height variability

Field data collected from automatic weather stations and annual mass balance stake surveys are used to estimate the potential impact of snowpack height variability on surface elevation validation and changes in this study. Precise measurements from two automatic weather stations, at DICS (1300 m a.s.l.) and DV1H (1800 m a.s.l.) indicate a mean difference in snow surface height of 6.5 cm (Std dev. of 6.7 cm) with maximum difference of 13 cm over the five- to six-month time periods between the TanDEM-X recordings and surface reference measurements (Table 4), which is well within the range of uncertainty of the TanDEM-X DEMs (Wessel and others, 2018) and reference data (Table 2) and can therefore be ignored in this study.

Table 4. Changes in snowpack depth corresponding to the time period between acquisition of the TanDEM-X and reference data used in this study

Snowpack height change data was collected from Automatic Weather Stations that span the entire range of the long-term accumulation zone within the SGB.

Negative snowpack height change value in 2011 likely due to wind scour. Negative value in 2018 due to reference data collected after TanDEM-X data.

4.2. Penetration bias

Figure 5 depicts the estimated penetration bias for each acquisition date. The mean values are in a similar range, with the 2010 dataset (Fig. 5a) showing the greatest absolute average penetration bias of about |−3.4| m and the 2012 dataset (Fig. 5d) the lowest absolute average penetration bias of about |−2.4| m. The spatial pattern of estimated penetration bias is different for each acquisition. In general, the lowest penetration bias is found at low elevations on the outlet glacier, increasing at higher elevations and decreasing again towards the ice cap summit. At all acquisition dates, except 2013, a zone of low penetration bias is found at altitudes between 1600 and 1800 m a.s.l. with values of between −5 and −1 m. This zone is particularly evident in 2010 and 2018. The adjacent zone at lower altitudes between 1500 and 1600 m a.s.l. with higher penetration biases is also particularly pronounced in 2010 and 2018. This pattern is reversed in 2013, where the zone of high penetration bias pronounced in the other years disappears and low penetration bias values ~ −2 m are observed between 1300 and 1500 m a.s.l.. Overall, it can be stated that although the penetration bias averaged over the basin of the individual acquisition dates is similar, the spatial pattern shows significant differences across the SGB.

Fig. 5. Estimated penetration bias from (a) December 2010, (b) November 2012, (c) December 2013, and (d) May 2018.

Variations in penetration bias among the acquisition dates shown in Figure 5 reflect the different subsurface properties formed by ablation and accumulation processes, as well as due to the different acquisition geometry (baseline, incidence angle, orbit, and viewing direction). The complex combination of these factors, some of which are constantly changing, leads to intra- and interannual variations in penetration bias of TanDEM-X DEMs. Interannual variations in the spatial pattern of penetration bias (Figs 5a–c) suggests that the magnitude of ablation and/or meltwater percolation that occurred the previous summer has important effect on these changes observed (König and others, 2001; Casey and Kelly, 2010).

In order to evaluate the accuracy of the penetration bias correction, the uncorrected and corrected TanDEM-X elevations were compared to the reference elevation datasets. As described above, the time lag between the TanDEM-X observations and the reference elevation data can be neglected, as the variability of snow surface height between winter and spring is negligible (Table 4). Figure 6 maps the spatial distribution of the reference data with respect to the individual penetration bias datasets. The corrected December 2010 TanDEM-X DEM and corresponding penetration bias were evaluated using large grid dGPS measurements and IceBridge data across the higher elevations of the ice cap (Fig. 6a). Temporally matching ArcticDEM strips encompassing nearly the entire range of estimated penetration bias (~300 to 1800 m a.s.l.) were used for validating the November 2012 and December 2013 datasets (Figs 6b and c), while the May 2018 dataset was validated with the localized (small grid) dGPS data along the Devon NW transect (Fig. 6d).

Fig. 6. Estimated penetration bias overlayed by the outline of the SGB and the footprints of the reference datasets. (a) Penetration bias of December 2010 TanDEM-X DEM with dGPS and IceBridge datasets from spring 2011, (b) Penetration bias of November 2012 and ArcticDEM strip from spring 2013, (c) Penetration bias of December 2013 with ArcticDEM strips from spring 2014, and (d) Penetration Bias of May 2018 with dGPS data from spring 2018.

Figure 7 shows the corresponding residual distributions, i.e. the differences between the uncorrected (dark blue) and corrected (turquoise) TanDEM-X DEM data, and the reference heights. Figures 7a1 and 7a2 show the comparison of the 2010 TanDEM-X DEM with dGPS (large grid) measurements and IceBridge data, respectively. Figures 7b and 7c show the comparison of the 2012 and 2013 TanDEM-X DEMs with ArcticDEM strips; and Figure 7d shows the comparison of the 2018 TanDEM-X DEM with dGPS (small grid) measurements. In all cases a significant improvement in height accuracy after penetration bias correction is found. While the mean deviation between the uncorrected TanDEM-X DEMs and the reference heights is about 3 m, the mean deviation between the corrected TanDEM-X DEMs and the reference heights is only up to 18 cm for the absolute mean values.

Fig. 7. Distributions of the residuals (i.e. the differences between TanDEM-X DEM and reference height) for the uncorrected (dark blue) and the corrected (turquoise) TanDEM-X elevations based on (a1) TanDEM-X DEM from December 2010 and dGPS measurements from spring 2011, (a2) TanDEM-X DEM from December 2010 and IceBridge measurements from 2011, (b) TanDEM-X DEM from November 2012 and ArcticDEM strip from spring 2013, (c) TanDEM-X DEM from December 2013 and ArcticDEM strips from spring 2014, and (d) TanDEM-X DEM from May 2018 and dGPS measurements from spring 2018.

Table 5 comprises four related summary measures characterizing the height accuracy based on uncorrected and corrected TanDEM-X DEMs in terms of the mean error (ME), the Std dev., median absolute deviation (MAD), and the root mean squared error of the residuals (RMSE) for each acquisition date.

Table 5. Summary statistics (mean error ME, Std dev. SD, median absolute deviation MAD, and root mean squared error RMSE of the residuals) characterizing the height accuracy of the uncorrected and the corrected TanDEM-X DEMs for each acquisition date

While the improved quality of the corrected TanDEM-X DEM relative to the uncorrected Data are reflected in the lower mean errors (ME), the corrected data show higher Std dev. relative to the uncorrected data. The RMSE values also indicate a significant improvement in height accuracy after correcting for signal penetration. While the mean RMSE between the uncorrected TanDEM-X DEM and the reference height averaged for all acquisitions is ~3.2 m, it is reduced by ~3-fold to 1.3 m after correction with residuals of 1 and 1.5 m. Although there is an increase in dispersion of the corrected compared to the uncorrected TanDEM-X DEMs (Table 5), the range of ME, SD, MAD and RMSE are consistent with the accuracy values for flat and unvegetated terrain reported by Wessel and others (2018).

In general, the validation results indicate that the modelled penetration bias corrections significantly improve the accuracy of TanDEM-X DEM data for elevation measurements over firn and ice surfaces.

4.3. Impact of TanDEM-X penetration bias on elevation change

Using the uncorrected and corrected TanDEM-X DEMs, the effect of penetration bias on elevation change can be evaluated. Differences between the uncorrected and corrected elevation changes are clearly evident in Figure 8 which shows the observed elevation changes for all possible combinations of TanDEM-X time series datasets based on the uncorrected and corrected DEMs. Similarly, discrepancies are evident in Figure 9 which depicts the corresponding elevation change per 100 m elevation bands based on uncorrected (dark blue) and corrected (turquoise) data.

Fig. 8. Elevation change based on the uncorrected and the corrected TanDEM-X DEMs and the corresponding difference between the uncorrected and the corrected elevation change for all possible combinations for the time series, i.e. (a1) 2010–2012, (a2) 2010–2013, (a3) 2010–2018, (b1) 2012–2013, (b2) 2012–2018, and (c) 2013–2018.

Fig. 9. Elevation change based on the uncorrected (dark blue) and the corrected (turquoise) TanDEM-X DEMs per 100 m elevation bands for all possible combination within the time series, i.e. (a1) between 2010 and 2012, (a2) between 2010 and 2013, (a3) between 2010 and 2018, (b1) between 2012 and 2013, (b2) between 2012 and 2018, and (c) between 2013 and 2018.

In all cases, the spatial pattern of the uncorrected elevation change across the SGB is significantly different from the spatial pattern of the corrected elevation change. In general, the more similar the magnitude and spatial pattern of the penetration biases of the DEMs, the less the influence of the penetration bias on the detection of elevation change. For example, the 2010 and 2018 TanDEM-X DEMs have very similar penetration biases (Figs 5a and d), which results in only small differences between the uncorrected vs corrected elevation changes measured (Figs 8a3 and 9a3) with, almost the same mean elevation loss for the SGB of about −2 m over the eight year period. By contrast, the magnitude and spatial pattern of the penetration biases from 2012 (Fig. 5b) and 2018 (Fig. 5d) are significantly different, causing considerable differences between the elevation changes derived from uncorrected vs corrected TanDEM-X DEMs (Fig. 8b2). For these years, an overestimation of the elevation decrease by ~60% using uncorrected TanDEM-X DEMs was found, with maximum discrepancies occurring between 700 and 1600 m a.s.l. (Fig. 9c). Overall, the results show a less negative surface height change estimate for time intervals 2010–2012, 2010–2013, and 2010–2018 (Figs 8a1–3 and 9a1–3) from differencing the corrected TanDEM-X DEMs, and more positive surface height changes for time intervals 2012–2013, 2012–2018, 2013–2018 (Figs 8b1–2 and c and 9b1–2b and c) compared to the uncorrected elevation change estimates.

4.4. Geodetic mass balance

Figure 10 shows the comparison of mass change rates for all possible combinations for the times series (2010–2012, 2010–2013, 2010–2018, 2012–2013, 2012–2018, 2013–2018) in the SGB. The comparison of the uncorrected geodetic mass change rates and the corrected geodetic mass change rates (Fig. 10a) illustrates the significant impact of signal penetration on the estimation of the mass balance using TanDEM-X DEMs. Depending on the combination of DEMs with their individual penetration biases, the use of uncorrected TanDEM-X DEMs results in significant over- and underestimates of mass change rates with an absolute deviation up to 0.29 ± 0.07 Gt a⁻¹. In general, the more similar the magnitude and spatial patterns of the penetration biases of the DEMs, the lower the distortion of the geodetic mass change rates will be when using uncorrected TanDEM-X DEMs. Compared to the glaciological mass change rate estimates, the corrected mass change rate (Fig. 10c) estimates fit much better than the uncorrected estimates (Fig. 10b) with mean absolute deviation of 0.10 ± 0.21 Gt a⁻¹.

Fig. 10. Comparison of mass change rates of the SGB on DIC between 2010 and 2018. (a) Uncorrected vs corrected geodetic mass change rate, (b) glaciological vs uncorrected mass change, and (c) glaciological vs geodetic mass change.

Considering the mass balance in mm w.e. a maximum absolute deviation between the uncorrected and the corrected mass balance is 780 ± 182 mm w.e. for the time period between December 2010 and November 2012, which corresponds to a misestimation of 409 ± 95 mm w.e. a⁻¹.

5. Discussion

In this study a multiple regression model based on interferometric coherence and backscatter intensity according to Abdullahi and others (2019) was applied to a time series of four TanDEM-X datasets to estimate the penetration bias and to quantify the impact of X-Band InSAR penetration bias on elevation and mass-change detection. We demonstrated that the effect of the penetration bias on elevation changes can vary tremendously depending on the combination of DEMs. In addition, small-scale variations and nonlinearity of the penetration bias with altitude influence the derived elevation and mass changes significantly. The results demonstrate that penetration bias corrections are necessary to minimize random errors in elevation changes measured and make reliable conclusions regarding mass gains and losses quantified. While comparisons between independent estimations of glacier mass balance demonstrates the applicability of geodetic mass-balance measurements from corrected TanDEM-X DEM elevation data, several potential complications persist.

5.1. Penetration bias

The accuracy at which TanDEM-X data detects surface heights across the SGB is entirely dependent on the accuracy of the modeled penetration bias. The estimated penetration biases show significant variations in magnitude and spatial pattern between acquisition dates and across the SGB. these variations in the penetration bias are not only due to the changing acquisition geometry from one acquisition date to the next (large baseline area and different viewing directions, Table 1), but also to different subsurface properties caused by varying climatic conditions. Since annual precipitation is extremely low on DIC (Braithwaite, 2005), interannual variations and longer-term trends of the ice cap are mainly caused by variations ins summer melt (Koerner, 2005). It can therefore be assumed that the interannual changes in penetration bias are mainly due to ablation processes in the summer months. 2012 was an exceptional warm summer season throughout the Arctic with high firn temperatures on DIC (Bezeau and others, 2013). The more pronounced melting season probably led to increased percolation and refreezing of meltwater into cold snow and firn at higher elevations forming ice lenses and layers near the surface, which acted as scatterers for the X-band signals and caused less signal penetration. In contrast, the summer season 2013 was characterized by outstanding cold surface temperatures (Mortimer and others, 2016) indicating a shorter and/or less intense melting season with less meltwater percolation which in turn could explain the increase in penetration bias. The variable melting conditions of the last decade on the DIC (Mortimer and others, 2016) and the associated shifts in the ablation zone led to varying subsurface structures across the SGB, which could explain the significant differences in the magnitude and spatial pattern of the penetration bias between 2010 and 2018

5.2 Validation and elevation change

Surface elevations across the ice cap can vary significantly over time due to changes in snowpack height, potentially affecting validation and/or height change measurements over intra- and interannual time frames. For all years in this study, intra-annual surface elevation measurements from TanDEM-X DEMs collected in November/December were compared with precision reference datasets (dGPS, NASA ATM, and ArcticDEM) collected in April/May. In-situ field measurements for these years indicate a maximum (average) change in snowpack height of 13 cm (6 cm) from December to May for three of the four years in which TanDEM-X data were validated. These results are consistent with long-term observations that >80% of the annual snowpack accumulates during the first 3 months of the autumn/winter (September 1 to June 1) period across this region (Koerner, 2005). Similarly, annual mass-balance field surveys along the Devon NW transect indicate an interannual variation of ±15 cm over the period when TanDEM-X data were collected for this study. Variability of this magnitude is ~12 times lower than the uncertainties associated with corrected TanDEM-X DEMs over the SGB. Therefore, the magnitude of intra- and inter-annual variability of snowpack height over the time periods required for validating TanDEM-X data (Table 2) and deriving multi-annual measurements of surface elevation change across the SGB, have negligible impact on the results of this study.

Estimating surface elevation change over glacierized regions from geodetic methods involve potential uncertainties related to ice dynamics and firn densification. Both factors can change the geometry across an ice cap or glacier surface while having no apparent impact on changes on the ice mass as a whole. Dynamic processes of submergence (emergence) that develop in response to long-term surface mass balance introduce a downward (upward) vertical ice motion across the accumulation (ablation) zone of a glacierized basin (Cuffey and Paterson, 2010). In addition, nonsteady horizontal ice motion associated with basal sliding can lead to localized thickness change through compressional thickening and extensional thinning. In all years which basin-wide coverage exceeded 98% of the total basin area (i.e. 2010, 2013 and 2018), biases associated with dynamic thickness changes were not considered as they cancel out at the basin-wide scale (Belart and others, 2017). In the 2012, a data gap due to missing acquisitions reduced coverage to ~78% of the total basin area. This data gap was however limited to elevations from 900 to 1600 m a.s.l. where ice is frozen to the bed (Burgess and Sharp, 2004) and average velocities are <20 m⁻¹ (Burgess and Sharp, 2004; Van Wychen and others, 2020). We therefore extrapolate the value of the changes measured for each 100 m elevation band across the areas for which coverage is missing in the elevation bands affected.

Comparisons of the derived uncorrected and corrected elevation changes clearly shows the high impact of the penetration bias on the detection of elevation changes. Uncertainties in the TanDEM-X DEMs associated with (i) sloping terrain (Table 5) and (ii) deviations in measured surface height relative to on-ice reference data (Table 3) indicate elevation differences measured in this study are only significant for changes exceeding ~1.8 m. In this study, annual thickness changes across the SGB measured average only ~40 cm per year.

5.3. Mass balance

Apart from mass changes at the ice cap surface due to accumulation and melt, iceberg calving and mass loss (or gain) at the glacier bed could also represent additional factors whereby ice is gained or lost from the SGB. Mass loss of ~0.01 Gt a⁻¹ due to ice-berg calving from the terminus of the Sverdrup Glacier however account for <−1 cm a⁻¹ thinning when averaged across the entire SGB and are therefore considered irrelevant in this study. Concerning basal processes, as more than 80% of the SGB is frozen to the bed (Burgess and Sharp, 2008), we assume that melt and accretion at the glacier bed also have a negligible impact on the thickness change observed.

Additional complications in deriving mass change from the TanDEM-X DEM differencing stems from using uniform density values for the entire regions below and above the equilibrium line. The use of single density values ignores the existence of different glacier facies characterized by different (sub-)surface structures, which are well reflected in the different penetration biases (Fig. 5). To minimize uncertainties related to using a single density value it is recommend that the time interval over which temporal change is measured is at least 5 years (Huss, 2013), which exceeds most of the time intervals for which elevation changes are measured in this study.

Finally, important uncertainties related to the effects of a warming firnpack across the SGB are often associated with internal accumulation, and densification which can result in significant thinning throughout the upper firn. Bezeau and others (2013) suggest that enhanced rates of firn densification in the near surface firn of DIC has accounted for up to 17 cm a⁻¹ surface lowering within the upper 1.25 m of the firn column alone. These findings imply that changes associated with ice cap thinning across the firn zone of the SGB (i.e. >1600 m a.s.l.) do not necessarily correspond to changes in ice mass.

6. Conclusion

The current study demonstrates the large impact of X-band InSAR penetration bias on elevation and mass-change detection over glacierized regions, and highlights the urgent need to correct for penetration bias to obtain reliable results. The penetration corrected TanDEM-X DEMs agreed to within ±14 cm of spatially and temporally coincident precise in situ kinematic dGPS data (±10 cm RMSE). In addition, geodetic mass change estimates derived from differencing TanDEM-X DEMs over multi-year periods between 2010 and 2018 showed good agreement with a mean deviation of 338 ± 166 mm w.e. with in situ mass change measurements collected over the same time-periods. These results demonstrate the applicability of TanDEM-X elevation data for cryosphere research and thus opens up the opportunity to use the enormous and growing database of area-wide high-resolution and highly accurate elevation information from the TanDEM-X mission in the context of glacier monitoring worldwide.