Ocean forced variability of Totten Glacier mass loss

Abstract A large volume of the East Antarctic Ice Sheet drains through the Totten Glacier (TG) and is thought to be a potential source of substantial global sea-level rise over the coming centuries. We show that the surface velocity and height of the floating part of the TG, which buttresses the grounded component, have varied substantially over two decades (1989–2011), with variations in surface height strongly anti-correlated with simulated basal melt rates (r = 0.70, p < 0.05). Coupled glacier–ice shelf simulations confirm that ice flow and thickness respond to both basal melting of the ice shelf and grounding on bed obstacles. We conclude the observed variability of the TG is primarily ocean-driven. Ocean warming in this region will lead to enhanced ice-sheet dynamism and loss of upstream grounded ice.

The East Antarctic Ice Sheet (EAIS) is the world's largest potential source of sea-level rise (53.3 m; an order of magnitude greater than either West Antarctica or Greenland), with the marine-based component (i.e. where the ice sheet is grounded below sea level) containing enough ice to raise sea levels by 19.2 m (Fretwell et al. 2013). The Totten Glacier (TG: see Fig. 1c for the location) has the largest outflow (70 ± 4 km 3 a −1 ) of any glacier in East Antarctica, and has the third highest ice flux of all Antarctic glaciers, behind only the Pine Island and Thwaites glaciers (Rignot & Thomas 2002), both in the Amundsen Sea (AS) region of West Antarctica. The Aurora Subglacial Basin is an extensive region of grounded ice that drains primarily through the main trunk of the TG (Roberts et al. 2011;Young et al. 2011;Wright et al. 2012) and contains enough marine-based ice to raise the global sea level by 3.5 m (Greenbaum et al. 2015). The TG has an extensive embayed floating ice shelf adjoining a weakly grounded 'ice plain' containing topographical features  that have previously been interpreted as 'ice rumples', with the ice surface slightly above the flotation limit (Greenbaum et al. 2015). It is thought that, by analogy with observed changes in the AS region of West Antarctica (Shepherd & Wingham 2007;Pritchard et al. 2009Pritchard et al. , 2012, these features make the TG susceptible to rapid grounding zone retreat via the marine ice-sheet instability mechanism (Schoof 2007) and potentially sensitive to changing ocean conditions (Li et al. 2016). Indeed, evidence of past large-scale retreat cycles has been inferred from erosion patterns beneath the ice sheet (Aitken et al. 2016), potentially driven by changing ocean conditions. Analysis of satellite laser-altimeter measurements (ICESat), of relatively short-duration (2003-08), showed a mass-loss trend for the TG, with surface lowering in some regions in excess of 1.9 m a −1 (Pritchard et al. 2009). Mass balance estimates from accumulation over the catchment basin, and velocities at the grounding line, show the TG catchment was losing ice at a rate of approximately 8 Gt a −1 over the period 1992-2006 (Rignot & Thomas 2002;Rignot et al. 2008), despite a period of anomalously high regional snowfall (van Ommen & Morgan 2010). Furthermore, analysis of satellite gravity data from the GRACE mission indicated accelerated mass loss of the grounded region from slightly over 2 Gt a −1 in mid-2002 to nearly 18 Gt a −1 in 2009 (Chen et al. 2009). Recent analysis of a longer time series (1994)(1995)(1996)(1997)(1998)(1999)(2000)(2001)(2002)(2003)(2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012) of surface-height data for the ice shelf revealed large interannual variability in elevation over the floating portion of the TG (Totten Ice Shelf; TIS), with an overall slight negative trend over the full time interval that becomes more pronounced during the ICESat period (Paolo et al. 2015).
Here we present further evidence showing that the variability of the TG is apparent in observations of both surface ice velocities and surface height. Numerical ocean modelling indicates that interannual differences in water mass properties in the region can lead to large variations in basal melting of the TIS (Khazendar et al. 2013;Gwyther et al. 2014), suggesting that both TIS elevations and ice dynamics may be modulated by variability in ocean forcing (Li et al. 2016). Thus, the aim of this paper is to provide a mechanism which explains that the variability of the TIS is in response to changes in ocean forcing. This aim will be addressed using both satellite observations and modelling.
Velocities were obtained from fast Fourier transform (FFT)-based correlations derived from surface feature displacements from Landsat4 or Landsat7 image pairs. The images are 'Systematic Terrain Correction' (Level 1GT) products in standard polar stereographic projection (reference latitude 71°S). Level 1GT images are corrected for terrain distortions using the Radarsat Antarctic Mapping Project (RAMP) digital elevation model (DEM), which may introduce potentially large errors. However, as our area of interest is an ice shelf of relatively uniform surface height, we expect any such effects to be small. To further minimize these effects, co-registration of the image pairs was based on persistent features near the edge of the ice sheet, where the ice is thin and surface features are a reflection of the underlying bedrock features. These regions tend to be at broadly similar elevations to the ice shelf.
To account for Scan Line Correction off (SLCoff ) problems in 2003-present images, the feature tracking software (Scambos et al. 1992) was modified to fill these gaps with random data selected from the intensity histogram for the valid region of the image (Warner & Roberts 2013). To reduce variability, displacements were calculated for 1920 × 1920 m reference 'chips' on a 300 m grid and then binned onto a 1 × 1 km grid based on a least trimmed squares (Rousseeuw & Hubert 1997), requiring at least three data points per bin. To further reduce noise, we only retained data where at least four of the nine neighbouring points had both displacement components within 150 m of each other.

