Vulnerability of Southeast Greenland Glaciers to Warm Atlantic Water From Operation IceBridge and Ocean Melting Greenland Data.

We employ National Aeronautics and Space Administration (NASA)'s Operation IceBridge high-resolution airborne gravity from 2016, NASA's Ocean Melting Greenland bathymetry from 2015, ice thickness from Operation IceBridge from 2010 to 2015, and BedMachine v3 to analyze 20 major southeast Greenland glaciers. The results reveal glacial fjords several hundreds of meters deeper than previously thought; the full extent of the marine-based portions of the glaciers; deep troughs enabling warm, salty Atlantic Water (AW) to reach the glacier fronts and melt them from below; and few shallow sills that limit the access of AW. The new oceanographic and topographic data help to fully resolve the complex pattern of historical ice front positions from the 1930s to 2017: glaciers exposed to AW and resting on retrograde beds have retreated rapidly, while glaciers perched on shallow sills or standing in colder waters or with major sills in the fjords have remained stable.

Introduction

The Greenland Ice Sheet has been losing mass rapidly, and the rate of ice mass loss has been accelerating (e.g., Enderlin et al., 2014; Rignot et al., 2008). Between 1996 and 2000, southeast Greenland (SEG) was a major contributor to the mass loss from the Greenland Ice Sheet. In 1994–1999, ice thinning extended to the ice divide with a total mass loss above 17 Gt/yr (Krabill et al., 1999). This could not be explained by an increase in ice melt alone and required a significant participation from ice dynamics (Abdalati et al., 2001; Rignot et al., 2004; Rignot & Kanagaratnam, 2006). During the last two decades, the glaciers experienced distinct periods of synchronous speedup and slowdown (Howat et al., 2008; Moon et al., 2012). Individual glaciers have, however, shown high variability in flow regime and frontal position, indicating that the glacier response to climate forcing is not uniform at the fjord level (Bjørk et al., 2012; Moon et al., 2012; Murray et al., 2010).

Off the coast of SEG, the Irminger Current branches off from the North Atlantic Current to bring warm, salty Atlantic Water (AW) in contact with Greenland glaciers to fuel high melt rates along coastal margins (Christoffersen et al., 2011; Murray et al., 2010). Interpreting the glaciological changes in this region, however, has been hampered by a dearth of data, including knowledge of glacier thickness, fjord depth, and ocean temperature. Most fjords have not been mapped, many glaciers have no name, and measurements of glacier thickness in terminal valleys are affected by uncertainties due to high snowfall and liquid water content that makes the radar depth sounder (RS) challenging (Miege et al., 2016). The presence of sills or deep fjords that could block or facilitate the access of AW is unknown or has remained hypothetical (Buch, 2002). Without these details, it is difficult to understand the impact of ice‐ocean interactions on the glaciers and resolve the pattern of ice retreat that affected the ice sheet mass balance in this region in the past decades (Enderlin et al., 2014).

Upstream of the glacier fronts, RS have only measured ice thickness at high elevation (≥2,000 m), tens of kilometers from the calving fronts (Gogineni, 2012). The radar results have been used to reconstruct the bed topography using a mass conservation (MC) approach (Morlighem et al., 2014, 2017). Downstream of the flux gates, few ice thickness data exist since bed echoes are rare in the echograms. The quality of the MC depends on the precision of the surface mass balance (SMB) reconstruction in a region with steep topography and narrow fjords (a few kilometers) versus the resolution of regional atmospheric climate models (11 km for RACMO2.3, Noël et al., 2016). The MC method also requires a precise knowledge of glacier thinning, which is obtained from sparse tracks of OIB or spaceborne altimetry (Khan et al., 2014) over complex topography. Hence, the uncertainty in ice thickness reconstruction accumulates toward the glacier fronts and is maximal at the ice margin, possibly exceeding hundreds of meters.

The 2016 Ocean Melting Greenland (OMG) mission collected multibeam echo sounding (MBES) data in the fjords of SEG. In some fjords, the data extended to the glacier fronts, providing precise constraints on calving front thicknesses. In other fjords, the glacier fronts could not be reached by boat due to the presence of heavy brash ice and iceberg debris in front of the glaciers. In 2016, we conducted an extensive airborne high‐resolution gravity survey of key glacier fronts to complement the OMG data and provide observational constraints on the most challenging sectors of Greenland. This gravity survey builds upon prior efforts at using high‐resolution gravity data to resolve bed topography in challenging sectors (An et al., 2017).

