Grounding Line Retreat of Denman Glacier, East Antarctica, Measured With COSMO‐SkyMed Radar Interferometry Data

Denman Glacier, East Antarctica, holds an ice volume equivalent to a 1.5 m rise in global sea level. Using satellite radar interferometry from the COSMO‐SkyMed constellation, we detect a 5.4 ± 0.3 km grounding line retreat between 1996 and 2017–2018. A novel reconstruction of the glacier bed topography indicates that the retreat proceeds on the western flank along a previously unknown 5 km wide, 1,800 m deep trough, deepening to 3,400 m below sea level. On the eastern flank, the grounding line is stabilized by a 10 km wide ridge. At tidal frequencies, the grounding line extends over a several kilometer‐wide grounding zone, enabling warm ocean water to melt ice at critical locations for glacier stability. If warm, modified Circumpolar Deep Water reaches the sub‐ice‐shelf cavity and continues to melt ice at a rate exceeding balance conditions, the potential exists for Denman Glacier to retreat irreversibly into the deepest, marine‐based basin in Antarctica.


Introduction
Over the last decades, the Antarctic Ice Sheet has been contributing to global sea level rise, and the marine-based Wilkes Land sector of East Antarctica has been an important contributor to the total mass loss . Denman Glacier (DG), one of the largest glaciers in East Antarctica in terms of ice discharge, holds an ice volume equivalent to a 1.5 m global sea level rise . Radar altimetry data from the European Envisat mission indicated that DG thinned at 0.4 m/year between 2002 and 2010 upstream of its grounding line (Flament & Rémy, 2012), where the glacier speed reaches 1.4 km/year (Young et al., 1989). DG has experienced a cumulative mass loss of 268 ± 19 Gt between 1979 and 2017, or 7.0 ± 0.5 Gt/year, or 13.2 ± 1% of its balance flux of 52.6 ± 4.2 Gt/year  ( Figure S1 in the supporting information). The 24,000 km 2 floating extension of the glacier, which includes Shackleton Ice Shelf and Denman Ice Tongue (DIT), has been melting at an area-average rate of 3.1 ± 0.7 m/year, which is above average among other ice shelves in East Antarctica .
The first and only delineation of the grounding line of DG was derived using Earth Remote Sensing (ERS-1/2) satellite data from year 1996 (Rignot, 2002). The bathymetry in front of and beneath DIT is not well known. Seismic traverses and airborne surveys over land have revealed the presence of a deep trough, more than 2,000 m below sea level, 100 km upstream of the 1996 grounding line (Young et al., 1989). The trough connects with the deep marine Aurora Subglacial Basin, about 400 km to the southeast (Young et al., 2011). There is no detail about the depth of DG near its grounding line.
Here, we present a new delineation of DG grounding line using a time series of Differential Interferometric Synthetic Aperture Radar (DInSAR) data from years 2016 to 2018. To quantify the thinning rate of the glacier and the melt rate of its floating ice shelf, we use a time series of high resolution Digital Elevation Models (DEMs) of the glacier surface topography from the TanDEM-X (TDX) radar mission and laser altimetry data from NASA's Operation IceBridge (OIB) mission. We present a new bed topography of the glacier that combines OIB ice thickness and vector ice velocity data on land. The map is extended offshore using a three-dimensional inversion of gravity data to reveal the bathymetry in front and beneath the ice shelf. We evaluate the short-term variability in grounding line position using the 2016-2018 data and its long-term variability in reference to the ERS-1/2 data from 1996. We compare the results with the newly reconstructed bed topography and bathymetry to conclude on the recent and future evolution of DG.

