1. Introduction
Ongoing global warming is causing accelerating melting of ice sheets, ice caps and glaciers worldwide (IPCC, Reference Pörtner2022). 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, Reference Wong and Field2014). 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, Reference Berthier and Brun2019; Howat and others, Reference Howat, Porter, Smith, Noh and Morin2019; Dehecq and others, Reference Dehecq2020), (2) radar altimetry (Nilsson and others, Reference Nilsson, Gardner, Sørensen and Forsberg2016; Sørensen and others, Reference Sørensen2018; Jakob and others, Reference Jakob, Gourmelen, Ewart and Plummer2021), (3) gravimetry (Tapley and others, Reference Tapley2019; Sasgen and others, Reference Sasgen2020; Velicogna and others, Reference Velicogna2020), and (4) laser altimetry (Pritchard and others, Reference Pritchard, Arthern, Vaughan and Edwards2009; Sørensen and others, Reference Sørensen2011; Smith and others, Reference Smith2023).
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, Reference Krieger2007) 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, Reference Krieger2013). 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, Reference Berthier2023). Especially for monitoring the narrow, fast-flowing outlet glaciers of the ablation zones with the highest rates of elevation change (Vaughan and others, Reference Vaughan and Stocker2013; Cooper and Smith, Reference Cooper and Smith2019), 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, Reference Dall2007). 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, Reference Weber Hoen and Zebker2000). 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, Reference Abdullahi2019; Fischer and others, Reference Fischer, Jäger, Papathanassiou and Hajnsek2019; Rott and others, Reference Rott2021). 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, Reference Rignot, Jezek, Drinkwater and Lou1993; Liu and others, Reference Liu, Wang and Jezek2006). 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, Reference Dall2007). 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, Reference Dehecq2016; Leinss and Bernhard, Reference Leinss and Bernhard2021; Wessel and others, Reference Wessel2021).
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, Reference Abdel Jaber, Rott, Floricioiu, Wuite and Miranda2019; Farías-Barahona and others, Reference Farías-Barahona2020) and/or consider the penetration bias as an error contribution in the uncertainty estimation of the mass balance calculations (Braun and others, Reference Braun2019; Seehaus and others, Reference Seehaus2019; Krieger and others, Reference Krieger, Strößenreuther, Helm, Floricioiu and Horwath2020). 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, Reference Weber Hoen and Zebker2000; Dall, Reference Dall2007; Fischer and others, Reference Fischer, Papathanassiou and Hajnsek2020; Rott and others, Reference Rott2021) or to determine the penetration bias more generally with respect to zones of different backscattering (Rignot and others, Reference Rignot, Echelmeyer and Krabill2001; Rizzoli and others, Reference Rizzoli, Martone, Rott and Moreira2017). 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 (Reference Abdullahi2019) 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 km2 (Randolph Glacier Inventory, Version 6.0, 2017) of the northwest sector of the Devon Ice Cap (DIC) (Fig. 1). The DIC covers ~14 000 km2 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, Reference Boon, Burgess, Koerner and Sharp2010). Excluding the stagnant southwest arm (Burgess and Sharp, Reference Burgess and Sharp2004), ~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, Reference Van Wychen2017), with 30% discharged from the Belcher Glacier alone (Burgess and Sharp, Reference Burgess and Sharp2004). 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, Reference Gardner and Sharp2007). The intensity and duration of summer melt represents the primary control on total mass balance of glaciers and ice caps in this region (Koerner, Reference Koerner2005; Millan and others, Reference Millan, Mouginot and Rignot2017).
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, Reference Apollonio1961). 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, Reference Koerner2005; Burgess, Reference Burgess2017). 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, Reference Burgess and Danielson2022). 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, Reference Sharp2011).
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, Reference Fritz2011; Rossi and others, Reference Rossi, Rodriguez Gonzalez, Fritz, Yague-Martinez and Eineder2012). 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, Reference Studinger, Koenig, Seelye and Sonntag2010), 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, Reference Porter2022). 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, Reference Martin2012; Studinger, Reference Studinger2013).
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-XtremeTM 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, Reference Noh and Howat2015). 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 (Reference Nuth and Kääb2011). 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, Reference Cogley2011) 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, Reference Burgess2017).
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 (Reference Abdullahi2019) 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, Reference Gruber, Wessel, Huber and Roth2012). 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 (Reference Wessel, Bertram, Gruber, Bemm and Dech2016), 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, Reference Moreira2013). 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, Reference Rizzoli, Dell'Amore, Bueso-Bello, Carcereri and Martone2022). 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, Reference Weber Hoen and Zebker2000; Martone and others, Reference Martone, Rizzoli and Krieger2016; Rizzoli and others, Reference Rizzoli, Dell'Amore, Bueso-Bello, Carcereri and Martone2022). For the TanDEM–X datasets acquired under acquisition geometries given in Table 1, the interferometric coherence layers were adjusted according to Abdullahi and others (Reference Abdullahi2019) 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 (Reference Ashcraft and Long2005), the backscatter coefficient σ 0 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
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, Reference Wessel2018), 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, Reference Fritz and Eineder2008; 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 (Reference Abdullahi2019) 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 (Reference Benson1962): 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 (Reference Abdullahi2019)).
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 2) (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, Reference Koerner2005) 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 (Reference Abdullahi2019) 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, Reference Wessel2018) 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
where ΔV is the total volume change over the basin in m3, 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−3 according to Huss (Reference Huss2013)
where B geo is mass change in kg and ρ is the density in kg m−3.
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, Reference Rolstad, Haug and Denby2009). The uncertainty of the spatial average of the elevation change rate σ (Δh/Δt) was assessed according to
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, Reference Abdel Jaber, Rott, Floricioiu, Wuite and Miranda2019).
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 (Reference Hugonnet2022).
Following the rule of error propagation, the uncertainty of mass change rates σ (ΔM/Δt) is calculated according to
where μ (ΔM/Δt) is the mass change rate in Gt a−1, δA the uncertainty related to the basin area in m2, and δρ the uncertainty related to the volume-to-mass conversion in kg m−3.
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−1 (Adams, Reference Adams1966; Cogley and others, Reference Cogley, Adams, Ecclestone, Jung-Rothenhäusler and Ommanney1996) was applied, consistent with other studies on Devon by Koerner (Reference Koerner1970) and on global data sets by Zemp and others (Reference Zemp, Hoelzle and Haeberli2009). Furthermore, an uncertainty of ±68 km2 of the basin area was taken into account for the conversion of the glaciological specific mass balance into Gt according to Burgess and Sharp (Reference Burgess and Sharp2004).
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, Reference Wessel2018) 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, Reference König, Winther and Isaksson2001; Casey and Kelly, Reference Casey and Kelly2010).
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 (Reference Wessel2018).
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−1. 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−1.
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−1.
5. Discussion
In this study a multiple regression model based on interferometric coherence and backscatter intensity according to Abdullahi and others (Reference Abdullahi2019) 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, Reference Braithwaite2005), interannual variations and longer-term trends of the ice cap are mainly caused by variations ins summer melt (Koerner, Reference Koerner2005). 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, Reference Bezeau, Sharp, Burgess and Gascon2013). 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, Reference Mortimer, Sharp and Wouters2016) 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, Reference Mortimer, Sharp and Wouters2016) 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, Reference Koerner2005). 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, Reference Cuffey and Paterson2010). 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, Reference Belart2017). 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, Reference Burgess and Sharp2004) and average velocities are <20 m−1 (Burgess and Sharp, Reference Burgess and Sharp2004; Van Wychen and others, Reference Van Wychen, Copland, Burgess, Slaymaker and Catto2020). 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−1 due to ice-berg calving from the terminus of the Sverdrup Glacier however account for <−1 cm a−1 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, Reference Burgess and Sharp2008), 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, Reference Huss2013), 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 (Reference Bezeau, Sharp, Burgess and Gascon2013) suggest that enhanced rates of firn densification in the near surface firn of DIC has accounted for up to 17 cm a−1 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.
Acknowledgements
Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) – 457342986. Funding for D.B. provided by the Climate Change Geoscience Program, Lands and Minerals Sector, Natural Resources Canada. Logistical support in the field was provided to D.B. by the Polar Continental Shelf Program, Natural Resources Canada. The authors would like to thank the anonymous reviewers for many constructive comments on the manuscript.
Open access