Here we combine the OIB airborne gravity data and OMG MBES data in 10 major fjords occupied by more than 20 outlet glaciers. We constrain a three‐dimensional inversion of the gravity observations using the MBES data in the fjord and MC on land to obtain a seamless bed bathymetry mapping across the ice margin. We estimate the uncertainty of the reconstruction and discuss the impact of the results on interpreting the detailed, recent deglaciation history of this part of Greenland.

Data and Methods

Southeast Greenland is characterized by a rough topography, with elevation ranging from sea level to 3,000 m elevation for the surrounding peaks. The region includes some of the fastest tidewater glaciers in Greenland, with ice velocity ranging from 2 to 8 km/yr (KøgeBugt C., Figure S1 in the supporting information). The drainage basin is 81,109 km2 in area (Mouginot et al., 2017) with a sea level equivalent of 0.3 m (Table S1). We adopt the glacier name nomenclature of Rignot and Mouginot (2012). The survey area is divided into nine blocks spreading from Narsarsuaq in the south to north of Sujunikajik (Figure 1). The nine blocks include B1 for the 4 Ikertivaq glaciers, B2 for the 3 KøgeBugt glaciers, B3 for Graulv and Gyldenløve, B4 for A. P. Bernstorff, B5 for Skinfaxe and Rimfaxe, B6 for Tingmiarmiut and the 3 Mogens glaciers, B7 for 2 glaciers south of Puisortoq, B8 for Anorituup, and B9 above Qajuutap and Eqalorutsit glaciers (Figure 1).

Figure 1

Free‐air gravity anomalies (red to blue) in southeast Greenland overlaid on a shaded relief of the 30 m resolution latest version of the Greenland Ice Mapping Project (GIMP) DEM. Ocean Melting Greenland multibeam echo soundings are in shaded relief on a color scale from blue (deep) to orange (shallow). Green diamonds are Ocean Melting Greenland conductivity‐temperature‐depth (CTD) measurements. Glacier symbols mark the stability of the present front (unstable = triangle and stable = circle), size of symbol is proportional to the balance flux, and color qualifies the retreat (red = retreat, blue = no retreat on a sill, and green = no retreat, not understood). AW = Atlantic Water.

OIB Gravity data were acquired on a AS350‐B3 helicopter by the Sanders Geophysics Airborne Inertially Referenced Gravimeter. Eight flights took place between 27 July 2016 and 9 August 2016 for a total 5,445 km between Narsarsuaq and Kulusuk airports, which were used as gravity references (Figure 1). The gravity lines were flown at a line spacing of 1 km, with a ground clearance of 80 m, and at a ground speed of 70 knots. The gravity field is sampled at a frequency of 128 Hz. GPS data were recorded at 10 Hz. We used a Riegl laser altimeter with a single optical laser beam to measure distance to the ground, with a range of 1,500 m, a resolution of 1 cm, and an accuracy of 5 cm at a 3.3 Hz sampling rate. Eotvos, normal gravity, free‐air, static, and level corrections were applied to the data to yield free‐air gravity anomalies, which we refer to as “gravity data” in the remainder of the paper. The final product we use was processed using a 20‐sec filter, gridded, and low‐pass filtered with a 750 m half wavelength. Data noise has been estimated at less than 0.2 mGal based on a crossover analysis.

Gravity inversion is performed using Geosoft GM‐SYS 3D, which implements the method of Parker (1973) on a 3‐D representation of the surveyed area with (1) a solid ice layer with a density of 0.917 g/cm3, (2) a sea water layer with a density of 1.028 g/cm3, and (3) a rock substrate layer with a density discussed below (An et al., 2017; Gourlet et al., 2016; Millan et al., 2017). A 3‐D forward model of the gravity anomalies is obtained by combining BedMachine version 3 (referred as BM3) over land, OMG at sea, and a linear interpolation in between to obtain gravity values at every point of the inversion domain (Figure 1). We merge the modeled and observed gravity by adjusting their mean value within a subregion where the OIB gravity and OMG MBES overlap. Where there is no overlap, we use an open land area to calculate this constant shift (B3, B4, and B9). We allow for a smooth transition between observed and modeled gravity over a length scale of 500 m to 1,500 m depending on the width of the fjord. The inversion is conducted everywhere on grounded ice and extends from the edges of the OMG data to the limit of the OIB gravity survey. The resolution of the gravity data is 750 m, but the results are gridded at a 150 m spacing for consistency with BM3 (Millan et al., 2017).

