Holocene thinning and grounding-line retreat of Darwin and Hatherton Glaciers, Antarctica

We present exposure ages of glacial deposits at three locations alongside Darwin Glacier and its tributary Hatherton Glacier that record several hundred meters of late Pleistocene to early Holocene thickening relative to present. As the 15 grounding-line of the Ross Sea Ice Sheet retreated rapidly southward, Hatherton Glacier thinned steadily from about 9 kyr BP until about 3 kyr BP. Our data are equivocal about the maximum thickness and mid-to-early Holocene history at the mouth of Darwin Glacier, allowing for two possible deglaciation scenarios: (1) ~500 m of steady thinning from 9 kyr BP to 3 kyr BP, similar to Hatherton Glacier, or (2) ~950 m of thinning, with a rapid pulse of ~600 m thinning at around 5 kyr BP. We test these two scenarios using a 1.5-dimensional flowband model of Darwin and Hatherton Glaciers, forced by ice-thickness 20 changes at the mouth of Darwin Glacier and evaluated by fit to the chronology of deposits at Hatherton Glacier. Our modeling shows that the constraints from Hatherton Glacier are consistent with the interpretation that the mouth of Darwin Glacier thinned steadily by ~500 m from 9 kyr BP to 3 kyr BP; rapid pulses of thinning at the mouth of Darwin Glacier are ruled out by the data at Hatherton Glacier. This contrasts with some of the available records from the mouths of other outlet glaciers in the Transantarctic Mountains, many of which thinned by hundreds of meters over roughly a one-thousand-year period in the 25 early Holocene. The deglaciation histories of Darwin and Hatherton Glaciers are best matched by a steady decrease in catchment area through the Holocene, suggesting that Byrd and/or Mulock glaciers may have captured roughly half of the catchment area of Darwin and Hatherton Glaciers during the last deglaciation. An ensemble of three-dimensional ice-sheet model simulations suggest that Darwin and Hatherton Glaciers are strongly buttressed by convergent flow with ice from neighboring Byrd and Mulock glaciers, and by lateral drag past Minna Bluff, which could have led to a pattern of retreat 30 distinct from other glaciers throughout the Transantarctic Mountains.


Introduction
Quaternary glacial deposits preserved in ice-free regions around the Darwin-Hatherton glacier system (DHGS) have been used to infer the glacial history of the region (Bockheim et al., 1989;Conway et al., 1999;Storey et al., 2010;Joy et al., 2014;King et al., 2020). In this paper, we present new geologic and exposure age results from Hatherton and Darwin Glaciers and use Bockheim et al. (1989) mapped lateral moraines and drift sheets in ice-free valleys alongside Hatherton Glacier, and dated 115 these deposits based on weathering, soil characteristics, and radiocarbon ages of freeze-dried algae. These algae grew in glacier-dammed ponds and died when the ponds disappeared as the glacier retreated downslope; the radiocarbon ages of the algae provide a proxy for the glacier margin position through time (Bockheim et al., 1989;King et al., 2020). Bockheim et al. (1989) distinguished five drift units: Hatherton, Britannia-I, Britannia-II, Danum, and Isca Drifts, in order from youngest to oldest. All subsequent work has adopted these names and confirmed the original mapping. Although we mapped and sampled 120 some material from the three oldest deposits, this paper focuses on the Britannia-I and Hatherton Drift.
The chronology and correlation of these deposits along the Hatherton-to-Darwin Glacier flowline remains uncertain, and the Britannia-II, Britannia-I and Hatherton Drifts have each been identified with the last maximum glaciation in this sector of the TAM (Bockheim et al., 1989;Storey et al., 2010;Joy et al., 2014;King et al., 2020). In addition to determining which deposits represent the last local maximum, projection of the associated ice surface downstream from sites on Hatherton Glacier 125 to the surface of the Ross Sea Ice Sheet has been complicated by the absence of comparable deposits on Diamond Hill at the glacier mouth. The ice surface elevation at the outlet of Darwin Glacier is key to understanding the dynamics of grounded LGM ice south of Ross Island, and is likely a pointer to the behavior of much larger Byrd and Mulock Glaciers, which dominate ice flow into this region.
Although Bockheim et al. (1989) tentatively attributed Britannia-II drift to the LGM, and Britannia-I deposits to a Holocene 130 still-stand or re-advance, Joy et al. (2014) have since dated Britannia-II drift to the penultimate glacial maximum. Younger and fresher Britannia I drift occurs up to ~350 m above the modern margin of Hatherton Glacier in Dubris and Bibra Valleys, Magnis Valley and at Lake Wellman (Bockheim et al., 1989;Storey et al., 2010;Joy et al., 2014;King et al., 2020). At these sites, the deposit is a boulder and cobble-rich ablation till, with clasts showing varying degrees of weathering. Large but discontinuous ice-cored moraines mark the limit of the deposit at Lake Wellman and in Magnis Valley, and separate it from 135 the older Britannia-II drift. At Dubris and Bibra Valleys the LGM limit varies from a thin scatter of fresh clasts to a low ridge of cobbles and boulders lacking an ice core. Clasts are distinctly less weathered than those in the underlying Britannia-II drift.
Similar recessional drift extends from the Britannia-I limit towards the present glacier margin, including low moraine ridges and numerous relict lake-ice boulder pavements (King et al., 2020).
The Hatherton drift is a minimally weathered deposit, commonly underlain by relict ice, and found within 50 -100 m 140 presence of an ice core and intermittent moraine along its margin suggest that the deposit marks a brief stillstand during glacial retreat, though there is no indication in the chronologies discussed below of a prolonged pause or readvance. 150 Given the absence of recognizable Britannia deposits at Diamond Hill, Bockheim et al. (1989) extrapolated their elevations from sites up-glacier to infer a thickening of 1100 m at the confluence of Darwin Glacier with the Ross Ice Sheet during the last glaciation; however, this was based on the assumption that the Britannia II drift represented the LGM limit, which Joy et al. (2014) showed to be incorrect. Subsequent numerical modeling experiments (Anderson et al., 2004) yielded a more modest 800 m of thickening relative to present, although they had to infer the bed topography by mass conservation. The inference by 155 Storey et al. (2010) that Hatherton rather than Britannia I drift marks the LGM high-stand of Hatherton Glacier led these authors to propose no more than ~50 m of LGM thickening. King et al. (2020) were able to date the advance and retreat of Hatherton Glacier in the Lake Wellman valley using radiocarbon dating of freeze-dried algae. They showed that the glacier advanced to the Britannia-I limit at 9.5 kyr BP, and argued that both the majority of their own exposure ages and those of Storey et al. (2010) were overestimates due to prior exposure of the samples. Joy et al. (2014)

dated deposits in Dubris and 160
Bibra Valleys alongside the upper Hatherton Glacier with cosmogenic 10 Be and 26 Al. They showed that Britannia II drift is in fact a deposit of the penultimate glaciation, and obtained Holocene ages similar to those of King et al. (2020)

