Mid-Holocene thinning of David Glacier, Antarctica: chronology and controls

Quantitative satellite observations only provide an assessment of ice sheet mass loss over the last four decades. To assess long-term drivers of ice sheet change, geological records are needed. Here we present the first millennial-scale reconstruction of David Glacier, the largest East Antarctic outlet glacier in Victoria Land. To reconstruct changes in ice thickness, we use surface exposure ages of glacial erratics deposited on nunataks adjacent to fast-flowing sections of David Glacier. We then use numerical modelling experiments to determine the drivers of glacial thinning.

Thinning profiles derived from 45 ¹⁰Be and ³He surface exposure ages show David Glacier experienced rapid thinning of up to 2 m/yr during the mid-Holocene (∼ 6.5 ka). Thinning slowed at 6 ka, suggesting the initial formation of the Drygalski Ice Tongue at this time. Our work, along with ice thinning records from adjacent glaciers, shows simultaneous glacier thinning in this sector of the Transantarctic Mountains occurred 4–7 kyr after the peak period of ice thinning indicated in a suite of published ice sheet models. The timing and rapidity of the reconstructed thinning at David Glacier is similar to reconstructions in the Amundsen and Weddell embayments.

To identify the drivers of glacier thinning along the David Glacier, we use a glacier flowline model designed for calving glaciers and compare modelled results against our geological data. We show that glacier thinning and marine-based grounding-line retreat are controlled by either enhanced sub-ice-shelf melting, reduced lateral buttressing or a combination of the two, leading to marine ice sheet instability. Such rapid glacier thinning events during the mid-Holocene are not fully captured in continental- or catchment-scale numerical modelling reconstructions. Together, our chronology and modelling identify and constrain the drivers of a ∼ 2000-year period of dynamic glacier thinning in the recent geological past.

1 Introduction

Since the Last Glacial Maximum (LGM), ice sheets retreated in both hemispheres, causing a sea-level rise of ∼ 130 m (Clark et al., 2009). This period, the last glacial termination, represents Earth's last major period of climate warming, between ∼ 20 and 11.7 ka (Denton et al., 2010). Ice sheet retreat continued into the Holocene (11.7 ka–present), and global sea level stabilised at near pre-industrial levels by 6–7 ka (Bentley et al., 2014; Lambeck et al., 2014). Ice sheet thinning and retreat was rapid at times, potentially contributing to periods of rapid sea-level rise (Weber et al., 2014). However, quantitative age control on this retreat behaviour is only available for a limited number of sites in Antarctica and is entirely absent from David Glacier, which is the focus of this study (Bentley et al., 2014; Small et al., 2019). Discovered by the British National Antarctic Expedition (1901–1904), David Glacier drains the East Antarctic Ice Sheet (EAIS), traverses and incises the Transantarctic Mountains (TAM), and discharges into the western Ross Sea as the floating Drygalski Ice Tongue (Figs. 1 and S1). Comprising an area ∼ 210 000 km², the glacier is the largest in Victoria Land, representing a significant element of the Antarctic cryosphere, draining from Dome C and Talos Dome. Improved understanding of the timing and the processes that caused deglaciation help us to better understand the processes driving observed mass loss in parts of Antarctica today.

.

Figure 1 Map of the greater David Glacier/Terra Nova Bay region including mapped submarine geomorphic features and geographic features mentioned in text. Geomorphic cross sections (A–A’, B–B’ and C–C’) highlighted in Fig. S1. Inset shows extent of main map along with drainage basin of David Glacier and position of two local ice domes: Dome C (DC) and Talos Dome (TD). Data sources: satellite imagery of Bindschadler et al. (2008), ice surface velocity of Rignot et al. (2011), surface elevation of Howat et al. (2019), bed elevation of Fretwell et al. (2013) and bathymetry of Arndt et al. (2013)

Geological reconstruction of the Antarctic Ice Sheet (AIS) through time provides critical insights into the history of ice sheet change and can narrow uncertainty on Antarctica’s contribution to global sea-level rise. A geologic perspective on ice sheet behaviour is particularly useful due to the following.

1. Quantitative satellite observations only extend back ∼ 40 years (IMBIE, 2018; IPCC, 2013; Miles et al., 2013; Rignot et al., 2019), and these observations do not fully capture natural variability of ice sheet behaviour.

2. Integration of the data constraining marine extent and terrestrial thickness of an ice mass through time provides a robust data set for evaluation of numerical ice sheet model simulations, which are used to predict future ice sheet contributions to sea-level rise.

3. Recent studies of the complex interactions between the solid earth, cryosphere and ocean highlight the need to understand short- and long-term changes in the AIS (Barletta et al., 2018; Whitehouse et al., 2019; Meredith et al., 2019). Modern continuous GPS measurements of bedrock deformation record the integrated response to both modern and ancient ice mass change. The elastic component due to modern ice mass fluctuations can be modelled using satellite observations constraining ice mass balance. The viscous component represents a “memory” of ice load history over centuries to millennia. Ice load reconstructions constrained by cosmogenic dating of glacial samples provide a powerful tool to robustly model glacial isostatic adjustment (GIA), a complex suite of ice-sheet–bedrock–sea-level feedbacks critical for projections of future global sea level.

Reconstructions of marine-terminating outlet glaciers provide an opportunity to constrain and understand the terrestrial and marine sectors of ice sheets. At the LGM, the David Glacier thickened and expanded as an ice stream, coalesced with grounded marine-based ice and extended hundreds of kilometres into the western Ross Sea (Anderson et al., 2014; Livingstone et al., 2012; Shipp et al., 1999). During deglaciation, grounding-line retreat from north of Coulman Island was initiated by ∼ 13 ka, and by 11 ka, the grounding line had retreated near the Terra Nova Bay (TNB) region (Licht et al., 1996; Domack et al., 1999a; McKay et al., 2008; Anderson et al., 2014; Yokoyama et al., 2016). By ∼ 10 ka, ¹⁴C dating of acid insoluble organic (AIO) matter from sediment cores suggests the grounding line in the Ross Sea had migrated south of the David Glacier (Licht and Andrews, 2002). However, AIO dates in the Ross Sea are known to be unreliable due to potential input of anomalously old carbon by reworking, and this 10 ka constraint is likely to represent a maximum age of this retreat (Andrews et al., 1999; Rosenheim et al., 2013). Marine sedimentary cores from the deepest parts of the inner Drygalski Trough provide evidence of a lingering ice shelf, but anomalously old surface ¹⁴C ages mean that the timing of this ice shelf presence is uncertain (Licht et al., 1996; Frignani et al., 1998; Domack et al., 1999a; Licht and Andrews, 2002; Prothro et al., 2020). North of TNB, meteoric ¹⁰Be and compound-specific radiocarbon ages capture the retreating calving line and onset of open marine conditions by 8 ka (Yokoyama et al., 2016). Grounding-line retreat north of Ross Island was achieved by ∼ 8.6 ka, and a modern configuration was established by ∼ 2–4 ka (Anderson et al., 2014; McKay et al., 2016).

Outlet (Reeves and Priestley glaciers) and valley glaciers (Tucker and Campbell glaciers) along the northern Victoria Land sector began thinning at ∼ 17 ka, and the majority of thinning ceased by ∼ 7.5 ka, broadly coincident with a linear rise in sea level and ocean temperatures throughout deglaciation (Baroni and Hall, 2004; Johnson et al., 2008; Smellie et al., 2018; Goehring et al., 2019; Rhee et al., 2019). In contrast, outlet glaciers covering a large swath of the TAM from Southern Victoria Land to Southern TAM experienced episodic thinning during the early-mid Holocene, likely due to local topographic effects associated with marine ice sheet instability (Jones et al., 2015; Spector et al., 2017). Along the Transantarctic Mountains, mid-Holocene outlet glacier thinning and retreat occurred during a relatively stable period of air temperature and sea level after the majority of post-glacial sea-level rise and rising atmospheric temperatures occurred (Lambeck et al., 2014; Jones et al., 2015; Spector et al., 2017; Jones et al., 2020). Overall, variation in the timing of glacier thinning suggests sea-level rise, ocean warming and overdeepened subglacial topography controlled the reconstructed glacier behaviour.