We calculate an optimal bedrock density as the density that minimizes the root‐mean‐square error between modeled and observed gravity (An et al., 2017; Gourlet et al., 2016; Millan et al., 2017). Where we do not find a minimum, we use a default rock density of 2.67 g/cm3. Because we have no independent information on sediment thickness (Overeem et al., 2017), our gravity inversion does not include a sediment layer. An et al. (2017) showed that a 100 m thick sediment layer only changes the bed elevation by 30 m, which is within the uncertainty of our inversion, defined as the misfit between modeled and observed gravity. We translate this misfit into uncertainty in bed elevation using the procedure in An et al. (2017) and Millan et al. (2017), which uses a conversion factor of 5.8 ± 0.5 mGal in gravity per 100 m of water. Where reliable RS data are available, we use them to verify the accuracy of our gravity‐derived bed elevation.

Ice fluxes are used to evaluate the quality of our gravity‐derived bed elevation and ice thickness. We calculate ice discharge at the ice front using ice velocities spanning from 1985 to 1995 derived from Landsat and synthetic aperture radar observations (Mouginot et al., 2012, 2017). Ice thickness is deduced from the gravity‐derived bed elevation and a digital elevation model for 1981–1985 (Korsgaard et al., 2016). Ice fluxes are compared with the average SMB accumulated in each basin for the years 1961 to 1990 using the RACMO2.3 model downscaled at 1 km (Noël et al., 2016). We define glacier basins using flow direction from a surface velocity reference (Mouginot et al., 2017) (fast flowing ice) and from a reference topographic slope smoothed over 10 thickness from Howat et al., 2014 (slow moving ice). Our hypothesis is that the glaciers were near balance in the late 1980s, so ice front fluxes should be close to the reference RACMO2.3 balance fluxes (Rignot et al., 2008).

Ice front positions are digitized for the last 80 years from aerial photos and historical satellite images. We choose aerial photos from middle to late summer from 1930 to 1985 (Bjørk et al., 2012) and Landsat 1‐8 images from summer 1990 to summer 2017. We measure the average ice front position each year (Figures 3 and S5) and changes in position following the methods of Bjørk et al. (2012) and Jiskoot et al. (2012).

Figure 3

GIMP v2.1 (from BM3) surface elevation along profiles in Figure 2 with bed elevation from BM3 (dashed red), OMG bathymetry (dashed green), and bed elevation from this study (solid black). Ocean is blue, ice is light blue, and bed is light brown. Ice front positions are color coded from blue to red and labeled by year. Ocean Melting Greenland (OMG) temperature from conductivity‐temperature‐depth (CTD) casts in 2016 are color coded from blue (cold) to yellow (warm), with CTD position as a diamond. Limits of the gravity inversion are black triangles.

Ocean temperature and salinity were obtained by OMG in September 2016 in each fjord using two conductivity‐temperature‐depth (CTD) sensors: an AML Oceanographic Minos X CTD in thick brash/sea ice conditions and a Valeport Rapid CTD in ice‐free waters (OMG Mission, 2016). The measurements were taken at a sampling rate of 16 Hz. We use the Gibbs‐SeaWater Oceanographic Toolbox to convert the conductivity profiles into absolute salinity and pressure into depth. We average CTD profiles in 1 m bin to smooth spurious sensor readings.

Results

The inversion combines three data sets in a seamless fashion: (1) OMG in the fjords, (2) BM3 upstream, and (3) OIB gravity data at the ice‐ocean transitions. We find minimum densities for B1, B2, B4, B5, B6, and B9 (Figures 1 and S2) that are ranging from 2.5 to 2.85 g/cm3, which are consistent with the density of granodiorite and orthogneiss found in this region (Henriksen et al., 2009; Pechnig et al., 2005). After the inversion, the gravity misfit drops from 8.1 to 1.7 mGal (Table S1). Translated into meters, the inversion reduces the uncertainty in bed elevation from ∼140 to ∼30 m. Differences between BM3 and the gravity inversion are >100 m hence reflect net improvements in bed elevation mapping. For five glaciers (Ikertivaq N. N., M., S., KøgeBugt N., and Eqalorutsit), we observe larger differences in areas where bed mapping with MC did not extend to the coastline and was completed using a krigging interpolation. Finally, in areas with quality RS, the gravity results are in better agreement with the RS observations (Figure S6).

