Inferring Mass Loss by Measuring Contemporaneous Deformation around the Helheim Glacier, Southeastern Greenland, Using Sentinel-1 InSAR

Abstract

The elastic response of solid earth to glacier and ice sheet melting, the most important consequences of climate change, is a contemporaneous uplift. Here, we use interferometric synthetic aperture radar (InSAR) measurements to detect crustal deformation and mass loss near the Helheim glacier, one of the largest glaciers in southeastern Greenland. The InSAR time series of Sentinel-1 data between April 2016 and July 2020 suggest that there is a maximum cumulative displacement of ~6 cm in the line of sight (LOS) direction from the satellite to the ground near Helheim. We use an exponentially decreasing model of the thinning rate, which assumes that the mass loss starts at the lower-elevation terminal region of the glacier and continues to the higher-elevation interior. A linear inversion of the derived crustal uplift in the vicinity of bedrock using this model for surface loading in an elastic half-space suggests a mass loss of 8.33 Gt/year, which agrees with the results from other studies.

1. Introduction

The balance of ice mass in glaciers, ice caps, and ice sheets is a very important issue associated with global climate change. Several studies have indicated rapid ice melting in Greenland over the last two decades [1,2,3,4]. The total volume of ice melting raises the global mean sea level at a rate of approximately 0.7–1.8 mm/year [5,6], which likely will increase in the future due to further global warming. The solid earth beneath glaciers is under the pressure of ice accumulation and ice melting. The response of the ground to current and past ice mass loss is an uplift signal, which consists of two components: (1) an instantaneous elastic component, and (2) a delayed viscous component (Figure 1) [7,8,9].

The elastic component indicates current uplift in response to the existing melting. The viscous component, on the other hand, is the response to deglaciation of the Pleistocene ice sheets that occupied a significant proportion of the Northern Hemisphere, and the advance and retreat cycles of the Antarctic and Greenland ice sheets since the Last Glacial Maximum (LGM) approximately 20,000 years ago [10]. The viscous component is related to flow in the mantle, more than 40 km below the surface, above which the lithosphere behaves like an elastically bending plate, resulting in displacements at wavelengths of several hundred kilometers [11]. To interpret ground deformation with respect to mass balance, it is essential to separate these spatial and temporal components (Figure 1).

Space geodetic techniques provide important observations for estimating ice mass balance. For example, laser and radar altimetry [2], space gravimetry [12], global navigation satellite system (GNSS) [13], and InSAR [14,15] have been used to resolve ice mass balance for glaciers, ice caps, and ice sheets in regions that cannot be easily accessed. The integration of such space geodetic techniques with ice thickness data and regional models can be used to estimate changes in mass balance [16,17].

Given the global availability of suitable SAR data derived from various national and international missions, spaceborne differential interferometric synthetic aperture radar (DInSAR) allows high-density measurements of surface deformation over large areas by combining the information from two SAR radar images acquired from approximately the same position in space over a length of time. DInSAR has been utilized for a variety of applications, such as volcanoes [18,19], earthquakes [20,21], landslides [22,23,24,25], surface subsidence due to groundwater pumping [26,27,28], glacier kinematics and surface unloading [29,30], and infrastructure monitoring [31,32,33,34].

The ability of DInSAR to measure large gradient displacements is limited due to problems arising from phase unwrapping [35], for which other techniques such as offset tracking can be used [36,37]. Furthermore, DInSAR fails if coherence, a measure of the radar reflection between two SAR observations, is lost due to, e.g., vegetation and snow cover [38]. Additionally, atmospheric artifacts can seriously affect the accuracy of DInSAR measurements, particularly if the rate of deformation is small (few cm/year—mm/year). To overcome such limitations, multiple interferograms using a stack of SAR acquisitions are generated and analyzed using InSAR time series techniques to obtain a time series of surface deformation. This approach is known as multi-temporal InSAR (MT-InSAR), which can be divided into three main categories: (1) persistent scatterer InSAR (PSI or PSInSAR), focusing on permanent/persistent scatterers (PS) defined as individual pixels with stable scattering properties over time, (2) small baseline subset, or SBAS, focusing on distributed scatterers (DS) defined as pixels whose neighboring pixels share similar statistical behavior [39], and (3) combined MT-InSAR methods that incorporate information from both PS and DS [40,41,42].

The InSAR-based measurements from older SAR missions such as ERS and Envisat presented a major challenge for assessing the elastic response to ice sheet melting due to their small temporal resolution, i.e., their sampling intervals of several weeks. This disadvantage has been improved with the launch of the Sentinel-1A satellite in 2014 and the Sentinel-1B in 2016, providing a repeat cycle of 6 days over Europe and 12 days for the rest of the world [43]. Sentinel-1 has four operational modes (SM, IW, EW, and WV). Interferometric wide swath (IW) is its main operational mode over land, featuring data products in single (HH or VV) or double (HH plus HV or VV plus VH) polarization at 5 $\times$ 20 m spatial resolution and a 250 km swath. Sentinel-1 SAR images are free and available for all users around the world [43].

