Greenland Mass Trends From Airborne and Satellite Altimetry During 2011-2020.

We use satellite and airborne altimetry to estimate annual mass changes of the Greenland Ice Sheet. We estimate ice loss corresponding to a sea-level rise of 6.9 ± 0.4 mm from April 2011 to April 2020, with a highest annual ice loss rate of 1.4 mm/yr sea-level equivalent from April 2019 to April 2020. On a regional scale, our annual mass loss timeseries reveals 10-15 m/yr dynamic thickening at the terminus of Jakobshavn Isbræ from April 2016 to April 2018, followed by a return to dynamic thinning. We observe contrasting patterns of mass loss acceleration in different basins across the ice sheet and suggest that these spatiotemporal trends could be useful for calibrating and validating prognostic ice sheet models. In addition to resolving the spatial and temporal fingerprint of Greenland's recent ice loss, these mass loss grids are key for partitioning contemporary elastic vertical land motion from longer-term glacial isostatic adjustment (GIA) trends at GPS stations around the ice sheet. Our ice-loss product results in a significantly different GIA interpretation from a previous ice-loss product.

Introduction

In recent decades, several methods have been used to estimate the Greenland Ice Sheet (GIS) mass change (Shepherd and Ivins, 2020). The Gravity Recovery and Climate Experiment (GRACE) and GRACE Follow‐On (FO) missions over Greenland show an ice loss of more than 4,550 ± 784 Gt (12.6 ± 2.2 mm sea level equivalent (SLE)) during 2002–2019 (Velicogna et al., 2020), with a record high mass loss rate in 2019 of 532 ± 58 Gt/yr (1.5 mm SLE; Sasgen et al., 2020). The input‐output method (IOM), which provides the longest continuous time series of mass change discussed here, suggests accelerated mass loss of the GIS over the last three decades (Mouginot et al., 2019). The IOM typically uses surface mass balance (SMB) estimates from a regional climate model, combined with grounding line discharge of glacier speed and ice thickness (Rignot et al., 2008). This allows the method to go further back in time and provide continuous annual ice sheet‐wide mass loss estimates from 1972 onward (Mouginot et al., 2019). The IOM reveals that ice discharge was relatively constant from 1972 to 2002, and increased thereafter as a result of the speed‐up of outlet glaciers. The IOM suggests that the GIS contributed a total of 13.7 ± 1.1 mm to sea level rise from 1972 to 2018 (King et al., 2020; Mouginot et al., 2019).

Satellite gravimetry provides high temporal resolution (monthly) mass loss estimates (Sasgen et al., 2020; Velicogna et al., 2020) of the GIS including peripheral glaciers, but the spatial resolution of ∼300 km does not allow changes of individual glaciers to be studied, or to separate the ice sheet from peripheral glaciers. The IOM provides mean glacier catchment mass changes, but cannot reveal the evolving spatial pattern of mass loss within a catchment (Mouginot et al., 2019; Rignot et al., 2008). Satellite and airborne altimetry, however, can provide direct measurements of changing ice‐sheet surface elevation at the spatial resolution of a few kilometers. Previous altimetry studies used observations from one satellite mission or fusion of multi‐sensors to estimate trends over various time intervals, typically between 5 and 10 years (Bamber & Dawson, 2020; Csatho et al., 2014; Gardner et al., 2013; Helm et al., 2014; Hurkmans et al., 2014; Shepherd and Ivins, 2020; Smith and Fricker, 2020; Smith et al., 2019; Sørensen et al., 2018). A small number of studies provide annual elevation change estimates, for example, Simonsen et al. (2021) and Sørensen et al. (2018). Several methods have been used to estimate ice‐sheet surface elevation changes, for example, orbit crossing points, along repeated ground tracks, or using plane‐fit solutions (Hurkmans et al., 2014; Khan et al., 2014; Moholdt et al., 2010; Sørensen et al., 2011).

Here, we estimate annual mass changes of the GIS (excluding peripheral glaciers) on a high resolution (1 × 1 km) grid. We make these annual grids available as an open‐access product. These mass loss grids are a valuable input to several branches of geodynamic and cryosphere research. For example, stations within the Greenland GPS Network (GNET) are uplifting in response to spatially and temporally variable combinations of long‐term Glacial Isostatic Adjustment (GIA) and present‐day elastic vertical land motion (VLM) of the bedrock (Bevis et al., 2012; Khan et al., 2016; Ludwigsen et al., 2021). Knowledge of the magnitude and spatial distribution of Greenland's GIA fingerprint has broad interdisciplinary value, such as to calculate the non‐uniform global sea‐level rise. To isolate the GIA signal, high‐resolution mass‐loss grids are needed to apply an elastic correction to partition elastic VLM from the elastic VLM + GIA measured at each GPS station. We estimate an elastic VLM correction for all GNET stations using our high‐resolution grid, and compare this solution with other available mass loss grids. High‐resolution mass loss grids are also needed for the IOM approach, which requires knowledge about ice thickness changes between the fluxgate and grounding line (Colgan et al., 2019; Hansen et al., 2021). Finally, the high spatial and temporal resolution grid that we present allows us to examine the complex evolving regional variations in mass loss, driven by both SMB and ice dynamics, that comprise Greenland's ice‐loss fingerprint.

Data and Methods

Radar and Laser Altimetry Data

We have used radar altimetry data from ESA's Earth Explorer CryoSat‐2 mission (Wingham et al., 2006) to estimate annual mass changes of the GIS from April 2011 to April 2020. To better resolve mass loss along the ice‐sheet margin, we use laser altimetry observations from NASA's Operation IceBridge Airborne Topographic Mapper (ATM) flights from April 2011 to April 2019 (Studinger, 2014). NASA ended its Operation IceBridge over Greenland in spring 2019. We fill the subsequent gap in laser altimetry data with Ice, Cloud, and land Elevation Satellite‐2 (ICESat‐2) data from April 2019 to April 2020 (Smith et al., 2021).