We compare our results with the 2‐D gravity‐based bathymetry from Boghosian et al. (2015). We find that the bed elevations from Boghosian et al. (2015) were on average 160 m higher. In several places, the frontal bed elevations from Boghosian et al. are hundreds of meters shallower than the OMG MBES data (Figure S7). This underestimation is explained by the lower resolution of the earlier gravity data (4.5 km for Boghosian et al., 2015 versus 750 m for our data) which smears the ocean floor topography in fjords only a few kilometers wide (Figures S7a, S7b, S7e, and S7g).

Our new bathymetry includes fjords with no prior depth information (Figure S8). In B1, we find a 500 m deep channel in front of Ikertivaq N. N. ( Figures 2a and 3) and a 800 m deep channel for Ikertivaq M. and Ikertivaq S. S. that extends offshore and inland. These deep troughs are not present in BM3. In B2, KøgeBugt C. has an 8 km wide, 600 m deep, relatively flat‐based fjord (Figures 2b and 3). In B3, our new bathymetry reveals a 800 m deep fjord in front of Gyldenløve (Figures 2c and S5). The gravity data extends the bathymetry mapping significantly in B4 (A. P. Bernstorff, Figure 2d), B6 (Mogens and Tingmiarmiut, Figure 2f), and B9 (Eqalorutsit and Qajuutaap, Figure 2i) for the first time (Figures 2 and S8). Conversely, in B5, B7 and B8, the OMG coverage already extended to the ice fronts.

Figure 2

Bed elevation from this study over (a) Ikertivaq (B1), (b) KøgeBugt (B2), (c) Gyldenløve and Graulv (B3), (d) A. P. Bernstorff (B4), (e) Rimfaxe and Skinfaxe (B5), (f) Tingmiarmiut and Mogens (B6), (g) Puisortoq (B7), (h) Anorituup (B8), and (i) Qajuuttap Se. and Eqalorutsit Killiit Se. (B9). Contour lines are shown at a 200 m interval. White is no data.

Inland, our new gravity‐derived bed elevation shows significant improvements compared to BM3. In B1, the trough below Ikertivaq N. N. extends 10 km farther upstream of the ice front than in BM3 (Figures 2a and S8). For all four glaciers, the troughs are on average 400 m deeper than BM3 but shallowing inland by ∼500 m over 2–3 km (Figures 3 and S5). In B2–B4, our bed elevation is deeper than BM3 for all glaciers (Figures 3 and S9). On Skinfaxe (B5), our bed elevation agrees with BM3 at the edges of the domain but is 200 m deeper at the center (Figures 2e and 3) and consistent with available RS data (Figure S6). In B6, the gravity‐derived bed elevation is 200 m deeper than BM3 in the ice front region for Tingmiarmiut glacier (Figure 3), but in agreement with RS data (9 ± 78 m in Figure S8). Beneath Mogens N in B6, we find a 12 km long, 700 m deep channel, 400 m deeper than BM3, rising to sea level between km 5 and 10 (Figures 2f and 3). In B8, the bed is shallower for Puisortoq N., but deeper for Puisortoq S., upstream of the ice front (Figures 2g and S5). In B9, the channels are 300–500 m deeper than in BM3 (Figures 2h and 3). Qajuuttap Se. and Eqalorutsit Killiit Se. are marine‐terminating glaciers instead of land terminating as shown in BM3 (Figures 2I, 3, and S5).

In terms of ice discharge, the revised fluxes are within 0.2 ± 0.5 Gt/yr of the 1960–1990 balance fluxes versus 0.8 ± 1.5 Gt/yr with BM3 (Table S1), hence reducing errors by a factor of 3. More importantly, the total ice discharge for all 20 glaciers is 64.7 Gt/yr, or +5% above the balance flux, versus 47.2 Gt/yr or 24% too low with BM3 (Table S1). The ice discharge are therefore in better agreement with the balance fluxes. Discrepancies are only found for Anorituup and Ikertivaq S., with a flux in excess of 1.7 Gt/yr and 1 Gt/yr, respectively. Explanations for the difference includes overestimation of the depth of the bed, the glaciers were not in a state of balance, or SMB is too low.

