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. 5 Thinning profiles derived from 45 10Be and 3He surface exposure ages show David Glacier experienced rapid thinning of up to 2 m/yr during the mid-Holocene (∼6.5 ka). Thinning stabilised 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 ∼3 ka after the retreat of marine-based grounded ice in the western Ross Embayment. The timing and rapidity of the reconstructed thinning at David Glacier is similar to reconstructions in the 10 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 is initiated by interactions between enhanced sub-ice shelf melting and reduced lateral buttressing, leading to Marine Ice Sheet Instability. Such rapid glacier thinning events are not captured in continental or sectorscale numerical modelling 15

Abstract. 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 10 Be and 3 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 simulta-neous 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 subice-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

Introduction
Since the Last Glacial Maximum (LGM), ice sheets retreated in both hemispheres, causing a sea-level rise of ∼ 130 m . 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 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 . Discovered by the British National Antarctic Expedition (1901)(1902)(1903)(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 2 , 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.
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.
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 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, 14 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 14 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 10 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 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;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 . 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) 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;. 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 Jones et al., 2015;Whitehouse et al., 2017).

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 longlived nuclides such as 10 Be, this inherited component can be significant (Sugden et al., 2005;Balco, 2011). Inheritance of cosmogenic nuclides is widespread in Antarctica as the coldbased 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;. We use cosmogenic nuclides 10 Be and 3 He to determine surface exposure ages on 45 samples. For 10 Be analysis, quartz is the target mineral and is present in granites and sandstones. For 3 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 10 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 9 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 3 He, we targeted pyroxene from the 125-250 mum grain size fraction. Pure pyroxene was obtained using an established HF etching method, and 3 He / 4 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 10 Be or 3 He using methods outlined in Oberholzer et al. (2003Oberholzer et al. ( , 2008 and Di Nicola et al. (2009Nicola et al. ( , 2012. Exposure ages are calculated by converting the AMS-derived 10 Be / 9 Be ratio to 10 Be concentration by sub-tracting a known amount of 9 Be carrier added during chemical processing. Final exposure ages are calculated using 10 Be or 3 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;. For 10 Be, we employed  systematics (http://ice-tea.org/en/, last access: 21 October 2021), and for 3 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  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.

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 groundingline 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 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 where A is the temperature-dependent rate factor (for −20 • C), ρ i is ice density (0.910 g/cm 3 ) and g is gravitational acceleration (9.8 m/s 2 ).
As ice reaches flotation, ice is assumed to spread unidirectionally along the flowline, and ice flow is balanced by where and v is the vertically average effective viscosity defined as where δu δ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 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: .
(6) 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 palaeoice streams and outlet glaciers and is further supported by geomorphic interpretations using high-resolution bathymetry data Livingstone et al., 2012Livingstone et al., , 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 ) 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 groundingline 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 . 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 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.

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 megascale 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).
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 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 sed- iments 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 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.

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 highelevation 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 oceandriven 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 groundingline 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 groundingline 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  (Domack et al., 1999a(Domack et al., , 2003 and East Antarctic margin 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 500year 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.

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.

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 3 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 3 He exposure age of 550 ± 82 ka. The highestelevation 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 ( 10 Be, n = 5) (Fig. 4b). An older population spans 51-25 ka ( 10 Be, n = 8 and 3 He, n = 4). The oldest erratic age population shows scattered evidence of glacial behaviour between 123-63 ka ( 10 Be, n = 2, 3 He, n = 2). The older populations are likely an artefact of repeated burial by overriding ice and multiple periods of exposure (i.e. inheritance).

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. 10 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.

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 . In areas 400 m above the ice, bedrock samples contain weathering rinds and deep weathering pits filled with erratics 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).

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/) . These reconstructed glacier thinning rates are compared to modern thinning rates derived from satellite data and continental-scale ice sheet models . 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 Bentley et al., 2017), Amundsen Sea sector (Johnson et al., 2008(Johnson et al., , 2017(Johnson et al., , 2020, Marie Byrd Land (Stone et al., 2003)    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 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 oceanice-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, 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.

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 thick-ening (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 re- constructions 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., 2016Pollard et al., , 2017Pollard et al., , 2018Kingslake 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 datamodel 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).

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  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. 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 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.
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 groundingline retreat from the sill approximately 5-6 kyr after the main phase of thinning and retreat.
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. 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.

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. Highresolution, 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 . 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.

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 presentday 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.

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. 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.
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.

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.
Knowledge of this complex lingering ice history is poorly constrained as the majority of research has focused on trough axes 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(Domack et al., , 2003Leventer 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.  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.
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 shortlived 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.

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-, regionaland 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 highresolution 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 oceanforced 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.
Code availability. The code used for flowline modelling is available by request from the corresponding author.
Sample availability. Samples collected during the 2016/17 austral summer are curated at Te Herenga Waka Victoria University of Wellington and are available from the corresponding author. Samples collected prior to 2016/17 austral summer are curated at the University of Pisa.
Author contributions. JS contributed to project design, fieldwork planning and implementation, sample analysis, modelling, and preparation of manuscript. AM contributed to original project design, fieldwork, modelling and manuscript preparation. KN contributed to original project design, sample analysis and manuscript preparation. RW contributed to field and lab work. CB and MCS led previous regional fieldwork, contributed data and helped prepare manuscript. SC contributed to field and lab work. SJ led modelling work and contributed to preparation of manuscript. RSJ contributed to project design and modelling work. GB contributed to project design and BGC-related lab work. PS contributed to data-model comparison (DMC) work. JL contributed unpublished bathymetry data and regional marine geologic observations. YBS contributed to regional synthesis discussion. TW, LD, MI, FS, SIO and MC conducted AMS analyses at respective laboratories. LV contributed through preparation of structure-from-motion photogrammetry (SfMP) models, modelling and manuscript preparation. DL contributed to DMC work. RM contributed to regional marine geology synthesis.
Competing interests. The authors declare that they have no conflict of interest.
Disclaimer. Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. a great deal of gratitude for the many experienced researchers in the David Glacier area for providing field photos from past expeditions.
Financial support. Field work, AMS analyses and a research visit to BGC was funded by NZARI (grant 2017-1-3). Antarctic Science Bursary provided funding for a research visit to Durham University. Andrew Mackintosh was supported by ARC SRIEAS grant SR200100005, Securing Antarctica's Environmental Future.
Review statement. This paper was edited by Pippa Whitehouse and reviewed by Keir Nichols and James Lea.