CryoSat‐2 Data Processing

To estimate elevation changes over ice we re‐processed CryoSat‐2 level_1B Baseline_D data. In our processing, we followed the approach of Helm et al. (2014) by applying a Threshold First Maximum Retracker Algorithm (TFMRA) with a low threshold (0.25) for Low Resolution Mode (LRM) and the relocation slope correction to estimate the Point‐Of‐Closest‐Approach (POCA). The POCA of the Synthetic Aperture Radar Interferometric (SARIn) data was estimated by using a slightly filtered phase and a coherency threshold of 0.75. Phase ambiguities were minimized by using the ArcticDEM as a reference, similar to Nilsson et al. (2016). Nilsson et al. (2016) used a similar low threshold retracker for LRM and found remarkable evidence of low sensitivity to seasonal changes in the scattering regime and hence reduced seasonal penetration affecting trends of elevation change of the upper snowpack. Consequently, the 2012 melt event does not appear as an erroneous surface increase in our reprocessed Level_2 product (see Figures 1, 2, 3, 4).

Figure 1

Effect of backscatter correction for Threshold First Maximum Retracker Algorithm (TFMRA) derived elevations shown for a 500 km2 area located in Northern Greenland: (a) TFMRA elevation anomaly (black) and backscatter corrected elevation anomaly (blue) (Note: for better visibility, the blue graph was shifted by 1m); (b) backscatter anomaly; (c) Linear correlation between elevation and backscatter anomaly. The regression coefficients were used to apply the backscatter correction.

Figure 2

Effect of backscatter correction for ESA‐ICE2 derived elevations shown for a 500 km2 area located in Northern Greenland: (a) ESA‐ICE2 elevation anomaly (black) and backscatter corrected elevation anomaly (blue) (Note: for better visibility the blue graph was shifted by 1m); (b) backscatter anomaly; (c) Linear correlation between elevation and backscatter anomaly. The regression coefficients were used to apply the backscatter correction.

Figure 3

Effect of leading edge correction for Threshold First Maximum Retracker Algorithm (TFMRA) derived elevations shown for a 500 km2 area located in Northern Greenland: (a) TFMRA elevation anomaly (black) and leading edge corrected elevation anomaly (blue) (Note: for better visibility the blue graph was shifted by 1m); (b) leading edge width; (c) Linear correlation between elevation and leading edge width. The regression coefficients were used to apply the leading edge width correction.

Figure 4

Effect of leading edge correction for ESA‐ICE2 derived elevations shown for a 500 km2 area located in Northern Greenland: (a) ESA‐ICE2 elevation anomaly (black) and leading edge corrected elevation anomaly (blue) (Note: for better visibility the blue graph was shifted by −1m); (b) leading edge width; (c) Linear correlation between elevation and leading edge width. The regression coefficients were used to apply the leading edge width correction.

We further evaluate whether there is a seasonal elevation anomaly and the effect of proposed corrections to reduce the penetration effect (Simonsen & Sørensen, 2017) for a 500 km2 area in North Greenland (76°N, 44°W). This area was chosen because it was strongly affected by the 2012 melt event. Figures 1 and 2 show the elevation anomaly of the TFMRA (Helm et al., 2014) and ESA‐ICE2 retracker (Brockley, 2019) and the effect of the applied backscatter correction, respectively. TFMRA shows only a small elevation change caused by the change in the scattering regime (from volume to surface scattering due to the melt event), while the ESA‐ICE2 elevation increased roughly by 1.5 m between June–August 2012. The backscatter correction doesn't affect the long‐term trend for TFMRA and ESA‐ICE2. However seasonal amplitudes are reduced specifically for ESA‐ICE2. Consequently, the backscatter correction is not able to correct for scattering regime shifts. In addition, we evaluated the leading edge correction. Figures 3 and 4 show, that for ESA‐ICE2 the elevation increase can be reduced to 0.7 m, whereas for TFMRA there is nearly no effect visible. Again, the seasonal amplitudes are reduced, specifically for ESA‐ICE2, but the long‐term trend for that specific position is not changing. In summary, it can be seen that TFMRA shows reduced seasonal elevation amplitudes and is less sensitive to time‐varying backscatter than ESA‐ICE2. Furthermore, a correlation with backscatter or leading‐edge width only slightly changes the estimated trend. In addition, linear correlations of 0.6 and less show that the dependencies are not striking.

Based on these sensitivity analyses, we conclude that an additional correlation with backscatter or leading‐edge width, as presented by in Simonsen & Sørensen, 2017, or the introduction of a step function as suggested by McMillan et al. (2016), is not required for the TFMRA retracked elevations. The same, or even better, reduction of sensitivity due to time‐varying penetration/scattering is clearly achieved using low thresholds for LRM data. SARIn processed data is less affected by time‐varying penetration as shown by Aublanc et al. (2018).

Annual Elevation Changes From CryoSat‐2

To estimate elevation changes over the ice surface, we use a regular grid with a resolution of 1 × 1 km that covers the entire GIS. We denote the center of each grid point with C (x , y ). For each grid point, we select CryoSat‐2 data with coordinates P (x , y ), with a maximum distance of 1,000 m from C. The CryoSat‐2 data points with coordinates P (x , y ), have elevation h measured at time t . The index i denote the ith data point.

Here, we use all available CryoSat‐2 data measured between July 2010 and July 2020 to create surface elevation time series at each grid point C. To describe surface changes over 10 years, and considering large surface lowering due to melting in 2011, 2012, and 2019, and cold periods with less elevation change in between, we fit seventh‐order polynomial to describe changes in elevation and a third‐order polynomial equation to describe the shape of the surface. In addition, we fit a seasonal term to account for the annual surface changes. For each grid point with center (x , y ), we find the nearest data within 1,000 m (x i, y i, h i, t i) and fit a seventh‐order polynomial , a third‐order surface topography H and an annual term :