Since the early 1990s, the Greenland ice sheet has experienced an increase in mass loss. Space geodetic techniques such as satellite altimetry [44,45,46], satellite gravity anomaly experiments [47,48,49], and flux-balance calculations [14,16] reveal an increase in the mass loss rate between 1992 and 2011. Recently, the observations from the CryoSat-2 altimeter have demonstrated that the volume loss rates for the period 2011–2014 were 2.5 times greater than those for the period 2003–2009 [50]. The rapid unloading of ice from southeastern Greenland has caused an elastic uplift of ~16 mm/year at a GPS site on the Helheim glacier between 2001 and 2006 [1].

This study is one of the few to use SAR data from the more recent Sentinel-1 satellite to assess ice-unloading-induced deformations at rocky outcrops near glaciers. Previous studies such as Liu et al. (2012), Auriac et al. (2013) [51], or Zhao et al. (2014) all predated the Sentinel-1 era: Liu et al. (2012) used RADARSAT-1 SAR data to measure elastic crustal uplift in west Greenland, Auriac et al. (2013) utilized ERS-1, ERS2 and Envisat SAR data to constrain uplift of the ground all the way up to the edge of the largest ice cap in Iceland, and Zhao et al. (2014) used C-band ERS imagery and the SBAS method for measuring the contemporaneous deformation in Iceland.

In this paper, we use the latest archive of Sentinel-1 to assess the contemporaneous displacement due to mass loss at the Helheim glacier, one of the largest glaciers in southeastern Greenland [66.5°N, 38°W] (Figure 2). Previous studies suggest that the rate of thinning of the Helheim glacier has not been constant in the last two decades. For example, Stearns et al. (2007) [52] analyzed sequential digital elevation models (DEMs) and showed that the thinning rate was 44 ± 21 mm/year between 2002–2003, 15 ± 21 mm/year between 2003–2004, 60 ± 13 mm/year between 2004–2005, and 7 ± 21 mm/year between 2005–2006. Using optical and SAR remote sensing data analysis, Bevan et al. (2012) [53] suggested a greater glacier discharge of the Helheim glacier since the mid-1990s. In this research, multi-temporal InSAR (MT-InSAR) analysis using the SBAS method is carried out to constrain the spatial and temporal pattern of the thinning rate for the period 2016–2020. The thinning rate is calculated by an exponential decreasing model which assumes the mass loss begins from the lower-elevation near the edge of glacier towards the higher-elevation interior and by a linear inversion-derived surface loading in an elastic half-space.

2. Materials and Methods

2.1. Interferometric Analysis of Crustal Uplift

For the MT-InSAR analysis, we used GMTSAR software and the SBAS technique. The interferometric processing in GMTSAR comprises three major elements: (1) a pre-processor for each type of satellite data to convert the native format and orbital data into a generic form; (2) an InSAR processor to coregister a stack of images and convert topography into phase; and (3) a post-processor built on GMT [54,55], which is used for filtering interferograms and creating interferometric products such as coherence and line-of-sight displacement in both geographic and radar coordinates. The small baseline subset (SBAS) technique [39] was utilized as implemented in GMTSAR software [56] to observe the temporal and spatial changes in crustal deformation around the Helheim glacier. The SBAS technique uses all interferometric pairs with short baselines, either in temporal or spatial terms, to estimate the unknown deformation at each acquisition time. Including a large number of interferometric pairs in SBAS time series analysis, on one hand, reduces temporary uncorrelated effects caused by the atmosphere or decorrelation noises; on the other hand, this approach increases the density of measurement points on the resultant displacement maps [39,57,58].

Our study area was well covered by SAR images from Sentinel-1A that were collected in the 5 years from April 2016 to July 2020. In glacial regions, acquisitions in seasons with less snowfall are more suited for interferometric processing than those in other seasons due to snow cover that causes significant signal decorrelation in C-band interferograms. We used 73 SLC images in IW mode that were acquired in a descending track (path 112 and frame 370). We did not use ascending track images as all of them were acquired in winter (the last days of December to the beginning of February), hindering us from making a reliable network for SBAS analysis. Excluding the winter images leads to missing the elastic response of the earth due to the snow load collected in winter, and makes it impossible to identify the seasonal signal.