The pattern of ice front retreat varied significantly from fjord to fjord (Figure 3). Fast retreats took place along retrograde slopes and retreats slowed or stopped where the glaciers met prograde slopes. We find that 9 of the 14 glaciers for which we have complete coverage advanced on prograde slopes by up to 1 km during the period 1930–1940 (Figures 2 and S5). In contrast, KøgeBugt C., S. (B2) and Puisortoq N. (B7) retreated by 3, 1.5, and 1 km, respectively, on retrograde slopes (Figures 3 and S5). Between 1970 and 1990, 13 out of 20 glaciers remained stable or readvanced on prograde slopes (Figures 2, 3, and S5).

Going from north to south, in B1, Ikertivaq N. has been stable on a 1.5 km wide and 170 m deep sill for 80 years, while Ikertivaq N. N. started retreating in 1990 along a retrograde bed and retreated by 4 km until it met a 450 m deep sill in 2010. The glacier sped up by 1 km/yr or 20 %, while Ikertivaq N., M., and S. maintained steady speeds (Mouginot et al., 2017). In this sector, glaciers with no retreat in the last 10 years have been standing on a sill or had retreated to higher ground by 2005. Since the 1990s, KøgeBugt N., S. have been stable (Figures 2, 3, and S5). For KøgeBugt C., the largest retreat of 3 km took place from 1990 to 2015 along a retrograde slope (Figure 3). Rimfaxe and Skinfaxe have been the most stable glaciers in the region, with an ice front migration of less than 500 m since the 1930s, which we explain by the presence of sills only 350 m deep (Figure 3). Farther south in B6, we detect the largest retreat (8 km) on a retrograde bed for Mogens N. after 1965 (Figures 2 and 3). This glacier doubled its speed since the 1990s (Mouginot et al., 2017). Between 1972 and 1981, the ice fronts of Mogens C. and N. split in half (Figure 2) and retreated 1.4 km in 10 years. During 2003–2016, Mogens C. has remained stable on a 2.5 km wide and 1 km long sill. In 2016–2017, glaciers in B6 experienced a widespread retreat along retrograde slopes of 1.3 km (Mogens N.), 2 km (Mogens C.), and 0.9 km (Mogens S.) (Figures 3 and S5). We note three anomalies: Graulv in B3, Anorituup in B8, and Qajuutaap Se. in B9 that retreated less than 2 km since the 1930s despite deep fjords (Figures 3 and S5).

Using the 2016 CTD data, we find warm, salty AW widespread in the northern part of the study region (B1, B2, and B5) with potential temperatures greater than 3°C below 350 m depth. In the southern part, for example, in Mogens fjord, the warm water layer (>3°C) starts at 250 m depth. In Puisortoq N., Puisortoq S., and Anorituup fjords, nearly entire water column is warm. We find no major sill at fjord entrances, except for Puisortoq N. and Anorituup, where we detect sills less than 300 and 200 m deep, respectively (Figure S10). Ocean temperature is high in these fjords and reaches 5°C below −200 m (Figure S10). The sills, however, seem to limit the access of warm AW since water temperature is at least 1° lower upstream (Figure S10).

Discussion

These results shed new light on the spatial pattern of glacier retreat in southeast Greenland. The stability of 70% of the glaciers before 1990 is corroborated by two factors: (1) most calving fronts rested on prograde slopes and (2) sea surface temperature were stable in 1930–1940 and cooled down by 1°C in 1960–1970 (Bjørk et al., 2012). The anomalously fast retreats observed in B2 in 1930–1940 are consistent with the presence of retrograde slopes. In B1, Ikertivaq N. N. is the only glacier that retreated due to a retrograde bed in a fjord with warm AW. In this area, BM3 has ice fronts at sea level, which is not compatible with the observed changes. In B4, the stability of A. P. Bernstorff since 2005 is explained by a fjord depth of 300 m (versus 500 m in BM3), which means low exposure to AW based on mapped fjords nearby. In B3, Graulv has been stable for 80 years despite its ice front being grounded at 400 m deep. We have no CTD measurements near the ice front and no data between the OIB gravity inversion and the OMG campaign (green circle in Figure 1). We posit that a sill may exist in the fjord that limits the access of AW. Farther south in B6, BM3 is not compatible with the evolution of Mogens N. because it displays a shallow (200 m depth) prograde bed elevation not exposed to AW. The widespread retreat observed in 2016–2017 in B6 is, however, consistent with our mapping: the ice fronts were dislodged from their stabilizing sills (Figures 3 and S5) to retreat along retrograde beds.