The seventh‐order polynomial is, where t is the time when the ith measurement was observed. For simplicity, we used 1 January 2010 as reference time t = 0. a to a are polynomial parameters. The parameter a represents surface elevation at time 1 January 2010. We assume the shape of the surface remains constant (Schenk et al., 2014) during the time interval considered in this study. While many studies use a planar surface (Hurkmans et al., 2014; Sørensen et al., 2011), here, we fit a third‐order polynomial to describe a surface topography of the area of 1 × 1 km (Csatho et al., 2014), where a to a are parameters that describe the slope and concavity/convexity of the surface. x and y are coordinates of the CryoSat‐2 data point, but in a system with x and y as the center. that is, x = x − x and y = y − y . The annual term is given by, where a denotes the amplitude, a is the phase, and ω is the angular frequency of the annual signal. For each grid point, we create a time series and use least squares adjustment to simultaneously estimate parameters a to a . Our procedure for deriving ice‐surface elevation changes is very similar to the method used by, for example, Schenk et al. (2014).

Figure 5 shows examples of time series of surface elevation corrected for third‐order surface topography and the annual term. The red curve shows the best‐fitting seventh‐order polynomial.

Figure 5

(a) Mean surface elevation change during 2011–2020. Thickness change time series derived from CryoSat‐2 data for a single point on (b) Jakobshavn Isbrae, (c) Kangerlussuaq Glacier (KG) and (d) Helheim Glacier (HG). The location of points is shown in (a) with symbols denoting Jakobshavn Isbræ (white square), KG (blue star), and HG (black circle). The solid red curves in (b) to (d) show the best‐fitting seventh‐order polynomial.

We use the parameters a to a for each grid point to estimate elevation changes over the entire GIS over consecutive 1‐year periods, starting with 15 April 2011–2015 April 2012, then from 15 April 2012–2015 April 2013, etc., until 15 April 2020.

Elevation Changes From NASA's Operation IceBridge ATM Flights

We estimate elevation changes using NASA's ATM surveys in Greenland from spring 2011 to spring 2019. The ATM flights are mainly concentrated along the margin of the GIS. To estimate elevation changes, we take the height difference between overlapping points from two different campaigns, that is, we take the height differences between the 2011 survey and 2012 survey, and between the 2012 survey and 2013 survey, etc. However, it should be noted that no survey was conducted in spring 2020. Annual elevation change rates from overlapping ATM surveys from spring 2011 to spring 2019 are shown in Figure 6.

Figure 6

Annual elevation change rates in m/yr from Airborne Topographic Mapper surveys from spring 2011 to spring 2019.

Elevation Changes From ICESat‐2 Data

Since NASA did not conduct an ATM survey in Greenland in spring 2020, we use ICESat‐2 observations from spring 2019 to spring 2020 and estimate elevation changes for overlapping ICESat‐2 points (Smith et al., 2019; Smith and Fricker, 2020). We use ICESat‐2 Algorithm Theoretical Basis Document for Land Ice Height (ATL06) Release 004 retrieved from https://nsidc.org/data/atl06 (Smith et al., 2021). Annual elevation change rates from overlapping ICESat‐2 points during spring 2019 to spring 2020 are shown in Figure 7a.

Figure 7

(a) Annual elevation change rates in m/yr from overlapping Ice, Cloud, and land Elevation Satellite‐2 (ICESat‐2) tracks during spring 2019 and spring 2020. (b) Annual elevation change rates from CryoSat‐2 during spring 2019 and spring 2020 using the method described in Section 2.1.2. (c) Difference between ICESat‐2 and CryoSat‐2 derived thinning rates.

The level of agreement between radar and laser altimetry data is shown in Figure 7. In general, CryoSat‐2 and ICESat‐2 capture large thinning rates in the same regions. We estimate a mean difference of 0.03 m/yr and standard deviation of 0.08 m/yr between ICESat‐2 and CryoSat‐2 thinning rates.

Elevation Changes Over 1‐Year and Their Uncertainty

It is important to note that we do not merge CryoSat‐2, ICESat‐2, and NASA's ATM surveys when we create elevation time series or elevation change rates (Figures 5, 6, 7) because the CryoSat‐2 radar signal can penetrate through the snowpack of the ice surface. For example, we do not use ICESat‐2 data from 2020 and ATM surveys from 2019 to estimate elevation change rates between these datasets. Instead, we estimate annual elevation change rates from April to April for each sensor's data set independently, and only afterward merge the estimated change rates from each data set when creating a multi‐sensor annual grid.

The estimated annual elevation change rates from CryoSat‐2, ICESat‐2, and NASA's ATM were used to interpolate elevation change rates onto a regular grid of 1 × 1 km. The interpolation was performed using the ordinary kriging method (Hurkmans et al., 2014; Nielsen et al., 2013). We first used the observed annual elevation change rates to estimate an empirical semi‐variogram. Next, we fit an exponential model variogram for each 1‐year interval (with a range of 65 km) to the empirical semi‐variogram to take the spatial correlation of elevation change rates into account in the error budget. For each grid point, we interpolate (using kriging) elevation change rate dhi,krig and the associated error σ .

Correction for Elastic VLM, Glacial Isostatic Adjustment, and Firn Compaction