We used HH polarization data, which have a higher signal-to-noise ratio. Unwrapping of the differential interferograms was performed using the Snaphu technique [59]. We adapted a multilooking factor of 16 $\times$ 4 in range and azimuth to further suppress the noise and improve the efficiency of phase unwrapping. We used 2 by 2 m elevation cells for topography phase correction and geocoding. These cells were designed by the National Science Geospatial-intelligence Agency (NGA) and the National Science Foundation (NSF) public-private initiative, and they automatically produce a high-resolution, high-quality digital surface model (DSM) of the Arctic (ArcticDEM).

For the SBAS analysis, we focused on interferograms with small perpendicular and temporal baselines (<150 m and 150 days) to ensure that errors imposed by topography and temporal decorrelation remained small, and we minimized orbital gradients across the interferograms by utilizing coherent data from the entire SAR frame to flatten the phase. A coherence threshold of 0.3 was used to select reliable pixels as measurement points after phase unwrapping. To minimize the effect of temporal decorrelation caused by interferograms covering the winter seasons, the time period has been divided into five phases, each covering a different summer season. The networks of each year are shown in Figure 3; interferograms with low quality were excluded from the networks before the least-squares inversion.

Usually, orbital errors appear as a large extended signal, formed like a ramp that does not affect local deformations and cannot be mistakenly interpreted as a deformation signal. Fattahi and Amelung (2014) [60] showed that for Sentinel-1, the range uncertainty depends on the magnitude of the orbital errors, about 0.2 mm/year per 100 km, and the azimuth uncertainty is about 0.5 mm/year per 100 km. We used precise orbital determination (POD) data from ESA (https://scihub.copernicus.eu/gnss) (accessed on 29 August 2020) to filter out orbital errors from the interferograms.

The stratified atmospheric delay is a major error source in the SBAS time series analysis which is computed in GMTSAR by a method that uses a common point stacking technique to estimate the delays in the radar phase induced by propagation through the atmosphere [61]. It is worth noting that at high latitudes, it is possible for InSAR measurements to also be influenced by ionospheric disturbances. However, this effect is usually insignificant for C-band data [50,62] and appears as a ramp in high latitudes such as Greenland (https://iono.jpl.nasa.gov) (accessed on 24 May 2020). Accordingly, we estimated this effect by fitting a first order polynomial to the interferograms during the SBAS analysis. The set of small baseline interferograms were finally translated into a range-change time series that describes the time-dependent surface deformation by a least-squares inversion [63].

2.2. Ice Mass Loss Modeling

Ice melting in glaciers results in unloading that causes crustal uplift due to the viscoelastic deformation of the Earth’s crust and mantle. The viscous component is a long-wavelength phase contribution to glacial rebound and needs many years of observations to be determined due to its rather small uplift rate. Thus, for a short period of time such as the one considered in this paper, the viscous effect is very small and was ignored in our analysis.

To model the surface displacements observed by the MT-InSAR measurements, we use a model developed by Zhao et al. [30]. The model assumes a unit point mass in elastic half-space cylindrical coordinates with axial symmetry [30,64], as follows:

$$V\left(r\right)=\frac{g}{2r\pi }\frac{1-{\vartheta }^2}{E}\frac{1}{2}$$

(1)

$$U\left(r\right)=\frac{g}{2r\pi }\frac{\left(1+\vartheta \right)\left(1-2\vartheta \right)}{E}$$

(2)

where $V$ is the vertical and $U$ is the radial displacement, ϑ is the Poisson’s ratio, $g$ is the gravitational acceleration, $r$ is the distance from the point source, and $E$ is the Young modulus. Under the assumption that the mass loss rate remains constant with time, it is possible to use velocities instead of displacements in Equations (2) and (3). As the trend of ice melting decreases from the edge of the glacier to the internal parts, an exponential decay function to describe the thinning rate in the margin of the glacier is considered [30]:

$$t\left(x\right)=b+\left(a-b\right)e^{-\frac{x}{h}}$$

(3)

where $t$ is the thinning rate, $h$ is the decay factor, $b$ indicates the thinning rate of the inland portion (away from the ice edge), $x$ denotes distance from the ice edge, and $a$ is the thinning rate of the ice edge ($x=0$). A negative $b$ means ice accumulation at an area away from the ice edge. In an attempt to extend this one-dimensional model (Equation (3)) into a two-dimensional model, the load change (Equations (1) and (2)) is split into a series of square blocks and the displacement rate is computed for each block. In the next step, solutions for all blocks are added as follows [65]:

$$d_l^v={\displaystyle \sum }_{i=1}^N\left(b+\left(a-b\right)e^{-\frac{x_i}{h}}\right)\frac{{\Delta }^2\rho g\left(1-{\vartheta }^2\right)}{\pi E{r}_{l_i}}$$

(4)

$$d_l^{e,n}={\displaystyle \sum }_{i=1}^N\left(b+\left(a-b\right)e^{-\frac{x_i}{h}}\right)\frac{{\Delta }^2\rho g\left(1+\vartheta \right)\left(1-2\vartheta \right)}{\pi E{r}_{l_i}}$$

(5)

where $d$ stands for displacement, and the vertical, north, and east directions are depicted by the superscripts $v$, $n$, and $e$, respectively. $x_i$ indicates the distance of an ice block from the ice edge for i = 1…, $N$ ($N$ is the number of blocks), $\rho$ is the ice density, $\Delta$ denotes the block dimension, and $r_{l_i}$ represents the distance of point load $i$ from observation point $l$ with $l$ = 1, 2, …, $L$ ($L$ indicates the number of observations). For a given $h$, a linear system of equations can be constructed as follows:

$$L=AM$$

(6)

where $L$ is the vector of observation velocities (LOS), $M$ is the vector of model parameters (the unknown parameters $a$ and $b$), and $A$ is the design matrix of the equation system, which acts as a mapping function. The $l$-th line of $A$ reveals mapping of the mass loss rate onto the $l$-th velocity field [30]. This system of equations has more equations than variables and is an overdetermined system of linear equations that can be solved using the least squares technique [66]:

$$M={\left(A^{T}A\right)}^{-1}{A}^{T}L$$

(7)

$${\sigma }_M^2={\displaystyle \sum }_l^{}{\left({(A^{T}A)}^{-1}{A}^T{\sigma }_l\right)}^2$$

(8)

where ${\sigma }_M^2$ is the model standard error of the $i$th ice block and ${\sigma }_l$ is the standard error of the $l$th velocity field. $A$ in Equation (6) is constructed from the two matrices of $B$ and $C$ as follows:

$$A=B\times C$$

(9)

$$\begin{array}{c}B=\frac{{\Delta }^2\rho g}{\pi E}\left[c_1\left(\frac{1-{\vartheta }^2}{r_{11}}\right)\dots s{c}_1\frac{(1+\vartheta )(1-2\vartheta )}{2r_{11}}⋮⋮\dots s{s}_1\frac{(1+\vartheta )(1-2\vartheta )}{2r_{11}}\dots ⋮\right.\hfill \\ \hfill \left.⋮c_1\left(\frac{1-{\vartheta }^2}{r_{L1}}\right)\dots s{c}_1\frac{(1+\vartheta )(1-2\vartheta )}{2r_{L1}}\dots s{s}_1\frac{(1+\vartheta )(1-2\vartheta )}{2r_{L1}}\dots \right]\end{array}$$

(10)

$$c_1=cos{\theta }_1$$

(11)

$$s{c}_1=-sin{\theta }_{1}cos{\alpha }_1$$

(12)

$$s{s}_1=sin{\theta }_{1}sin{\alpha }_1$$

(13)

$$\begin{array}{c}C=[1-e^{-\frac{x_1}{h}}⋮1-e^{-\frac{x_N}{h}}{e}^{-\frac{x_1}{h}}⋮e^{-\frac{x_N}{h}}1-e^{-\frac{x_1}{h}}⋮1-e^{-\frac{x_N}{h}}{e}^{-\frac{x_1}{h}}⋮e^{-\frac{x_N}{h}}1-e^{-\frac{x_1}{h}}⋮1-e^{-\frac{x_N}{h}}{e}^{-\frac{x_1}{h}}⋮e^{-\frac{x_N}{h}}]\end{array}$$

(14)

where B is a $2L\times 3N$ matrix generated from Equations (4) and (5) in which the first $N$ columns correspond to the vertical component of displacement, the second $N$ columns correspond to the east component, and the third $N$ columns correspond to the north component. $c_1$, $s{c}_1$, and $s{s}_1$ contain the mapping parameters from the three-dimensional velocities to LOS direction, $\alpha$ is the looking angle of the satellite, $\theta$ is the incidence angle, and $C$ is a $3N\times 2$ matrix for the thinning model parameters. The mass loss estimated by this approach is reasonable for glaciers with exponential thinning from the edge of glaciers towards the interior of glaciers. Therefore, a thickness reduction of the load points that are estimated by the model is multiplied the area and ice density to find the mass loss volume and ice mass loss, respectively. The modeling approach discussed above was implemented in an in-house Matlab code.

The overall flowchart of the InSAR processing and model development is illustrated in Figure 4.

2.3. Auxiliary Data

To investigate the relation between the local land surface temperature (LST) on the ice mass loss in the study area, we used LST data from Moderate Resolution Imaging Spectroradiometer (MODIS). MODIS is a sensor operating on the Terra and Aqua satellites, which were launched by NASA in December 1999 and May 2002, respectively. It helps scientists to determine the amount of water vapor in a column of the atmosphere and the vertical distribution of temperature and water vapor, measurements that are crucial to understanding Earth’s climate system. The MODIS Land Surface Temperature and Emissivity weekly data (MOD11C2) in the study area are collected for the time span of InSAR observations, and the LST measurements were extracted at a 1 km resolution to use them as local land surface temperature data.

3. Results

3.1. InSAR Velocity Field

Figure 5 shows the LOS displacements resulting from the SBAS analysis for the five summer seasons between 2016 and 2020. Although the area is glacial and rugged and snow and ice cover reduce the coherence of the interferograms, InSAR has been able to show a clear signal of deformation. All five plots in Figure 5a–e show a 15–25 ${\mathrm{km}}^2$ wide yellow-red area in the northern flank of the Helheim glacier, illustrating uplift in an LOS direction. The rate of uplift, however, is not constant. The maximum rates of approximately 18 mm/year occurred in the summer 2019 at the edge of the ice sheet in the northern flank of the Helheim glacier. The average standard error of the LOS velocity fields is 0.59 mm/year, following Zhao et al. [30], which was calculated in a non-deforming area by using the velocity of the points in that area, as shown by the cyan rectangle in Figure 5a–e, and then was averaging it as a standard error for all the LOS velocity fields.

To better illustrate the temporal behavior of the displacement field, we selected a point (P) near the ice edge and plotted its LOS time series for all five periods. The velocity of the point was estimated using a linear regression (a constant velocity at each period). As illustrated in Figure 6a, point P shows uplift with a rate of 15 mm/year in summer 2016, decreasing by up to 46% and 26% in 2017 and 2018, respectively, and increasing again by 20% to 18 mm/year in summer 2019. The rate of uplift at point P in summer 2020 is reduced by more than 60% as compared to summer 2019.

In order to provide an assessment of the overall cumulative displacement in the summertime from 2016 to 2020, we assume that there is not any significant deformation in the period between November to March, and we merge the individual time series in Figure 6a by shifting the beginning of each time series to the end of the previous one in order to derive an estimation for the cumulative uplift during 2016–2020 (Figure 6b). The results show a cumulative displacement of up to 58 mm in the summer of the years 2016–2020, a long-term LOS uplift rate of 12.7 mm/year (assuming a constant velocity). The fluctuations with respect to the linear trend observed in Figure 6a,b are caused either by transient melting during the summer months, the remaining atmospheric artefacts from InSAR processing, or the temporary effect of precipitation in the elastic response of the crust in this area. To study the mass loss that caused the observed uplift in Figure 5, an area containing the uplift on the northern flank of Helheim is selected (black rectangle in Figure 5a–c) and the LOS velocities are used as observations in the model equations.

3.2. Mass Loss Estimation

The method described in Section 2 is applied to model displacements and to estimate the mass loss rate for the five observation periods. An area around the northern flank of Helheim (the black rectangle in Figure 5; area ~1500 km² (30 km × 50 km)) was selected for the modeling. To construct the model, we chose a snow block with a radius of 30 km around the velocity field as a load in the model. The load change was applied via a gridded point source based on square blocks with a size of 1 km. The decay factor h in Equation (4) was changed from 1 to 30 km with a step size of 1 km (Figure 7) to obtain the best decay factor in the model.

The best fit between observations and models occurred at h = 7500 m. The cumulative estimated ice volume loss of all periods of observation is 8.65–9.53 ${\mathrm{km}}^3$/year, equivalent to an ice mass loss of 7.93–8.75 Gt/year, assuming a snow density of 916.7 kg/m³. This range of 0.81 Gt/year between the minimum and maximum suggests that variations in h are largely compensated by variations in the parameters a and b. The fit of the model to the data is presented by the root mean square error: $\mathrm{RMSE}=\sqrt{\frac{{{\displaystyle \sum }}_{i=1}^L{\left(l_i-m_i\right)}^2}{L}}$, where $l_i$ is the observation, $m_i$ is the model prediction, and L is the number of observations. We assume unit variance for all observation points and observe that as $h$ changes, the RMSE does not vary significantly (h = 1–30 km, RMSE = 1.58–1.67 mm/year). Therefore, the data are not sufficient to resolve h. According to the residuals of the data and model, the best fit of the model occurs at h = 7500 m, which we selected to find the parameters of the model to estimate the ice mass loss for each observation period. The average of the model parameters for all observation time spans is a = 5.07 m/year and b = −2.42 m/year. The model parameters and the estimated ice mass loss for each period of observation are listed in

Table 1. Model parameters and ice mass loss and its averages estimated for the five periods of observation.

Model Parameters	2016/Apr–2016/Oct	2017/May–2017/Oct	2018/May–2018/Oct	2019/Apr–2019/Oct	2020/Apr–2020/Jul
a (m/year)	5.32	4.5	5.14	6.03	4.36
b (m/year)	−2.78	−2.14	−2.11	−2.90	−2.17
Ice mass loss (Gt)	7.38–7.875	5.625–6.39	6.48–7.02	10.67–11.56	5.63–6.30

.

By comparing the average of the ice mass loss in the years 2016–2020, we found that the maximum ice mass loss occurred in 2019. The increase in ice mass loss in 2019 in comparison to 2020, 2018, 2017, and 2016 is about 46%, 39%, 45%, and 31%, respectively.

In order to assess how our model illustrates ice sheet/block melting, we plotted two-dimensional average thinning rates of the ice sheet near the edge of the Helheim glacier in Figure 8. The figure shows how the thinning rate decreases from the front of the ice sheet to the interior. On average, the maximum thinning rate of 5.07 m/year occurs at the edge of the ice sheet at the northern flank of the Helheim glacier. By moving away from the edge, the ice sheet height increases while the thinning rate decreases.

4. Discussion

4.1. InSAR and Modeling

In this study, we investigated how InSAR-based measurements of spatial and temporal variations in crustal uplift along the northern flank of the Helheim glacier in southeastern Greenland could be used to model and estimate the amount of ice melting in the region as an unloading problem. The model that we used to detect the thinning rate in the study area is based on the direct effect of ice melting near the edge of the ice sheet. As we assumed that ice melting decreases from the edge of a glacier to its interior, we selected a portion of the InSAR results for modeling, which was located close to the edge and clearly shows the rebound effect.

The InSAR displacement time series results showed a maximum cumulative LOS uplift of 58 mm for 2016–2020 and 40.86 ${\mathrm{km}}^3$ ice volume loss in an area of approximately 1500 ${\mathrm{km}}^2$ in the time span from 2016 April to 2020 July. Simpson et al. (2011) [67] showed that the viscous uplift rates for their best fitting earth model are about ± 0.2 mm/year all around Greenland, although Khan et al. (2016) [68] argued that the viscous deformation is larger in this region and estimated a viscous uplift rate of about 3–5 mm/year in southeast Greenland. Therefore, we also accounted for the viscous component in our results and removed the modeled viscous deformation (4 mm/year) from the InSAR measurements before inferring the present-day mass change. Our findings suggest that ignoring the viscous component leads to an estimation of 7.92 Gt/year for the ice mass loss, and so it will change our results by only 5%, indicating the result was not sensitive to the viscous deformation.

It is worth noting that for homogeneous elastic half-space modeling, it is important to justify our Young modulus, E, because the surface displacement due to loading is inversely proportional to the Young modulus (Equations (1) and (2)). Auriac et al. (2013) showed that one method to select the appropriate Young modulus is by minimizing the residuals between the observed and modeled values of the InSAR phase. Therefore, we applied this method and varied E between 1 GPa and 90 GPa. The minimum value for the residuals occurred at E = 48 GPa, which we used in our model.

The Poisson’s ratio of 0.23 in our model was adapted from the study by Dahl-Jensen et al. (2003) [69]. To consider the sensitivity of our model to this factor, we changed it from 0.2 to 0.4 by 0.01 steps and considered the results. The RMSE of the residuals did not change significantly, indicating that the model is not sensitive to this factor.

Figure 9 shows the comparison between observations (a, d, g, j, and m) and the models (b, e, h, k, and n) with h = 7.5 km. In this area, the model predictions closely resemble the observations (see the residuals in Figure 9). The standard deviation of the residuals is 2.52–3.77 mm/year and the RMSE between the observation and model varies from 1.493 to 1.789 mm/year for the observation periods (

Table 2. Standard deviation of the residuals and the RMSE between the observations and models of the five periods.

Statistic Parameters	2016	2017	2018	2019	2020
Standard deviation of residual (mm/year)	2.52	2.83	3.02	3.77	3.26
RMSE (mm/year)	1.493	1.612	1.657	1.789	1.665

).

Despite the strongly scattered measurement points in this region, the model aligns well with the observations. To show the consistency of the observations and model, we considered two profiles, A′A and B′B, at the edge of the ice sheet for each observation period, and we compared the discrepancy between the observation and model along the profiles. The plots in Figure 10 illustrate how the model fits to the observation. The RMSE values are between 1.467 and 1.875 mm/year for A′A and between 1.152 and 1.698 mm/year for B′B.

In order to consider the absolute difference between the observation and model in Figure 10, we used the average of absolute percentage deviation (AAPD) [70] as follows:

$$\mathrm{AAPD}=100\times \frac{{{\displaystyle \sum }}_{i=1}^n\left(\left|\frac{Ob{s}_i-Mode{l}_i}{Ob{s}_i}\right|\right)}{n}$$

(15)

where $|.|$ computes the absolute values. The minimum (maximum) values of the AAPD correspond to the best (worst) performance of the model to the observation.

Table 3. AAPD between the observations and models of the five periods in Figure 10.

AAPD (Percentage)	2016	2017	2018	2019	2020
A′A	3.4	3.7	3.1	4.9	3.9
B′B	1.8	2.8	4.2	4.7	4.3

shows that in Figure 10, the best performance of the model is 1.8%, while the worst is 4.9%. Therefore, the model is fitted to the observations with a maximum five percent difference and the RMSE values are acceptable.

4.2. Inferring Elastic Rebound from GPS

In order to infer the elastic rebound from the GPS, we used a GPS time series of displacement recorded at the permanent GPS station, Hel2, which is the only GPS station in the study area, from http://geodesy.unr.edu/gps_timeseries/txyz/IGS14/HEL2.txyz2 (accessed on 14 January 2021). The data were analyzed in the Nevada Geodetic Laboratory using GipsyX software that was developed at JPL. The precision of the GPS components is reported 0.9, 0.8, and 3 mm in the north, east, and up directions. As we had only one GPS station in the study area, we could not directly compare the SBAS and GPS measurements with each other; as DInSAR is a relative measuring system, we would need at least two GPS stations for comparison. Therefore, in the following, we do not focus on the direct comparison between the GPS and InSAR, but rather on the interpretation that we infer from the GPS. To do this, we first converted the GPS time-series displacement into the LOS direction by using the following equations:

$$d_{LOS}=d_u\left({\theta }_{inc}\right)-d_{ALD}\left({\theta }_{inc}\right)$$

(16)

$$d_{ALD}=d_{e}coscos\left({\alpha }_h\right)-d_{n}sinsin\left({\alpha }_h\right)$$

(17)

where $d_n$, $d_e$, and $d_u$ are the three components of the displacement vector $d$ in the north, east, and up directions, respectively. ${\theta }_{inc}$ is the incidence angle and ${\alpha }_h$ is the azimuth of the satellite heading. $d_{ALD}$ is the projection of the horizontal components on the azimuth look direction $ALD$ [71].

Figure 11 shows InSAR and GPS time series for four periods (due to the gap in the GPS data for 2020, the years 2019–2020 are assumed as one period), where the former was calculated based on the average of points within a radius of 100 m around the permanent GPS station Hel2. Both the GPS and InSAR time series show an increasing trend in the LOS displacement time series. In general, the GPS time series data are noisier than the InSAR, and they have much more fluctuations. The higher fluctuations in the GPS time series might be due to the continuous record of data that can affected by environmental effects such as rain and snowfall. In contrast, the InSAR data are discontinuous and are collected every 12 days, and so they may miss such effects. In addition, due to the employment of temporal filters in the SBAS analysis, the InSAR displacement time series are much smoother than the GPS, which may mask part of the ground displacement related to ice melting. As seen in Figure 11, the maximum LOS displacement in the GPS is more than the InSAR. It may be because the GPS time series is in an absolute reference frame, while the InSAR is relative to a reference area. The reference area experiences many of the same geophysical phenomena as the area close to the GPS, e.g., the seasonal elastic deformation due to snow cover would be similar for the area around the GPS station and for the reference area. The differencing will consequently remove a part of the displacement signal in the InSAR, while the signal remains in the un-differenced GPS. This causes the differences in the amount of uplift between the GPS and InSAR time series. The absolute precision of the GPS is ~3 mm and the relative precision of the InSAR is ~0.6 mm. Therefore, the integration of the InSAR measurements from the missions with shorter revisit times (1–3 days) with more continuous GPS stations in the region would improve the capability for a better interpretation of non-linearity related to ice melting in different seasons and glacier areas.

4.3. Comparison with Previous Studies

As the observation period covered by our SAR data is relatively short and the detected average uplift of ~13 mm/year around the Helheim glacier in our study can be assigned to the elastic component of the crust, triggered by the very quick response of crust to atmospheric and oceanic fluctuations in short time scales. Khan et al. (2007) used the dem-extraction method to estimate the volume changes in Helheim and determined that there was 16.6 mm of uplift in this region and 180 ${\mathrm{km}}^3$ volume loss for 2002–2005 in an area larger than 48,000 ${\mathrm{km}}^2$. Our study has revealed a maximum cumulative line-of-sight uplift of 58 mm during 2016–2020, corresponding to a cumulative volume loss of 40.56 ${\mathrm{km}}^3$ in an area of about 1500 ${\mathrm{km}}^2$. If we divide the volume loss by the area in these two studies, the results show that the peak rates of ice melting near the Helheim glacier increased by almost four times between 2016 and 2020. This temporal acceleration of ice melting has also been mentioned in previous studies. For example, Bevis et al. (2012) used GPS data recorded at the Hel2 station and showed that the velocity of the vertical displacement was 11.8 mm/year between 2007 and April 2010, but it increased to 15.3 mm/year between April 2010 and March 2011, which the authors attributed to the summer of 2010 being the hottest and longest on record over most of southern Greenland (until 2011). Helm et al. (2014), using Cryosat-2 altimeter observations, showed that the volume loss rate at the Helheim glacier between 2011 and 2014 was 2–5 times greater than that between 2003 and 2009. Wake et al. (2016), using GPS data from Hel2 for 2008 to 2013, suggested that Hel2 was uplifting at a rate of 16 ± 0.36 mm/year, with significant acceleration (3.7 ± 0.24 mm/year from April 2006 to September 2012). Mouginot et al. (2019) [72], using an extensive survey of surface elevation, thickness, velocity, and surface mass balance, suggested that the mass loss acceleration shifted from positive in the 2000–2010 period to negative in the 2010–2018 period. They suggested a 20 Gt/year mass loss for all glaciers in southeastern Greenland over 46 years. Sasgen et al. (2020) [73] observed a mass loss in 2017 and 2018 by GRACE and GRACE-FO and reported that these were lower than in any other two-year period between 2003 and 2019. In June and July 2019, the ice sheet experienced high rates of melt water runoff close to the maximum rates observed over the full 2003–2018 period. Williams et al. (2021) [74] used a range of remotely sensed datasets and showed that during 2014–2020, the Helheim glacier had continued to lose mass and thin, and it had retreated inland, up to 100 m thinner than it was during 2005. Our results also show that in 2017 and 2018, the LOS displacements related to ice mass loss were lower than other year periods between 2016 and 2019, and in 2019, the most LOS displacement took placed. Our InSAR-based analysis indicates that the Helheim glacier plays a significant role in this process, and it shows a large area of the regional mass loss in southeastern Greenland.

4.4. Influence of Land Surface Temperature

Snow et al. (2021) [75] examined the influence of the sea surface temperature (SST) and ocean currents near the Helheim glacier and showed that fluctuations in subsurface Atlantic water temperatures coincided with some of the variability in discharge rates previously found at Helheim, although this relationship is not straightforward. For example, the SST was increasing in 2000 when Helheim experienced discharge rates. However, the SST was relatively warm in 2010, although Helheim did not experience any substantial increase in ice discharge. These findings suggest that there is no straight relationship between the SST and discharge rate of the glacier. In addition to the SST and ocean currents, other environmental factors such as land surface temperature, glacier configuration, and mélange rigidity may influence glacier discharge rates.

As analyzed above, the InSAR time series suggest faster melting in 2019, in comparison to periods before and after. Reports from the World Meteorological Organization (WMO) indicate that 2016–2020 were the five warmest years on record [76], with 2019 being warmer than the record-level temperatures seen in 2016 and 2018, making 2019 a fast melt year for the ice sheets.

We compared the average of the local LST in our study area for each year by using MOD11C2 data from MODIS (Figure 12). For the 2020 local LST, the average of the first 9 months, for which we had InSAR data, is considered. As seen in Figure 12, the regional analysis from the LST shows that the LST rate increased from 2016 to 2019 but decreased in 2019–2020. However, as illustrated in Figure 6a, displacement rates decreased from 2016 to 2017 and 2018 and increased again in 2019, before decreasing in 2020. Thus, there is no direct link between elastic rebound of the crust and land surface temperature, and as stressed above, other parameters may also influence ice melting and the resulting unloading at the edges of the glaciers.

5. Conclusions

InSAR time series interferograms are useful for assessing the elastic surface deformation on the edges of large glaciers subject to ice melting, providing insight into the rheological characteristics of the region. Based on Sentinel-1 InSAR measurements covering summer seasons between 2016 and 2020 over the Helheim glacier, we could infer the uplift on the northern flank of the Helheim margin up to 58 mm in the LOS direction, attributed to the elastic response of solid earth to ice melting. InSAR time series analysis suggests faster melting in 2019 than in 2016, 2017, 2018, and 2020. Modeling the thinning rate of ice thickness and the mass loss rate due to unloading in this region results in an average thinning rate of 5.07 m/year at the edge of ice sheet and a mass loss of 8.33 Gt/year on the northern flank of the Helheim margin. Although the distribution of InSAR measurement points in this region is rather sparse, this study shows the ability of the Sentinel-1 SAR system to provide a reliable assessment of the surface deformation on the edges of a large glacier, which makes it promising for applying to other areas to evaluate the elastic response of solid earth to glacier and ice sheet melting.