Geochronological methods
We measured 10 Be and 26 Al in erratics we collected from Diamond Hill, as well as from Bibra, Dubris and Magnis Valleys.
We selected the freshest available samples, targeting rocks with differentially weathered top and bottom surfaces, interpreted as indicating that a sample had remained stable since deposition. Collection of isolated erratics rather than from moraines avoids samples which may have been disturbed, overturned, or temporarily covered during sublimation and collapse of ice in 175 the moraine core. This process can take thousands of years in Antarctica, where some glacial maximum moraines remain icecored up to the present day (including those in Magnis Valley and at Lake Wellman). On the floors of Dubris, Bibra and Magnis Valleys we collected samples from Britannia and Hatherton Drift. Where it was possible, we attempted to sample areas of thin or discontinuous drift, targeting rocks resting on underlying older, consolidated drift, though the thickness of Britannia-complication for these samples is coverage by proglacial lakes at the ice margin during retreat. The pervasive occurrence of sub-fossil algae (desiccated fragments of lacustrine cyanobacterial mats), discontinuous shorelines, and lake-floor pavements in both valleys indicates that retreating ice was fronted by one or more shallow proglacial lakes throughout deglaciation.
Although we sampled to avoid cover by lakes, some valley floor exposure ages may have been affected and therefore underestimate depositional ages. We also noted the relative height of nearby boulders and the dominant orientation of local 185 snow tails so as not to sample rocks that may have been covered by drift snow. In these wind-swept areas adjacent to blue ice glacier margins, it is unlikely that deep snow cover would persist for long enough to strongly skew the exposure ages.
We also measured cosmogenic 14 C produced in situ in quartz along an elevation transect of granitic bedrock from Diamond Hill to constrain the LGM ice thickness. We preferentially sampled very weathered bedrock to ensure minimal erosion over the last glacial cycle. Because any rock exposed during the few tens of kyrs prior to the LGM may contain inherited 14 C, 190 apparent 14 C exposure ages provide an upper bound on the timing of exposure since the LGM. However, the short half-life of 14 C makes these exposure ages less sensitive to exposure prior the LGM than 26 Al or 10 Be ages. Within the range of subaerial erosion rates typical of Antarctic bedrock (<1 m/Myr), 14 C concentrations reach secular equilibrium in ~30 kyr (Balco et al., 2016). Because of the short half-life of 14 C, a few thousand years of ice cover during the LGM will create a detectable signal of burial, while rock that was not covered by ice in the last 30 kyr will remain at its equilibrium concentration. Our goal at 195 Diamond Hill, where there are no deposits that indicate the upper limit of glacial maximum ice cover, was to bracket LGM ice elevation between the lowest sample that is saturated with respect to 14 C (i.e., not ice-covered during the LGM) and the highest unsaturated sample (i.e., ice-covered during the LGM).
We used conventional purification methods including heavy liquid separation, surfactant separation, and dilute HF etching at ~60°C (Kohl and Nishiizumi, 1992) to isolate pure quartz for exposure dating. Samples analyzed for in-situ 14 C were 200 prepared and measured twice. For the first set, quartz purified at UW was further treated at the Tulane lab with 1:1 concentrated HNO3:de-ionized water for 30 minutes at room temperature prior to loading into the carbon-extraction system. Following a 500°C pre-heating in vacuo, 14 C was extracted using methods described by Goehring et al. (2019). Five of the six samples in this first batch gave 14 C concentrations exceeding theoretical saturation concentrations (by up to ~75%). Realizing that this was due to contamination by modern carbon from residual surfactant-separation reagents that survived the procedures 205 described above, we re-analyzed all six samples using a more aggressive pre-treatment protocol, as described by Nichols and Goehring (2019). Results reported in this paper are from this second set of analyses; the results from the first set are reported by Nichols and Goehring (2019). As a further check on the validity of the second set of measurements, we also analyzed two erratic samples from Diamond Hill and Danum Platform, which gave exposure ages in agreement with 10 Be and 26 Al results.
Complete analytical results are given in the supplementary information. We isolated Be and Al for isotopic analysis at the 210 University of Washington Cosmogenic Nuclide Laboratory following the standard ion-exchange chromatography method of Ditchburn and Whitehead (1994). 14 C cathodes were prepared at Tulane University using the new fully-automated carbon extraction and graphitization system (Goehring et al., 2019a). We analyzed samples for in situ 10 Be at the Center for Accelerator Mass Spectrometry (CAMS) at Lawrence Livermore National Laboratory. As a further check, we analyzed 26 Al from the same aliquots of quartz at the Purdue Rare Isotope 215 Measurement Laboratory (PRIME). 14 C cathodes were analyzed at the National Ocean Sciences Accelerator Mass Spectrometry Laboratory and the Woods Hole Oceanographic Institution. The exposure ages presented in this paper have been calculated using the latitude-altitude scaling scheme of Stone (2000). While the choice of another scaling scheme changes the individual exposure ages, there is no major impact on the overall results, and we have therefore chosen to present results based on this scheme for simplicity. 220 Sample data and calculated ages for all samples in this study are found in the Supplemental Material and are archived at ice-d.org.

Dubris and Bibra Valleys
The Britannia I limit across Dubris and Bibra Valleys ( Figure 2) represents at least 370 m of thickening relative to the present 225 Hatherton Glacier surface, and dates to 8-10 kyr BP (Figure 3). Slight variations in the age of the Britannia I limit between different sites is to be expected, as changes in local meteorology could allow small-scale fluctuations of the ice margin. After 8 kyr BP, the glacier margin began to retreat steadily towards its present position. Only one exposure age exhibits significant prior exposure to cosmic rays (13-HAT-133-BV; 79 kyr). This sample was the closest to the glacier margin, and thus prevents us from better constraining the last 150 m of glacier thinning. The thinning profile on Danum Platform is steeper than in Dubris 230 and Bibra Valleys, which could be an effect of either the complex topography or lake cover in the valleys.  . The blue curve represents the Britannia I limit; erratic samples outboard of that limit are taken from the limit of the older Britannia II deposit. Contours and satellite imagery from Land Information New Zealand (LINZ, 2010). 235 Inset imagery from Bindschadler et al. (2008).

Magnis Valley
The Britannia I limit in Magnis Valley (Figure 4) predates the limit at Dubris and Bibra Valleys, with ages spanning 7.8 -245 13.9 kyr on the valley walls, and 8.3 -12.4 kyr on the valley floor ( Figure 5). A single pre-exposed sample at the drift limit dates to ~32 kyr BP. Ignoring this one date, if the ages at the upper end of these distributions are valid and not due to small degrees of prior exposure, they imply that the glacier margin was close to its limit in Magnis Valley for ~ 6 kyr. At Lake Wellman, 14 C ages of subfossil algae suggest that the glacier margin there advanced by ~ 500 m (50 m in altitude) from 13 -11 kyr BP (King et al., 2020). Differences between the apparent ages and durations of maximum advance at Bibra, Dubris and 250 Magnis Valleys and at Lake Wellman may reflect small-scale differences in microclimate between the sites, or may be an artifact of having dated too few samples at each site to constrain the full span of each maximum advance. Based on flowline modeling in Section 3 that suggests there should not be significant time lags between advances at these sites along Hatherton Glacier, we view the data cumulatively and infer that Hatherton Glacier was close to maximum thickness from at least ~14 kyr BP until ~8 kyr BP, although there is considerable variation in timing between sites. 255 Two erratics from a recessional deposit on the upglacier valley wall record 120 m of glacier thinning by 5.8 kyr BP. Erratics on the valley floor largely agree with this history, showing 100 m of thinning by 6.1 kyr BP. Two of the erratics on the valley floor evidently have had prior exposure to cosmic rays (14-HAT-045-MV and 14-HAT-059-MV). Based on the lack of dated erratics close to the ice margin, we are unable to determine the time at which the glacier margin stabilized at its present position in Magnis Valley. However, if thinning continued at the same rate of ~5 cm/yr, the glacier would have established its modern 260 margin position ~2 kyr BP, which is in agreement with the time at which the glacier margin stabilized at Lake Wellman.  Panel c shows ages on a log scale in order to show samples with prior exposure to cosmic rays.

Lake Wellman
The exposure age chronology at Lake Wellman suffers from numerous pre-exposed erratics, resulting in a wide spread of 275 ages (King et al., 2020;Storey et al., 2010) (). However, radiocarbon ages of algae from lakes and ponds show a strong dependence of age on elevation, with a maximum position at 9.5 kyr BP, followed by steady thinning to modern elevations by the late Holocene (King et al., 2020). Figure 6: Simplified map of Lake Wellman and sample locations from King et al. (2020). Contours and satellite imagery from Land Information New Zealand (LINZ, 2010). Inset imagery from Bindschadler et al. (2008).

