Glaciological Setting.

ZI is 20 km wide. Satellite imagery shows that the ice shelf area decreased from 706 $$ {\mathrm{k}\mathrm{m}}^2$$ in 1985 to 616 $$ {\mathrm{k}\mathrm{m}}^2$$ in 2002, and, subsequently, to 37 $$ {\mathrm{k}\mathrm{m}}^2$$ in 2014; 377 $$ {\mathrm{k}\mathrm{m}}^2$$ of the detached ice shelf remains in the fjord to the north of the glacier (Fig. 1). Between 2014 and 2019, we find that the ice front retreated 1.6 km to lose another 32 $$ {\mathrm{k}\mathrm{m}}^2$$, leaving virtually no floating section. In ref. 21, we reported that the glacier speed increased from 1.2 km/y in 1979–1992 to 2.1 km/y in 2019, or 81% (22); the grounding line flux increased from 9.7$$ \pm $$1 Gt/y in 1979 to 16.2$$ \pm $$2 Gt/y in 2019, or 81%, versus a balance flux of 8.9 Gt/y; the mass loss averaged 5 Gt/y in 2009–2019 (SI Appendix, Table S1).

Glacier 79N is 20 km wide, with a speed that increased from 1.3 km/y in 1979 to 1.5 km/y in 2019, or 12%, and a grounding line flux that increased from 11.7 km/y to 13.1 Gt/y versus a balance flux of 9.6 Gt/y (21). The mass loss of 79N averaged 2 Gt/y in 1979–1990 and 4 Gt/y in 2009–2019.

ZI lost mass from enhanced flow but also from increased surface melt, which produces runoff. Runoff production reconstructed by the Regional Atmospheric Climate Model (RACMO)2.3p2 regional atmospheric climate model (23) increased from 1.5$$ \pm $$0.3 Gt/y in the 1980s to 3.3$$ \pm $$0.7 Gt/y in the 2010s, or 126%. A time series of surface elevation updated from ref. 21 indicates that ZI thinned at 0.8 m/y in the 1980s and 1990s and peaked at 3 m/y in 2012–2013, before turning into slight thickening in the last couple of years (SI Appendix, Table S1). For 79N, runoff production increased from 1.4 Gt/y to 4.8 Gt/y from 1979–1989 to 2009–2019, or 240%. The glacier thinned at 0.1 m/y in the 1970s, peaking at 1.0 m/y in 2010, with lower thinning in 2011–2017 and slight thickening in 2018–2019.

Bathymetry.

The glacier bed topography reconstructed from BedMachine Greenland (24) reveals a grounding line depth of 600 m in 1996 at the center of ZI and 79N and lower values toward the margin. Operation IceBridge (OIB) airborne gravity in front of the glacier revealed that Jøkelbugten Fjord is overdeepened, with 900 m depth about 75 km south of the grounding line. From 28 May to 30 June 2016, the OMG mission collected higher-resolution, airborne gravity data over the ocean between ZI and NT using the Sander Geophysics Ltd. AIRGrav deployed on a Cessna Grand Caravans 208B from the airports of Kulusuk and Station Nord (25). The survey was conducted at 110 knots, with a ground clearance of 150 m, 2-km spacing inshore, and 4-km spacing offshore (Fig. 1). Onboard receivers were NovAtel OEMV-3, while the reference station used a NovAtel OEM4 GPS receiver. Outer lines were gridded at 1-km grid cell size versus 500 m for inner lines. The data were filtered with a 1-km half-wavelength filter with an rms error of 1.5 mGal.