The observed ice surface elevation changes were corrected for bedrock movement caused by elastic VLM from present‐day mass changes and long‐term past ice changes GIA. To correct for GIA, we use the GNET‐GIA empirical model of Khan et al. (2016). For each grid point on a 1 × 1 km grid, we adopt the GIA rate dhGIA and the associated uncertainty σGIA directly from Khan et al. (2016) (see their Figure 3). We correct for elastic VLM of the bedrock by convolving mass loss estimates from CryoSat‐2, ATM, and ICESat‐2 with the Green's functions derived by Wang et al. (2012) for elastic Earth model iasp91 with refined crustal structure from Crust 2.0 (Bassin et al., 2000). For each grid point, we estimate the elastic VLM rate dhelas and the associated uncertainty σelas. To estimate σelas, we convolve uncertainties of mass loss with the Green's functions. We converted the volume loss rate into a mass loss rate, taking firn compaction into account, as described by Ligtenberg et al. (2018). The firn compaction model is forced with surface mass fluxes and temperature from Regional Atmospheric Climate Model (RACMO2.3p2) as described in Ligtenberg et al. (2018). We estimate firn compaction uncertainties as described in Kuipers Munneke et al. (2015) (see their Equations 8 and 9) using input fields from RACMO2.3p2 until April 2020. For each grid point, we estimate the firn compaction rate dhi,firn and the associated uncertainty σi,firn. The total elevation change rate is for each grid point i is,

Assuming the uncertainty in each of these terms is independent, we summarize the total associated uncertainty by summing each uncertainty term in quadrature

Table 1 shows Ice‐sheet wide annual corrections due to elastic VLM of the bedrock, GIA, and firn compaction.

Table 1

Annual Rates of Corrections Due To Elastic Vertical Land Motion (VLM) of the Bedrock, Glacial Isostatic Adjustment, and Firn Compaction

Time interval (year)	Elastic VLM (km3/yr)	Firn compaction (km3/yr)	GIA (km3/yr)
2011–2012	9.4 ± 1.1	−57.0 ± 9.2	1.1 ± 0.1
2012–2013	12.2 ± 1.4	−46.7 ± 9.2	1.1 ± 0.1
2013–2014	4.4 ± 1.5	−69.9 ± 9.2	1.1 ± 0.1
2014–2015	6.4 ± 1.5	−58.3 ± 9.2	1.1 ± 0.1
2015–2016	4.2 ± 1.4	−71.0 ± 9.2	1.1 ± 0.1
2016–2017	5.4 ± 1.5	−40.7 ± 9.2	1.1 ± 0.1
2017–2018	4.0 ± 1.6	−27.7 ± 9.2	1.1 ± 0.1
2018–2019	5.6 ± 1.1	−15.8 ± 9.2	1.1 ± 0.1
2019–2020	14.3 ± 1.1	−8.4 ± 9.2	1.1 ± 0.1

Ice Mass Change

To estimate ice ice‐sheet wide mass change rates, we used the ice mask of Pedersen et al. (2013). The mass change dM of the GIS or a drainage area was estimated as the sum of elevation changes dhi over the considered area, where n is an index of grid nodes within the given area and ρice is the density of ice (917 kg/m3). The associated uncertainty was similarly estimated as,

To estimate SMB induced elevation change and mass change, we used a SMB product from the Regional Atmospheric Climate Model (RACMO2.3p2) at 5.5 km horizontal resolution statistically downscaled to 1 km resolution (Noël et al., 2018, 2019). We use RACMO2.3p2 extended to April 2020. SMB uncertainty has been previously estimated in this product as the bias between modeled and observed SMB. We adopt this previous uncertainty estimate of ±17 mm w. e./yr in the accumulation zone, and ±70 mm w. e./yr in the ablation zone (Noël et al., 2019).

To estimate elevation changes due to changes in ice dynamics for each grid point, dhi,dyn, we subtracted SMB derived elevation changes dhi,SMB from RACMO2.3p2 (Noël et al., 2019) from airborne and satellite altimetry dhi.

The associated uncertainty of changes in ice dynamics is

Results

The observed elevation changes from April 2011 to April 2020 are displayed in Figure 8 (top row of panels) and Table 2, partitioned into SMB anomalies (middle row) and annual elevation change rates due to ice dynamics (bottom row). From April 2011 to April 2020, the GIS lost 2459 ±154 Gt of ice (see Figure 9 and Table 2). In line with recent GRACE studies (Sasgen et al., 2020; Velicogna et al., 2020), we observed a peak mass loss of 498 ± 45 Gt froms April 2019 to April 2020, which exceeds the previous peak mass loss anomaly of 453 ± 45 Gt observed from April 2012 to April 2013 (see Figure 9).

Figure 8

(Top row) Annual (April to April) elevation change rates of the Greenland Ice Sheet from April 2011 to April 2020 from CryoSat‐2, ICESat‐2, and NASA's Airborne Topographic Mapper flights (Middle row) annual elevation change rates due to surface mass balance anomalies (Bottom row) annual elevation change rates due to ice dynamics.

Table 2

Total Greenland Ice Sheet (GIS) Mass Change (Negative Values Indicate Loss) in Gigatons (Gt = 1012 kg) for the Period From April 2011 to April 2020 for the GIS Partitioned Between Surface Mass Balance Processes and Ice Dynamics

	Total GIS
Time	Total	SMB	Dynamic
Year	Gt/yr	Gt/yr	Gt/yr
2011.3–2012.3	−352.7 ± 42.0	−257.8 ± 52.2	−94.9 ± 67.0
2012.3–2013.3	−452.8 ± 44.7	−298.6 ± 52.3	−154.2 ± 68.8
2013.3–2014.3	−163.1 ± 49.9	−22.2 ± 35.2	−140.9 ± 61.1
2014.3–2015.3	−246.2 ± 51.9	−116.2 ± 52.2	−130.0 ± 73.6
2015.3–2016.3	−165.3 ± 49.1	−126.3 ± 52.2	−39.0 ± 71.7
2016.3–2017.3	−222.5 ± 48.6	−183.5 ± 52.2	−39.0 ± 71.3
2017.3–2018.3	−167.0 ± 50.9	−15.3 ± 37.4	−151.7 ± 63.2
2018.3–2019.3	−191.2 ± 48.2	2.0 ± 35.3	−189.2 ± 59.7
2019.3–2020.3	−498.0 ± 45.0	−380.5 ± 52.4	−117.5 ± 69.1
2011.3–2020.3	−273.2 ± 17.1	−155.3 ± 33.6	−117.4 ± 47.5