Brown Hills and Diamond Glacier
The Brown Hills lie north of Diamond Hill adjacent to Diamond Glacier-a distributary lobe of Darwin Glacier-and are bordered by the Ross Ice Shelf to the east ( Figure 8). During glacial periods, Diamond Glacier would likely have crossed the 285 Brown Hills and connected with the Ross Ice Sheet. Exposure ages of five erratics from the Brown Hills range from 7-205 kyr BP (Figure 9). Despite sampling the freshest, differentially-weathered erratics available, four of these ages are older than the maximum advance at sites on Hatherton Glacier, and they are not ordered by elevation or distance from the Diamond Glacier tongue. We interpret these samples as previously exposed erratics deposited by ice receding from its last maximum position, though it is also possible that they are undisturbed remnants of one or more older deposits. The youngest date from 290 the Brown Hills also may be affected by inherited 10 Be, so we can only conclude that ice retreated from the saddle separating Diamond Glacier from the Ross Ice Sheet no earlier than 7 kyr BP.
https://doi.org/10.5194/tc-2020-356 Preprint. Discussion started: 15 December 2020 c Author(s) 2020. CC BY 4.0 License. Figure 7: a Surface exposure ages and algae radiocarbon ages at Lake Wellman after King et al. (2020). Cyan arrow indicates elevation of sill that would have had to be overridden by ice to create the pond that hosted the samples in the cyan circle. 295 Horizontal blue line shows elevation of present-day glacier margin. b All ages are shown on a logarithmic scale, in order to include a significant number of erratics with prior exposure to cosmic rays.

Diamond Hill
At Diamond Hill, at the mouth of Darwin Glacier, deposits are sparse and cannot be correlated with the drift succession described above. On the south slope of the mountain, above Darwin Glacier, fresh glacially-worked erratics are scattered over 300 polished and striated bedrock outcrops. These clasts, stranded by the final stages of recession of Darwin Glacier, extend up to ~ 135 m above its present-day surface. Above this, bedrock becomes significantly more weathered and surfaces are largely devoid of transported boulders and cobbles. The geomorphic transition from surfaces shaped by wet-based glaciation at low altitude to heavily weathered upper slopes resembles transects on other peaks in the TAMs and Marie Byrd Land (e.g. Sugden et al., 2005). We attribute the geomorphic transition to differential erosion during the LGM and earlier glaciations. Ice in 305 contact with the low-elevation bedrock was warm-based and erosive, but graded into thinner, cold-based ice cover at higher altitudes. A similar geomorphic gradient occurs in the Brown Hills and on the north slope of Diamond Hill, though bedrock is everywhere more weathered and there are very few lightly weathered erratics. Till deposits conceal much of the floor of the valley north of Diamond Hill, and fill the cirque below and west of the summit at ~800 m. Exhumation of lightly weathered clasts from these deposits may be a source of erratics with inherited nuclides in this area. The only relatively fresh, unweathered erratics that we found on Diamond Hill were perched on glacially sculpted bedrock domes overlooking Darwin Glacier, ~10 km upglacier from the modern grounding-line ( Figure 8). There was no clear limit of 320 deposition, but we did not find any relatively unweathered erratics higher than 135 m above the modern glacier margin. We analyzed eight of these erratics (Figure 9), which yield 10 Be ages spanning the latter half of the Holocene, from 5.2 ± 0.2 kyr BP at 135 m above the current glacier margin, to 0.3 ± 0.03 kyr BP at the current glacier margin. The rate of thinning matches the record from sub-fossil algae 14 C dates at Lake Wellman (King et al., 2020) (Figure 7), corroborating our inference that Hatherton Drift represents the last stages of the glacier's recession, separated from the youngest Britannia-I deposits by no more than a few hundred years of still-stand or readvance. Thinning during these final stages appears to have been relatively constant between 5.2 and 3.1 kyr BP, after which it slowed. and side of Diamond Hill (including summit samples) closer to Diamond Glacier. c All samples adjacent to Diamond Glacier on a logarithmic scale to show outliers. Bedrock 14 C ages give a maximum time of last exposure, because inherited 14 C from exposure before the LGM cannot be quantified. The top of Diamond Hill, >900 m above the modern glacier surface, was covered during the LGM and was exposed sometime after 11 kyr BP. In contrast, the 14 C-saturated sample ~500 m above 335 present-day Darwin Glacier suggests more modest thickening at the LGM relative to present. Therefore, it is not straightforward to interpret these data without a glacier flow model, and not all ages can conceivably represent the glacier surface.
The lack of fresh deposits above 400 m on Diamond Hill does not necessarily imply ice-free conditions during the last 340 glaciation. Based on the preservation of heavily weathered bedrock we presume that ice above this elevation was cold-based and debris-free, and thus did not deposit erratics as it retreated. Steep terrain and the extent of Diamond Hill prevented us from examining all locations, but we climbed to the summit of the mountain via the east (ice-shelf proximal) ridge, and descended on the northern (Brown Hills proximal) side, and we did not find evidence of a glacial maximum lateral moraine or recessional deposits above 400 m. Given both potential burial of surfaces by cold-based ice and the lack of erratics that would be suitable for dating with 10 Be, we measured in-situ-produced 14 C in bedrock to determine the maximum thickening and retreat history of Darwin Glacier at its downstream end at Diamond Hill.
On the slopes of Diamond Hill immediately above Darwin Glacier, our highest bedrock sample (14-HAT-026-DH; 472 m a.s.l.) gives an apparent 14 C exposure age of 6.7 ± 0.7 kyr BP, extending the retreat history from 10 Be data to 200 m above the modern glacier margin (Figure 9). Similarly, bedrock sampled <2 m above the current ice margin (14-HAT-033-DH; 280 m 350 a.s.l.), and adjacent to the 300-year-old ( 10 Be age) erratic (14-HAT-032-DH) gives an apparent 14 C exposure age of 500 ± 200 years BP. These ages confirm and extend the thinning chronology given by the 10 Be exposure ages. The exposure ages in this transect are too young to record, but do not preclude, a large early-Holocene deglaciation event, as found at Beardmore (Spector et al., 2017) and Mackay glaciers (Jones et al., 2015). However, those glaciers only thinned by tens of meters near the glacier mouth following a large drawdown of several hundred meters in the early Holocene, whereas Darwin Glacier thinned by 200 355 m since 6.7 kyr BP, a history more similar to that at Reedy Glacier . The contrasting thinning histories may be a result of the ice sheet bed topography near the mouths of these glaciers, which is currently poorly known.
About 500 m above the glacier margin (593 m a.s.l.) on the east ridge of Diamond Hill, the bedrock is at or near saturation concentration with respect to in situ 14 C production-decay systematics (14-HAT-006-DH). This sample was either not covered by ice during the last glaciation or was covered for only a brief period (< ~2 kyr). The maximum ice thickness at other locations 360 in the Ross Sea lasted for 3-5 kyr Spector et al., 2017), so this sample may be an upper bound on LGM ice surface elevation near the modern grounding line. This is a surprising result, given that the LGM ice surface was ~700-900 m above the modern ice surface at the mouths of other TAM outlet glaciers Bromley et al., 2010Bromley et al., , 2012Spector et al., 2017), and 600-700 m above the modern ice shelf at Minna Bluff and on Ross Island Denton and Marchant, 2000). It is at odds with both Bockheim et al.'s (1989) extrapolated LGM ice surface 1000 m 365 above present at the mouth of Darwin Glacier, and with Storey et al.'s (2010) assertion that Darwin Glacier did not thicken during the LGM. Anderson et al.'s (2004) estimate of an LGM ice surface 800 m above the modern is closer to our potential LGM surface, but still appears to be an overestimate.
In contrast, on the summit and north flank of Diamond Hill -above Diamond Glacier -in situ 14 C concentrations are well below saturation. Apparent exposure ages decrease with elevation from 10.8 ± 1.4 kyr at the summit of Diamond Hill 370 (1287 m a.s.l.) to 5.3 ± 0.5 kyr at 1134 m a.s.l. and 4.3 ± 0.4 kyr at 813 m a.s.l. (Figure 9b). In contrast to the ages from the Darwin Glacier side of Diamond Hill, these ages suggest that the summit was covered by ice or snow for a long period during the last glaciation. Today, Diamond Glacier terminates in a bedrock saddle between Diamond Hill and the Brown Hills below these samples at ~350 m a.s.l. While Darwin Glacier had thinned to within 135 m of its modern thickness by ~5.1 kyr BP, these ages seem to suggest that ice north of Diamond Hill was still at least 785 m thicker than present at this time, which is not 375 glaciologically plausible.
We hypothesize that the 14 C-saturated sample (14-HAT-006-DH; 598 m a.s.l.) and the other samples near Darwin Glacier on the south side of Diamond Hill (Figure 9a) are more representative of large-scale ice fluctuations than the undersaturated samples from higher elevations on Diamond Hill (Figure 9b,c) for three reasons: (i) there is no evidence of glacial erosion during the last glacial cycle for the high elevation samples, and therefore there are no mechanisms for producing older 380 cosmogenic 14 C ages at low elevation than at high elevation, except for local ice or snow cover at those higher elevations; (ii) two of the three highest-elevation cosmogenic 14 C ages on Diamond Hill are younger than the bedrock and erratics ages nearest Darwin Glacier, which cannot be explained if all samples are interpreted to represent changes of the large-scale glacier surface; (iii) the summit of Diamond Hill (1287 m a.s.l.) is higher than the local LGM deposits reported by King et al. (2020) at Lake Wellman (1130 -1270 m a.s.l.; Figure 6; Figure 7), ~50 km upglacier. This interpretation implies that the exposure ages higher 385 on Diamond Hill are the result of the disappearance of local ice or snow fields. However, given the limited number of in-situ 14 C exposure ages at Diamond Hill, these results are equivocal.