We employ the Geosoft GM-SYS 3-D with Parker’s method (26) to minimize the misfit between calculated and observed gravity. The model domain has three horizontal layers: 1) a solid ice layer with a density of 0.917 g/$$ {\mathrm{c}\mathrm{m}}^3$$, 2) an ocean water layer with a density of 1.028 g/$$ {\mathrm{c}\mathrm{m}}^3$$, and 3) a rock layer with a density of 2.67 g/$$ {\mathrm{c}\mathrm{m}}^3$$. A forward model of the gravity is calculated using the bed elevation from BedMachine v3 (BMv3) (24) and new bathymetry measurement in front of 79N (10, 27) (SI Appendix, Fig. S2). We calculate the direct current (DC) shift, or offset between modeled and observed gravity in areas where bed elevation is known from conductivity, temperature, and depth (CTD), BMv3, or multibeam soundings. We interpolate the DC shift onto a regular grid using a minimum curvature algorithm, correct the observed gravity with the interpolated DC shift, fill the data gaps with the model results, and invert the resulting gravity field where bed elevation is not known (28). The interpolation of the DC shift accounts for natural variations in underlying geology across the model domain caused by variations in crustal thickness, sedimentary basins, or intrusions. The inversion stops when the misfit between observed and modeled gravity is less than 0.1 mGal (SI Appendix, Fig. S3).

The results reveal an 800-m-deep trough between Schnauders Ø Island and Norske ØER. Jøkelbugten Trough shallows to the south near a set of islands. The new trough, named here Schnauders Ø Trough (ST), connects with NT and the channels leading to 79N. ST has a maximum depth of 800 m versus 1,100 m for Jøkelbugten Fjord and a minimum depth of 350 m at two sills (Fig. 3B).

Fig. 3.

Longitudinal profiles of ocean temperature color coded from blue (–$$ {\mathrm{2}}^{○}$$C) to red (+$$ {\mathrm{2}}^{○}$$C) with 0.$$ {\mathrm{25}}^{○}$$C contour level for (A) 79N (profile A-C, purple in Fig. 1) and (B) ZI (profile B-C, blue in Fig. 1) with sea floor in black, and areas with no water temperature data in white. AXCTD 1 to 5 locations are indicated with a triangle and number (red). Temperature data (C) from AXCTD 1 to 5 with Gades melt water line and runoff line, and (D) from historical CTD data (SI Appendix, Fig. S1) in NT for years 1965–2019 with (E) salinity in practical salinity units (psu).

Oceanographic Data.

In the summer of 2019, NASA’s OMG mission deployed Airborne Expendable CTD (AXCTD) ocean probes in front of ZI (29). Additional AXCTD data were collected on a NASA Gulfstream III aircraft in September to October 2016, C-130 Hercules in October 2017, a Basler DC-3 Turbo Prop in September 2018, and a Basler BT-67 in September 2019. AXCTD accuracy is $$ \pm 0.{\mathrm{1}}^{○}$$C for temperature and $$ \pm $$2% for depth. The probe in front of ZI reached 627 m depth, with a bottom temperature of 1.$$ {\mathrm{4}}^{○}$$C (Fig. 3). The probe in the deepest part of 79N reached 453 m depth, with a bottom temperature of 1.$$ {\mathrm{36}}^{○}$$C.

At AXCTD 1 to 5, the temperature/salinity relationship indicates a mixture of warm AIW at 1.$$ {\mathrm{5}}^{○}$$C and 35 psu with ice shelf melt water along the Gades meltwater mixing line. At AXCTD 4 at the front of ZI, we detect a strong runoff signal consistent with the presence of a plume of subglacial water (Fig. 3). We complement the CTD data with Estimating the Climate and Circulation of the Ocean (ECCO) model annual averages between 1992 and 2011 (see Materials and Methods) (30, 31). The rms error between CTD and ocean model is 0.55 $$ °$$C.