We note three low retreat anomalies in B3, B8, and B9. In B8, the fjord has a 200 m deep sill that limits the access of warm AW to Anorituup (Figure S9). This glacier also calves at 8 m/d, which is 1 order of magnitude higher than typical ocean‐induced melt rates at ice fronts (Rignot et al., 2016). In B3 and B9, we have no CTD or bathymetry data outside of the OIB survey and cannot determine the cause of the stability. We speculate that a sill may be present in the unmapped part of the fjord that limits the access of warm AW to the ice.

Overall, except for two glaciers (Qajuuttap Se. and Graulv) for which we need additional data, the new bed map reveals pathways for warm AW to reach glacier fronts. For 17 out of 20 glaciers, we find bed elevations several hundreds of meters deeper than BM3, which makes the glaciers more vulnerable to AW (Figure S9). We have confidence in the results because the derived ice fluxes for 1980s to 1990s are in better agreement with the balance fluxes, our results match available RS data, and the misfit between observed and modeled gravity is lower (Table S1). The 20 glaciers that drain this sector of Greenland occupy 5% of the ice sheet in area, with a balance flux of 64 Gt/yr, or 17% of Greenland's ice mass flux in the 1980s to 1990s.

During the 1930–1990 period, the ice front positions were several kilometers downstream of the present positions for 16 out of 20 glaciers (Figure 1). Around the year 2000, 12 glaciers started a significant retreat, while 8 glaciers remained stable (Figure 1). Among the glaciers that retreated, the new CTD data show that more than 60% were exposed to warm AW greater than 3°C (Figure 1). For the remaining 40% without CTD data, the fjords are deep enough that exposure to AW is highly likely. Among the glaciers that have remained stable, >60% have ice fronts with low exposure to warm AW and/or standing above 300 m depth (Figure 1). For two glaciers (Graulv and Qajuuttaap Se.), we need additional oceanographic data to understand the stability of the glaciers (Figure 1).

The simultaneous retreat of six glaciers in summer 2017 (KøgeBugt C., Mogens N.‐C.‐S., and Puisortoq N.‐S.) is consistent with the retrograde bed upstream of the 2016 front positions. We expect that these glaciers will continue to retreat by several kilometers before reaching a more stable bed position (Figures 3 and S5). This new evolution provides additional confidence in our results because earlier maps did not indicate the presence of deep retrograde bed and could not explain this retreat.

Conclusions

We present a new bathymetry and bed mapping of southeast Greenland from a 3‐D inversion of OIB gravity data over 20 glaciers in 10 major fjords. Our multisensor approach provides the first reliable mapping of these inner fjords and terminal valleys where glaciers have been retreating during the last 80 years. The results reveal deep channels hundreds of meters deeper than estimated previously that provide natural pathways for AW to reach the glacier fronts. The uncertainty in bed elevation has dropped to 30–50 m, or only about 10% of the ice front thickness (versus >50% in BM3). Where we have reliable RS data, we confirm the precision of the inversion. In addition, the ice fluxes from the 1980s are much closer to the balance fluxes, as expected. We also find that the pattern of retreat or advance of nearly all glaciers is consistent with three major factors: (1) the presence of warm AW in the fjords, (2) the existence of retrograde bed, and (3) the presence of stabilizing sills or prograde slopes. These results provide multiple, independent lines of evidence of the reliability of the gravity‐based bed mapping. We recommend that other glaciers be mapped using a similar method to obtain a reliable glacier and fjord depths. Surveys should extend from the glacier fronts to the mouth of the fjords to fully understand the connection with the surrounding ocean.

Supporting Information S1

Click here for additional data file.

Table S1

Click here for additional data file.