Summary of geochronologic constraints
We have presented geochronological data that constrain ice thickness along a flowline from near the head of Hatherton Glacier to near the modern Darwin Glacier grounding line, based on: (i) 10 Be and 26 Al exposure ages of deposits from Bibra, Dubris 390 and Magnis Valleys, (ii) radiocarbon dates of algae and 10 Be and 26 Al exposure ages from Lake Wellman (King et al., 2020), and (iii) 10 Be ages of erratics and in-situ 14 C exposure ages of bedrock from the south side of Diamond Hill adjacent to Darwin Glacier. Collectively these data suggest glacier advance and/or fluctuations at close to maximum thickness from ~14-8 kyr BP. At Lake Wellman, Magnis Valley, and Bibra/Dubris Valleys ice was at its maximum limit at 9.5 kyr BP, a time when large outlet glaciers in the southern TAM were thinning at their mouths Spector et al., 2017), and the 395 grounding line of the Ross Sea Ice Sheet was retreating southwards (Hall et al., 2013;Goehring et al., 2019b). Dates of maximum thickness are not demonstrably different between the three sites on Hatherton Glacier, which also differs from the diachronous behavior of upper and middle Reedy Glacier , although the sizes of these glaciers and the distance between sites are very different. Our data from Diamond Hill show that Darwin Glacier thinned steadily by 200 m since 6.7 kyr BP. However, the data are equivocal about the magnitude and rate of ice thickness changes near the grounding 400 line in the early-to mid-Holocene. Data from nearest Darwin Glacier on Diamond Hill suggest slow and steady thinning of ~500 m through the Holocene, while data at the summit and north side of Diamond Hill suggest ~900 m of rapid thinning in the early-to-mid Holocene. We prefer the former interpretation, but we do not have enough data to rule out the latter. To examine these aspects of the glaciers' history and evolution we use a flowband model to simulate how these deglaciation scenarios at the mouth of Darwin Glacier would be reflected in the chronologies from upstream on Hatherton Glacier. We also 405 test the alternative hypothesis that 500 m of LGM-to-present thinning included a rapid drawdown in the early Holocene, consistent with records from other TAM outlet glaciers (e.g. Jones et al., 2015;Spector et al., 2017), followed by slow and steady thinning recorded in the ages of the erratics. These simulations are presented in the following section. We modeled Darwin and Hatherton Glaciers since the LGM using a 1.5-D shallow ice glacier flowband model to evaluate possible deglaciation scenarios consistent with our geochronological data. Applying a simple ice-flow model to constrain iceflow conditions and deglaciation behavior has been done for other TAM outlets (e.g., Anderson et al., 2004;Golledge et al., 2014;Jones et al., 2016). The model we apply solves the mass conservation equation to calculate ice thickness evolution using the finite volume method. This model is computationally inexpensive and reduces the number of poorly constrained ice-flow 415 parameters and boundary conditions. The model domain starts near Diamond Hill at 10 km upstream from the modern grounding-line where we have geochronological data. We do not model grounding-line evolution, but we are able to use geochronological data to prescribe surface-elevation change at a location that is always upstream of the grounding line over the past 20 kyr. This avoids the use of the shallow ice approximation near the grounding line, where the inherent assumptions in the model do not apply. We are thus able to use data to constrain aspects of the model that would otherwise require more 420 complexity, but where a more complex model would introduce additional parameters that are poorly known or completely unknown.
Ice flow in a flowband (1.5-D) is described by the time-evolving mass conservation equation (Cuffey and Paterson, 2010): where H(x,t) is the ice thickness, q(x,t) is the volumetric ice flux, W(x) is the glacier width, and ̇( , ) is the surface mass balance. Mass balance effects of melting or freezing at the glacier bed are neglected due to a lack of data. The total ice velocity 425 (U) is the sum of contributions from internal deformation (Ud) and basal sliding (Us) taken to be of this form: where fd is the deformation factor, fs is the sliding factor, n=3 is the flow-law exponent, 1/3 ≤ m ≤ 4 is the sliding exponent, τd is the driving stress, and τb is the basal shear stress. Cross-sectional area is accounted for at each grid point. The modern fluxes from small tributary glaciers contributing to Darwin and Hatherton Glaciers are estimated using a flux gate calculation and we make the assumption that these tributary flux contributions are constant through time. 430

Model parameters, tuning, and boundary conditions
Surface mass balance and basal conditions of Darwin and Hatherton glaciers are poorly known. We use the RACMO2.1 5.5-km resolution model of the modern surface mass balance (Lenaerts et al., 2012a), as this is the only data product to include the significant surface ablation that is observed. While these ablation areas do not occur in the same locations observed from satellite imagery, we argue that this is the best choice of surface mass balance because it matches reasonably well the overall 435 flux out of the glacier system (Gillespie et al., 2017). Furthermore, the model is a function of the integral of the surface mass  Figures S1 and S2). 440 We tune a steady-state version of the flowband model to estimate poorly constrained ice-flow parameter values by minimizing the mismatch between the modeled and observed modern glacier surface elevation and the modern surface velocity at each model grid point. This strategy has been applied in related problems (e.g., Golledge et al., 2014). We vary the basal sliding factor (fs) and the ice-flow deformation factor (fd) defined at each model grid point within a plausible range of values.
The deformation factor and basal sliding parameter vary spatially over the model grid but we assume that the spatial patterns 445 that optimize the fit to modern data have not changed over the past 20 kyr because there is not enough information available to constrain a different choice for the parameter values in the past. In general, we have limited information available to constrain the values of the sliding and deformation factors, and while our inferred patterns of these parameters generate a surfaceelevation profile and surface-velocity profile that match modern values within their uncertainties, it is likely that it is not a unique solution. However, we are using this simple model to constrain the time evolution of the glacier system that is consistent 450 with our geochronological observations, and not to interpret inferred parameter values and patterns of controls on ice flow. Hatherton Glaciers, these catchment boundaries may have changed in the past. Thus, we evaluate the influence of time-varying flux that enters the domain from the upstream boundary, as required by a flux-balance calculation. We are also able to modulate flux entering the glacier through the upstream boundary as a way of investigating the role of catchment geometry changes 460 (including ice divide migration) on the evolution of the glacier profiles.

