Extensive inland thinning and speed-up of Northeast Greenland Ice Stream – Nature

Calving front positions

We map the positions of the calving fronts of ZI and NG using aerial and Landsat 5–8 optical satellite imagery^51,52 from 2000 to 2021 (Extended Data Fig. 1). For ZI, the most notable changes occurred between 2011 and 2013, when a large portion of its floating extension collapsed, resulting in a retreat of 10 km (ref. ⁵³). This collapse has been followed by a steady retreat of about 700 m year⁻¹. NG lost large sections of its floating tongue from 2002 to 2004; after this, calving events of the main floating tongue have been less frequent. However, in 2020, the northern section of NG completely collapsed, releasing more than 120 km² of floating shelf ice into the ocean (Extended Data Fig. 1).

GNSS data processing

We process the GNSS data using the GIPSY-OASIS software package with high-precision kinematic data processing methods²⁷ and with ambiguity resolution using the orbit and clock products of the Jet Propulsion Laboratory (JPL). We use GIPSY-OASIS version 6.4, which was developed at the JPL³⁰. We use JPL final orbit products, which include satellite orbits, satellite clock parameters and Earth orientation parameters. The orbit products take the satellite antenna phase centre offsets into account. The atmospheric delay parameters are modelled using the Vienna Mapping Function 1 (VMF1) with VMF1 grid nominals⁵⁴. Corrections are applied to remove the solid Earth tide and ocean tidal loading. The amplitudes and phases of the main ocean tidal loading terms are calculated using the Automatic Loading Provider (http://holt.oso.chalmers.se/loading/), which is applied to the FES2014b ocean tide model including correction for the centre of mass motion of the Earth owing to the ocean tides. The site coordinates are computed in the IGS14 frame⁵⁵. We convert the Cartesian coordinates at 15-s intervals into local up, north and east coordinates for each GNSS site monitored at the NEGIS surface. An example of a 15-s solution is shown in Extended Data Fig. 2a–c.

GNSS-derived surface speeds and their uncertainties

We use the 15-s solution to estimate daily solutions of the ice surface speed (blue dots, Extended Data Fig. 2d–f). The time series have been screened for outliers. To remove outliers, we fit and remove a trend to each time series of speed, latitude and longitude. We estimate the mean of detrended speed, latitude and longitude and define outliers as values greater than three standard deviations from the mean. For NEG1, we removed in total eight data points or (8/387 = 0.02) 2% of data. For NEG2, we removed in total four data points or (4/195 = 0.02) 2% of data and 0% for NEG3. We use daily solutions of the ice surface speed (screened for outliers) to estimate the monthly mean ice speed (red error bars). We estimate the root-mean-square of each monthly mean to assign uncertainties. The black lines denote the trends (using least square adjustment) that best fit the observed ice speeds and represent mean flow accelerations from 2016 to 2019 of 2.7 ± 0.2 m year⁻² at NEG1, 4.5 ± 0.3 m year⁻² at NEG2 and 4.9 ± 0.3 m year⁻² at NEG3 (corrected for downhill acceleration).

Surface speed and acceleration from mosaics based on ESA Sentinel-1 SAR offset tracking

We derive ice speeds from mosaics based on ESA Sentinel-1 SAR offset tracking obtained from https://dataverse01.geus.dk/dataverse/Ice_velocity. The ice velocity maps of the NEGIS are derived from the intensity tracking of the ESA Sentinel-1 data with a 12-day repeat; the operational interferometric post-processing chain is applied for the analysis²⁹. We use all available speed mosaics (provided on a grid with a spatial resolution of 500 m) and associated standard deviation of the underlying shift maps generated by the offset tracking to create a time series of speed for each grid point (Extended Data Fig. 3b–d). We then remove outliers. To remove outliers, we fit and remove a trend to each time series of speed and estimate the mean. We define outliers as values greater than three standard deviations from the mean. For each grid point, we remove from 0 to 1.5% of data corresponding to 0 to 3 data points. Next, we use the screened time series at each grid point and fit a trend that represents the flow acceleration (the black line in Extended Data Fig. 3b,c). For each grid point, we use the least-squares fit to estimate the acceleration and associated uncertainty. Extended Data Fig. 3a shows the mean acceleration from 2016 to 2022 in m year⁻² at each grid point. The uncertainty is ± 0.7 m year⁻².

Downhill correction for GNSS

We use mosaics based on ESA Sentinel-1 to estimate the downhill correction of the flow acceleration at GNSS stations. We estimate the time series of the ice speed at two pixels on a velocity map on a 0.5 × 0.5-km grid. Pixel 1 is located at the GNSS starting position (when the station was deployed) and pixel 2 is located at the GNSS end-point position. The difference between the two time series is used to estimate the downhill correction for the flow acceleration. The uncertainty level in Extended Data Fig. 3a is ±0.7 m year⁻². However, two pixels close to each other (such as 1,000 m apart) experience almost the same noise. Thus, when estimating the difference between two neighbouring time series, the noise is reduced to ±0.3 m year⁻². To estimate the downhill correction of the flow acceleration, we fit a trend to the time series of the differences between two pixels.

Our corrections of the flow acceleration, which compensate for the downhill movement of the GNSS stations, are as follows: 0.06 ± 0.18 m year⁻² at NEG1, 0.14 ± 0.29 m year⁻² at NEG2 and 0.12 ± 0.30 m year⁻² at NEG3.

Elevation changes from CryoSat-2

To estimate elevation changes over the ice surface, we use a regular grid with a resolution of 500 × 500 m that covers the NEGIS. We denote the centre of each grid point as C(x₀, y₀). For each grid point, we select CryoSat-2 data with coordinates P(x_i, y_i), with a maximum distance of 1,000 m from C. The CryoSat-2 data points with coordinates P(x_i, y_i) have an elevation h_i measured at time t_i. The index i denotes the ith data point. We use all available CryoSat-2 data measured between July 2010 and July 2021 to create surface elevation time series at each grid point C. Previous studies have used a third-order polynomial equation to describe changes in the elevation and a third-order polynomial equation to describe the shape of the surface^2,56,57. However, to describe surface changes over 10 years, we fit a polynomial with a degree of 7 to describe changes in the elevation and we use 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 a centre C(x₀, y₀), we find the nearest data point within 1,000 m and fit a seventh-order polynomial H(t_i)_poly, a third-order surface topography polynomial H_topo and an annual term H(t_i)_Annual:

Extended Data Fig. 4a shows all the grid points on the NEGIS where we have successfully estimated the time series of H(t).

Extended Data Fig. 4b,c shows two examples of elevation time series H(t). P1 is a point closest to the ice margin and is located at about 590 m elevation. Extended Data Fig. 4b shows CryoSat-2 elevations corrected for the topography H_topo (black error bars) and a combination of our best-fitting seventh-order polynomial H(t)_poly and the annual term H(t)_Annual (red curve). The annual term at this elevation has an amplitude of 1.41 ± 0.10 m. Extended Data Fig. 4d shows CryoSat-2 elevations corrected for topography and the annual term (black error bars), and our best-fitting seventh-order polynomial (red curve) for point P1.

Extended Data Fig. 4c shows the same information as Extended Data Fig. 4b but for the point P2, which is located at about 1,078 m elevation. Here the amplitude of the annual signal is 0.16 ± 0.08 m.

We use the best-fitting seventh-order polynomial (such as the red curve in Extended Data Fig. 4d,e) for each location shown in Extended Data Fig. 4a to estimate the annual rates of elevation changes from April to April, for example, from April 2011 to April 2012, from April 2012 to April 2013 etc.

Elevation changes from ICESat-2

We use ICESat-2 data from October 2018 to June 2021 and estimate elevation changes over the ice surface³³. We use the ICESat-2 Algorithm Theoretical Basis Document for Land Ice Along-Track Height (ATL06) Release 004, which was retrieved from https://nsidc.org/data/atl06 (ref. ⁵⁸).

We estimate elevation changes using the same method described in the previous section. However, to describe surface changes over about 2.5 years, we fit a third-order polynomial (as opposed to the seventh-order polynomial used for CryoSat-2 data) and we also use a third-order polynomial equation to describe the shape of the surface and a seasonal term to account for the annual surface changes. Extended Data Fig. 5a shows the grid points on the NEGIS where we have successfully estimated the ICESat-2 elevation time series.

Extended Data Fig. 5 shows two examples of ICESat-2 elevation time series H(t). P3 is a point close to the ice margin and is located at about 267 m elevation (Extended Data Fig. 5b). Extended Data Fig. 5b shows ICESat-2 elevations corrected for topography (black error bars) and a combination of our best-fitting third-order polynomial and the annual term (red curve). Extended Data Fig. 5d shows ICESat-2 elevations corrected for topography and the annual term (black error bars), and our best-fitting third-order polynomial (red curve) at P3.

Extended Data Fig. 5c shows the same information as Extended Data Fig. 5b but for the point P4, which is located at about 1,071 m elevation.

We use the best-fitting third-order polynomial (red curves in Extended Data Fig. 5d,e) for each location shown in Extended Data Fig. 5a to estimate the annual rates of elevation changes from April 2019 to April 2020 and from April 2020 to April 2021.

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 (ref. ³²). The ATM flights are mainly concentrated along the margin of the GrIS. 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 the 2012 survey, between the 2012 survey and the 2013 survey etc. However, it should be noted that no survey was conducted in spring 2020. Extended Data Fig. 6 shows ATM flight lines over the NEGIS.

Extended Data Fig. 6 (lower-right panel) shows an example of elevation changes at the overlapping points from 2016 to 2017. The largest elevation change is observed near the glacier margin.

Gridded elevation changes and their uncertainties

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. Instead, we estimate the annual elevation change rates from April to April for each dataset independently and then merge the rates.

The observed annual elevation change rates from CryoSat-2, ICESat-2 and NASA’s ATM surveys are used to interpolate elevation change rates onto a regular grid of 1 × 1 km. The interpolation is performed using the ordinary kriging method^59,60. We use the observed annual elevation change rates to estimate an empirical semi-variogram. We fit a model variogram to the empirical semi-variogram to take the spatial correlation of elevation change rates into account. For each grid point, we estimate the elevation change rate dh_i,krig and the associated error σ_i,krig.

We correct the observed ice surface elevations for bedrock movement caused by elastic uplift owing to present-day mass changes and long-term past ice changes (glacial isostatic adjustment (GIA)). We correct for GIA using the GNET-GIA empirical model⁶¹. For each grid point, we estimate the GIA uplift rate dh_GIA and the associated uncertainty σ_GIA. We correct for the elastic uplift of the bedrock by convolving ice loss estimates (from CryoSat-2, ATM and ICESat-2) with the Green’s functions derived by Wang et al.⁶² for the elastic Earth model iasp91 with a refined crustal structure taken from CRUST 2.0. For each grid point, we estimate the elastic uplift rate dh_elas and the associated uncertainty σ_elas. The elevation change rate for each grid point is

and the associated uncertainty is

Extended Data Fig. 7 shows annual elevation change rates from April to April for the period from April 2011 to April 2021 that are corrected for GIA and elastic uplift. Extended Data Fig. 7 (lower panels) shows the uncertainties associated with the annual elevation change rates from April to April.

Ice flow modelling

We use the Ice-sheet and Sea-level System Model (ISSM)⁶³, a finite-element ice flow model, to model the NEGIS. The horizontal mesh resolution varies from 200 m near the ice front to 20 km inland, and it is vertically extruded into four layers. The model is based on a 3D higher-order approximation of the full Stokes model to include vertical shear for the stress balance; this is a good approximation for both fast-moving and slow-moving regions^64,65. We use the surface and bed geometry from BedMachine Greenland⁵⁰ (version 3). We infer the ice viscosity parameter over floating ice and basal conditions under grounded ice using inversions. To infer the initial basal conditions, we use the Budd friction law and invert for the friction coefficient using surface velocities from 2007 to 2008 (ref. ⁶⁶). We then analytically calculate the friction coefficient for a regularized Coulomb friction law that produces the same basal stress. The friction coefficient is kept constant in time during the simulations. Previous studies have shown that the choice of friction law can have a considerable impact on simulations^67,68,69,70. However, here we select the friction law based on observations.

The model uses the level-set-based moving boundaries to track ice front positions. We use the von Mises tensile stress calving law to calculate the calving rate^71,72. We calibrate the stress threshold of the calving law to match the observed ice front retreat from 2007 to 2017. We use the same stress threshold values, 1 MPa for grounded ice and 150 kPa for floating ice, used by Choi et al.²⁴. The calibrated stress thresholds for two glaciers (NG and ZI) are assumed to be constant for all simulations.

Using these nearly plastic friction laws, and the same atmospheric and oceanic forcing as Choi et al.²⁴, we calibrate the calving law to qualitatively match the observed front changes from 2007 to 2017 (ref. ⁷²). For hindcast simulation, the model is forced by monthly SMB data from the Regional Atmospheric Climate Model (RACMO) version 2.3 (ref. ⁷³). We apply ocean forcing that considers melting under floating ice for ZI and NG⁷⁴, along with the undercutting at the ice front of the terminus of ZI once it becomes grounded^75,76. The details of the ocean forcing parameterizations can be found in Choi et al.⁸. We keep the ice temperature constant, as it is not affected by the surface temperature on the timescales considered in this study⁷⁷. We find that the model with a regularized Coulomb law yields much better agreement with observed accelerations (Fig. 2d,e) and the observed mass loss from 2011 to 2021 (Fig. 4).

Friction law selection

Extended Data Fig. 8 shows maps of the observed and modelled ice flow acceleration from 2016 to 2022. The observed accelerations are identical to those shown in Extended Data Fig. 3a. Here we show the Budd friction law³⁵ with different exponents, the regularized Coulomb friction law^36,69, the Schoof friction law³⁷ and the Weertman friction law³⁸ with different exponents. We tested a wide range of models with different exponents; however, our model results indicate that only the regularized Coulomb friction law or the Budd friction law with exponents of 1/5 or 1/6 can produce the deep inland acceleration observed in the satellite data. For all models in Extended Data Fig. 8, we used NorESM1 and RCP4.5. The choice of the GCM and RCP has a much smaller impact than the choice of the friction law.

In Extended Data Fig. 9, we consider models that are able to reproduce deep inland acceleration. Therefore, the choice of friction laws is reduced to the regularized Coulomb friction law and the Budd friction law with exponents of 1/5 and 1/6, as all the other models were unable to generate deep inland acceleration.

Extended Data Fig. 9 shows the modelled cumulative changes in the ice mass from 2007 to 2100 for the regularized Coulomb friction law and the Budd friction law with exponents of 1/5 and 1/6. For each friction law, we model the mass change using RCP4.5 and RCP8.5, respectively. For all models, NorESM1 was used. For comparison, we include results from ref. ⁸ that were calculated using the Budd friction law (linear viscous). Extended Data Fig. 9 suggests that only the mass change from the regularized Coulomb friction law is within the uncertainty level of the observed mass change from 2007 to 2021. The Budd friction law with exponents of 1/5 or 1/6 slightly overestimates the mass change from 2007 to 2021. However, we note that, by 2100, the mass loss from the Budd friction law with exponents of 1/5 or 1/6 is almost the same as that from the regularized Coulomb friction law. The retreat of the terminus from 2007 to 2100 using NorESM1, RCP4.5 and regularized Coulomb friction law is shown in Supplementary Video 1.