Grounding Line Mapping
To delineate the grounding line of DG, we form SAR interferograms with a 1-day time interval by coherently combining SAR acquisitions from the ERS-1/2 satellites in year 1996 and from the Agenzia Spaziale Italiana's COSMO-SkyMed (CSK) constellation in years 2016, 2017, and 2018 (Table S1). We multilook each CSK SAR interferogram using an 8 × 8 averaging window in range and azimuth, equivalent to 24 m × 24 m on the ground and use a TDX DEM at 30 m spacing to remove the topographic signal (Milillo et al., 2017;Rignot, 2002). We form 24 CSK differential pairs (DInSAR) by differencing consecutive 1-day CSK SAR interferograms ( Figure S2). The differencing eliminates the horizontal motion of the glacier to reveal its differential motion, which includes the vertical motion of floating ice in response to changes in oceanic tides (Brunt et al., 2011;Rignot et al., 2011), changes in surface elevation caused by subglacial lake drainage or filling on land (Gray et al., 2005), and residual data noise (e.g., turbulent water vapor and thermal noise). We reprocess the two ERS DInSAR interferograms from year 1996 improving the interferometric baseline determination using OIB altimetry data and the topographic correction using a TDX DEM.
The differential signal associated with tidal motion is characterized by a step function immediately seaward of the grounding line, over a region typically 5-9 km wide called the flexure zone, where ice adjusts to hydrostatic equilibrium (Brunt et al., 2011). In the flexure zone, the deformation signal results from a quadruple difference in tidal state between the four epochs used to form the DInSAR interferogram (Rignot et al., 2011). The small amplitude of the step function makes arduous the delineation of the grounding line. Therefore, we use only DInSAR data with a strong differential tidal signal, that is, typically more than 10 fringes (or 360 • variation in interferometric phase) across the flexure zone ( Figure S2). The grounding line is identified manually at the inward limit of the vertical displacement zone of the ice surface. The most upstream interferometric fringe represents the farthest point inland where seawater intruded beneath the glacier among the four different tidal epochs and modified its ice surface elevation ( Figure S2) (Milillo et al., 2017;Rignot, 2002).
The tidal signal is easily separated from the signal associated with the drainage/filling of subglacial lakes, which has a circular to oval shape, a different signal amplitude than the tidal signal, and occurs upstream of the grounding line. Conversely, the number of fringes between the grounded ice (inland area with no residual fringes) and freely floating ice (seaward of the zone of tidal flexure) must always be the same in each DInSAR interferogram, that is, equal to the difference in tidal height between the four different epochs projected along the radar line of sight. This eliminates the possibility of placing the grounding line on the wrong set of fringes. This guideline is crucial to correctly map the grounding line across the entire glacier since the intrusion of seawater beneath grounded ice varies spatially across the grounding line as a function of tidal state, bed slope, and ice thickness. In some places, for example, the western flank of DG, the tidal deformation extends farther inland than on average. The fringe spacing in the flexure zone increases relative to the average fringe spacing elsewhere, yet the number of tidal fringes across the flexure zone remains the same, which helps reduce the risk of misidentifying the grounding line. For each DInSAR data, we digitize grounding line segments with an uncertainty in grounding line position based on the amplitude of the tidal signal and phase noise.
The CSK data do not provide a complete coverage of the grounding line of Scott and Northcliffe glaciers. For those glaciers, we complement the CSK data with a 6-day DInSAR interferogram from the European Space Agency (ESA) Sentinel-1A/B (S1-A/B) in year 2017 processed as reported in Scheuchl et al. (2016) 10.1029/2019GL086291 ( Figure S3). However, the S1-A/B DInSAR data cannot be used to detect the grounding line of DG since the repeat cycle of S1-A/B is too long to maintain phase coherence over DG.
To reconstruct the tidal state at the time of passage of the satellites (i.e., 4 epochs for ERS-1/2, 4 epochs for CSK), we use the CircumAntarctic Tide Simulation (CATS2008A) model at a center location of 64 • 43.2831 ′ S, 100 • 0.5712 ′ E, just off the DIT, in the open ocean (Padman et al., 2002). We assume that the floating ice responds elastically to changes in oceanic tides. To evaluate the tidal model, we compare the differential tidal signal measured from the DInSAR data (i.e., step function in height between grounded and freely floating ice) with the tidal differential height estimated from the tidal model ( Figure S4). For each DInSAR data, we calculate the maximum amplitude of the tidal signal among the four tidal states used to form the DInSAR interferogram, or h max .