Model evaluation
In order to compare the model output with our chronologies from each location, we must account for the fact that our geochronology data record the elevation of the glacier margin through time, rather than the glacier centerline that is calculated in the flowband model. Thus, the exposure and radiocarbon ages represent a minimum elevation for the glacier centerline 465 through time. While the glacier centerline generally lies 100 m above its margin today alongside our study sites, this would not necessarily have been the case at the LGM. We can determine the possible range of glacier centerline elevations at the LGM by calculating a linear best-fit curve to the LGM deposits in the valley floor (far from the modern glacier) and on the valley walls (close to the modern glacier). By extrapolating this linear fit out to the center of the glacier, we obtain a maximum https://doi.org/10.5194/tc-2020-356 Preprint. Discussion started: 15 December 2020 c Author(s) 2020. CC BY 4.0 License. constraint on the LGM glacier centerline elevation because the glacier should not have a concave-up profile in the transverse 470 direction. We use the elevation of the highest LGM erratics as a minimum constraint because the glacier margin can never be higher than the centerline. In the absence of further constraints, we use the mean of these elevations as the estimated LGM centerline elevation, and the span as the 2-sigma uncertainty. Because the age of the LGM deposit on the valley walls does not necessarily match the age of the corresponding limit on the valley floor, these are slightly time-transgressive estimates.
However, we consider this the best estimate possible without the use of a higher-dimensional glacier model that requires 475 additional assumptions. We obtain glacier centerline elevations and 2-sigma uncertainties of 1355 ± 67 m at Lake Wellman, 1402 ± 56 m at Magnis Valley, and 1550 ± 45 m at Dubris-Bibra Valleys. Elevations from recessional deposits were then projected to estimate the centerline elevation by determining the percentage of total thinning each sample elevation represents and multiplying this by the height of the LGM surface above the modern glacier elevation. Model results are compared against these projected elevations in Figure 10, Figure 11, and Figure 12. 480 At Diamond Hill, we lack the clear limit of LGM deposits necessary to project samples to a glacier centerline elevation.
Instead, we convert the elevations of bedrock and glacial erratics samples to a height above the nearest ice margin and use those values as constraints for glacier fluctuations and for the experiments discussed in Section 3.4. Likewise, we do not have a reliable means of assigning uncertainty to these glacier centerline elevations, which may be large. However, this uncertainty is greatly outweighed by the uncertainty in what samples represent the elevation of Darwin Glacier versus local snow and ice 485 fields.

Transient Experiments
Inputs for all experiments in this section are summarized in Table 1.

Modern surface mass balance, constant catchment size 490
Since we have no constraint on changes in accumulation and ablation patterns through time for the DHGS, we first set up a transient run in which the modern surface mass balance pattern is kept constant in time. The glacier surface at the outlet of Darwin Glacier is held at 500 m above the modern surface until 9 kyr BP to match the onset of thinning at Hatherton Glacier, and then lowered according to the geochronologic constraints from Diamond Hill (Figure 9), using the assumption that the 14 C-saturated bedrock sample represents an upper bound on the LGM ice surface of the main trunk of Darwin Glacier. This is 495 based on the interpretation that the Holocene 14 C exposure ages higher up on Diamond Hill represent cover by local ice or snow fields. The results of this experiment are shown in Figure 10   There is essentially no lag time between the application of the elevation-change at the mouth of Darwin and the response of Hatherton Glacier in this scenario. Thus, the onset of thinning at each location along the profile of Hatherton Glacier occurs 505 almost simultaneously with the onset of thinning at Diamond Hill. This is consistent with our chronologies from the valleys alongside Hatherton Glacier, each of which indicates that retreat began ~8-9 kyr BP. However, using the modern surface mass balance at the LGM leads to a thinner modeled Hatherton Glacier than is indicated by the glacial deposits. This discrepancy is more pronounced further upglacier: at Lake Wellman (close to the confluence of Darwin and Hatherton), the model predicts an LGM ice surface ~75 m below the LGM erratics at the entrance to the valley; at Magnis Valley, the model underpredicts 510 LGM thickness by ~100 m; at Dubris and Bibra Valleys, the modeled glacier is ~170 m too thin.
https://doi.org/10.5194/tc-2020-356 Preprint. Discussion started: 15 December 2020 c Author(s) 2020. CC BY 4.0 License. Figure 10: Flow-band model results from Experiment 1 compared with surface exposure ages from this study and surface exposure and radiocarbon ages from King et al. (2020), which have been projected to centerline elevations using the method described in Section 3.3. The two best-fitting scenarios require either a larger catchment for both glaciers at the LGM (green 515 curve; DH_run_1l), or an accumulation rate far higher than modern (purple curve; DH_run_1d). The black curve in the Diamond Hill panel represents the prescribed ice surface history at the downstream boundary.
There are three possible explanations for the underprediction of LGM ice thickness in this experiment. First, it is possible that the current size of the glacier catchments is smaller than it was during the LGM. Observations of striated bedrock near the 520 head of Hatherton Glacier indicate that at some unknown time in the past, ice likely flowed over into the Hatherton valley from what is now the catchment of Byrd Glacier (Bockheim et al., 1989). Likewise, there are no constraints on the former size of https://doi.org/10.5194/tc-2020-356 Preprint. Discussion started: 15 December 2020 c Author(s) 2020. CC BY 4.0 License.
the Darwin Glacier catchment. It is thus possible that more ice fed into the DHGS due to larger catchment size during the LGM. Second, our assumption of a constant pattern of surface mass balance through time is too simple. The large ablation areas found on both glaciers today are likely near their maximum extent, and a small climate perturbation would decrease their 525 size, leading to a more positive surface mass balance (Brown and Scambos, 2004). However, Hatherton Glacier must have had ablation areas (both wind scouring and surface melt) at its margins during the local LGM in order to provide the meltwater necessary to create the ponds that hosted algae, and in order to deposit erratics in the valleys. Furthermore, accumulation rates tend to increase as the atmosphere warms during a glacial termination; the Taylor Dome ice core record shows a roughly 100% increase in accumulation rate between 12 kyr BP and 0.7 kyr BP (Monnin et al., 2004). The presence of blue ice areas due to 530 scouring by katabatic winds would likely lead to a more complicated change in surface mass balance, but the reduction in snow accumulation during glacial periods would likely have been a common feature throughout the Transantarctic Mountains.
Third, our records at the mouth of Darwin Glacier (Figure 9) are equivocal as discussed above. The fact that the complicated and spatially sparse bedrock 14 C record is difficult to interpret leaves open the possibility that the average LGM ice thickness at Diamond Hill was higher than we have tentatively concluded. We explore the effects of time-evolving surface mass balance, 535 variable flux into the glacier canyons, and different LGM ice thickness at the mouth of Darwin Glacier in the following model experiments.