To visualize the influence of bathymetry on ocean temperature, we use the line of lowest elevation through the trough systems starting from the grounding line to the middle of NT (Fig. 3). Bed elevation rises to 350 m depth about 170 km from ZI at the junction between ST and NT, followed by a second sill at the same depth common to ZI and 79N at the junction with NT. For 79N, the bed elevation is deepest at the center of the calving front (480 m) but is limited by an upstream sill with a depth of 325 m (10). In contrast, ST offers a broader, 25-m-deeper pathway for AIW to reach ZI. Water temperature of 1.$$ \mathrm{5}\hspace{0.17em}°$$C in NT below 300 m depth is blocked by the first sill at AXCTD 3. Water at 1.$$ \mathrm{25}\hspace{0.17em}°$$C reaches only the second sill in front of 79N but reaches the ZI grounding line. The temperature of AIW changes little from AXCTD 1 to 5; that is, subsurface ocean waters are transferred with relatively low heat loss toward the glacier fronts.

Historical CTD data spanning 1965 to 2019 collected in the inner shelf of NT where AIW circulates (13) document the evolution of ocean temperature below the surface and above the main sills (Fig. 3 and SI Appendix, Table S2). We calculate ocean thermal forcing, $$ TF$$, as the depth-integrated difference between the in situ temperature of the ocean and the pressure-dependent and salinity-dependent freezing point of seawater. We integrate $$ TF$$ at two ocean flux gates, one located at the mouth of 79N ice shelf and one located at the northern entrance of ST (Fig. 1). We estimate thermal forcing in the lower 40% of the water column to have increased by 1.3$$ \pm $$0.5 $$ °$$C for both glaciers in 1965–2019, with significant interannual variability in between. The thermocline rose 100$$ \pm $$50 m from 1965–1979 to 2016–2019.

Ice–Ocean Interaction.

From $$ TF$$, water depth, and subglacial discharge, we calculate a rate of grounded ice undercutting, $$ q_m$$, using a parameterization derived from high-resolution (1 m) modeling of ice melt along a vertical wall by the Massachusetts Institute of Technology General Circulation Model (MITgcm) ocean model in three dimensions (32) (see Materials and Methods). Subglacial water discharged is imposed at the grounding line as a source of freshwater. Boundary conditions for the ocean model include temperature and salinity versus depth several kilometers away from the grounding line. Ice melt is predicted by the model to be uneven, with peak rates 30 m to 50 m above the sea floor and lower rates toward the surface (Fig. 2). The uneven distribution of melt yields glacier undercutting which has been confirmed in situ with multibeam echo sounder data in west Greenland (20, 33) and Alaska (34). Undercutting has not been observed at ice shelf grounding lines, but we note that radar echograms reveal sharp transitions in ice thickness at that location, on multiple years, and at multiple locations (Fig. 2). Radar echograms cannot reveal cavity shapes, however, because radar signals do not penetrate seawater.

Undercutting of grounded ice by the ocean directly affects the glacier force balance. Ice above the area of removal will not be supported by basal friction, hence will offer no resistance to flow. Ice below the area of removal will be too thin to provide basal resistance. Undercutting of grounded ice therefore directly reduces basal resistance to flow. To illustrate its importance, let us assume a basal drag of 100 kPa or 100 kN/$$ {\mathrm{m}}^2$$ at the grounding line and, similarly, a lateral drag of 100 kPa along the sides of an ice shelf. A 1-km retreat of a 20-km-wide glacier will reduce the buttressing force by 2,000 GN. For comparison, a 20-m thinning of a 50-km-long ice shelf will reduce the buttressing force by only 200 MN. The total removal of a 50-km-long, 300-m-thick ice shelf will be necessary to reduce the buttressing force by 3,000 GN. In terms of force balance, a 1-km retreat of the grounding line is therefore equivalent to the complete loss of a long ice shelf.

Using our ocean model, we calculate that $$ q_m$$ increased from 108$$ \pm $$28 m/y in 1979–1989 to 185$$ \pm 48$$ m/y in 2009–2019 for ZI. For 79N, $$ q_m$$ increased from 52$$ \pm 13$$ m/y in 1979–1989 to 100$$ \pm $$26 m/y in 2009–2019. These values, which are averages across the glacier width, are comparable to that reported near the grounding zone of 79N (17).