Surface exposure dating using in situ terrestrial cosmogenic nuclides has transformed the ability to reconstruct the AIS through time (e.g. Stone et al., 2003; Mackintosh et al., 2007; Todd et al., 2010; Balco, 2011; White et al., 2011; Jones et al., 2015; Spector et al., 2017; Small et al., 2019). Through surface exposure dating of glacial deposits, this study aims to constrain the LGM to present behaviour of the David Glacier. Guided by our resulting ice thinning histories, we investigate the forcings and controls on glacier thinning using a glacier flowline model. We use geologic data and ice sheet model outputs to frame and inform the first detailed deglacial sensitivity experiments for David Glacier, building on previous flowline modelling approaches in Antarctica (Jamieson et al., 2012; Jones et al., 2015; Whitehouse et al., 2017).

2 Methods

2.1 Field and laboratory methods

Glacial erratics and glacially striated bedrock were collected adjacent to the David Glacier both upstream and downstream of the present-day grounding line. Downstream sites, Hughes Bluff and D'Urville Wall, are expected to show the most recent and dramatic changes in ice thickness, while smaller and longer-term changes are expected at the upstream site, Mt Kring. (Fig. 1) (Anderson et al., 2004; Bockheim et al., 1989). Using the structure-from-motion photogrammetry technique (Vargo et al., 2017), we use helicopter-based photography to construct a high-resolution digital elevation model of each field site, allowing integration of glacial geological field surveys with the regional and local geomorphology (Baroni et al., 2004) (Figs. S3 and S4).

We primarily sampled glacial erratic cobbles from local glacial deposits and glacially moulded bedrock surfaces along vertical transects perpendicular to the modern glacial flow direction. The cosmogenic nuclide inventories in cobbles and bedrock comprise signals of both exposure since the most recent deglaciation and nuclides accumulated during any prior exposure (known as “inheritance”). For long-lived nuclides such as ¹⁰Be, this inherited component can be significant (Sugden et al., 2005; Balco, 2011). Inheritance of cosmogenic nuclides is widespread in Antarctica as the cold-based nature of some ice acts to cover (preventing build-up of cosmogenic nuclides) without causing significant erosion (Sugden et al., 2005; Atkins, 2013). However, since cosmic rays attenuate rapidly with depth into rock, subglacial erosion of more than a few metres can reset the cosmic-ray signal. As such, we primarily focused on sampling faceted and striated glacial erratics because they indicate a subglacial origin which reduces the likelihood of a prior exposure history, and these intact features indicate minimal post-depositional erosion. Exposure ages from these samples are interpreted to track glacier thinning at the site because such samples would have been brought to the glacier surface by upward-flowing ice before being deposited at the ice margin (Stone et al., 2003; Mackintosh et al., 2007; Fogwill et al., 2012; Jones et al., 2015; Small et al., 2019).

We use cosmogenic nuclides ¹⁰Be and ³He to determine surface exposure ages on 45 samples. For ¹⁰Be analysis, quartz is the target mineral and is present in granites and sandstones. For ³He, pyroxene is the target mineral and present in dolerites. In the field, we preferentially sampled on isolated, local topographic high points, distant from areas of blue ice, snow drifts or local till deposits. While it is the focus of this study to use glacial erratics to track glacier thinning, we also collected bedrock samples to better understand complex exposure histories due to burial and the non-erosive nature of cold-based ice (Atkins, 2013; Joy et al., 2014).

Samples were processed at the cosmogenic nuclide facilities at Te Herenga Waka Victoria University of Wellington. For ¹⁰Be, we separated quartz by crushing and sieving to extract the 250–500 mum grain size fraction. Magnetic minerals were separated using a Frantz isodynamic magnetic separator. For granitic rocks, feldspars were removed by froth flotation. We etched samples for 1 d in 10 % hydrochloric acid, for 1 d in 2.5 % hydrofluoric acid (HF) and for 3 d in 1 % HF to further purify the quartz separates. Beryllium was extracted from quartz following an established method including the addition of a ⁹Be carrier, enhanced quartz etching and digestion followed by ion exchange chromatography and BeOH precipitation (Norton et al., 2008). Targets of BeO were packed and sent to the Purdue Rare Isotope Measurement Laboratory (PRIME lab) for analysis using accelerated mass spectroscopy (AMS). For ³He, we targeted pyroxene from the 125–250 mum grain size fraction. Pure pyroxene was obtained using an established HF etching method, and ³He $/$ ⁴He ratios were measured at the Berkeley Geochronology Center noble gas spectrometer (Bromley et al., 2014; Blard et al., 2015; Balter-Kennedy et al., 2020). We supplemented this data set with eight samples collected prior to the 2016–2017 austral summer, and those samples were processed for¹⁰Be or ³He using methods outlined in Oberholzer et al. (2003, 2008) and Di Nicola et al. (2009, 2012). Exposure ages are calculated by converting the AMS-derived ¹⁰Be $/$ ⁹Be ratio to ¹⁰Be concentration by subtracting a known amount of ⁹Be carrier added during chemical processing. Final exposure ages are calculated using ¹⁰Be or ³He concentration, field information (elevation, shielding and sample thickness) and production rate scaling method (LSDn) (Balco et al., 2008; Lifton et al., 2014; Marrero et al., 2016; Jones et al., 2019). For ¹⁰Be, we employed Jones et al. (2019) systematics (http://ice-tea.org/en/, last access: 21 October 2021), and for ³He, we use Balco et al. (2008) systematics (https://hess.ess.washington.edu/, last access: 21 October 2021). All relevant data used to calculate exposure ages are served on the ICE-D online database (Balco, 2020) within the David Glacier region (http://antarctica.ice-d.org/, last access: 21 October 2021) and in Supplement Tables S1–S5. The resulting glacier thinning chronologies derived from surface exposure ages inform our glacier modelling approach by providing geometric targets for evaluating the simulated maximum thickness of ice and the subsequent magnitude and duration of glacier thinning.

2.2 Glacier modelling approach

To understand potential controls on grounding-line migration and onshore thinning, we apply a glacier flowline model to the David Glacier. Originally designed to track grounding-line migration using a moving grid, the flowline model we employ has been described fully elsewhere (e.g. Vieli and Payne, 2005; Nick et al., 2009; Jamieson et al., 2012; Enderlin et al., 2013; Jamieson et al., 2014; Clason et al., 2016; Whitehouse et al., 2017). The underlying theory and governing equations are outlined in Vieli and Payne (2005), and here we describe the components pertinent to our research questions.

Our glacier modelling focuses on a ∼ 1600 km long flowline from Dome C to the western Ross Sea shelf edge. The domain consists of varying flow regimes and extends to areas where the ice stream converged with grounded ice in the western Ross Embayment during the LGM. Ice flow is governed by fundamental Eqs. (1)–(6).

Ice flow along the grounded portion of the model domain, u, is calculated using the shallow ice approximation, which assumes ice flow only from internal ice deformation and horizontal velocity defined as

$$\begin{array}{}\text{(1)}& u=C\left({\displaystyle \frac{\mathit{\delta }s}{\mathit{\delta }x}}\right)H^{n+\mathrm{1}},\end{array}$$

where n(=3) is Glen's flow law exponent, s is ice surface elevation, H is ice thickness and C is a constant given by

$$\begin{array}{}\text{(2)}& C={\displaystyle \frac{\mathrm{2}A({\mathit{\rho }}_{i}g{)}^n}{n+\mathrm{2}}},\end{array}$$

where A is the temperature-dependent rate factor (for −20 ^∘C), ρ_i is ice density (0.910 g/cm³) and g is gravitational acceleration (9.8 m/s²).

As ice reaches flotation, ice is assumed to spread unidirectionally along the flowline, and ice flow is balanced by

$$\begin{array}{}\text{(3)}& \mathrm{2}{\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }x}}Hv{\displaystyle \frac{\mathit{\delta }u}{\mathit{\delta }x}}={\mathit{\rho }}_{\mathrm{i}}gh{\displaystyle \frac{\mathit{\delta }s}{\mathit{\delta }x}},\end{array}$$