Time-evolving surface mass balance
We first examined the effect of scaling the modern surface mass balance to the normalized accumulation rate history 540 recorded in the Taylor Dome ice core (Steig et al., 2000;Monnin et al., 2004). As expected, this also leads to a modeled Hatherton Glacier that is too thin at the LGM, with 150-250 m of thickening relative to modern. This gives a very poor fit to all of our glacial geologic data.
Next, we used a simpler time-varying surface mass balance pattern in an attempt to determine the magnitude of glacialinterglacial change needed to match the records from Hatherton Glacier using only variations in SMB. The modern surface 545 mass balance is taken from RACMO 2.3. We defined the LGM surface mass balance by multiplying the RACMO 2.3 surface mass balance -which does not include ablation areas -by a scaling factor. We use values of 60%, 100%, and 200%. The choice of 60% is made to match the accumulation rate change in the Vostok ice core between the LGM and present (Petit et al., 1999). This is of course a crude estimation, but it results in reasonable LGM accumulation rates of 2-10 cm/yr. The surface mass balance is varied linearly in time between the modern and LGM states. After a 5 kyr spin-up to allow the model to 550 equilibrate with LGM climate, the LGM SMB is held constant until 15 kyr, and then varied linearly to the modern SMB.
In this scenario, 60% scaling leads to underprediction of LGM ice thickness by about 100 m at all three Hatherton Glacier sites (DH_run_1b in Figure 10). Using a 200% scaling of RACMO2.3 SMB at the LGM leads to excellent agreement between the glacier model and most of the records from Hatherton Glacier (DH_run_1d in Figure 10). LGM surface elevations are matched to within a few tens of meters at all three locations, the timing and rate of thinning agree very well, and the modern 555 surface elevations are reproduced to within 40 m at all three locations.
While the 200% scaling results in a reasonable match between the model and surface-exposure age data, it produces an unlikely surface mass balance history. This requires LGM accumulation rates of 7-30 cm/yr, compared to ~3 cm/yr accumulation at Taylor Dome during the same time period. The fact that at least some surface ablation had to occur to create ice-marginal ponds means total surface mass balance may have been even lower than at Taylor Dome. This scenario also 560 requires the overall surface mass balance to decrease during the termination of the glacial period, which is unlikely. Although this scenario fits our data quite well, we search for a more reasonable explanation.
Added flux to account for changing catchment area When additional flux is added only to Hatherton Glacier, the glacier profile tends to steepen, leading to overprediction of 575 ice thickness by >100 m at Dubris-Bibra and Magnis valleys, while ice thickness is still ~75 m too low at Lake Wellman.
Conversely, when a proportional amount of flux is added to the upstream boundary of the Darwin Glacier domain, LGM ice thickness is underpredicted at all locations on Hatherton Glacier by 75-100 m, while the head of Darwin Glacier is ~200 m thicker than Bockheim et al. (1989) found at Darwin Nunatak. However, by adding a more modest amount of flux to both Darwin and Hatherton Glaciers, and by using the scaled SMB of 60% modern for the LGM, we achieve a better fit to the LGM 580 limits and retreat histories at all three locations along Hatherton Glacier, without unduly thickening the head of Darwin Glacier (DH_run_1l in Figure 10). Given the uncertainties in bed properties, accumulation rate, and tributary fluxes through time we consider this to be a satisfactory fit to the data. If the catchment boundaries stabilized before the grounding line of Darwin Glacier stopped retreating, this could also explain the slowdown of thinning at Dubris-Bibra Valleys several kyr prior to the slowdown of thinning at Lake Wellman, ~30 km closer to the grounding line. 585 We now consider whether the amount of additional incoming flux required to match the LGM limits is a physically reasonable quantity. For the additional 6.4 x 10 7 m 3 /yr added to Hatherton and the additional 1.2 x 10 8 m 3 /yr added to Darwin and a reasonable LGM accumulation rate of 3 cm/yr over the East Antarctic Plateau, this would require 2,100 km 2 and 4,100 km 2 of additional LGM catchment area for Hatherton and Darwin glaciers, respectively. This represents a ~75% increase over the modern catchment area (Gillespie et al., 2017), but only a 5% decrease in the Mulock Glacier catchment or a 0.5% decrease in the Byrd Glacier catchment areas based on values calculated by Stearns (2011). Because this total areal change is roughly equivalent to the reported uncertainty in the Mulock Glacier catchment area, and only about 11% of the uncertainty in the Byrd Glacier catchment area, we consider this amount of divide migration over a glacial-interglacial cycle to be reasonable.

Experiment 2: Rapid early-Holocene thinning at the mouth of Darwin Glacier
This experiment, like Experiment 1, assumes the 14 C-saturated bedrock at Diamond Hill represents the upper limit of the 595 LGM ice surface. Other records from the western Ross Embayment show a rapid drawdown event 9-8 kyr BP, presumably indicating widespread deglaciation of the region in the early Holocene (Spector et al., 2017). There is no record of such abrupt thinning at Hatherton Glacier, but such an event could have occurred in the data gap between the 14 C-saturated and the 6.7 kyr bedrock (separated by 300 m of elevation) at the mouth of Darwin Glacier. We explore this possibility by imposing a rapid thinning event of 275 m (similar to thinning at Beardmore and Mackay glaciers (Jones et al., 2015;Spector et al., 2017)) from 600 9-8 kyr BP, followed by gradual thinning consistent with the glacial geologic constraints. Surface mass balance and additional flux are the same as the best-fit scenario from Experiment 1 (i.e., 75% larger catchment, SMB varies linearly through time).
The results of this experiment are shown Figure 11.
The rapid thinning imposed at the mouth of Darwin Glacier propagates upglacier in the model with no significant lag.
While the amplitude of the signal decays with distance upglacier, it is still readily detectable in the modeled glacier changes 605 at Dubris-Bibra Valleys. This leads to a much poorer fit to the data than the situation with gradual thinning at the Darwin Glacier mouth in Experiment 1. This shows that the records from Hatherton Glacier are not consistent with an episode of rapid thinning at the mouth of Darwin Glacier, which implies that the DHGS did not respond to the last deglaciation in the same way as other outlet glaciers to the north and south (e.g., Jones et al., 2015;Spector et al., 2017). https://doi.org/10.5194/tc-2020-356 Preprint. Discussion started: 15 December 2020 c Author(s) 2020. CC BY 4.0 License. Figure 11: Flow-band model results for Experiment 2. All inputs are the same as for Experiment 1, except for the downstream boundary condition, which contains a period of rapid deglaciation 9-8 kyr BP. It is evident that this event would have been recorded in the deposits at Hatherton Glacier, but there is no clear evidence of this in our chronologies. Thus, the case of slow thinning in Experiment 1 provides a better fit to the data.

Experiment 3: 950 m of thinning at the mouth of Darwin Glacier, with a rapid pulse of thinning in the mid-Holocene
This experiment, unlike Experiments 1 and 2, assumes the 14 C-saturated bedrock at Diamond Hill is merely an outlier, and that the exposure ages at higher elevations represent the fluctuations of Darwin Glacier. While we have argued that the higherelevation ages are not as representative of the surface of Darwin Glacier as they are of more local ice configuration, we do not 620 have enough data to prove this assertion. Thus, in this experiment we evaluate a scenario in which Darwin Glacier was 1 km thicker than present at the LGM and covered the top of Diamond Hill, after which it thinned by ~800 m from 6-4 kyr BP. This thinning history is defined by the bedrock samples in Figure 9b. Results are shown in Figure 12. Figure 12: Flow-band model results for Experiment 3, in which we assume that the Holocene exposure ages high on Diamond 625 Hill represent thinning of the main trunk of Darwin Glacier. Fitting the LGM thickness of Hatherton Glacier requires a much smaller Darwin Glacier catchment at the LGM, but we cannot fit the shape of the chronologies with any combination of simple assumptions about surface mass balance or catchment area. Thus, we rule out this scenario as unlikely and conclude that the young exposure ages at high elevation on Diamond Hill may reflect more local ice fluctuations.

630
For RACMO2.1 accumulation rates and modern catchment boundaries, this deglaciation scenario leads to a modeled Hatherton Glacier that is ~160 -200 m too thick at all locations on Hatherton Glacier (DH_run_3a in Figure 12). Scaling the RACMO2.1 accumulation rate to the Taylor Dome record, which is our lowest LGM accumulation rate, produces essentially no change in the glacier profiles, while scaling the RACMO 2.3 SMB by 60% leads to a worse fit due to the lack of ablation https://doi.org/10.5194/tc-2020-356 Preprint. Discussion started: 15 December 2020 c Author(s) 2020. CC BY 4.0 License. zones (DH_run_3c in Figure 12). The fit to Hatherton Glacier LGM ice thickness can be improved by moving the Darwin 635 Glacier catchment boundary ~25 km into the model domain (i.e., decreasing the catchment area; DH_run_3i in Figure 12).
There is no geologic evidence for this, but it may not be an unreasonable amount of change, given the enormous catchment areas of Byrd and Mulock glaciers. However, the rapid thinning during deglaciation propagates upglacier in all three model runs, leading to a very poor fit to the records of Hatherton Glacier fluctuations. Thus, the pattern of thinning cannot be matched using simple assumptions about catchment size or surface mass balance, and thus we conclude that this is an unlikely scenario. 640 This supports our earlier hypothesis that the elevation transect near the modern margin of Darwin Glacier at Diamond Hill