Figure 9

Left axis: Cumulative anomalies in total mass (gray) in gigatons (Gt = 1012 kg) for the period from April 2011 to April 2020 for the Greenland Ice Sheet (GIS). Right axis: Dashed blue curve denotes annual (April to April) mass change rate of the total mass of the GIS in Gt/yr. The thickness of the curves denotes the error estimates.

Regional Mass Loss

To investigate regional scale ice mass variability, we divided the GIS into eight major drainage basins denoted as D1‐D8 (Zwally et al., 2012) and provide basin‐wide cumulative total mass change, along with SMB and dynamic ice losses, in each basin (see Figure 10 and Table 3). D1 consists of north Greenland characterized by many marine‐terminating outlet glaciers, including Petermann Glacier, with the second largest floating shelf of Greenland (Rignot et al., 2001) and Humboldt Glacier. In this sector, we observe a total mass loss of 229 ± 17 Gt from April 2011 to April 2020 with large inter‐annual variability (Figure 10b). The northeast sector, D2, in total lost 97 ± 21 Gt during 2011–2020 (Figure 10c). This sector is home to several major outlet glaciers for example, Nioghalvfjerdsfjorden Glacier, Zachariae Isstrøm, and Storstrømmen Glacier that define the northeast Greenland ice stream, which extends >600 km inland (Hvidberg et al., 2020). Continued dynamic mass loss of 5–10 Gt/yr dominates this sector (Mouginot et al., 2019).

Figure 10

Cumulative anomalies in total mass (gray) in gigatons partitioned between surface mass balance processes (light red) and ice dynamics (purple) for (a) Drainage D8, (b) Drainage D1, (c) Drainage D2, (d) Drainage D7, (f) Drainage D3, (g) Drainage D6, (h) Drainage D5, (i) Drainage D4. The thickness of the curves denotes the error estimates. (e) Map of average elevation changes in m/yr from April 2011 to April 2020. Location names are shown for Jakobshavn Isbræ (JI), Helheim Glacier, Kangerlussuaq Glacier, Nioghalvfjerdsfjorden Glacier, the Zachariae Isstrøm, Storstrømmen Glacier, Petermann Glacier, Humboldt Glacier, and northeast Greenland ice stream.

Table 3

Annual Drainage Wide Mass Change Rate (in Gigatons/Yr) for the Period From April 2011 to April 2020 Partitioned Between Surface Mass Balance (SMB) Processes and Ice Dynamics

	Drainage D1			Drainage D2
Time	Total	SMB	Dynamic	Total	SMB	Dynamic
Year	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr
2011.3–2012.3	−32.1 ± 4.7	−29.7 ± 7.9	−2.4 ± 9.2	−22.2 ± 5.8	−0.1 ± 9.6	−22.1 ± 11.2
2012.3–2013.3	−33.7 ± 4.8	−39.8 ± 7.9	6.1 ± 9.3	−7.3 ± 6.5	−11.9 ± 9.6	4.6 ± 11.6
2013.3–2014.3	−17.8 ± 5.4	9.8 ± 5.3	−27.6 ± 7.6	5.0 ± 6.9	7.7 ± 6.4	−2.7 ± 9.4
2014.3–2015.3	−29.9 ± 5.9	−21.3 ± 7.9	−8.6 ± 9.8	−7.1 ± 7.6	2.6 ± 9.6	−9.7 ± 12.2
2015.3–2016.3	−29.3 ± 5.7	−38.3 ± 7.9	9.0 ± 9.7	−10.2 ± 6.9	5.3 ± 9.6	−15.5 ± 11.8
2016.3–2017.3	−13.9 ± 5.5	−21.7 ± 7.9	7.8 ± 9.6	−10.0 ± 7.2	−13.2 ± 9.6	3.2 ± 12.0
2017.3–2018.3	−5.4 ± 5.5	−3.3 ± 5.8	−2.1 ± 8.0	−6.0 ± 7.4	−9.5 ± 7.0	3.5 ± 10.2
2018.3–2019.3	−18.4 ± 5.1	−5.2 ± 5.3	−13.2 ± 7.4	−4.5 ± 6.7	25.5 ± 6.4	−30.0 ± 9.2
2019.3–2020.3	−48.5 ± 4.7	−54.9 ± 7.9	6.4 ± 9.2	−34.9 ± 6.2	−31.5 ± 9.6	−3.4 ± 11.4
	Drainage D3			Drainage D4
Time	Total	SMB	Dynamic	Total	SMB	Dynamic
Year	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr
2011.3–2012.3	−51.4 ± 8.7	1.1 ± 7.8	−50.3 ± 11.7	−49.8 ± 4.8	−37.5 ± 4.5	−12.3 ± 6.5
2012.3–2013.3	−94.9 ± 8.3	−65.0 ± 7.8	−29.9 ± 11.4	−60.5 ± 5.0	−38.9 ± 4.5	−21.6 ± 6.7
2013.3–2014.3	−14.5 ± 9.6	−11.4 ± 5.3	−3.1 ± 11.0	−14.7 ± 5.5	1.3 ± 3.0	−16.0 ± 6.2
2014.3–2015.3	27.2 ± 10.6	13.5 ± 7.8	13.7 ± 13.1	−41.7 ± 5.7	3.3 ± 4.5	−45.0 ± 7.2
2015.3–2016.3	28.7 ± 9.4	18.3 ± 7.8	10.4 ± 12.2	−2.5 ± 5.3	−9.3 ± 4.5	6.8 ± 7.0
2016.3–2017.3	−1.3 ± 9.4	−9.2 ± 7.8	7.9 ± 12.2	−33.7 ± 5.3	−2.3 ± 4.5	−31.4 ± 7.0
2017.3–2018.3	−31.8 ± 9.7	18.0 ± 5.8	−49.8 ± 11.3	−38.5 ± 5.6	−13.3 ± 3.1	−25.2 ± 6.4
2018.3–2019.3	−37.1 ± 9.1	28.2 ± 5.3	−65.1 ± 10.6	−3.9 ± 5.7	−0.5 ± 3.0	−3.4 ± 6.4
2019.3–2020.3	−105.8 ± 8.8	−53.5 ± 7.8	−52.3 ± 11.8	−26.3 ± 5.4	4.9 ± 4.4	−31.2 ± 7.0
	Drainage D5			Drainage D6