where and v is the vertically average effective viscosity defined as

$$\begin{array}{}\text{(4)}& v=A^{\frac{-\mathrm{1}}{n}}{\left[{\left({\displaystyle \frac{\mathit{\delta }u}{\mathit{\delta }x}}\right)}^{\mathrm{2}}\right]}^{(\mathrm{1}-n)/\left(\mathrm{2}n\right)},\end{array}$$

where $\frac{\mathit{\delta }u}{\mathit{\delta }x}$ is the effective strain rate.

For an ice stream, Eq. (3) is modified by including a basal traction coefficient β which is linearly related to the horizontal velocity u and is given by

$$\begin{array}{}\text{(5)}& \mathrm{2}{\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }x}}hv{\displaystyle \frac{\mathit{\delta }u}{\mathit{\delta }x}}-\mathit{\beta }u={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\mathit{\delta }s}{\mathit{\delta }x}}.\end{array}$$

We also incorporate a width term, including f_lat, a lateral buttressing factor, to account for lateral buttressing along a coupled ice stream shelf. Thus, as in previous applications of this model (Nick et al., 2010; Jamieson et al., 2012; Whitehouse et al., 2017) the final ice thickness (h) and width (W) averaged ice flow (u) is computed using the following equation:

$$\begin{array}{}\text{(6)}& \begin{array}{rl}\mathrm{2}{\displaystyle \frac{\mathit{\delta }}{\mathit{\delta }x}}\left(Hv{\displaystyle \frac{\mathit{\delta }u}{\mathit{\delta }x}}\right)& -u\mathit{\beta }|u{|}^{m-\mathrm{1}}-{\displaystyle \frac{H}{W}}{\left({\displaystyle \frac{\mathrm{5}}{\mathrm{2}AW{f}_{\mathrm{lat}}}}\right)}^{\frac{\mathrm{1}}{n}}|u{|}^{\frac{\mathrm{1}}{n}-\mathrm{1}}u\\ & ={\mathit{\rho }}_{\mathrm{i}}gH{\displaystyle \frac{\mathit{\delta }s}{\mathit{\delta }x}}.\end{array}\end{array}$$

f_lat controls the strength of lateral drag exerted at the ice stream margins. For all experiments in this study, f_lat is only altered where the ice stream is floating. When f_lat equals 1, full effects of lateral drag are felt; when values of f_lat are greater than 1, the magnitude of buttressing applied to the margins of the ice stream is reduced non-linearly.