Ice sheet model ensemble
Our flowband model is only applied to the grounded portions of Darwin and Hatherton Glaciers, and thus cannot be used 645 to directly examine the effect of grounding-line retreat on the ice thickness of Darwin and Hatherton glaciers. To investigate regional controls we use the Pennsylvania State University 3-D ice sheet model (PSUICE) (e.g., Pollard and DeConto, 2012a) to examine the LGM-to-present retreat of Byrd, Darwin-Hatherton, and Mulock glaciers. PSUICE uses a combination of the shallow ice and shallow shelf approximations (SIA and SSA, respectively) along with a parameterization of grounding-line flux (Schoof, 2007) to allow for full continent-scale simulations on kyr to Myr timescales, as well as an optimized model-650 parameter set that is tuned to match glacial geologic data from around West Antarctica over the past 20 kyr (Pollard et al., 2016). However, because we are only interested here in the ice sheet in the Ross Sea, this parameter set tuned at the continent scale may not be the most suitable for our regional investigation. So, we setup an ensemble of model runs to investigate the sensitivity to the choice of key parameter values. First, we run the model at 20-km resolution over the whole continent from 25 ka to present, and use this to establish boundary conditions for an ensemble of 48 nested model simulations at 10 km 655 resolution since 20 ka in which we vary model parameters around the values optimized from Pollard et al. (2016). While this resolution may not be sufficient to accurately represent the dynamics of Darwin and Hatherton glaciers, we use these model runs to investigate regional deglaciation. The modern ice discharge from those glaciers (~0.2 Gt/yr at present; Gillespie et al., 2017) is small compared to that from Byrd and Mulock glaciers (~27.5 Gt/yr at present; Stearns, 2011), so DHGS does not affect large-scale and long-time evolution of the Ross Sea Ice Sheet. We are interested in the relationship between ice thickness 660 at the TAM front and the position of the Ross Ice Sheet grounding-line, and thus do not need to model the individual glaciers at high resolution. The 10-km resolution was chosen as a trade-off between computing time and the ability to resolve pinning points that may have a large effect on grounding-line evolution. We vary four model parameters that were examined by Pollard et al. (2016): basal traction on the modern seafloor, isostatic rebound rate, ice shelf melt sensitivity to ocean temperatures, and a calving rate factor. We also use two different sea-level curves to explore the dependence of the model on time resolution in 665 proxy records (Lisiecki and Raymo, 2005;Spratt and Lisiecki, 2016). Parameter values used in the ensemble are listed in Table  2. Because our results show more gradual and recent deglaciation of the DHGS than at other locations in the TAM (Spector et al., 2017), our choice of parameter values is intended to slow down grounding-line retreat relative to the optimized parameter set of Pollard et al. (2016). Sub-shelf melting factor 0.5, 1 1, 3

670
Basal sliding coefficient on seafloor 10 -5 , 10 -6 , 10 -7 10 -5 Calving factor 0.7, 1 1 Sea-level curve Lisiecki and Raymo, 2005;Spratt and Lisiecki, 2016 Only used Lisiecki and Raymo, 2005 We investigate the relationship between rates of grounding-line retreat and ice thickness change at the mouth of Darwin Glacier. Ideally, we seek a combination of model parameters that enable the grounding line to slowly approach the mouth of Darwin Glacier, despite the steeply back-sloping bed topography south of Minna Bluff, while still leading to early Holocene 675 deglaciation of much of the western Ross Embayment (Spector et al., 2017). However, we also recognize that the parameters could vary in space and time, so this solution could be non-unique. Furthermore, higher-order ice physics and adaptive mesh refinement to resolve grounding-line evolution may be required to fully capture the processes involved in the retreat but require too much computational expense to be included in an ensemble of model runs over timescales of tens of thousands of years.
Results from our ice sheet model ensemble are shown in Figure 13. We find that no combination of the parameters we 680 explored reproduces the slow and steady drawdown through the Holocene that we observe in the glacial geologic data at Diamond Hill or at other locations along Hatherton Glacier. The basal sliding coefficient leads to a difference of ~600 m ice thickness between the fastest sliding value (10 -5 ) and the slowest (10 -7 ). While the ice at the mouth of Darwin Glacier begins thinning ~1 kyr earlier for the faster sliding scenarios, modern ice thickness is achieved by 5-6 kyr BP in all model runs, which is much earlier than our glacial chronologies show (~2-3 kyr BP). While none of the model runs in this ensemble achieve what 685 we consider to be a good fit to our data from Diamond Hill, the three clusters of model runs defined by the three basal sliding coefficients display internally consistent behavior characteristics that are worth examining in order to better understand the sensitivity of the Ross Sea Ice Sheet near the mouth of Darwin Glacier. rate of grounding-line retreat, we refrain from interpreting our glacier chronologies in terms of the rate of grounding-line retreat. However, it is still valid to interpret the cessation of major changes in glacier thickness ~3 kyr BP as a sign that the grounding line had reached close to its modern position at the mouth of Darwin Glacier.