Time	Total	SMB	Dynamic	Total	SMB	Dynamic
Year	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr
2011.3–2012.3	−25.0 ± 2.1	−28.9 ± 1.6	3.9 ± 2.6	−77.9 ± 5.5	−76.6 ± 5.7	−1.3 ± 7.9
2012.3–2013.3	−26.9 ± 2.0	−24.5 ± 1.7	−2.4 ± 2.6	−96.7 ± 5.7	−75.4 ± 5.7	−21.3 ± 8.0
2013.3–2014.3	−9.5 ± 2.7	−1.6 ± 1.1	−7.9 ± 2.9	−15.3 ± 6.3	−9.7 ± 3.9	−5.6 ± 7.4
2014.3–2015.3	−11.6 ± 2.7	−17.8 ± 1.6	6.2 ± 3.2	−32.5 ± 5.3	−48.8 ± 5.7	16.3 ± 7.8
2015.3–2016.3	−23.2 ± 2.4	−17.8 ± 1.6	−5.4 ± 2.9	−19.6 ± 6.0	−7.1 ± 5.7	−12.5 ± 8.3
2016.3–2017.3	−24.8 ± 2.1	−26.9 ± 1.7	−2.1 ± 2.6	−68.4 ± 5.2	−75.3 ± 5.7	6.9 ± 7.7
2017.3–2018.3	−17.4 ± 2.1	−15.4 ± 1.2	−2.0 ± 2.4	−10.6 ± 6.3	−3.6 ± 3.9	−7.0 ± 7.4
2018.3–2019.3	−1.90 ± 2.5	−7.0 ± 1.2	5.1 ± 2.7	−29.6 ± 5.9	−20.6 ± 3.9	−9.0 ± 7.1
2019.3–2020.3	−12.8 ± 2.0	−13.1 ± 1.7	−0.3 ± 2.7	−88.3 ± 5.6	−113.9 ± 5.8	25.6 ± 8.0
	Drainage D7			Drainage D8
Time	Total	SMB	Dynamic	Total	SMB	Dynamic
Year	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr	Gt/yr
2011.3–2012.3	−43.4 ± 5.5	−42.4 ± 6.8	−1.0 ± 8.7	−50.9 ± 5.8	−43.6 ± 8.4	−7.3 ± 10.2
2012.3–2013.3	−69.1 ± 6.1	−12.7 ± 6.8	−56. ± 9.2	−63.5 ± 6.2	−30.0 ± 8.4	−33.5 ± 10.5
2013.3–2014.3	−39.9 ± 6.4	−12.6 ± 4.5	−27.3 ± 7.8	−56.3 ± 7.1	−5.7 ± 5.6	−50.6 ± 9.1
2014.3–2015.3	−47.6 ± 6.8	−21.4 ± 6.8	−26.2 ± 9.6	−103.0 ± 7.4	−26.2 ± 8.4	−76.8 ± 11.2
2015.3–2016.3	−32.3 ± 6.5	−20.3 ± 6.8	−12.0 ± 9.4	−76.9 ± 6.9	−56.6 ± 8.3	−20.3 ± 10.8
2016.3–2017.3	−21.8 ± 7.0	−18.8 ± 6.8	−3.0 ± 9.8	−48.6 ± 7.0	−15.8 ± 8.4	−32.8 ± 10.9
2017.3–2018.3	−22.6 ± 7.3	5.5 ± 4.7	−28.1 ± 8.6	−34.5 ± 7.0	6.5 ± 5.9	−41.0 ± 9.2
2018.3–2019.3	−41.1 ± 6.1	−16.4 ± 4.5	−24.7 ± 7.6	−54.5 ± 7.1	−2.8 ± 5.6	−51.7 ± 9.0
2019.3–2020.3	−93.9 ± 6.1	−65.3 ± 6.8	−28.6 ± 9.1	−87.3 ± 6.0	−52.6 ± 8.4	−34.7 ± 10.3

Central East Greenland (D3) consists of numerous outlet glaciers, with Kangerlussuaq Glacier being the largest contributor to dynamic ice loss over the past decade (Mouginot et al., 2019). In this sector, we observe a modest mass gain from April 2014 to April 2017. However, during 2011–2020, this sector in total lost 281 ± 29 Gt of ice (Figure 10f). Like D3, the southeast basin (D4, Figure 10i) is dominated by numerous outlet glaciers that show accelerated flow speed, and have thinned and retreated over the past three decades (King et al., 2020; Mouginot et al., 2019).

In the sectors that make up south (D5) and southwest (D6) Greenland, the ice‐sheet margin is located between 50 and 100 km inland from the coast, and only a few minor outlet glaciers can contribute to dynamic ice loss (King et al., 2020; Mouginot et al., 2019). Mass changes are dominated by SMB fluctuations. The increased loss in 2011–2012, and 2020 (see Figure 10g) coincide with extensive summer melting (Bevis et al., 2019; Sasgen et al., 2020; Velicogna et al., 2020).