The inclusion of the width term allows modelling of changing offshore trough width, which is common in palaeo-ice streams and outlet glaciers and is further supported by geomorphic interpretations using high-resolution bathymetry data (Jamieson et al., 2012; Livingstone et al., 2012, 2016). To determine offshore trough width, we map all glacial geomorphic features using the global multi-resolution topography data set within GeoMapApp (http://www.geomapapp.org, last access: 21 October 2021) based on analogues (Dowdeswell et al., 2016) and prior regional mapping (Stutz, 2012; Shipp et al., 1999; Anderson et al., 2014; Halberstadt et al., 2016).

Whitehouse et al. (2017) improved the ice shelf dynamics to incorporate the treatment of large horizontal grounding-line movements that are expected for the palaeo-David Glacier based on the long-term presence of the Drygalski Ice Tongue and evidence of extensive sub-ice-shelf conditions inferred from marine sedimentary cores (Domack et al., 1999a; McKay et al., 2008). In the model, the ice shelf geometry evolves with time, and variations in ice shelf thickness and extent are fed into lateral drag calculations (Whitehouse et al., 2017). Improved treatment of ice shelf dynamics allows exploration of the grounding-line sensitivity to ice shelf feedbacks. Importantly for this study, a reduction in lateral buttressing from adjacent, coalesced ice is expected as the expanded David Glacier and grounded ice in the Ross Sea decouple, which has been suggested from interpretation of submarine geomorphic features in the western Ross Sea (Shipp et al., 1999; Halberstadt et al., 2016).

Accumulation in the model is applied using modelled modern rates at the central flowline and multiplied to account for the ice stream width (W). Accumulation is then scaled to represent warmer or cooler conditions in our experiments (see Sect. 2.2.2). Where an ice shelf is present, sub-ice-shelf melt rates (SIMRs) are applied as a linear function of the depth of the ice shelf draught. From a minimum rate of 0.1 m/yr at 0 m (the ocean surface), SIMR increases to a maximum value at a depth of 500 m, and where ice shelf draught is deeper than 500 m, the maximum value is applied.

2.2.1 Boundary conditions

The model domain, sampled at 5 km resolution, runs along a ∼ 1600 km flowline from Dome C to the continental shelf break (Figs. 2 and S6). Figures 1 and S1 provide map extent of cross section views detailing the subglacial topography surrounding the modern-day grounding line. To simplify the flowline domain, we focus on the dominant ice stream emanating from Dome C. Using standard basin delineation tools in ArcGIS, we construct a drainage basin outline from surface elevation and ice velocity (Rignot et al., 2011; Howat et al., 2019). We then construct a flowline that follows the centre of the ice stream and thereafter calculate ice stream width perpendicular to this flowline (W) to create a defined lateral ice stream margin across the entire domain. Following Jamieson et al. (2012), the lateral ice stream margins in the offshore parts of the domain are determined using geophysical data to map the distribution of trough parallel mega-scale glacial lineations (MSGLs) that are indicators of past ice stream flow and thus show the width of the ice stream. Onshore, the lateral margins are defined by the valley width. The ice stream width perpendicular to the flowline is thereafter used to control the lateral stress applied by coalesced ice along the flowline (Eq. 6).

Figure 2 Map of model domain with modern ice surface velocity of Rignot et al. (2011) and bathymetry of Arndt et al. (2013). Flowline follows red line; ice stream width in black. RIS: Ross Ice Shelf; BG: Byrd Glacier; MUG: Mulock Glacier; SG: Skelton Glacier; MKG: Mackay Glacier; MAG: Mawson Glacier; RG: Reeves Glacier; PG: Priestley Glacier; MG: Mariner Glacier; AG: Aviator Glacier; TG: Tucker Glacier; CI: Coulman Island; DT: Drygalski Trough; HB: Hughes Bluff; MK: Mt Kring.

To evaluate model performance and the appropriate boundary conditions for deglacial experiments, we tune basal traction (β), accumulation, sub-ice-shelf melt rate (SIMR) and rate factor (A) to reproduce the modern geometry and flow speed of the David Glacier as closely as possible. Modern ice surface elevation, surface mass balance and ice velocity are well constrained from satellite and in situ measurements (Rignot et al., 2011; van Wessem et al., 2018; Howat et al., 2019). We employ an ice mass accumulation scheme whereby tributary ice mass is added along the length of the ice stream (Jamieson et al., 2014; Whitehouse et al., 2017) to account for confluent ice from tributary glaciers. This tributary injection from a secondary ice stream originating from Talos Dome is included in the surface mass budget for the David Glacier model domain (Fig. 2).

Ongoing aerogeophysical and remote sensing techniques continue to reveal new detail of the subglacial environment of the AIS (Schroeder et al., 2020; MacGregor et al., 2021). Common among many TAM outlet glaciers, a prominent bedrock high immediately upstream of the modern grounding line provides a significant stabilising impact on the David Glacier (Fretwell et al., 2013; Morlighem et al., 2019). The bed geometry along the flowline is derived solely from Bedmap2 and the International Bathymetric Chart of the Southern Ocean (Fretwell et al., 2013; Arndt et al., 2013).

Geologic data from former glaciated terrains and satellite observations can be used to approximate basal conditions, yet a general lack of in situ data from the ice sheet basal environment leads to enduring uncertainties in the basal stress regime (Joughin et al., 2006; Stokes et al., 2015). Based on existing knowledge of sea floor morphology and consistent with previous modelling experiments, the basal traction parameter used in this model attempts a first-order approximation of two subglacial basal environments: (1) relatively high basal traction for onshore ice flow over bedrock and (2) relatively low basal traction for offshore ice flow over soft sediments associated with basal till deposits (Dowdeswell et al., 2004; O Cofaigh et al., 2005).

We use modelled surface mass balance (van Wessem et al., 2018), calculated SIMR (Wuite et al., 2009) and modern geometry (Howat et al., 2019; Fretwell et al., 2013) to set the model up so that it replicates modern flow conditions and geometry at a steady state (Table 1). We then tune the basal traction parameter to produce ice thickness and velocities which do not change significantly. The model is unable to reproduce stable conditions along the steep ice surface profile as the David Glacier descends from the ice sheet interior. In this zone, the modelled ice surface is steep and undergoes thinning throughout a 2000-year initiation period. In an effort to stabilise the ice surface upstream of the grounding line, we tune the basal traction parameter, effectively stiffening the bed, to reduce the modelled instability. We note here that according to Eq. (5), as the ice surface evolves, basal traction evolves naturally; for example, as the ice profile gets steeper or thicker, the basal traction evolves. Resultant modelled estimates of basal shear stress approach 100 kPa in this setting. While difficult to constrain with in situ measurements, Zoet et al. (2012) suggest higher stresses, such as these, should be expected near the modern-day grounding zone, which is consistent with the modelled stress distribution in this study. Overall, the modern grounding-line position, ice shelf thickness and ice shelf length do not vary significantly over a 2000-year modelled initiation period.

Table 1 Parameter values used for sensitivity experiments indicating choices for model tuning to modern conditions over a spin-up period. Glacial conditions are parameter values used to simulate deglaciation. SIMR: sub-ice-shelf melt rate. Accumulation values reported as percentage of RACMO2 (van Wessem et al., 2018).

2.2.2 Deglaciation approach

LGM geometric boundary conditions incorporate modelled ice thickness and surface velocity from W12, an ice sheet reconstruction based on available geologic data (Fig. S6) (Whitehouse et al., 2012). Our ice flow model does not adjust for isostatic deformation as it evolves. Therefore, we use the W12 modelled ice thickness and modern bed elevation (Fretwell et al., 2013), which allows for a calculated estimate of upper ice surface without introducing uncertainty involving variable along-flow isostatic response and dynamic topography associated with the long-term evolution of the Antarctic subglacial topography (Stern et al., 2005; Whitehouse et al., 2019; Paxman et al., 2019). We do not compare our sensitivity experiment's results with those of the W12 ice reconstruction; we purely use a calculated upper ice surface as a starting point for our model during spin-up. The resulting ice thickness geometrically fits to geologic data (i.e. covers highest-elevation Holocene aged erratics) and serves as an initial ice profile for deglaciation sensitivity experiments.

To account for environmental changes during deglaciation, accumulation and internal ice temperature are tuned over the initial 7500 model years to ensure an LGM configuration that is consistent with geological constraints and has a grounding line that does not move significantly. Knowledge of past accumulation over glacial–interglacial cycles is restricted to ice core data and internal ice sheet layer mapping near high-elevation ice domes (Siegert, 2003; Frezzotti et al., 2005; Buiron et al., 2011; Cavitte et al., 2018). Generally, modern accumulation patterns show relatively low accumulation (0.02 m/yr) over the East Antarctic interior and high accumulation (0.2 m/yr) near the coastline (Arthern et al., 2006; Lenaerts et al., 2012; van Wessem et al., 2018). Interpretations from ice core records suggest LGM accumulation rates were lower than modern accumulation rates and generally correlate well with temperature proxy records throughout the Holocene (Siegert, 2003; Veres et al., 2013). We use a scaling relationship between modern and palaeo-accumulation patterns and estimate that accumulation at the start of the model run was roughly 75 % of modern accumulation (Veres et al., 2013). Internal ice temperature is increased through time to represent the increase in temperature that occurred during deglaciation. The internal ice temperature during deglaciation is not known, and for this study, we used values of −25 ^∘C for the first 7500 model years and −20 ^∘C during the remaining 7500 model years as a way to represent an appropriate amount of warming in the ice column. Deglacial sensitivity experiments use a range of accumulation and internal ice temperature forcings representing the potential scale of change experienced through a deglaciation (Table 1), allowing exploration of the relative influence of lateral buttressing reduction (LBR) and sub-ice-shelf melt rates (SIMR) on glacier thinning and retreat.

During the LGM, the David Glacier coalesced with grounded ice in the Ross Sea and was susceptible to ocean-driven basal melt and changes in lateral drag (buttressing) from surrounding grounded ice. We therefore run a suite of individual experiments that allow us to initiate grounding-line retreat by either linearly increasing SIMR or decreasing lateral buttressing over a 500-year period, with each model run applying different perturbed values for SIMR or LBR (Table 2). For combined forcing experiments (MS1–3), we alter both SIMR and LBR as above until the grounding-line retreats to a near modern configuration. We introduce a 1  kyr timing lag between SIMR forcing and LBR. This staggered forcing seeks to simulate ocean-driven melt occurring first followed by a reduction in lateral buttressing. This scenario is designed to represent a “calving bay” environment, potentially led by ocean-heat-driven sub-ice-shelf melt, as recognised in the western Ross Sea (Domack et al., 1999a, 2003, 2006) and East Antarctic margin (Leventer et al., 2006; Mackintosh et al., 2014).

All sensitivity experiments run for 15 kyr with an initial spin-up period lasting 7.5 kyr for SIMR forcing and 8.5 kyr for LBR forcing, at which point the forcing perturbation relating to increased SIMR or reduction in lateral buttressing is applied for the remaining modelled period. We run our experiments for 15 kyr to provide ample spin-up time at glacial conditions and force our model from 7.5 kyr to explore forcings that, by the end of the modelled period, result in a configuration consistent with modern observations. When forcing the model, we linearly increase the forcing over a 500-year period. To simulate the natural variability that might be expected in an ocean forcing record, we apply a fluctuation of up to 0.5 m/yr magnitude using a random noise generator, and this variation is then added on top of the 500-year increase in forcing that is already applied. Further, it is not the intention of these sensitivity experiments to reproduce the exact timing of marine-based grounding-line retreat along the David Glacier, but rather we aim to explore the range of possible drivers that would enable the overall scale of retreat observed between the LGM and present day. All idealised scenarios are presented in Table 2.

Table 2 Idealised deglacial sensitivity experiments. Sub-ice-shelf melt rate (SIMR) varies within a 0.5 m/yr window around the value reported. SIMR is forced at 7500 model years, and lateral buttressing reduction (LBR) is forced at 8500 model years. X/Y indicates initial conditions and forced conditions, respectively. If held constant, one value is reported.

3 Results – chronology

Of the 45 samples analysed in this study, 21 yield mid-Holocene exposure ages. Holocene aged samples are interpreted to have received minimal prior exposure (i.e. inheritance) and suggest a simple exposure history. Focusing on two sites, we derive a high-resolution chronology of glacier thinning from Mt Kring and Hughes Bluff.

3.1 Mt Kring

Mt Kring, situated in the interior of the EAIS, flanks the streaming ice of David Glacier draining Dome C. The nunatak, composed of layered dolerite, stands ∼ 300 m above the modern ice surface, henceforth “above the ice”. Of the 24 samples processed, 3 are dolerite bedrock and 21 are glacial erratics of mixed lithologies (Figs. 3 and S3). The bedrock surface at the peak of Mt Kring lacks glacial striations or erratics and has an apparent ³He bedrock exposure age of 554 ± 91 ka. A discontinuous, steep ridge line between ∼ 200–300 m above the ice was not sampled due to inaccessibility. At ∼ 200 m above the ice, the bedrock is striated parallel to the modern ice flow direction and has an apparent ³He exposure age of 550 ± 82 ka. The highest-elevation erratic is found at 171 m above the ice with increasing erratic abundance culminating in a drift sheet covering nearly all bedrock lower than ∼ 110 m above the ice. Striated, bulleted cobbles and boulders of dolerite, basalt and sandstone are common in the thin, patchy drift covering the bedrock. Based on our erratic exposure dating, the glacial drift here is a composite of three age populations (Fig. 4a). The younger population spans 7.4–5.5 ka (¹⁰Be, n=5) (Fig. 4b). An older population spans 51–25 ka (¹⁰Be, n=8 and ³He, n=4). The oldest erratic age population shows scattered evidence of glacial behaviour between 123–63 ka (¹⁰Be, n=2, ³He, n=2). The older populations are likely an artefact of repeated burial by overriding ice and multiple periods of exposure (i.e. inheritance).

Figure 3 Orthomosaic map of Mt Kring with all exposure ages, ¹⁰Be (red) and ³He (grey), and total errors. Erratics and bedrock ages plotted as circles and squares, respectively. Bold ages indicate Holocene age erratics, and italicised ages indicate bedrock exposure. Large inset displays ice surface velocity of Rignot et al. (2011), and small inset displays LIMA data of Bindschadler et al. (2008). DG: David Glacier; DIT: Drygalski Ice Tongue; TNB: Terra Nova Bay; MK: Mt Kring.

Figure 4 Relative elevation, above the ice, vs. calculated surface exposure ages for Mt Kring. (a) All exposure ages. Erratics plotted as circles, bedrock as squares(¹⁰Be, red, and ³He, grey). Errors bars show total uncertainty. Panel (b) shows only Holocene ages.

3.2 Hughes Bluff

Hughes Bluff, situated on the Scott Coast, is a granite outcrop along the southern flank of the David Glacier. The outcrop is glacially scoured exhibiting spectacular rochés moutonnées, crescentic gouge marks and striations parallel to modern flow directions (Figs. 5 and S4). Scattered glacial erratics blanket the entire outcrop. ¹⁰Be exposure ages from 15 erratics span 6.7–4.3 ka (Fig. 6). The majority of glacial erratics are dated to between 6.7 and 6.2 ka. Two bedrock surface exposure ages sampled from rochés moutonnées at the highest and lowest outcrops (20.55 ka ± 2.10, 5.5 ± 0.47 ka) suggest significant wet-based glacial erosion since the LGM. At 20 m above the ice, similar bedrock and erratic ages provide evidence of recent emergence of this low-lying outcrop since 5.5 ka. The fact that well-developed landforms of glacial erosion occur at the highest outcrops at Hughes Bluff, including evidence of abundant basal sliding and plucking, indicate that ice thickness during the LGM was considerably greater than 230 m, the maximum elevation at Hughes Bluff (Frisia et al., 2017). While the onset and magnitude of thinning prior to 6.7 ka is not constrained, the Hughes Bluff chronology clearly indicates a period of rapid ice surface lowering during the mid-Holocene.

Figure 5 Orthomosaic map of Hughes Bluff with all exposure ages, ¹⁰Be (red) and ³He (grey), and total errors listed. Erratics and bedrock ages plotted as circles and squares, respectively. Bold ages indicate Holocene age erratics, and italicised ages indicate bedrock exposure. Large inset shows ice surface velocity of Rignot et al. (2011), and small inset shows LIMA data of Bindschadler et al. (2008). DG: David Glacier; DIT: Drygalski Ice Tongue; TNB: Terra Nova Bay; HB: Hughes Bluff; CY: Camp Yellow.

Figure 6 Relative elevation, above the ice, vs. calculated surface exposure ages for Hughes Bluff. (a) All exposure ages. Erratics plotted as circles and bedrock as squares (¹⁰Be, red, and ³He, grey). Errors bars show total uncertainty. Panel (b) shows only Holocene ages.

3.3 High-elevation constraints

In an effort to identify higher-elevation glacial activity and long-term erosion history, fieldwork was undertaken along the northern flank of David Glacier from the D'Urville Wall area (Mt Neumayer to Cape Phillipi) (Fig. S5). The D'Urville Wall is a steep, high-elevation (> 400 m above the ice) granite outcrop delimiting the northern flank of David Fjord. The geomorphology and exposure ages from two field sites situated high above the David Glacier (Mt Neumayer and Cape Phillipi) limit the reconstruction of the glacier surface during the LGM. Mt Neumayer extending above the D'Urville Wall forms a rounded summit with faint striations sub-parallel to modern ice flow and few scattered erratics (Baroni et al., 2004). In areas 400 m above the ice, bedrock samples contain weathering rinds and deep weathering pits filled with erratics (Baroni et al., 2004; Giorgetti and Baroni, 2007). Bedrock exposure ages from Mt Neumayer (649 m above the ice, 642 ± 61 ka), a terrace on top of the D'Urville Wall (418 above the ice, 116 ± 10 ka) and Cape Phillipi (∼ 300 m above the ice, 532 ± 52 and 957 ± 98 ka) suggest either a thin cover of cold-based ice or ice-free conditions through the LGM. Bedrock exposure ages from the D'Urville Wall area (including Mt Neumayer and Cape Phillipi) do not allow a precise estimate of the past ice surface along the northern flank of David Glacier. Supported by geomorphic evidence from Hughes Bluff, which indicates ice thicker than 230 m, and the LGM limit of ∼ 400 m a.s.l. derived from drift deposits in TNB, we suggest the past ice surface was between 300 and 649 m higher than today in this area (Stuiver et al., 1981; Orombelli et al., 1990; Di Nicola et al., 2009) (Fig. S2).

3.4 Palaeo-thinning rates

Using the high-resolution chronology from Hughes Bluff and Mt Kring, we derive a mean estimate of palaeo-thinning rates along David Glacier using a weighted least squares regression scheme within the iceTEA plotting tools (http://ice-tea.org/en/) (Jones et al., 2019). These reconstructed glacier thinning rates are compared to modern thinning rates derived from satellite data and continental-scale ice sheet models (Small et al., 2019). At Mt Kring, we reconstruct a 2 ka thinning event from 7.5–5.5 ka of up to 0.19 m/yr (0.06–0.19, 2σ) thinning rate (Fig. 7A). At Hughes Bluff, we reconstruct a period of rapid thinning from 6.7–6.2 ka of up to 2 m/yr (0.19–2.06 m/yr, 2σ) followed by ∼ 4 kyr period of minimal thinning (Fig. 7b). The reconstructed palaeo-thinning along the David Glacier during the mid-Holocene is synchronous with rapid thinning reconstructed at a number of sites in the Weddell embayment (Nichols et al., 2019; Johnson et al., 2019; Small et al., 2019; Bentley et al., 2017), Amundsen Sea sector (Johnson et al., 2008, 2014, 2017, 2020), Marie Byrd Land (Stone et al., 2003) and the Transantarctic Mountains (Todd et al., 2010; Jones et al., 2015; Spector et al., 2017; Jones et al., 2020). However, the rate of palaeo-thinning reconstructed at Hughes Bluff is one of the highest rates recorded in Antarctica, comparable to a reconstruction from nearby Mackay and Mawson glaciers (Small et al., 2019; Jones et al., 2015, 2020).

Figure 7 Linear thinning rates derived from relative elevation, above the ice, and mid-Holocene exposure ages estimated using http://ice-tea.org/en/ from Jones et al. (2019) for (a) Mt Kring Holocene samples and (b) Hughes Bluff samples with outliers identified and omitted using outlier detection in http://ice-tea.org/en/.

Previous work from inland sites in a similar setting as Mt Kring suggested that such sites changed little or may have thickened during the Holocene due to accumulation increases (Bockheim et al., 1989; Denton et al., 1989). In contrast, the synchronous timing between thinning at Hughes Bluff and Mt Kring but lower thinning rate observed at Mt Kring suggests that this thinning event was driven by ocean–ice-sheet interactions at the coast that rapidly propagated inland. Such dynamic thinning, where increased ice shelf basal melting and grounding-line retreat leads to accelerated flow and inland thinning, is well documented along modern marine-terminating margins of the Greenland and Antarctic ice sheets (Pritchard et al., 2009, 2012), yet identification in the geologic record is rare, primarily due to a lack of exposed bedrock in the upper reaches of glacier catchments. The thinning history from David Glacier allows for a unique comparison with the broad pattern of dynamic thinning derived from the modern satellite record and suggests dynamic thinning can occur > 100 km into the interior of the EAIS and can persist over multiple millennia.

3.5 Exposure age data – ice sheet model comparison

During the LGM, ice core records and numerical model outputs suggest that the EAIS experienced widespread thinning in its interior, while coastal sites experienced extensive thickening (Verleyen et al., 2011; Mackintosh et al., 2014). Our data from Mt Kring show that the ice sheet was thicker than present at this site during the LGM and recorded at least 171 m of thinning during the Holocene. This provides a critical tie point between the high-elevation, low-accumulation ice domes where ice cores are drilled and the low-elevation, high-accumulation coastal sites with more abundant geologic data. Mt Kring is the first site along the high-elevation margin of the EAIS to constrain LGM ice sheet behaviour, and it represents a unique site to compare with other ice sheet reconstructions and continental-scale ice sheet modelling experiments.

The chronologies from David Glacier provide critical insights to the rate, timing, style and magnitude of thinning from a marine-based outlet glacier. The thinning histories presented here provide context to modelled continental-scale ice sheet reconstructions. Using a data–model comparison software, http://dmc.ice-d.org/, last access: 21 October 2021, we extracted a modelled elevation history for David Glacier from five different ice sheet models (four continental-scale and one regional-catchment scale) and one ice sheet reconstruction created to drive a GIA model (Argus et al., 2014; Pollard et al., 2016, 2017, 2018; Kingslake et al., 2018; Lowry et al., 2019). To overcome differences in spatial resolution, the scheme extracts an interpolated ice elevation history from the model grid cell containing the field site and its neighbouring model grid cells. The suite of ice sheet models in which we compare against our geologic data represent a variation in model grid size, flow approximations, model physics used, forcing and boundary conditions. We use them to provide a first-order approximation for expected results from the ice sheet modelling community.

The resulting data–model comparisons for Hughes Bluff and Mt Kring reveal a noticeable mismatch in time during the main phase of thinning (Fig. 8). For Mt Kring, ice sheet models indicate a phase of thinning that precedes our thinning history by ∼ 4–7 kyr. For Hughes Bluff, the timing lag is comparable to Mt Kring with one notable exception being the deglacial model of ICE-6G (Argus et al., 2014). This improved match is likely because ICE-6G is constrained by multiple relative-sea-level curves along the Scott Coast (Baroni and Hall, 2004; Hall, 2009; Argus et al., 2014). The style of thinning is variable between models with noticeable short-lived pauses in thinning, mainly during the Antarctic Cold Reversal (∼ 15–13 ka) (Pedro et al., 2016). Relative to Hughes Bluff, the magnitude of modelled elevation changes at Mt Kring geometrically fit better with the surface exposure age data. At Mt Kring, the average modelled thickness change throughout the deglacial phase for all models in Fig. 8a is 190 m ± 117 m, which compares well with the 173 m thickness change derived from our ice thinning chronology. In contrast, at Hughes Bluff, we capture only 171 m of thickness change over the deglacial period, and the average modelled thickness change for all models in Fig. 8b is 623 m ± 142 m. These comparisons suggest that our surface exposure ages capture more of the overall ice thinning at Mt Kring than at Hughes Bluff. Overall, the data–model mismatch may be related to (1) individual topographic features/outcrops not being spatially resolved in the models, (2) limited constraints on ocean/atmosphere forcing (e.g. Lowry et al., 2019) and (3) poorly constrained model parameters that influence basal sliding, isostatic adjustment and ice flow/rheology, which impact the rate of ice sheet response to a climate forcing (Lowry et al., 2020; Kingslake et al., 2018).

Figure 8 Surface exposure age data – ice sheet model comparisons for (a) Mt Kring and (b) Hughes Bluff. Ice sheet modelled reconstruction data and supporting information from http://dmc.ice-d.org/ and Lowry et al. (2019). Kingslake2018_WDC uses WAIS Divide accumulation forcing, and Lowry2019_EMV is a model in which enhanced mantle viscosity is applied.

4 Results – deglacial sensitivity experiments

Based on the data-constrained thinning histories presented in this study, the main episode of glacier thinning occurred during the mid-Holocene and captures > 200 m of glacier change along the David Glacier. This time period does not correlate with significant increases in atmospheric temperature or global mean sea level (e.g. Lambeck et al., 2014; Menviel et al., 2011; Liu et al., 2009); therefore, we ask the following questions. (1) What is the role of ocean heat in driving the observed glacier thinning and retreat? And (2) given that ice from David Glacier coalesced with grounded ice in the Ross Sea during the LGM, what impact does ice sheet buttressing have on the timing and style of glacier thinning and retreat? The style and rate of modelled thinning and retreat from our flowline experiments are compiled in Table 3. We note we carried out many more experiments than reported here, and only those that show notable changes are presented. While we do not force our simulations by any date-specific reconstructions or proxy data, we co-view our model results against our geologic data to allow generalised time-varying geometric relationships to be compared.

Focusing solely on SIMR, a set of sensitivity experiments (M1–3) simulate the impacts of enhanced ocean heat on grounding-line retreat (Fig. 9). Modern SIMR have been calculated along the Drygalski Ice Tongue with an average of 4.89 m/yr ± 3.38 m/yr (Wuite et al., 2009). For our experiments, we progressively increased SIMR until retreat is initiated. Threshold values represent this stepwise increase in SIMR. After a 7.5 kyr spin-up period and for melt rates between 2 and 10 m/yr (represented by experiment M2), grounding-line retreat is initially rapid, but the final grounding line remains pinned to the prominent sill at the mouth of the David Fjord. A threshold SIMR of 11 m/yr achieves rapid grounding-line retreat behaviour and is consistent with the modern grounding-line position. Final modelled ice surface reconstructions for experiment M2 place the upper ice surface ∼ 300 m above the Hughes Bluff site. At Mt Kring, the final ice surface reconstruction for experiments M2 and M3 lies below the modern ice surface. In the high melt case (M3), the grounding-line position pauses at the sill for approximately 5 kyr. The final retreat phase from the sill to a modern position correlates with a final ice surface broadly consistent with modern observations.

Figure 9 (a) Profile view: 15 kyr evolution of modelled upper ice surface along flowline for sub-ice-shelf melt rate (SIMR) experiments. Black bars represent elevation range from three sites in this study. Orange triangles represent offshore features marking observed grounding zone wedges. (b) Onshore elevation and grounding-line evolution for SIMR experiments. Top subpanel: forcings applied (no relationship between LBR and SIMR). In SIMR cases, we varied the SIMR perturbation within a 0.5 m/yr window at 500-year temporal frequency to simulate short-lived pulses of relatively warmer or cooler water mass changes. Top middle subpanel: time-transgressive elevation profile for Mt Kring with exposure ages; bottom middle subpanel: time-transgressive elevation profile for Hughes Bluff surface exposure ages; and lower subpanel: evolution of grounding-line position with modern grounding-line position, bathymetric sill at the mouth of David Fjord and mapped grounding zone features. We co-view exposure age data in these plots to allow generalised time-varying geometric relationships to be compared.

Focusing on LBR experiments (S1–3), we simulate the impacts of glacier–ice-sheet decoupling on grounding-line retreat, a scenario suggested from glacial geomorphic features in the western Ross Sea (Shipp et al., 1999; Halberstadt et al., 2016) (Fig. 10). After a 8.5 kyr model spin-up period, ice shelf buttressing is incrementally reduced in successive simulations until retreat is initiated. Retreat occurs when lateral buttressing is reduced by 4 %. Ice shelf debuttressing by at least 40 % retreats the modelled grounding line to near-modern configuration, deep in the David Fjord. In both cases, the reduced buttressing forces rapid grounding-line retreat to a prominent sill at the mouth of the David Fjord. In these scenarios, the resulting modelled upper ice surface remains above the Hughes Bluff site when the grounding line is pinned to the sill. At Mt Kring, modelled rapid thinning is synchronous with Hughes Bluff, yet the 40 % LBR simulation (S3) results in an unrealistic final ice surface elevation hundreds of metres below the observed modern ice surface elevation.

Figure 10 (a) Profile view: 15 kyr evolution of modelled upper ice surface along flowline for lateral buttressing reduction (LBR) experiments. Black bars represent elevation range from three sites in this study. Orange triangles represent offshore features marking observed grounding zone wedges. (b) Onshore elevation and grounding-line evolution for LBR experiments. Top subpanel: forcings applied (no relationship between LBR and SIMR). Top middle subpanel: time-transgressive elevation profile for Mt Kring with exposure ages; bottom middle subpanel: time-transgressive elevation profile for Hughes Bluff surface exposure ages; and lower subpanel: evolution of grounding-line position with modern grounding-line position, bathymetric sill at the mouth of David Fjord and mapped grounding zone features. We co-view exposure age data in these plots to allow generalised time-varying geometric relationships to be compared.

Finally, we force our model with a combination of SIMR and LBR. Our approach is similar to that used for individual forcings, except that this time we use a melt rate (MS1: 5 m/yr) and progressively reduce lateral buttressing until the grounding line retreats to near-modern configuration. This approach is repeated for two additional cases of enhanced SIMR (MS2: 7 m/yr and MS3: 9 m/yr). Overall, when forcings are combined, lower threshold values are required to initiate thinning and retreat (Fig. 11). For experiments MS1–MS3, the final modelled upper ice surface reconstruction at Mt Kring lies below the modern ice surface. For experiments MS2–3, Hughes Bluff remains ice covered until grounding-line retreat from the sill approximately 5–6 kyr after the main phase of thinning and retreat.

Figure 11 (a) Profile view: 15 kyr evolution of modelled upper ice surface along flowline for combined forcing experiments. Black bars represent elevation range from three sites in this study. Orange triangles represent offshore features marking observed grounding zone wedges. (b) Onshore elevation and grounding-line evolution. Top subpanel: forcings applied (no relationship between LBR and SIMR). In SIMR cases, we varied the SIMR perturbation within a 0.5 m/yr window at 500-year temporal frequency to simulate short-lived pulses of relatively warmer or cooler water mass changes. Top middle subpanel: time-transgressive elevation profile for Mt Kring with exposure ages; bottom middle subpanel: time-transgressive elevation profile for Hughes Bluff surface exposure ages; and lower subpanel: evolution of grounding-line position with modern grounding-line position, bathymetric sill at the mouth of David Fjord and mapped grounding zone features. We co-view exposure age data in these plots to allow generalised time-varying geometric relationships to be compared.

We are confident that our modelling results reasonably reconstruct the period of multi-millennial glacial change during the mid-Holocene. However, our approach does have limitations including imperfect knowledge of past boundary conditions, treating ice as isothermal, lack of isostasy and subglacial hydrology parameterization, and no use of a calving law to control ice shelf length. Through modern sensitivity experiments, geologic control on geometry, and a general fit to modern SIMR and basal stress, we remain confident our results support a first-order approximation for the dominant controls on the David Glacier's thinning history.

Table 3 Results for deglacial sensitivity experiments. We compare onshore and offshore modelled glacial behaviour, as well as thinning and grounding-line retreat rates. ¹ Modelled retreat rate is calculated during the deglacial phase of experiments. ² Data constrained maximum thinning rate at Mount Kring (MK) and Hughes Bluff (HB) is 0.19 and 2.06 m/yr, respectively. CI: Coulman Island.

5 Discussion

The thinning history of David Glacier places modern observations in a long-term perspective and allows for local, regional and continent-wide comparisons with other glacial histories and modelled ice sheet reconstructions. High-resolution, low-inheritance exposure ages obtained from Hughes Bluff and Mt Kring overlap in time reveal different thinning styles with Hughes Bluff recording rapid ice surface lowering and Mt Kring revealing a much slower thinning signal. The chronology from Mt Kring implies glacier thinning reached far inland along zones of streaming ice and provides rare constraints on ice behaviour from the margins of the EAIS.

The results of flowline modelling experiments along the expanded David Glacier reveal a threshold-driven sensitivity to both SIMR and ice stream lateral buttressing. Combined forcing experiments show a similar rate and duration of onshore thinning compared with our onshore geologic records constrained by surface exposure studies at two locations along the flowline. The results show the glacier response occurs at lower threshold values when the processes act in combination, consistent with previous applications of this model (Jamieson et al., 2014). Although the sensitivity experiments do not resolve which process or combination of processes forced the observed onshore thinning, we discuss potential explanations for differences in our observed vs. modelled onshore glacier thinning.

5.1 Coastal thinning and impacts on local oceanography

Continental-scale ice sheet model reconstructions and available geological data suggest that the glacier thinning profile at Hughes Bluff likely only captures the final 200 m of the approximately 400+ m thinning since the LGM (Fig. S2). In total, marine geological evidence, ice modelling and new results in this study suggest that at the LGM the expanded David glacier had an ice surface profile above Hughes Bluff and a grounding-line location near Coulman Island (Licht et al., 1996; Domack et al., 1999a; Shipp et al., 1999; McKay et al., 2008). Flowline modelling experiments starting from this ice extent and thickness capture both the magnitude and rate of onshore thinning derived from surface exposure data.

The new geological reconstruction of ice surface elevation changes at Hughes Bluff shows a rapid lowering of the David Glacier at 6.5 ka and a period of slow thinning from ∼ 6–4 ka. Given the marked slowdown in thinning rate from ∼ 6 ka at Hughes Bluff and its location downstream of the present-day grounding line (Fig. 6b), we suggest that stabilisation of the Drygalski Ice Tongue occurred after ∼ 6 ka. This finding is consistent with Orombelli et al. (1990) and Baroni and Hall (2004), who mapped a series of raised beaches along the TNB coastline that mark beach depositional processes in an open ocean setting (i.e. no grounded ice) initiating at 7.2 ka. Stevens et al. (2017) show that the modern Drygalski Ice Tongue is essential for the development of the modern TNB polynya. Thus, if the Drygalski Ice Tongue formed at ∼ 6 ka as our exposure chronology from Hughes Bluff and the raised beach chronology from TNB suggest, it is likely that the TNB polynya has also existed since this time.

5.2 Palaeo-thinning rate comparison

General agreement between palaeo-thinning rates derived from surface exposure ages of glacial deposits and those from our modelling experiments provides confidence in the approach used in this study. In Table 3, we show that in the combined forcing experiments we achieve a close match between data-driven and glacier modelled thinning rates. For Hughes Bluff, a maximum thinning rate of 2.06 m/yr compares well with the average modelled thinning rate of 2.27 ± 0.29 m/yr (experiments MS1–3). At Mt Kring, the maximum thinning rate is 0.19 m/yr, and, for experiments MS1–3, the average modelled thinning rate is 0.36 m/yr ± 0.02 m/yr.

Our modelling demonstrates that glacier thinning at Mt Kring is sensitive to grounding-line migration during ice retreat from the outer to the inner continental shelf. Mt Kring is part of a broad bedrock platform along the northern flank of the major glacial trough dissecting the TAM (Fig. S1). The platform is comprised of three high-elevation outcrops protruding ∼ 200 m above the ice, with Mt Kring nearest to the zone of streaming ice (e.g. ice surface velocities > 100 m/yr). Geophysical characterisation of the northern TAM suggests the presence of individual tectonic blocks bounded by faults which likely serve as zones of relatively weak rock strength allowing preferential ice flow and glacial erosion (Salvini and Storti, 1999; Dubbini et al., 2010). In West Antarctica, dynamic thinning of the inland ice sheet has been linked to underlying tectonic controls (Bingham et al., 2012). At David Glacier, similar tectonic controls may have conditioned the spatial pattern of dynamic thinning during the Holocene.

Along the upper reaches of the David Fjord, there are no outcrops in the zone of highest ice surface velocity (> 100 m/yr). Mount Kring lies 40 km away from the modelled flowline, in an area of relatively slow flowing ice. Modern ice surface velocities near Mt Kring are ∼ 15 % of the ice surface velocities at the projected flowline position (Rignot et al., 2011). On modern ice sheets experiencing dynamic thinning, satellite-derived thinning estimates are largest in the centre of the ice streams or outlet glaciers and become progressively smaller at lower velocity sites further from the central flowline. For example, the central parts of Greenland’s outlet glaciers are currently thinning at rates of ∼ 0.84 m/yr, while marginal areas with slower ice velocities are thinning at 0.12 m/yr (Pritchard et al., 2009). Therefore, the Mt Kring site likely reflects ice stream marginal thinning (maximum of 0.19 m/yr, derived from surface exposure data) rather than the greater thinning rate that was experienced in the centre of an ice stream (modelled average maximum 0.36 m/yr ± 0.02 m/yr). Taken together, the Hughes Bluff and Mt Kring chronologies suggest that ∼ 2 kyr of dynamic thinning occurred at David Glacier and that this thinning propagated significantly into the ice sheet interior.

5.3 Controls on thinning and grounding-line migration

Overall, the style of modelled grounding-line retreat is controlled by marine ice sheet instability, a positive feedback where grounding-line retreat into overdeepened subglacial basins leads to progressively enhanced ice discharge and ice sheet thinning (Weertman et al., 1974; Mercer, 1978; Schoof, 2007). Experiments designed to simulate glacier–ice-sheet decoupling show that once independent, the David Glacier grounding line rapidly retreats through the Drygalski Trough and stabilises on a bathymetric sill at the mouth of the David Fjord.

This study provides insights into the processes that occurred as a large grounded section of an ice sheet retreated into discrete outlet glaciers. Previous descriptions of this retreat have focused on grounded ice in the Ross Sea as a whole (Licht et al., 1996; Shipp et al., 1999; Domack et al., 1999a; McKay et al., 2008; Anderson et al., 2014). However, at the scale of David Glacier, retreat was likely influenced by local processes including interactions with adjacent glaciers and other ice bodies. For example, a recent ice sheet retreat reconstruction of Halberstadt et al. (2016) suggests that ice lingered on higher-elevation banks as the grounding line retreated in adjacent, large bathymetric troughs. Modelled grounding-line retreat may have been slower in reality due to lateral buttressing from outlet glaciers on the TAM coast and stagnant ice on bathymetric highs immediately seaward of the expanded David Glacier, and such lingering ice is not accounted for in our simplified modelling approach.

Figure 12 (a) Map of David Glacier, Terra Nova Bay and surrounding areas. Synthesis map focused on two-phase retreat between > 8 and 5.5 ka, including geographic and geomorphic features mentioned in text. Ice surface velocity of Rignot et al. (2011) and LIMA satellite data of Bindschadler et al. (2008). Inset with red box for main extent shows bathymetry of Arndt et al. (2013) and ice surface velocity of Rignot et al. (2011). (b) Synthesis profile focused on two-phase retreat between 11 and 5.5 ka. Grounding zone features: triangles with age constraints (orange), triangles without age constraints (green) of Lee (2019). GL: grounding line.

Knowledge of this complex lingering ice history is poorly constrained as the majority of research has focused on trough axes (Anderson et al., 2014; Halberstadt et al., 2016), whereas offshore Mawson Glacier, south of David Glacier, submarine geomorphology reveals a complex record of lingering ice (Stutz, 2012; Greenwood et al., 2018; Prothro et al., 2020). Analysis of marine sediments, primarily from trough axes, indicates that a calving bay environment formed during grounding-line retreat (Domack et al., 1999a, 2003; Leventer et al., 2006; Mackintosh et al., 2014). In this scenario, grounding-line retreat is rapid along the trough axis, while lingering ice remains grounded along the lateral margins of the ice stream/glacier. In addition to sedimentary evidence, the abundant large iceberg keel marks seaward of the Coulman Island grounding zone wedge (GZW) and smaller keel marks within TNB provide evidence for a calving bay during deglaciation. Short-term grounding-line stagnation may have been facilitated by complex interaction with other outlet glaciers (Reeves and Priestley glaciers together forming proto-Nansen Ice Sheet/Shelf) of TNB or grounded ice lingering on the banks surrounding the Drygalski Trough, or both.

During the first ∼ 500 model years, the modelled grounding line initially retreats to the location of a large grounding zone wedge (Lee, 2019). A lack of modelled grounding-line stability at other, smaller GZWs suggests that the grounding line may only have experienced short pauses in these locations. Further, in the case of the much smaller GZWs observed in the deepest portion of the Drygalski Trough, the elongate morphology of these smaller GZWs may reflect a point source versus a line or zone source associated with more classically defined sheet-like GZWs. Regardless, the small mounds in the trough axis are likely to reflect short-lived grounding line pauses during overall retreat.

Generally, once initiated, the modelled David Glacier grounding line retreats rapidly and pauses at a prominent sill at the mouth of David Fjord for up to 5 kyr before subsequent grounding-line retreat to its modern configuration. This simulated two-phase grounding-line retreat aligns with phases of enhanced (at ∼ 6.5 ka) and reduced thinning (from 6 ka) identified in our onshore reconstructions at Mt Kring and Hughes Bluff, both in terms of timing and rates of past glacier thinning. Models forced by moderate SIMR and LBR (MS1–3) fit best with onshore geologic constraints.

Figure 12 synthesises the results of the terrestrial thinning chronology, modelled glacier flowline behaviour and the existing regional marine retreat chronologies. This synthesis suggests that beginning at 7.5 ka, the grounding line unpins from the sill at the mouth of the David Fjord, the David Glacier and proto-Nansen Ice Shelf decouple, and widespread onshore thinning is initiated. Between 6–5.5 ka, the grounding line retreats to near its modern position, thinning slows significantly and open marine conditions prevail regionally. In summary, geological evidence from offshore David Glacier indicates that retreat of grounded ice through the Drygalski Trough and the formation of open marine conditions similar to today occurred immediately prior to the dynamic thinning of David Glacier recorded in this study (Licht et al., 1996; Domack et al., 1999b; McKay et al., 2008; Yokoyama et al., 2016). Our onshore surface exposure data record the final stages of glacial thinning and retreat along the David Glacier.

6 Conclusions

Ice surface reconstructions from the David Glacier reveal a period of glacier thinning along a large swath of the drainage basin during the mid-Holocene (∼ 6.5 ka). The reconstructed thinning style between two sites separated by ∼ 130 km reveals a dynamic thinning event that endured for two millennia. This chronology is synchronous with local-, regional- and continental-scale geological records of ice sheet change yet is not fully captured in continental-scale ice sheet models. Our flowline modelling results suggest that thinning and grounding-line retreat were driven by increased sub-ice-shelf melt (SIMR) and decreased ice stream lateral buttressing and that the combination of these two processes reduces the individual forcing thresholds required to initiate retreat. Modelled patterns of grounding-line retreat correlate well with patterns of onshore thinning, constrained by our high-resolution surface exposure ages. Data–model mismatches likely highlight processes or feedbacks not represented in our simplified modelling approach related to the relative role of local topographic pinning points and the nature of ocean-forced dynamic thinning.

Through careful collection of glacial deposits from numerous sites along the David Glacier, we have improved the spatial and temporal resolution of the onshore reconstruction of the AIS. Our data constrain a 130 km portion of one of the largest outlet glaciers in the world which carries regional significance for Victoria Land and the western Ross Sea, as well as offering clues about processes currently underway in rapidly changing sectors of Antarctica and Greenland. While we acknowledge the data and modelling presented in this work may not apply to other settings, we hope that our study may serve as a template for future work aiming to extend the observations of ice sheet change beyond the last ∼ 40 years of satellite data.