Ice Surface and Bed Elevation
We generate a time series of DEMs at 12 m spacing, with meter-scale vertical precision, using bistatic SAR acquisitions from the TDX mission and a standard processing chain described in Rizzoli et al. (2017). We use data from the Ice Cloud and land Elevation Satellite (ICESat) mission, (GLAH14, Version 34; Zwally et al., 2014) and OIB (i.e., Riegl laser altimeter; Blankenship et al., 2012) to calibrate the TDX DEMs absolute height over rocky outcrops (Milillo et al., 2019). We estimate the height precision and absolute height accuracy of the TDX DEMs to be better than 2 m (Rizzoli et al., 2017).
We employ a new reconstruction of DG bed topography from BedMachine Antarctica using a mass conservation approach, with a grid spacing of 500 m and a nominal accuracy of 100 m . The mass conservation method over the study region combines ice surface velocity vectors from years 2007-2008 ( Mouginot et al., 2017), ice thickness data from the HiCARS sounder (Hi-CApability Radar Sounder) (Young et al., 2011), surface elevation from the TDX DEM from year 2017, and surface mass balance averaged for the years 1979 to 2008 from Regional Atmospheric Climate Model Version 2.3p1 (RACMO2.3p1) van Wessem et al., 2014). In this region, RACMO2.3p1 is superior to RACMO2.3p2 (Mohajerani et al., 2018(Mohajerani et al., , 2019. Temporal changes in ice thickness (≤1 m/year) and basal melt beneath grounded ice (≤cm/year) are negligible. Full radar sounder-derived thickness crossings of DG are only available 120 km upstream of the grounding line.

Elevation Changes and Ice Shelf Melt
We quantify changes in surface elevation by differencing TDX DEMs separated by 2 to 3 years using both an Eulerian (i.e., fixed reference grid) and a Lagrangian (i.e., reference grid moving with the ice) approach. These methods provide distinct information over the grounded and floating portion of DG. In an Eulerian framework, changes in surface elevation over grounded ice reflect variation in ice thickness, H, caused by ice dynamics and snowfall accumulation since the elevation of the ice base does not change (Shean et al., 2019). Under the assumption of incompressibility and assuming a column-average velocity v, mass conservation applied to an ice column of thickness a is the surface accumulation rate, and . b is the basal melt rate. Due to the advection of heterogeneities in ice thickness, this methodology is not useful on floating ice. To calibrate dH∕dt, we use ice-free areas near Burger Hills where changes in surface elevation should be zero. We smooth the results using a 15 × 15 pixels (420 m × 420 m) moving window and grid the results at 1.5 km spacing. We estimate a mean error of 0.05 m/year and a standard deviation less than 0.1 m/year. On floating ice, temporal changes in surface elevation obtained in a Lagrangian framework (Moholdt et al., 2014;Shean et al., 2019) yield estimates of ice shelf melt rates assuming a uniform rate of horizontal motion and ice to be in hydrostatic equilibrium as Dh∕dt where Dh∕dt is the material derivative, ∇ the gradient operator, v the vector velocity, . a the surface mass balance, m b the ocean-induced ice shelf melt rate with positive (negative) values indicating melting (freezing), i = 0.917 g/cm 3 is the density of ice, and sw = 1.028 g/cm 3 is the density of sea water. We use a chip size of 30 × 30 pixels (840 m × 840 m) with step size of 16 pixels (448 m) in both directions and smooth the results with a 15 × 15 pixels (420 m × 420 m) moving average window to obtain melt rates at 1.5 km spacing. We follow the trajectory of each pixel (on a monthly basis) using the velocity map and correct for tidal motion at the time of the satellite overflights using the CATS2008A model (Table S2).
The main sources of uncertainty of the estimation of ice shelf melt are the measurement errors (e.g., DEM, ice velocity), correction of ocean tides, changes in surface elevation due to snowfall accumulation, and ice Geophysical Research Letters 10.1029/2019GL086291 dynamic thinning. We assume these error sources to be unbiased. We assume uncorrelated error sources for the TDX surface elevation (2 m) (Rizzoli et al., 2017) and yearly-averaged ice velocity maps (2-3 m/year) . The error in tide correction is 6 cm. We use a depth-averaged density for ice and ocean water, both assumed spatially and temporally constant. We extract the surface component of the transient melt from RACMO2.3p1 model (van Wessem et al., 2014), which is less than 1 m/year, yielding an uncertainty of few cm per year. These results must be interpreted with caution over grounded ice and in the transition zone between grounded and floating ice wherein the assumption of hydrostatic equilibrium is not strictly valid due to bending stresses of ice and adjustments of ice to floatation over several km. We correct our estimates of ice shelf melt rate for dynamic thinning (0.4 m/year) measured on grounded ice, assuming that the same dynamic thinning applies to floating ice. The corresponding error is 0.1 m/year in height or 1 m/year in ice thickness. Considering all the error sources, the total uncertainty of the melt rate estimates is 4 m/year.

Oceanographic data
We use hydrographic data collected by instrumented marine mammals in 2004-2005 and 2011-2015 distributed by the Marine Mammals Exploring the Oceans Pole to Pole (MEOP) consortium (Treasure et al., 2017). Each mammal is equipped with a Conductivity Temperature Density Satellite Relay Data Logger (CTD-SRDL). Before deployment, each sensor is calibrated with a precision of 0.005 • C and 0.005 mS/cm, respectively, for temperature and conductivity. Tags recovered after deployments (i.e., 72-136 days) show no change in calibration within the manufacturer's specification (Roquet et al., 2014;Treasure et al., 2017). Each time the mammal surfaces, the data collected throughout a dive are processed by the CTD-SRDL sensors and telemetered through the Argo satellite system. The geolocation of each dive is determined by the Argo satellite triangulation with an accuracy of 4 km. To ensure data quality, the hydrographic data are postprocessed as in Roquet et al. (2014). We find many CTD casts near the edge of the continental shelf but few in the vicinity of DIT (Figures S5 and S6; Table S3).

Bathymetry
The International Bathymetric Chart of the Southern Ocean (IBCSO) (Arndt et al., 2013) includes no depth data on the continental shelf and beneath DIT. The water column beneath DIT was arbitrarily set at 40 m. The ice shelf thickness was inferred from a DEM of Antarctica assuming ice in hydrostatic equilibrium. Here, we employ a three-dimensional inversion of airborne gravity data to infer the bathymetry beneath DIT using constraints from BedMachine Antarctica on land , CTD maximum depth data at sea, islands, and the presence of ice rises (zero water column) ( Figure S6). To invert the gravity data, we use the Geosoft GM-SYS 3-D software, which uses the Parker method (Parker, 1973). This method calculates the gravity anomaly caused by an uneven, uniform layer of material by means of a series of Fourier transforms and modifies it to minimize the misfit between calculated and observed gravity. We model the domain as three horizontal layers: (1) a solid ice layer with a density of 0.917 g/cm 3 ; (2) an ocean water layer with a density of 1.028 g/cm 3 ; and (3) a rock layer with a uniform density of 2.67 g/cm 3 following An et al. (2019). Airborne gravity data are from OIB, Investigating the Cryospheric Evolution of the Central Antarctic Plate (ICECAP) (Young et al., 2011) and from the Gravity and Geoid in Antarctica (AntGG) international compilation (Scheinert et al., 2016) (Figure S7).
We run a forward model of the gravity field using an initial bed solution. In areas with bed elevation constraints, we calculate the DC shift (Direct Current) or mean difference between modeled and observed gravity. We then interpolate the DC shift over the entire domain using a minimum curvature algorithm . The DC shift captures the long wavelength (100 km) variations in geology below the model domain and removes it from the observations. As a result of the inversion, the misfit between modeled and observed gravity is reduced from 15.5 mGal with the initial bed solution to 5.8 mGal with the final bed solution. We estimate that the nominal precision of the bathymetry to be 150 to 200 m, but local errors larger than 200 m cannot be excluded .

Results
The availability of multiple grounding line measurements with different tidal states on the same calendar year allows to identify the grounding zone (Figures 1, 2, and S8) of DG, that is, the region where the grounding line migrates back and forth with changes in oceanic tides (Milillo et al., 2017;Sayag & Worster, 2011) and also to determine the average rate of grounding line retreat with more confidence. The grounding line  (Mazloff et al., 2010) at 310 m depth color-coded between −2 to +2 • C. Color-coded rectangular markers represent MEOP-CTD in situ temperatures (Roquet et al., 2014) at the maximum depth of the cast. (c) At sea, bathymetry from a three-dimensional gravity inversion (this study) and IBCSO (Arndt et al., 2013)   On neighboring glaciers, we find no retreat on Scott Glacier, a 2 to 3 km retreat on Reid Glacier, and no conclusion possible for Northcliffe Glacier, which is incompletely covered by CSK data.
The reconstructed bed topography ( Figure 1c) reveals a 10 km wide subglacial ridge, 400 to 700 m below sea level, 4 km upstream of the present-day grounding line on the eastern flank, and a 5 km wide, 1,800 m below sea level trough on the west flank, with no ridge. The western trough deepens for another 18 km in the inland direction to 1,900 m below sea level, before widening to a total width of 12 km and deepening further to 3,400 m below sea level over the next 30 km inland, that is, a steep retrograde bed slope of 5%. This maximum depth makes this region one the deepest locations in Antarctica (Fretwell et al., 2013;Young et al., 2011). Beyond 50 km inland, the bed elevation rises but remains well below sea level (Figure 1).
From the Eulerian approach, we find that the main trunk of DG thinned at 0.4 ± 0.1 m/year between 2011 and 2013, with areas surrounding the glacier thickening at 0.35 ± 0.1 m/year (Figure 3c). This magnitude thinning is consistent with prior altimetry records (Flament & Rémy, 2012;McMillan et al., 2018) where thinning was also confined in the areas of fast flow. This signal suggests that ice flows faster than the speed required to maintain a state of mass balance. Indeed, DG has been speeding up by 16 ± 2% since the 1970s and is estimated to be 10 ± 2% out of balance at present . Most glacier thinning is confined in the last 50-100 km of grounded ice, coincident with the area of speed up (Figure 3b). On DIT, the speed up was larger, 33% between 1957 and 1996, and another 10% until 2016 (Dolgushin, 1965;Rignot et al., 2019), implying that the ice shelf is thinning dynamically as well.  the methodology is affected by uncertainties in ice floatation (Wilson et al., 2017), but our results indicate melt rates in the range of 10-20 m/year. Such ice shelf melt rates near the grounding line are high compared to other East Antarctic ice shelves , but comparable in magnitude to those observed in warm parts of West Antarctica, for example, the Amundsen Sea Embayment (ASE). In the grounding zone, a narrow region that is not always in contact with seawater, there are no prior estimates of ice shelf melt rates, but the rates we estimate are significant.
The 45 m/year melt on the ice shelf requires the presence of relatively warm water, typically in the range of 4 to 5 • C above the in situ freezing point of seawater (Rignot & Jacobs, 2002). The five MEOP casts from year 2011 near DIT reveal in situ ocean temperatures ranging from −1.31 to −0.26 • C, at depths of 550 to 850 m, which correspond to a thermal forcing (difference between the in situ temperature and the depth-dependent, salinity-dependent freezing point of seawater) of 2.6 to 3.2 • C at 1,000 m depth (Fujino et al., 1974).
The new reconstructed bathymetry reveals a sinuous, deep trough beneath DIT intersected by three ridges at about 500 m depth (Figure 1f). The first ridge, approximately 60 km north of the grounding line, yields a water column thickness (difference between the seafloor depth and the depth of the ice shelf draft) less than 200 m. At this location, the DInSAR data do not suggest the presence of any ephemeral grounding points, that is, locations at which the ice shelf touches the sea floor at low tides (Schmeltz et al., 2002). Therefore, it is likely that warm water, possibly of Circumpolar Deep origin, circulates freely beneath the shelf (more than 400 m deep) and reaches DG grounding line. The two ridges farther north may also affect the heat transfer to the cavity, but the water column in that location is thicker since the ice shelf is thinner. To quantify the heat transfer, it would be essential to better constrain the depth of the ridge and obtain observations of the water structure in the ice shelf cavity.

Discussion
Our analysis of the DInSAR data spanning from 1996 to 2018 indicates that the grounding line of DG has been retreating asymmetrically, with a maximum retreat on the western flank along a previously unknown deep trough with a retrograde slope, and almost no retreat on the eastern flank, which is blocked by a prominent ridge. The 5 to 6 km retreat in 21 years along the western flank is consistent with the presence of a retrograde bed and high ice shelf melt rates along a grounding line more than 2 km below sea level. To the east, Scott Glacier is thinning upstream of its grounding line, but the grounding line is stabilized on a ridge. We do not have enough data on Northcliffe Glacier to quantify its retreat, but its grounding zone was exceptionally wide (≥10 km) in 1996, which suggests the presence of a nearly flat bed at the grounding line (Sayag & Worster, 2011). Since the glacier does not extend below sea level far inland, it does not hold much potential for sea level rise. Similarly, Reid Glacier is retreating but does not drain a large marine-based sector.
On the west flank of DG, the grounding zone is several km wide, that is, seawater infiltrate beneath grounded ice over considerable distances at high tide (Brunt et al., 2011). This wide grounding zone cannot be explained by the interaction of solid ice with a rigid bed but can be explained by the presence of a viscoelastic bed as in Sayag and Worster (2011). One important consequence is that ice melts in contact with ocean waters over considerable distances in a region that exerts a critical control on glacier stability (Sayag & Worster, 2013). If warm water present on the continental shelf reaches the deeper part of DG sub-ice-shelf cavity, the thermal forcing increases to 4.2 • C, in principle sufficient to justify an average melt rate of 45 m/year (Rignot & Jacobs, 2002). Vigorous ocean-induced ice melt within the grounding zone will directly reduce basal resistance to flow and entrain glacier speed up and further retreat. Such type of ice-ocean interaction, which is not included in current ice sheet models (Nowicki et al., 2016;Seroussi & Morlighem, 2018), could significantly increase projections of sea level rise from this region.
To the best of our knowledge, the bed slope configuration near the grounding line of DG is unique in this sector of East Antarctica. Totten Glacier, which is retreating at a comparable rate, must retreat for 50 km along prograde bed slopes before reaching retrograde bed slopes (Li et al., 2015). Moscow University Glacier is also flanked by prograde slopes that should provide more stability .
The grounding line retreat of DG is one order of magnitude smaller than that observed in the ASE (Scheuchl et al., 2016), However, the sector containing DG holds an ice volume equivalent larger than for the ASE sector (1.2 m); hence, any signal of retreat along DG is of considerable importance for the future. The new bed topography data suggest the presence of a direct retrograde pathway for the glacier, that is, no bump in bed topography above the current grounding line that could slow down the retreat. One factor that could slow the retreat is the narrowness of the deep channel, which may provide significant lateral drag, similar to Jakobshavn Isbrae in Greenland which retreats at 0.5 km/year (Bondzio et al., 2017;Habermann et al., 2013;Thomas, 2004). A second factor would be the efficiency of ocean heat transfer to the DG grounding line. In the future, it will be critical to collect oceanographic data in that region and continue observations of the grounding line dynamics, which is currently only possible with the CSK constellation.

Conclusions
A detailed mapping of the grounding line of Denman Glacier, East Antarctica, with satellite radar interferometry, reveals a complex and asymmetric pattern of retreat. By comparing satellite data acquired 21 years apart, we detect a 5.4 ± 0.3 km retreat on the western flank, proceeding along a previously unknown, smooth, retrograde bed dropping to more than 3,400 m depth over the next 50 km. The retreat is negligible on the eastern part of the glacier, which is protected by a subglacial ridge. Existing ocean data and observations of ice shelf melt rates suggest the presence of vigorous ice-ocean interaction beneath the Denman Ice Tongue, which must play a central role in the current evolution of the glacier, similar to what is observed in the ASE sector of West Antarctica. We conclude that the potential exists for a rapid retreat along the Denman trough in the future. A strengthening of the westerlies (Spence et al., 2014), caused by rapid midlatitude climate warming and cooling of the stratosphere by the ozone depletion (Swart et al., 2018), may bring more warm Circumpolar Deep Water toward the coast of Antarctica and explain the present-day evolution of this sector. These observations challenge the view of glacier stability in East Antarctica.