Discussion
Darwin and Hatherton glaciers continued to adjust to regional deglaciation until ~3 kyr BP, which means that the grounding 720 line likely did not arrive at its present location until about that time. Darwin and Hatherton Glaciers lie roughly halfway between McMurdo Sound and Beardmore Glacier, which both deglaciated in the early Holocene Spector et al., 2017), yet the timing of grounding-line arrival at Darwin Glacier is apparently much more recent. This could imply that a large region extending across Byrd, Darwin, and Mulock glaciers remained grounded for about 4 kyr after the grounding-line in the central Ross Embayment had retreated further south. It remains an open question how far along the TAM front this 725 grounded ice persisted, and if it comprised a single, grounded ice mass or local piedmont lobes (Lee et al., 2017).
It is unclear if the chronology from Skelton Glacier supports this interpretation, as there are no erratics near the grounding line that date to the last deglaciation (Jones et al., 2017). A small number of erratics at locations far from the modern grounding line suggest that most thinning in the Skelton Névé had occurred by 6 kyr BP (Jones et al., 2017;Anderson et al., 2020), but the number of samples is too small to determine the rate of thinning or the timing of its onset. More data from Skelton Glacier 730 are needed for thorough comparison with the chronology at Darwin and Hatherton glacier.
Steady thinning of Hatherton Glacier through the Holocene supports a later arrival of the grounding line relative to glaciers to the north and south (Hall et al., 2015;Jones et al., 2015;Spector et al., 2017). In our flowband model, thinning at the mouth of Darwin Glacier propagates rapidly upglacier; thus, the deposits alongside Hatherton Glacier can be interpreted as recording the relative rates of ice thinning at the mouth of Darwin Glacier. There is no record of an exceptionally fast period of thinning 735 at Hatherton Glacier, in agreement with our interpretation of the sparse data near the modern grounding line. However, our ice sheet model ensemble analysis suggests that periods of exceptionally fast grounding-line retreat -such as across the Discovery Deep region -would not necessarily correspond to noticeably faster ice thinning rates at the mouth of Darwin Glacier. A two-to four-fold increase in the modeled rate of grounding-line retreat as ice in Discovery Deep goes afloat does not cause faster drawdown at Darwin Glacier for most of the parameter values we investigated. In the scenarios with the fastest 740 basal sliding, there is an increase in the rate of thinning when the grounding-line retreats into Discovery Deep. However, this lags the change in the grounding-line retreat rate and reaches a maximum thinning rate after the grounding line has already slowed. This suggests that the time at which the grounding-line arrived at the glacier mouth may be deducible from geologic data, but that there are significant complications in drawing conclusions about relative rates of grounding-line retreat from glacier-elevation proxies. 745 While we are not able to determine the rate of grounding-line retreat, we can make a first-order estimate of the sensitivity of Darwin and Hatherton glaciers to grounding-line position using records from elsewhere in the Ross Sea. The groundingline retreated past Ross Island >8.6 kyr BP (McKay et al., 2016), and McMurdo Sound became free of grounded ice between 7.5 and 9 kyr BP Jones et al., 2015;Anderson et al., 2017;Christ and Bierman, 2020). This is similar to the time at which Hatherton Glacier began to retreat from its last high-stand (8 -9 kyr BP), and our glacier modeling suggests 750 Darwin likely began thinning around this same time as well. Thus, Darwin and Hatherton glaciers may be strongly dependent on the ice configuration near Ross Island. The results of our ice sheet model ensemble suggest that this relationship is complex, as different parameter choices lead to opposing behavior as the grounding-line retreat past Cape Crozier. This discrepancy should be addressed with further modeling studies of outlet glacier response to grounding-line retreat and/or more data from near the modern grounding line of Darwin Glacier. The more rapid thinning of Mackay and Beardmore glaciers could be due 755 to the steeply back-sloping bed topography in the immediate vicinity of the glacier mouths (Jones et al., 2015;Spector et al., 2017).
Our estimated LGM ice surface elevation of ~600 m asl at the mouth of Darwin Glacier from our flowband modeling needs to be reconciled with the elevations of the LGM limit of 637 m asl at Minna Bluff (Denton and Marchant, 2000) and 520 m asl at Mt. Discovery . The flowband of Darwin Glacier likely would have passed between these sites at 760 the LGM (Denton and Hughes, 2000), and with our inferred ice surface elevation the glacier surface profile would have been very flat at the LGM, resulting in lower driving stresses than we would expect in comparison to modern glaciers and ice streams (c.f. Cuffey and Paterson, 2010, pg 297). We identify three factors that could reconcile these elevation differences and lead to a sufficient LGM surface slope between Darwin Glacier and Minna Bluff-to-Mt. Discovery.
First, there is significant uncertainty in the large-scale ice surface elevation that corresponds to the 14 C-saturated bedrock 765 elevation at Diamond Hill. Ablation by katabatic winds can cause large difference in surface elevations around a nunatak (Bintanja, 1999). The downglacier side of Diamond Hill may have been subject to especially vigorous wind-driven ablation at the LGM, and thus could have been significantly lower than the large-scale ice surface. As we have shown, the fit of our flowband model results to the data at Hatherton Glacier is more sensitive to rates of thickness change at the glacier mouth than to the absolute LGM ice thickness. Therefore, the possibility that the glacier centerline LGM ice surface could have been 770 significantly higher (perhaps ~100-200 m) than the 14 C-saturated bedrock sample at Diamond Hill is not in necessarily in conflict with our preferred deglaciation scenario, as long as thinning was smooth, steady, and initiated around 9 kyr BP.
Second, none of the reported LGM surface elevations include the effects of glacial isostatic adjustment, which could vary widely between Darwin Glacier and Minna Bluff. Geothermal flux south of Ross Island may be 58-125% higher than at Diamond Hill (Morin et al., 2010;An et al., 2015). If the crustal densities around Minna Bluff and Ross Island are 775 correspondingly lower than at the mouth of Darwin Glacier, the loading of the crust by grounded ice would have caused greater isostatic depression at Minna Bluff than at Darwin Glacier, potentially reconciling these elevation differences.
Third, the majority of the Ross Sea Drift on Minna Bluff is found significantly lower than the reported LGM maximum of 637 m asl reported by Denton and Marchant (2000). Only one small patch occurs higher than 500 m asl; this defines the reported 637 m asl LGM limit (see Plate 1 in Christ and Bierman, 2020). As there are no numeric age constraints at this 780 location, it is possible that this represents a short-lived (~1 kyr) transient maximum that could have also occurred at Diamond Hill, but would not be detectable by our in-situ 14 C measurement. The short-lived maximum at Lake Wellman reported by King et al. (2020)  Glacier around 10 kyr BP as the grounding-line retreats between Minna Bluff and the TAM front. The models generally agree poorly with glacial geologic data in the McMurdo Sound region as well, causing it to deglaciate thousands of years earlier than glacial geologic data suggest. This is not surprising, given that ice-sheet model sensitivity to climate forcing is strongly sensitive to a number of poorly known parameters (Pollard et al., 2016;Lowry et al., 2020), climate forcing itself is poorly 790 known, and the Bedmap2 topography beneath the modern day Ross Ice Shelf used by these models is coarsely resolved (Fretwell et al., 2013). Including new constraints on bed topography (Tinto et al., 2019;Morlighem et al., 2020) could lead to different model behavior. Additionally, the large ensemble experiments that are tuned to glacial geologic data (Pollard et al., 2016;Albrecht et al., 2020) attempt to fit data from around West Antarctica or the whole continent and thus might sacrifice predictive skill in the Ross Embayment to improve the overall ensemble score when tuning spatially homogeneous parameter 795 values. Tuning parameters using a large ensemble over the Ross Embayment alone could lead to a better fit between models and data.
We suggest that the relatively recent and slow deglaciation of the DHGS is likely due to the convergence of Byrd and Mulock Glaciers near the mouth of Darwin Glacier and/or lateral drag past Minna Bluff and Cape Crozier, which would both have led to dynamic ice thickening that opposed grounding-line retreat. Numerical investigations of grounding-line dynamics 800 have shown that convergent flow can counteract the acceleration of dynamic thinning as the grounding line retreats down a reverse bed slope (Gudmundsson, 2013); positive strain rates measured south of Minna Bluff show that convergent flow and compression are causing dynamic thickening of the ice shelf in this region today (Thomas et al., 1984). Whillans and Merry (2001) identified Cape Crozier and Minna Bluff as controlling obstacles to the flow of the modern Ross Ice Shelf, and we expect they would have had a similar effect during the last deglaciation. Together, the effects of convergent flow and lateral 805 drag could have created a protected embayment that resisted grounding-line retreat for longer than glaciers farther to the south, even given the steeply back-sloping bed topography of Discovery Deep. This lends support to the hypothesis that Byrd Glacier is of fundamental importance to the stability of ice in the Ross Sea sector (Hughes et al., 2017).

Conclusions
We have dated deposits of the DHGS in order to help constrain the timing and pattern of grounding-line retreat in the Ross 810 Embayment since the Last Glacial Maximum. The data suggest a later and slower deglaciation than that experienced by glaciers both farther south (e.g. Spector et al., 2017) and farther north (Jones et al., 2015). The scarcity of glacial deposits near the modern grounding line of Darwin Glacier and absence of deposits immediately adjacent to the grounding line make it difficult to directly interpret the thickness and timing of the LGM at the mouth of Darwin Glacier. We used a 1.5-D glacier flowband model to evaluate possible deglaciation scenarios for Darwin Glacier that are consistent with our new data and we used a 3-D 815 ice sheet model to evaluate the sensitivity of model reconstructions of regional deglaciation to poorly known model parameters.
Our key findings are: • LGM glacial deposits in ice-free valleys alongside Hatherton Glacier record up to 450 m of thickening relative to present at Lake Wellman, 350 m at Magnis Valley, and 300 m at Dubris Valley. The glacier margin extended several kilometers into each valley during its maximum before receding slowly and steadily through the Holocene. It did not reach its present thickness 820 until ≤2.8 kyr BP.
• Erratics perched stably on granitic bedrock at Diamond Hill, which is 10 km upstream of the modern grounding-line, record 135 m of thinning of Darwin Glacier between 5.1 kyr BP and 300 yr BP, with most of that thinning complete by 3 kyr BP.
• Maximum bedrock 14 C exposure ages indicate two possible, but conflicting, thinning histories at the mouth of Darwin Glacier: the LGM surface of Darwin Glacier was either >190m but <~500 m above the modern glacier and thinned steadily through 825 the mid-Holocene, or it was ~950 m thicker than present with ~600-800 m of rapid thinning in the mid Holocene.
• We used a glacier flowband model to test the two conflicting deglaciation scenarios suggested by the data at the mouth of Darwin Glacier. Flowband model results show that our data are most consistent with ~500 m slow and steady thinning at the mouth of Darwin Glacier between ~9 kyr BP and ~3 kyr BP, accompanied by a large decrease in catchment area through the Holocene. Rapid deglaciation at the glacier mouth -like that exhibited by other major glaciers in the Ross Embayment -is 830 not consistent with our geochronologic data. The fit of the model results to the data from Hatherton Glacier is more sensitive to rates of ice thickness change than to the absolute LGM ice thickness.
• An ensemble of 48 runs using a 3-D ice-sheet model indicates that using glacial geologic data of ice-elevation change to constrain even relative rates of grounding-line retreat is not straightforward. Counter to expectation, periods of more rapid grounding-line migration may correspond to relatively constant thinning at outlet glacier mouths. The presence of pinning 835 points such as Cape Crozier reduces the rate of modeled grounding-line retreat, but not necessarily the rate of thinning at Darwin Glacier. Thus, rates of thinning from glacial geologic records at Darwin and Hatherton glaciers should not be directly interpreted as reflecting rates of grounding-line retreat. This may also be true of other glaciers that we have not examined here.
• We suggest that the relatively slow thinning of Darwin and Hatherton Glaciers through the last deglaciation could be the result of convergent flow of Darwin Glacier with Byrd, Mulock, and Skelton glaciers, along with lateral drag due to the flow 840 past Minna Bluff.