To quantify the impact of undercutting on the glacier retreat, we use a reference state, $$ q_m^{ref}$$, when the glacier was in steady state. For 79N, thinning was small in 1979 (0.1 m/y), and the glacier was not speeding up, so we use year 1979 as a reference. We calculate a cumulative $$ Q_m$$ of 1.7 km in 1979–2019. For ZI, ice thinning was 0.6 m/y in 1979, which is significant. The water temperature was 0.$$ \mathrm{2}\hspace{0.17em}°$$C cooler in 1965 compared to 1979, which could be due to seasonal differences, since the CTD was collected in winter, or express colder ocean conditions prior to 1979 (13). Using the 1965 reference for ZI, we calculate $$ q_m^{ref}$$ of 59 m/y (vs. 80 m/y in 1979–1980) and a cumulative $$ Q_m$$ of 3.5 km in 1979–2019.

Thinning-Induced Retreat.

Thinning-induced retreat, $$ q_s$$, is deduced from the observed glacier thinning and calculated, time-dependent bed and surface slopes from BMv3 (SI Appendix, Fig. S4). The $$ q_s$$ increased from 104$$ \pm $$27 m/y in 1979–1989 to 217$$ \pm $$56 m/y in 2009–2019 for ZI, yielding a $$ Q_s$$ of 7.6 km for 1979–2019 ($$ q_s^{ref}$$ = 0). For 79N, $$ q_s$$ increased from 27$$ \pm $$7 m/y in 1979–1989 to 52$$ \pm $$15 m/y in 2009–2019, yielding a $$ Q_s$$ of only 2.1 km for 1979–2019, because the glacier retreated along prograde slopes (bed elevation rises in the inland direction) instead of retrograde slopes for ZI.

Ice thinning includes enhanced surface melt from warmer air temperature and dynamic thinning from flow acceleration. In 41 y, ZI thinned by 49$$ \pm $$7 m at the thickness line used for ice fluxes (Fig. 1). At that location, changes in surface mass balance relative to the reference period 1961–1990 resulted in a cumulative surface lowering of 5$$ \pm $$1 m, or 10 times lower than observed (SI Appendix, Fig. S4). Most of the thinning is therefore of dynamic origin. For 79N, total thinning averaged 12$$ \pm $$2 m versus 8$$ \pm $$1 m from surface mass balance; that is, surface mass balance dominates, consistent with its more steady flow.

Grounding Line Retreat.

Grounding line retreat, $$ q_{gl}$$, has been observed using differential radar interferometry (InSAR) since 1992 on 79N, 1996 on ZI (35), and 2000, 2011, and 2014 on both glaciers to reveal a rapid retreat for ZI versus a slow retreat for 79N (21). To complement this record, we use a digital elevation model (DEM) from year 2001 (36) to calculate where ice first reached hydrostatic equilibrium, which is a proxy for the grounding line position. We select a firn depth correction (7 m) to transform ice surface elevation into solid-ice equivalent thickness that best fits the DEM-derived grounding line and the InSAR-derived grounding line in 2001. We verify the reliability of the 7-m firn correction by comparing the results of a DEM from 2015 with the SAR-derived grounding line from 2014. We then apply the same firn correction to a DEM from 1978 (37) and to ICESat-2 data from 2019 to derive the grounding line locations in 1978 and 2019, respectively. The grounded ice loss estimated from these data totals 217 $$ {\mathrm{k}\mathrm{m}}^2$$ for ZI and 130 $$ {\mathrm{k}\mathrm{m}}^2$$ for 79N from 1979 to 2019, which translates into a linear grounding line retreat $$ Q_{gl}$$ of 13 km for ZI and 4.4 km for 79N, with an uncertainty of 500 m (Fig. 4).

Fig. 4.