Surface heights
We constructed an 18 year record of ice-shelf height using data from three European Space Agency (ESA) satellite radar altimeter (RA) missions: the European Remote Sensing Satellite-1 and Satellite-2 (ERS-1, 1991-96;andERS-2, 1995-2011), and the Environmental Satellite (Envisat, 2002-12). We obtained Level-2 RA data as Version 5 Ice Data Records (IDRs) for each mission from the NASA/ GSFC (Goddard Space Flight Center) Ice Altimetry Group (http://icesat4.gsfc.nasa.gov/). At GSFC, the RA waveforms were retracked using a rangeretracking algorithm that fits a multi-parameter function to each one (the β-retracker) (Zwally & Brenner 2001); one of these parameters is the retracking point, which is the mid-point on the waveform leading edge. The following corrections were applied by GSFC: atmospheric range corrections; instrument corrections; slope corrections; ocean and solid Earth tides (Brenner et al. 1983;Zwally & Brenner 2001;Zwally et al. 2005); (for ERS) removal of a 0.41 m bias from ERS-1 heights to account for a change in instrument parameter used for ERS-2 (Femenias 1996); corrections for drifts in the ultrastable oscillator and bias changes in the scanning point target response that are obtained from ESA; and upgraded orbits (DGM-E04 orbits for ERS), which have a radial orbit precision of 0.05-0.06 m (Scharroo & Visser 1998).
Our determination of height changes over the ice shelves from multi-mission satellite RA data is based on 'crossover analysis' (Zwally et al. 1989(Zwally et al. , 2005Davis & Ferguson 2004;Wingham et al. 2009), which estimates change in surface height at intersections between time-separated ascending and descending satellite tracks. The following processing steps were undertaken (Paolo et al. 2016): subsetting data over floating ice; data editing and additional corrections required for floating ice and for radar signal interactions with the surface layer on the ice shelves; crossover analysis; methods for merging data records from the different missions; and the averaging scheme required to achieve satisfactory data coverage and accuracy. The output is 18 year-long records of ice-shelf surface height, which we then analysed using statistical techniques to derive long-term trends, acceleration and uncertainties. Ice sheet/shelf modelling Ice sheet and shelf modelling was undertaken (Donnelly 2014) with a drainage basin regional model using the BISICLES model (Cornford et al. 2013(Cornford et al. , 2015. The 'full Stokes' equations describing ice dynamics are approximated with a so-called higherorder approximation, with vertically integrated treatment of momentum for computational efficiency. Horizontal grid resolution varies (1-16 km) throughout the computational domain as a function of the ice speed and proximity to the grounding line. The geometry is based on the BEDMAP2 (Fretwell et al. 2013) compilation, and an initial thermal state is specified (Pattyn 2010). The applied pattern of ocean-driven melt at the base of the ice shelf is diagnosed from observed ice-flux divergence, with the temporal modulation of melt specified as sinusoidal with a period of 10 years and an amplitude of 0.5 of the temporal average (Donnelly 2014).

Ocean modelling and basal melt
We simulated the ice shelf-ocean interaction with the Regional Ocean Modeling System (ROMS) ( (Tamura et al. 2008), which is described fully in Gwyther et al. (2014). This model is spun-up by 16 years of repeated 1992 forcing, at which point the ocean heat content has approximately plateaued. The other model run simulated 1992-2012 conditions with the same finite-difference code, but employed a more modern ECCO2 version, the reanalysis product ERA-Interim (Dee et al. 2011), for surface forcing and Fig. 2. Change in speed of the TG ice shelf. Difference in speed between image pairs: (a) (28 March 1989, 25 December 1989) and (1999, 2002; (b) (1999,2002) and (2008,2010); and (c) (2008, 2010) and (2010, 2011). Also shown is the ASAID grounding line (Bindschadler et al. 2011a). updated surface sea-ice production rates. This model was run with two periods of 1992-2012 forcing (42 years in total), with the first period spinning the model up and the second period being used for analysis.
For comparison with the ocean-model solution, basal melt for the TIS was also calculated from icepenetrating radar thickness profiles (Young et al. 2011) and surface ice velocities (Rignot et al. 2011). Both datasets were averaged onto the same 1 × 1 km polar stereographic grid. Local mass flux estimates from coincident thickness and velocity data points were advected downstream until reaching the next local mass flux point using a modified version of a the Lagrangian streamline code (Roberts et al. 2011), accounting for across-stream convergence/divergence. The mass flux difference between these points was assigned uniformly along the streamline as the average change in mass flux along the streamline. All such estimates for each grid cell where averaged and a 15 km Gaussian filter applied. This represents an estimate of the local basal melt rate, in general an underestimate as surface snowfall will have contributed to the along-streamline mass flux changes.

Results
We analysed a sequence of ice-surface velocity measurements from automated correlation between pairs of co-registered visible-band light satellite images (Scambos et al. 1992;Warner & Roberts 2013) over the period 1989-2012. The results show temporal variability in the flow of the TIS during these 24 years, in agreement with mass-budget estimates (Li et al. 2016). Between 1989 and 2001, the TIS slowed by over 7%, followed by a 4% speed-up between 2002 and 2009 and another, more tightly constrained, slowdown during 2010-11 of similar magnitude (Fig. 1a). While the lack of suitable satellite imagery leads to sparse temporal coverage, our velocity measurements are of sufficient quality to capture interannual variability (Fig. 2), as it is significantly larger than the uncertainty in each velocity measurement. Ice-shelf speeds in the vicinity of the rumple near the grounding line are in the range of 900-1000 m a −1 (Rignot et al. 2011). The observed variability of velocity appears linked with an 18 year (1994-2011) record of surface-height change, averaged over the central TIS area (Paolo et al. 2015), with an overall trend not significantly different from zero (Fig. 1b). An approximately 6 year period (2004-09) is characterized by a relatively lower elevation corresponding to the period of faster surface velocities.
Hindcast area-averaged TIS basal melt rates from an ocean model (Gwyther et al. 2014), which shows high variability over a 20 year period , are anti-correlated with surface-height variations (r = 0.70 with a 95% confidence interval of 0.69-0.71, calculated using a bootstrap method taking into account signal auto-correlation: Ólafsdóttir & Mudelsee 2014). Increased melting leads to a lowering of surface height and vice versa (Fig. 1b), driven primarily by changes in ocean temperature. This result was supported by a second ocean-model simulation, although of a shorter 16 year period, using a different source of hindcast forcing (not shown). The average melt rates for the TIS are 9.1 ± 2.7 m a −1 for the 16 year run and 12.1 ± 3.0 m a −1 for the 20 year run, in agreement with steady-state mass-budget estimates of 9.2 ± 0.7 m a −1 .
A sequence of three Landsat visible light images separated by approximately a decade and spanning the epoch 1989-2010 (Fig. 3) shows two stationary features (white arrows) downstream of the grounding line at around 67°15′S, which we interpret to be ice rumples (Matsuoka et al. 2015), areas where the underlying bedrock shoals, re-grounds the ice shelf and forces the ice to locally rise above hydrostatic equilibrium (but, as distinct from an ice rise, does not result in significant changes in the local flow direction, only its magnitude). This is confirmed by the along-flow strain rate (Fig. 4) calculated by differentiation of surface velocity data (Rignot et al. 2011). This shows a region of flow retardation upstream of the rumples, associated with the basal drag there. Furthermore, time-varying basal melt will produce waves of thickness variations (undulations shown by black arrows in Fig. 3) that propagate at the ice-shelf velocity, which is equal to the product of the wavelength of the thickness undulations and their frequency of generation (Fig. 5) (Donnelly 2014). The calculated period of the thickness variation is 6-7 years (Donnelly 2014), based on an estimate of the wavelength from Landsat7 imagery, which will also correspond to the period of variation in the basal melt rate. High melt regions are near the deep grounding line and around the two rumples (Fig. 6a). For comparison, the spatial distribution of average basal melt rate from the 1992-2012 oceanmodel simulation is shown in Figure 6b, showing that the regions of highest melt correspond with the deepest ice near the southern grounding zone.
The observed significant anti-correlation between variations in elevation and modelled basal melt rates (Fig. 1b) supports the realism of the timing of the basal melt rates as driven by the time series of forcing fields representing the surrounding atmosphere and ocean. Furthermore, comparing the relative magnitudes of the variations indicates that the melt rates are quite sensitive to the range of natural variability described by the forcing. This basal melt rate variability also influences the ice flow of the TIS. Fig. 4. The TG ice-shelf along-flow strain rate from surface velocity data (Rignot et al. 2011); flow compression (negative area, dashed contours around 114°E) upstream of the ice-shelf rumples (white arrows) shows flow retardation associated with basal drag at these rumples. In contrast, the surface undulations (black arrows) are not associated with flow retardation. Also shown is the ASAID grounding line (Bindschadler et al. 2011a) (thick black).

Discussion
We suggest that the observed velocity and height variations are controlled by ice-shelf back-stress generated around two prominent ice rumples near the southern grounding zone. The reduced backstress arises from a combination of basal melting around the rumples directly reducing the area of contact, and an indirect, lagged response to thinner ice that is transported over the rumples from increased melting upstream. This back-stress from the rumples has two consequences for the grounded ice stream. First, the drag force of the ice flowing over the rumples provides a back-stress, retarding the flow of the TIS, as indicated by the observed compressive longitudinal strain rate. Secondly, reduction in the back-stress induces dynamic thinning of the ice shelf, with subsequent feedbacks on the grounded flow. Modelling shows the thickness of ice flowing over these rumples is highly sensitive to melt rates. Stronger melting thins the ice and reduces backstress, leading to an increase in ice velocity (Donnelly 2014) and probable thickness, which is likely to impact the location of the grounding zone. The observed thinning (thickening) appears coincident with faster (slower) velocities, which when considered along with the modelling results suggests a dynamic ice response to the variations in ice-shelf thickness and flow, driven by changes to the retarding (back-stress) forces present in the system due in turn to ocean-driven variability in the basal melt rates.
Coupled ice sheet-shelf modelling (Cornford et al. 2015) of the TG system (Donnelly 2014) forced with periodic varying basal melt rates produced associated modulations in ice velocities (differences in the range 90-150 m a −1 ) controlled by changing back-stress of the ice shelf at rumples (Matsuoka et al. 2015), illustrating the connection between basal melt rates, contact topography and flow modulation. The temporal modulation of basal melting also led to train undulations in TIS thickness (Donnelly 2014), analogous to observed TIS surface features in the along-flow direction. The observed spacing suggests variations in melt rates with a period of 6-7 years (Donnelly 2014). Similar undulations have been observed on the Pine Island Glacier ice shelf, where time-varying (in that case, seasonal) basal melt rates have been suggested as the cause (Bindschadler et al. 2011b).
Interactions between heat exchange across the continental shelf break, atmosphere-ocean processes over the continental shelf and melt-water advection from the Dalton Ice Shelf to TIS combine to produce a range of regimes for the supply of oceanic heat to drive ice-shelf melting (Gwyther et al. 2014). For example, a few percent decrease from the average in heat lost to the atmosphere above coastal polynyas can result in TIS melt rates that are double the average. This result is consistent with previous modelling studies of East Antarctic continental shelf seas, showing the high sensitivity of basal melting to relatively small changes in the amount of oceanic heat that is lost to the atmosphere (Cougnon et al. 2013;Khazendar et al. 2013). Warm modified Circumpolar Deep Water (MCDW) has been observed on the continental shelf near the TIS (Williams et al. 2011), and recent airborne geophysical surveys (including gravity, radar and magnetics) have revealed bathymetric pathways of suitable scale and depth to deliver this water into the ocean-filled cavity beneath the TIS (Greenbaum et al. 2015). As a result, modulation of the intrusions of the MCDW to the deepest part of the TIS cavity, near the southern grounding zone, is the most likely scenario driving the strong basal melting and its large variability.
Other plausible mechanisms that could explain the observed variations in velocity include: (1) natural changes in the surface accumulation of ice; or (2) changes in the subglacial hydrology. We doubt either  can account for the observations. Increased accumulation could act to make the ice plain more firmly grounded, which could slow both the TG and TIS. However, we expect the adjustment in ice flow in response to accumulation changes to be much slower than observed. Alternatively, there is mounting evidence for an extensive hydrological system in the Aurora Subglacial Basin (Wright et al. 2012), including subglacial lake discharge events (Smith et al. 2009). Variable basal hydrology is thought to influence glacier dynamics (Stearns et al. 2008;Stearns 2011). If the primary mechanism is changing basal hydrology, then this requires relatively less influence of the basal hydrology during the periods 1999-2002 and post-2010. Opportunities for temporary water storage (such as lakes that drain and fill) in the region are limited (Wright et al. 2012), suggesting that obvious hydrological-induced ice-flow changes are unlikely to be responsible for the observations reported here (Fig. 7). In addition, no substantial basal water transfer, as has been inferred from satellite altimetry in other regions, has been observed in the TG.
We therefore hypothesis the two main modes of variation of this system are as follows: (A) during periods of cool ocean forcing, basal melt rates are lower and the ice shelf thickens, leading to increasing back-stress at rumples, and the ice flow decreases (Fig. 7a); and (B) during periods of warm ocean forcing, basal melt rates are higher and the ice shelf thins, decreasing the back-stress at rumples, and the ice flow increases (Fig. 7b). Observational evidence does not exist to quantify how promptly the response should be between enhanced melting, which may occur in a region upstream of the rumple, and the advection of ice thickness into the region controlling back-stress that is sufficient to change the ice flow. However, we suggest that a time delay between melting and velocity changes is likely to be small given the proximity of the rumples to the region of deepest ice and highest melting near the grounding line.

Conclusions
From observations and modelling results presented here, we conclude that the variability of the TG over the past two decades is primarily ocean-driven. In particular, the anti-correlation of observed broadscale elevation changes and modelled melt rates support the timing and overall magnitude of the modelled melt rates, and indicates a high sensitivity of the TIS melt rates to ocean forcing. We argue that the strong variation in melt rate modulates the ice dynamics via the rumples that lie near the grounding line. It is likely that thinning, acceleration and grounding zone retreat of the TG will occur if subice-shelf ocean temperatures, and associated basal melt rates, rise in the future (Sun et al. 2016). Our results also demonstrate that the TIS is sensitive to relatively small changes in the supply of oceanic heat into the ice-shelf ocean cavity. In particular, variation in basal melting leads to periodic changes in back-stress at topographical rumples, which act to control ice-shelf velocity. This interplay between ocean forcing, bedrock topography and ice dynamics is analogous to changes observed for the Pine Island region in West Antarctica (Jacobs et al. 2011), and if TG responds in the same way, the melting may exceed the ice inflow, leading to sub-ice-shelf cavity enlargement and retreat of the grounding line.