The D7 basin contains Jakobshavn Isbræ (JI, also known as Sermeq Kujalleq), which has been losing mass since the collapse of its ice tongue in the early 2000s (Joughin et al., 2008). This sector lost 412 ± 20 Gt of ice during 2011–2020 (see Figure 10d) with large negative mass loss anomalies in 2012 and 2019 with respect to 2011–2020 mean. The northwest sector (D8) is known for its many marine terminating glaciers. The area is characterized by little unglacierized land, and almost no fjord systems, with the ice sheet instead terminating directly in the open ocean. This drainage basin is dominated by significant dynamic ice loss (Mouginot et al., 2019). This basin is the largest contributor to mass loss among basins D1–D8, with a total mass loss of 576 ± 21 Gt of ice during 2011–2020 (Table 3).

Complex Evolution of Mass Loss at Greenland's Largest Outlet Glacier: Jakobshavn Isbræ

Dynamic changes cause a complex pattern of thinning and thickening in space and time (Moon et al., 2012) and must be resolved to improve projections of mass loss. Jakobshavn Isbræ has been retreating, thinning, and speeding up since 2000, with record flow speeds measured in 2012 (Joughin et al., 2020). Several studies have shown that subsurface ocean warming and cooling is the dominant factor controlling JI ice dynamics (Holland et al., 2008; Khazendar et al., 2019; Wood and Rignot, 2021).

Elevation changes from airborne and satellite altimetry show annual thinning at rates >10 m/yr (near the glacier front) between April 2011 and April 2016 (Figure 11). This thinning is dominated by ice dynamics (Khazendar et al., 2019). From April 2016 to April 2019, we observed episodic dynamic thickening of about 5–15 m/yr, from near the glacier terminus to about 50 km upstream. Despite thickening near the glacier terminus in 2016–2019, we estimated a JI basin‐wide net dynamic mass loss (Figure 11k). This is due to persistent thinning of the upstream part of the glacier. Interestingly, our results showed that the thickening was short lived and subsequently replaced by thinning of about 5 m/yr near the terminus after April 2019 (Figure 11j).

Figure 11

Elevation change rates of Jakobshavn Isbræ in m/yr during (a) April 2011 ‐ April 2012, (b) April 2012 ‐ April 2013, (c) April 2013 ‐ April 2014, (d) April 2014 ‐ April 2015, (e) April 2015 ‐ April 2016, (f) April 2016 ‐ April 2017, (g) April 2017 ‐ April 2018, (h) April 2018 ‐ April 2019, (i) April 2019 ‐ April 2020. The black dashed line in (a) denotes the main flowline. Landsat image of Jakobshavn Isbræ from 3 July 2020 is used as background (j) Elevation change rates from April to April each year from 2011 to 2020 along the main flowline shown in (a) (k) Left axis: Cumulative anomalies in total mass (gray), dynamic loss (blue), and surface mass balance (pink) in gigatons for the time period from April 2011 to April 2020 for the Jakobshavn Isbræ drainage basin. Right axis: Cumulative mass change converted to equivalent sea‐level rise. JI drainage basin is marked as red area in the Greenland map.

Comparison With Previous Studies