Evaluation of the components of the grounding line retreat of (A and C) ZI and (B and D) 79N, Greenland. (A and B) Time series of runoff production, $$ q_{sg}$$, (cubic kilometers per year, blue, left scale), thermal forcing, $$ TF$$ (degrees Celsius, green, right scale) and change in surface elevation, $$ dh/dt$$, (meters per year, red, right scale) for (A) ZI and (B) 79N. (C and D) Cumulative anomaly in glacier undercutting, $$ Q_m$$, (kilometers, blue), cumulative thinning-induced retreat, $$ Q_s$$, (kilometers, red) versus the observed grounding line retreat, $$ Q_{gl}$$, (kilometers, black) with SE in light color for (C) ZI and (D) 79N for the years 1979–2019. If the components are correct, $$ Q_m$$ + $$ Q_s$$ should balance $$ Q_{gl}$$.

Mass Conservation.

Mass conservation at the grounding line dictates that the observed grounding line retreat, $$ q_{gl}$$, results from three physical processes: 1) grounded ice undercutting, $$ q_m$$, above equilibrium conditions; 2) thinning-induced retreat, $$ q_s$$, caused by surface melt and increases in flow speed; and 3) calving of grounded ice blocks, $$ q_c$$, if applicable (Fig. 2). The cumulative grounding line retreat, $$ Q_{gl}$$, balances the cumulative anomaly in undercutting, $$ Q_m$$, cumulative thinning-induced retreat, $$ Q_s$$, and—if applicable—the cumulative calving of grounded blocks, $$ Q_c$$. $$ Q_s$$ and $$ q_s$$ are in meters and meters per day, respectively, and are positive when forcing a grounding line retreat.

On ZI, we observe a grounding line retreat of 200 m/y in 1979–1989 and 550 m/y in 2009–2019, for a $$ Q_{gl}$$ of 13.2 km. Compared with the estimates of undercutting and thinning, 30% of the retreat is caused by $$ Q_m$$, 58% by $$ Q_s$$, and 12% is an error or $$ Q_c$$ after 2014. For 79N, the grounding line retreated at 131 m/y in 1979–1989 and 149 m/y in 2009–2019 for a total $$ Q_{gl}$$ of 4.4 km. Compared with our calculations, 39% of the retreat is due to $$ Q_m$$, 47% to $$ Q_s$$, and 14% to errors. The agreement between calculations and observations is therefore excellent for both glaciers. If undercutting by the ocean were ignored, the effect of ice thinning would only explain half of the observed retreat on 79N and ZI. In the absence of grounded ice removal by the ocean, the grounding line retreat would be significantly underestimated.

Sea Level Equivalent.

Using BMv3, the volume of ice, $$ V$$, in the drainage basins of each glacier above flotation is converted into a sea level equivalent after a change of density from $$ {\rho }_i$$ = 917 kg/$$ {\mathrm{m}}^3$$ (ice density) to $$ {\rho }_w$$ = 1,028 kg/$$ {\mathrm{m}}^3$$ (seawater density) (i.e., $$ V$$ × $$ {\rho }_i/{\rho }_w$$) divided by the mass of ocean water equivalent to 1 mm of sea level rise (362 Gt), minus the volume of ice below flotation diluted from ice to seawater ($$ V$$ × ($$ {\rho }_w-{\rho }_i$$)/$$ {\rho }_i$$)), divided by 362 Gt, a smaller number. We find 0.54 m for ZI (−0.002 m for dilution) and 0.57 m for 79N (−0.004 for dilution). Flotation is from the freeboard height, $$ f_b=\left(1-{\rho }_w/{\rho }_i\right)\hspace{0.17em}b$$, where $$ b$$ is bed elevation (42).

Runoff.

Runoff, $$ q_{sg}$$, is reconstructed with the RACMO2.3p2 downscaled to 1 km (43) from 1958 to 2019. We integrate runoff production over the drainage basin of each glacier to obtain a flux in cubic meters per day and annual production in gigatons per year. We divide the annual flux by the cross-section of the ocean flux gate (8,986,662 $$ {\mathrm{m}}^2$$ for ZI and 6,081,095 $$ {\mathrm{m}}^2$$ for 79N) to yield an area average speed of subglacial discharge, $$ q_{sg}$$, in meters per day. The $$ q_{sg}$$ varies from 0.60 m/d to 1.64 m/d for ZI during 1979–2019 and from 0.98 m/d to 3.74 m/d for 79N (SI Appendix, Table S1).

Ocean Thermal Forcing.

Ocean thermal forcing, $$ TF$$, is from 19 CTDs from 1965 to 2019 in NT (SI Appendix, Fig. S1 and Table S2). In 2019, we estimate the fraction of ocean temperature below 300 m transmitted from NT to the ice fronts to be 91.6% for 79N and 90% for ZI. We use this transfer coefficient to scale the NT temperature to a temperature equivalent at the ice fronts. The ocean gate is located at the mouth of the ice shelf front for 79N (10) and at the northern entrance of ST for ZI. TF and $$ q_m$$ are calculated across the ocean flux gates over the lower 40% of the water column, because high-resolution modeling of ice–ocean interactions indicates that it is the lower 40% of the water column that controls thermal forcing. To complement the CTD data, we use two ECCO project model outputs: 1) a high-resolution (4 km) forward model of the Arctic from an initial state for the period 1992–2011 (30) and 2) a medium-resolution (13.5 km) global domain solution for 2001–2017 (LLC270) (31). The two solutions are calibrated using CTD data to remove an absolute bias and a linear trend in temperature.

Grounded Ice Undercutting.

Grounded ice undercutting, $$ q_m$$, is calculated based on high-resolution (1 m), three-dimensional simulations of a melt plume using the MITgcm ocean model with varying water depth, b, subglacial discharge, $$ q_{sg}$$, and thermal forcing, $$ TF$$, as $$ q_m$$ = (0.0003 b $$ q_{sg}^{0.33}$$ + 0.15) $$ {\mathrm{T}\mathrm{F}}^{1.18}$$ (19). We integrate $$ q_m$$ across the ocean flux gates using the same temperature data. The uncertainty in $$ q_{sg}$$ is 20%, $$ b$$ is 5%, $$ TF$$ is 0.55 $$ °$$C, and the uncertainty associated with the spatial distribution of subglacial discharge across the ice face is 15% (19), yielding a nominal uncertainty in width-averaged $$ q_m$$ of 26%. Error in $$ Q_m$$ progresses as the square root of the number of years. We derive a time series of cumulative undercutting, $$ Q_m$$, which adds up to 3.5 km for ZI and 1.7 km for 79N.

Thinning-Induced Retreat.

Thinning-induced retreat, $$ q_s$$, is deduced from a time series of elevation based on DEMs from 1978 (37), 2001 and 2006 (36), and 2015–2016 (44), combined with airborne altimetry from NASA’s OIB spanning 1993 to 2019 and a time series of ice velocity (21) (SI Appendix, Fig. S4). We calculate ice thinning averaged across the glacier width, bed slope, $$ \beta $$, and surface slope, $$ \alpha $$, at the glacier center over several ice thicknesses to convert the rate of ice thinning, $$ dh/dt$$, into a rate of flotation retreat, $$ q_s$$ = $$ dh/dt$$ / ([1 – $$ {\rho }_w$$/$$ {\rho }_i$$] $$ \beta $$ – $$ \alpha $$) (45) in meters per year. The cumulative thinning-induced retreat, $$ Q_s$$, totals 7.6 km for ZI and 2.1 km for 79N from 1979 to 2019. DEM have a 1- to 2-m noise level over 80-m peak thinning, or 2%. Average thinning is 10 m/y, hence 15%. The noise in bed and surface slopes is 10%, yielding an error in Q_s of 18% for ZI and 25% for 79N.