We compare our annual mass change rates with previous studies. Figure 12 shows ice‐sheet wide mass loss rates from this study (black curve), compared with IMBIE (Shepherd and Ivins, 2020), Colgan et al. (2019), and Simonsen et al. (2021). In general, our study agrees (within the uncertainty level) with the previous studies. However, on a regional scale, we observe large differences. Figure 13 shows elevation change rates of Jakobshavn Isbræ in m/yr from Simonsen et al. (2021) (downloaded from https://data.dtu.dk/). For consistency, we use the same color scale and area coverage as in Figures 11a–11i. The elevation change rates of Simonsen et al. (2021) do not show 10–15 m/yr thickening during 2016–2018 or large thinning during 2011–2016 along the main flow line of Jakobshavn Isbræ. However, our results are consistent with Khazendar et al. (2019) and Joughin et al. (2020), which both capture the substantial thickening during 2016–2018. The disagreement between this study and Simonsen et al. (2021) is likely due to a combination of coarse resolution used by Simonsen et al. (2021) and lack of observation close to the ice margin when using CryoSat‐2 only.

Figure 12

Annual mass change rates of the Greenland Ice Sheet. The black curve shows this study. The orange curve shows IMBIE (Shepherd and Ivins, 2020), the purple curve shows Colgan et al. (2019), and the blue curve shows Simonsen et al. (2021).

Figure 13

Elevation change rates of Jakobshavn Isbræ in m/yr from Simonsen et al. (2021). The grid resolution is 5 × 5 km grid. The elevation change rates are show during (a) 2011, (b) 2012, (c) 2013, (d) 2014, (e) 2015, (f) 2016, (g) 2017, (h) 2018, (i) 2019. The background shows a Landsat Image from 3 July 2020.

Implications for Solid Earth Deformation

Crustal VLM rates measured with GPS can be used to validate or constrain GIA models (Adhikari et al., 2021; Kappelsberger et al., 2021; Khan et al., 2016; King et al., 2010; Milne et al., 2018; Peltier et al., 2015). However, GPS measures the sum of GIA and elastic crustal VLM. Therefore, an estimate of elastic VLM is needed to isolate the GIA signal from the GPS time series. Figure 14a shows modeled mean elastic VLM rates during 2011–2020 using our high‐resolution grid as presented in this study. To model elastic VLM, we convolve our mass loss product with the Green's functions (see Section 2.3). Near the major outlet glaciers, we model elastic VLM rates of about 10–15 mm/yr. To highlight the striking influence of mass loss product on estimated elastic VLM, and or residual GIA, we model comparable elastic VLM rates for the same period and the same Green's functions using a comparable mass loss product derived from mass loss estimated of Simonsen et al. (2021) (see also Figure 14b). The differences between the two elastic VLM estimates are shown in Figure 14c. In general, our grids yield higher elastic rates: the difference is up to 10 mm/yr and, for some sites, is larger than the signal in the altimetry‐derived data (see Figure 14c). Our high‐resolution grid has a significant impact on modeled elastic crustal VLM.

Figure 14

(a) Mean elastic vertical land motion (VLM) rates in mm/yr based on 2011–2020 ice mass loss grids from this study, (b) mean elastic VLM rates based on 2011–2020 ice mass loss grids from Simonsen et al. (2021). (c) Difference between (a) and (b) in mm/yr.

Large differences in the magnitude and spatial distribution of ice elevation change rates between Figures 11 and Figure 13 yield, for example, significant differences in modeled elastic VLM rates at the GNET site KAGA (see Figure 15a). KAGA is a GPS site located on bedrock at the northern margin of Jakobshavn Isbræ. Figure 15b shows observed elastic VLM previously estimated at KAGA by Hansen et al. (2021) using GPS observations. The time series is corrected for GIA using the Khan et al. (2016) GNET‐GIA empirical model. We compare this previously observed elastic VLM with modeled elastic VLM derived from the annual elevation change rates presented in this study (red curve), as well as the modeled elastic VLM derived from annual elevation change rates presented in Simonsen et al., 2021 (blue curve). Due to the relatively coarse spatial and temporal resolution of the altimetry‐derived ice‐loss trends, both altimetry series do not adequately capture the patterns of mass loss needed to accurately estimate elastic VLM rates (Figure 15). To capture the seasonal cycle in melt and, therefore, elastic VLM, requires mass trend estimates at, ideally, about quarterly (i.e., 3 months) or finer temporal resolution. Additionally, as the uplift time series measured by the GPS is also sensitive to the spatial distribution of the ice loss load center it measures, and the effective resolution of marginal ice geometry is highly dependent on the spatial resolution of an altimetry product, the influence of marginal ice changes on uplift would be better captured with an altimetry product of better than the 1 km resolution presented here.

Figure 15

(top) Landsat Image (from 3 July 2020) of Jakobshav Isbræ. Red dot represent location of the KAGA GPS station (botton) black errorbars denote GPS time series of cumulative vertical land motion (VLM) at KAGA corrected for glacial isostatic adjustment. Red curve denotes cumulative elastic VLM based on mass loss from this study. Blue curve denotes cumulative elastic VLM based on mass loss from Simonsen et al. (2021).

Discussion and Conclusions

The average ice loss rate of the GIS was 273.2 ± 17.1 Gt/yr from April 2011 to April 2020. However, the last decade has been characterized by the progressive record‐breaking peak mass loss years. A record‐high mass loss occurred in 2011, followed by an even higher record mass loss in 2012, and finally the exceptionally high record mass loss of 2019 (Sasgen et al., 2020; Velicogna et al., 2020). We observed a total of 498 ± 45 Gt/yr of ice loss from April 2019 to April 2020, or a sea‐level equivalent of 1.4 mm during a single year. This rate of mass loss is equivalent to 15,850 tonnes of ice per second, continuously, for 12 months.

A striking feature of the eight basins considered, each of which includes tens of glaciers, is that there is a relatively weak or no correlation in the phase and/or pattern of mass change between sectors. Most, but not all, have a minimum in mass loss rate around 2017–2018. For basin D2, however, the minimum in mass loss rate is in 2015, followed by a monotonic increase in mass loss thereafter. No other basin, including those adjacent, shows the same behavior. This might be surprising given that the driver of changes in SMB is a relatively smooth function in space and changes in dynamics are influenced by some combination of SMB, ocean forcing, basal slipperiness/basal hydrology, along with changes in the driving stress. However, previous studies have shown that the dynamic changes of adjacent outlet glaciers have a complex response to a common ocean forcing, and are sensitive to highly local factors such as the subglacial topography close to the grounding line (Moon et al., 2012, 2014; Vijay et al., 2019; Wood & Rignot, 2021). In addition, Figure 8 indicates that, over a period of a decade, interannual variability in runoff and snowfall between basins can be large compared to the ice sheet wide trend.

Ice‐sheet model projections for the GIS, using a common climate forcing, have a large spread. For the RCP8.5 forcing scenario, the sea‐level contribution varies between 4 and 33 cm by 2100 (Aschwanden et al., 2019; Goelzer et al., 2020). The models also vary in the respective contribution from SMB and discharge that they estimate (Goelzer et al., 2020). It is, however, important that the models correctly partition trends between SMB and discharge and, as indicated in Figure 8, that they also correctly reproduce the spatial patterns of mass loss at a basin scale. We suggest, therefore, that our high resolution, annual time‐series is a valuable tool for calibrating prognostic ice sheet models (Choi et al., 2021).

Studies employing GNET GPS stations in Greenland should treat elevation change rates from altimetry with caution when using them to estimate the elastic component of VLM. The apparent agreement of ice sheet wide mass change estimates can be underlain by large regional scale differences, especially in terms of the resolution of marginal ice loss, that could result in significantly different estimates of associated elastic VLM rates GNET sites. Elastic VLM estimated using mass loss grids as presented in this study differs by up to 10 mm/yr from a recent comparable satellite altimetry study (Simonsen et al., 2021). Such large differences in elastic VLM correction can result in misinterpretations of past ice load history or earth's rheology when using elastic VLM corrected GPS observations to infer GIA. In addition to resolving recent mass loss to evaluate ice flow models, constraining longer term GIA signals provides a strong motivation to create improved community‐available altimetry mass‐loss products of high spatial and temporal resolution. Ideally, these mass‐loss products should be derived by aggregating data from the maximum number of altimetry sensors, and their ice‐sheet and sector‐scale mass‐loss rates should be cross‐evaluated with coincident independent gravimetry of input‐output solutions. This work represents an imperfect first step in this direction.