A topographic hinge-zone divides coastal and inland ice dynamic regimes in East Antarctica

The impact of late Cenozoic climate on the East Antarctic Ice Sheet is uncertain. Poorly constrained patterns of relative ice thinning and thickening impair the reconstruction of past ice-sheet dynamics and global sea-level budgets. Here we quantify long-term ice cover of mountains protruding the ice-sheet surface in western Dronning Maud Land, using cosmogenic Chlorine-36, Aluminium-26, Beryllium-10, and Neon-21 from bedrock in an inverse modeling approach. We find that near-coastal sites experienced ice burial up to 75–97% of time since 1 Ma, while interior sites only experienced brief periods of ice burial, generally <20% of time since 1 Ma. Based on these results, we suggest that the escarpment in Dronning Maud Land acts as a hinge-zone, where ice-dynamic changes driven by grounding-line migration are attenuated inland from the coastal portions of the East Antarctic Ice Sheet, and where precipitation-controlled ice-thickness variations on the polar plateau taper off towards the coast. The Heimefrontfjella escarpment in Dronning Maud Land, East Antarctica divides inland precipitation-dominated ice dynamics from coastal dynamics dominated by grounding line migration, suggest inverse model simulations constrained by cosmogenic nuclide measurements.

T he East Antarctic Ice Sheet (EAIS) is assumed to have been stable since the Mid Miocene climate transition (14.8-13.8 Ma 1,2 ). However, owing to its large ice volume (22.4 million km 3 /52.2 m sea-level equivalent 3 ) even relatively small changes in EAIS extent and/or thickness have potentially large impacts on global sea level. Constraining the configuration of the EAIS during past climates is critical to improve our understanding of past sea-level changes and make predictions for the future. Inverse modeling of ice-core records indicates interior EAIS thinning of~80-140 m during glacial maxima in the late Pleistocene 4,5 , while the ice sheet concurrently extended close to the shelf edge along most of its perimeter 6,7 . But beyond these first-order trends, the EAIS response to climatic changes remains poorly understood. Past configurations of the EAIS are difficult to constrain empirically, in part because the present-day ice sheet covers 99.8% of the continent 8 .
Cosmogenic nuclide exposure dating in Antarctica is often applied to constrain the ice-sheet demise since the Last Glacial Maximum (LGM) when rocks on coastal islands and on nunataks became exposed. However, the low glacial erosion rates associated with the hyper-arid polar climate present a challenge for dating the last ice retreat 2,[9][10][11] , owing to inheritance of long-lived and stable cosmogenic nuclides (such as 26 Al, 10 Be, and 21 Ne) from previous exposure periods [12][13][14] . On the other hand, contributions from previous periods of exposure also provide an opportunity to explore the long-term ice-sheet history, because long-lived cosmogenic nuclides retain a memory of past exposure, burial, and erosion [15][16][17][18] . Indeed, complex exposure and burial histories can be determined utilizing multiple nuclides with different halflives 18 . In this study, we set out to resolve the late Cenozoic iceburial history of nunataks in western Dronning Maud Land (DML), East Antarctica, by analyzing cosmogenic 36 Cl, 26 Al, 10 Be, and 21 Ne in exposed bedrock surfaces of the Heimefrontfjella escarpment and along the Jutulstraumen Ice Stream (Fig. 1). Glacial and subaerial erosion rates in Antarctica are low but need to be constrained to accurately assess the long-term ice history 18 . We, therefore, apply a Markov Chain Monte Carlo (MCMC) inversion procedure to this dataset to simultaneously quantify erosion and ice-burial histories. Our data inversions document a change in the duration of ice burial with distance to the coast, where near-coastal sites along the Jutulstraumen Ice Stream experienced ice burial up to 79-97% of time since 1 Ma, while inland sites on the Heimefrontfjella escarpment experienced <20% ice burial since 1 Ma. This difference indicates that the escarpment acts as a hinge-zone for ice-dynamic changes in DML, separating inland precipitation-dominated ice dynamics from coastal dynamics dominated by grounding-line migration.
Geomorphological and glacial setting of DML study sites. We investigate two regions in western DML. One is a coast-parallel transect along the inland escarpment in Heimefrontfjella, while the other is a coast-to-inland transect within the Jutulstraumen Ice Stream drainage basin (Fig. 1a).
Heimefrontfjella is a~120 km long segment of a NE-SW trending mountain range separating ice on the high polar plateau from the lower coastal portion of the EAIS, approximately 200 km inland of the present-day grounding line in western DML (Fig. 1). Heimefrontfjella has a passive margin escarpment morphology, with near-vertical visible bedrock cliffs up to 500 m high, facing N-NW and extending subglacially to the north. It is sub-divided into four massifs (from NE to SW: Milorgfjella, XU-fjella, Sivorgfjella, and Tottanfjella) that are separated by minor outlet glaciers which coalesce into the Veststraumen Ice Stream and the Riiser-Larsen Ice Shelf (Fig. 1b).
Ice on the polar plateau abuts the summits from the southeast. Heimefrontfjella forms a topographic barrier to ice flow and an orographic barrier for precipitation. Samples analysed in this study were collected from gneiss, quartzite, and granite bedrock outcrops in Milorgfjella, Sivorgfjella, and Tottanfjella (n = 13) during the austral summer of 2016-17, and span elevations from 1561 to 2215 m a.s.l., 20-530 m above the ice (Supplementary  Table S1, Fig. S1; See Methods section for details on Heimefrontfjella samples).
Cosmogenic nuclide chronometry and inversion approach. To investigate ice-burial duration and erosion depths during the late Cenozoic, we employ multiple cosmogenic nuclides to constrain model selection, including one stable nuclide ( 21 Ne) and three radionuclides with different half-lives ( 10 Be: T 1/2 = 1.39 Myr, 26 Al: T 1/2 = 0.705 Myr, 36 Cl: T 1/2 = 0.301 Myr; Supplementary Tables S2-S4, Fig. S2). We analyzed 2-4 nuclides in each of our samples depending on mineral availability (Supplementary  Table S1). We utilize a MCMC modeling approach similar to one used for Greenland 21,22 but adapted to the longer timescale of glaciation in East Antarctica and including 36 Cl and 21 Ne. Exposure and erosion histories are defined by a set of parameters, which are constrained through the MCMC inversions ( Fig. 2; Table 1). Each forward model represents a complex exposure history with periods of ice burial, when cosmogenic-nuclide production is inhibited, and periods of cosmogenic-nuclide production during subaerial exposure (Fig. 2). Two parameters define the exposure history: a global marine benthic δ 18 O ref. 23 threshold value, δ 18 O th , and the timing of last deglaciation, T dg , following the LGM. To characterize the sample trajectory from depth within bedrock to the surface, and to allow for (glacial and subaerial) erosion to vary through time, an erosion history is parameterized using four time-depth tie points and an initial steady-state erosion rate E 0 (Fig. 2). Inherent trade-offs between the time-depth parameters describing the exhumation trajectories entail that the overall exhumation path may be well constrained although single parameters are not. Due to these trade-offs, the distributions of single parameters are not necessarily informative. We, therefore, report the time when a sample reaches within 1 m of the present-day surface ('time of 1 m erosion') and the sample depth at 1 Ma ('erosion since 1 Ma'), to illustrate the exhumation paths to the present-day surface consistent with nuclide inventories for each sample site (Fig. 3a-d). The 'time of 1 m erosion' parameter is an indicator of the (minimum) memory of the cosmogenic nuclide inventories in the sample-within 1 m of the surface, the nuclide inventory becomes increasingly sensitive to surface production. This memory effect entails that for samples with low cosmogenic-nuclide abundances and where nuclide ratios converge on surface-production ratios, only the most recent history can be constrained, which becomes directly evident by the wider range of accepted parameter values. Finally, for each sample, we report the duration of ice burial since 1 Ma (%) derived from the modeled δ 18 O th distribution. Unless specified otherwise, we cite the full range of accepted models that fit the maximum number of measured nuclides for each sample within 1σ; for the few samples where the inversions do not find solutions within 1σ, we refer to the 2σ range, and for the BOR01 sample that has no solutions within 2σ of all measured nuclides we refer to the full range of all accepted models (Supplementary Table S5). Results and discussion Exposure and erosion histories derived from the inversion procedure. Rather than offering unique solutions, the MCMC inversion approach returns a range of exposure and erosion scenarios consistent with the measured nuclide concentrations. Figure 3 displays examples of inversion results for four different samples, each with at least two combinations of nuclide constraints ([ 10 Be, 26 Al] vs. [ 10 Be, 26 Al plus 21 Ne and/or 36 Cl]). These four samples were selected to illustrate the full span of exposure and erosion histories recorded by our samples. Particularly, the four samples differ in their cosmogenic-nuclide memory, i.e., how long they have been within 1 m of the presentday surface. In the following, they are presented in order from longest to shortest memory. The summit at Ristinghortane, the highest-elevation sample site in Heimefrontfjella (RH01; 2215 m a.s.l.), has likely experienced less than 1 m of erosion since the Late Miocene to Early Pliocene and has been mostly subaerially exposed (<14% ice burial; [ 10 Be-26 Al-21 Ne]-constrained inversion) since 1 Ma ( Fig. 3a; Supplementary Table S5). MID02, the highest sample from Midbresrabben, a summit in-between the Penck Trough and Jutulstraumen (1660 m a.s.l.), was likely exhumed to within 1 m of the present-day surface in the Pliocene and has been fully shielded by the EAIS for 17-75% of the last 1 Myr ([ 10 Be- 26 Al- 36 Cl]-constrained inversion; Fig. 3b; Supplementary Table S5). The Pleistocene erosion history of MID02 remains unconstrained, ranging from exhumation in the early Pleistocene to episodic plucking of up to~0.8 m within the last glacial cycle (Fig. 3b). The [ 10 Be, 26 Al]-constrained inversion of sample VK02 from Vardeklettane (1561 m a.s.l.), that crops out at the foot of the Heimefrontfjella escarpment, is comparable to the [ 10 Be, 26 Al]-constrained inversion of MID02. However, when adding 36 Cl to the inversion constraints, VK02 agrees with a recent erosion event and limited ice burial (<20% since 1 Ma; Fig. 3c; Supplementary Table S5). Like RH01, including an additional cosmogenic nuclide in the inversion leads to fewer accepted models that provide a good fit to 26 Al. A sample collected from Cottontoppen Jr (CJ01; 1609 m a.s.l.),~4 km from VK02 was likewise rapidly eroded since 1 Ma and has experienced limited ice burial during exhumation (<12% since 1 Ma; Fig. 3d; Supplementary Table S5). Cross plots between 'Erosion since 1 Ma' and the 'Ice burial since 1 Ma (%)' (Fig. 3e-h) highlight the trade-off between ice burial and recent rapid erosion for [ 10 Be, 26 Al]-constrained inversions of samples MID02 (Fig. 3f) and VK02 (Fig. 3g), where high degrees of ice burial correspond with low erosion rates, and vice versa. Importantly, including 36 Cl in the inversions facilitates the discrimination between these scenarios.  Table S5). Inversions with and without 21 Ne lead to similar ice-burial histories for this group of samples (Fig. 4b, gray vs. green tones), but the addition of 21 Ne pushes the accepted models towards more recent erosion ( Fig. 4c). At lower elevations the inversion results indicate more recent erosion with samples arriving within 1 m of the presentday surface in the Late Miocene to Middle Pleistocene (Fig. 4c). Of the five lowest samples (1561-1646 m a.s.l.), one sample stands out immediately with >1.7 m erosion and a short ice burial duration (<12%) since 1 Ma (CJ01; Fig. 3d; Supplementary  Table S5). For the four other samples (MAB01, 05, VK01, 02) inversions of [ 10 Be, 26 Al] indicate longer ice burial (interquartile range: 37-59% for all accepted models), although less ice burial is also possible as indicated by the tail of solutions towards 0% (Fig. 4b). Interestingly, including 36 Cl in the inversion for two of these samples (VK01, 02) markedly changes the modeled parameter distributions towards shorter ice-burial duration (<10-20%; Figs. 3c, 4b) and more erosion (0.5-0.9 m) since 1 Ma (Figs. 3c, 4d; Supplementary Table S5). It is noteworthy that the  Timing of last deglaciation 0 to 25 ka (varies between samples, see Table S1) Zdg Depth at time of last deglaciation 0-0.1 m dZ 2 , dZ 3, dZ 4 Depth change from previous tie point of exhumation path 0.05-10 m dT 2 , dT 3, dT 4 Time change since previous tie point of exhumation path 0-4 Myr E 0 Constant erosion rate from model start until T 4 0.01-100 m Myr −1 inversions including 36 Cl are more compatible with results from nearby samples CJ01, BRA01, and BRA02 regarding ice burial and a preference for high recent erosion rates like sample CJ01 (Fig. 4c, d). However, the addition of 36 Cl also entails that the inversions find fewer models that fall within 1-2σ of all measured nuclide concentrations as indicated by lighter colors in Fig. 4b-d.  Table S5).
Sensitivity of inversions to cosmogenic nuclide production. Inverse modeling of multi-nuclide data from bedrock samples in the Heimefrontfjella and Jutulstraumen regions provides novel constraints on long-term ice cover and erosion histories in the region. However, not all inversions yield models that fit all nuclide concentrations within 2σ (e.g., BOR01, VK01 ( 10 Be- 26 Al- 36 Cl-21 Ne); Fig. 4). This may (i) result from errors that are unaccounted for in the nuclide measurements or the applied production rates or (ii) reveal that the parameterization of the forward model does not fully reflect the exposure and erosion histories of the samples. Although uncertainties related to calibration and scaling of the production parameters are not incorporated into the inversion algorithm, sensitivity tests show that our results are largely insensitive to the choice of scaling scheme (Supplementary Discussion SD1; Supplementary Fig. S3) as well as potential nucleogenic 21 Ne in our samples (Supplementary Discussion SD2; Supplementary Fig. S4). Uncertainties related to muogenic production are also of minor importance for our samples that generally have long residence times in the near-surface zone (Fig. 4), which is dominated by spallation reactions. In the following, we, therefore, focus on the sensitivity of our inversion results to 36 Cl production. Inversions performed in this study demonstrate that measuring 36 Cl in combination with 26 Al-10 Be may enable the differentiation between exposure-erosion scenarios with-and-without episodic plucking (Fig. 3b, c, f, g). This is important since long periods of exposure followed by recent bedrock plucking can lead to similar 26 Al-10 Be inventories as prolonged burial under ice 18 . Introducing 36 Cl to the inversion procedure may constitute a promising tool for constraining long-term exposure in polar regions. However, 36 Cl has a broader range of production pathways for different elements compared to 10 Be and 26 Al, which may complicate the applicability of this method depending on bedrock composition. The ten samples with 36 Cl measurements presented in this study have widely varying compositions (Supplementary Tables S1, S3-S4) and are consequently dominated by different production pathways ( Supplementary Fig. S5). Generally, spallation and muon production on K and Ca are better constrained than low-energy neutron capture on native Cl 24 and uncertainties in the low-energy neutron capture production parameters, therefore, dominate the uncertainty of our results. For 36 Cl samples with a large production component stemming from low-energy neutron capture, an additional complication arises from the presence of snow cover. This is because the hydrogen in water both in and above the ground surface acts as a strong moderator of neutron fluxes [25][26][27] . Although we expect a generally thin or absent snow cover at our sites due to their exposed locations, low precipitation, and strong katabatic winds, the long-term evolution of snow thickness at our sites is not well constrained.
To test the sensitivity of our results to the production parameters for neutron capture, we ran three additional experiments. In the first experiment, we remove the production component that accounts for leaking of low-energy neutrons from the shallow subsurface to the atmosphere. This approach is designed to mimic the most extreme case of hydrogen moderation through snow shielding [25][26][27] and effectively removes the sub-surface maxima in the 36 Cl production profiles ( Supplementary Fig. S5) so all production pathways decline exponentially with depth. In the second and third experiments, we decrease or increase all parameters related to thermal and epithermal 36 Cl production by 20% ( Supplementary Fig. S5), corresponding to the estimated uncertainty in the neutron capture pathway 24 . These sensitivity tests demonstrate that changing the low-energy neutron 36 Cl production parameterization shifts the parameter distributions of accepted models somewhat but does not fundamentally change our results ( Supplementary Fig. S6). The largest effect is seen in experiment 1 ("maximum Hydrogen moderation"), which produces a wider range of potential ice burial since 1 Ma for all samples, and more erosion since 1 Ma for the Jutulstraumen samples. However, for the Jutulstraumen samples, this experiment also yields worse fits to measured nuclide concentrations than our standard model, with few or no accepted models fitting within 2σ of all nuclides for six out of seven samples. In contrast, decreasing (or increasing) all parameters related to thermal and epithermal 36 Cl production by 20% leads to a <10 percentage point decrease (or increase) in ice burial since 1 Ma and minor shifts in derived erosion histories, but allows for solutions that fit all measured nuclides within 1-2σ in most instances (Supplementary Fig. S6).
Sensitivity to ice-burial parameterizations. The results presented in this study are based on the assumption that ice burial at our sites fluctuates in sync with the global δ 18 O curve. This assumption is simplistic since the δ 18 O values measured in benthic foraminifera reflect a composite of global ice volume and ocean temperature and salinity 1,23,28 . Furthermore, the long-term increase in δ 18 O values towards the present combined with a single δ 18 O threshold value effectively enforces scenarios that go from more to less exposure through time (Fig. 2), which may not be realistic for all our sites. In support of the use of the δ 18 O-curve as a proxy for ice thickness, ice-sheet models simulating EAIS surface elevations over the past 5 Myr yield a longterm increase in ice thickness in coastal sectors of DML [29][30][31] ( Supplementary Fig. S7c-e) as would be expected if the coastal ice responds to grounding line migration in sync with global sea level 1 (Supplementary Fig. S7b). On the other hand, given current EAIS margin surface gradients in DML, our highest-elevation sites are unlikely to be covered by the EAIS even when the ice sheet advances to the continental shelf break during glacial maxima. It has been suggested that nunataks at high elevations in DML were last ice covered in the Pliocene when ice-sheet surface gradients at the margins were potentially steeper due to a stronger hydrological cycle resulting from a warmer Southern Ocean combined with diminished sea ice and ice-shelf extent [32][33][34] . Additionally, mountain peaks will have experienced gradual uplift due to flexural isostasy following subglacial erosion focused beneath warm-based, fast-flowing ice in adjacent troughs 35 . Both effects would result in samples at high elevations experiencing a gradually decreasing ice cover through timethe opposite of what is imposed by our δ 18 O-record approach.
To test the sensitivity of the exposure parameterization, we ran alternative inversions using mean daily insolation at 65°N on summer solstice 36 as a proxy for ice thickness ( Supplementary  Fig. S7f). Although this implementation does not yield realistic ice-cover scenarios, it effectively removes the long-term increasing trend in ice cover, and we thereby avoid imposing a Cenozoic cooling trend altogether. The results are comparable to the δ 18 O implementation regarding ice-burial histories for all samples, and erosion histories for high-elevation samples with long nearsurface residence times ( Supplementary Figs. S8-S9). The main difference lies in the long-term exhumation history for lowerelevation samples: without the forced change from exposure to ice burial at the advent of the Pleistocene, these samples no longer need to be situated at depths >1 m in the Pliocene and Miocene to avoid too much exposure ( Supplementary Fig. S8). This difference is reflected in the 'Time of 1 m erosion' distributions that are much broader for the insolation implementation ( Supplementary Fig. S9).

Regional variability of ice-sheet burial, escarpment hinge-zone.
Our results point to large spatial variations in ice-sheet burial and erosion histories of presently exposed nunataks in western DML. Indeed, our inversion results indicate that the ice-sheet surface in the Jutulstraumen drainage basin was above the present elevation for most of the Late Pleistocene, covering samples located 35-295 m above the present-day ice-sheet surface 13 (Supplementary Table S1) and 30-135 km upstream of the grounding line for up to 75-97% of the last 1 Myr (Fig. 4; Supplementary  Table S5). Only the two highest elevation samples near Jutulstraumen (>2100 m a.s.l.; 183-416 m above ice) remained ice free for the majority of the last 1 Myr (BOR01, 03; <8-18% ice burial). In contrast, most samples from Heimefrontfjella are consistent with dominantly ice-free conditions since 1 Ma within 1σ. Samples at lower elevation (1500-1700 m a.s.l.; 20-225 m above ice) require <10-71% ice burial during the last 1 Myr, whereas samples at higher elevation (1700-2200 m a.s.l.; 20-530 m above ice) require <1-18% ice burial.
Patterns of long-term exposure reconstructed in this study are consistent with the pattern of exposure since the LGM as inferred from the youngest [ 10 Be, 26 Al, 36 Cl] apparent exposure ages from Jutulstraumen 13,37 . These patterns demonstrate that the EAIS surface in the Jutulstraumen drainage basin was >200 m higher than present within the last glacial cycle near the ice-stream trunk 13 and 120 m higher~100 km to the east 37 , and that deglaciation occurred in the Middle to Late Holocene. In comparison, LGM ice-surface change at Heimefrontfjella is estimated to be <50 m above present 38,39 .
Differences in the long-term ice-cover histories between the Heimefrontfjella escarpment and the near-coastal Jutulstraumen Ice Stream are likely a result of their different settings. Heimefrontfjella forms part of the escarpment in a region with only few, relatively narrow (6-9 km) and shallow (<850 m) outlet glaciers 40 (Fig. 1). The escarpment forms an efficient barrier to ice flow and the discharge of ice from the polar plateau, advected through the escarpment, likely decreased during global glacial maxima when the interior ice sheet thinned due to moisture starvation from a more distal and cooler Southern Ocean 4,5 (Fig. 5). We suggest that the escarpment acts as a hinge-zone for ice-dynamic changes. Below the hinge-zone, the coastal portion of the EAIS experiences dynamical ice-thickness fluctuations, whereas the plateau ice above the hinge-zone is mainly controlled by changes in precipitation (Fig. 5). In contrast, the samples in the Jutulstraumen drainage basin are situated along the largest ice stream in DML (up to 50 km wide and grounded more than 1.6 km below present sea level), with current ice-flow velocities up to 750 m/a at the grounding line 41 . Such dynamic components of the margin respond mainly to the interplay between sea-level forcing, precipitation cycles, and glacial isostatic adjustments (Fig. 5). We do not discuss glacial isostatic adjustments in detail here as it is covered elsewhere 37,42,43 . However, in general, the outward migration of the grounding line in response to global sea-level fall leads to dynamic thickening of the coastal portions of the ice sheet and a propagation of thickening upstream (Fig. 5). These changes attenuate towards the escarpment, including the Heimefrontfjella sites that are located 110-155 km from the nearest present-day grounding line. The precipitation barrier imposed by the escarpment currently leads to a steep decline in snow accumulation rates from 0.2-0.7 m a −1 ice equivalents in the coastal sector to <0.1 m a −1 on the polar plateau 44,45 (Fig. 5). Samples collected only 20-50 m above the present-day ice surface on the edge of the polar plateau at Ristinghortane, have cosmogenic-nuclide inventories consistent with <14-18% ice burial since 1 Ma, supporting that thickening across this upstream sector of the ice sheet was of limited duration and/or frequency. That mountain ranges can act as barriers to ice-dynamic changes is in line with studies from the Transantarctic Mountains that document inland stability of the EAIS since >200 ka despite large fluctuations in marginal ice thickness in the Ross Sea Sector 46 .
Erosion and escarpment stability. In addition to an improved ice-burial history, the inversion results also unravel the erosional histories at our sample sites. Our results imply that erosion rates are generally low, with samples from both areas often residing within 1 m of the present-day bedrock surface since the Late Miocene, Pliocene, or Early Pleistocene (Fig. 4c, g). At the broad scale, we also resolve a trend of decreasing amounts of erosion with increasing elevation (Fig. 4c, d, g, h). Assuming erosion is tied to glacial processes, the decreasing erosion with elevation can be explained by shorter duration of ice burial for higher-elevation sites or lower erosivity of the ice sheet at higher elevation where the ice is thinner and colder [47][48][49] . Several studies document low erosion rates beyond major troughs in DML: (i) in bedrock During global interglacials (a) sea-surface temperatures (SST) and sea level increase while sea-ice extent shrinks, and the grounding line retreats. As a result, increased evaporation leads to increased precipitation and inland ice thickness. In contrast, icecore records (illustrated in a) indicate interior ice thinning during global glacial periods (b), when SST and sea level are lower, the grounding line advances, and sea ice is more extensive resulting in lower evaporation and precipitation. Cosmogenic nuclide inversions in this paper indicate that the escarpment acts as a hinge-zone where ice-dynamical changes resulting from grounding line migration taper off inland, and, inversely, inland precipitation-driven changes taper off towards the coast.
samples above the present-day ice-sheet surface using cosmogenic nuclides 10,33 and thermochronometry 50 , and (ii) below the present-day ice-sheet surface by ground-based 51,52 and airborne 20 radar surveys, or inversion of ice-surface imagery 53 revealing subglacially-preserved alpine and fluvial landscapes.
By inverse modeling 10 Be, 26 Al, and 36 Cl inventories in samples from the base of the Heimefrontfjella escarpment, we show that these sites experience episodic plucking during the last 1 Myr. Inversions of the three lowest-elevation samples in Heimefrontfjella indicate that >1.7 m (CJ01) and 0.5-0.9 m (VK01, VK02) of erosion occurred since 1 Ma (Fig. 3c, d). At Cottontoppen Jr, the erosion was of sufficient magnitude to remove previously accumulated nuclides, while for the Vardeklettane-samples, our inversions indicate that the erosion event(s) followed a long period of (mostly subaerial) exposure in the near subsurface (Fig. 3c, d; 4b-d). If the episodic erosion event(s) is (are) related to subglacial plucking, we can speculate that one or several glaciations since the mid-Pleistocene transition became sufficiently thick to cause patchy subglacial sliding and glacial erosion at these sites. However, we cannot rule out that erosion was due to subaerial slope processes. In the field, we often observed that bedrock is highly weathered and/or covered by a layer of regolith 54 (in-situ weathered rock). To ensure our samples recorded continuous exposure since they were last deglaciated, they were generally collected from spatially restricted outcrops of the locally most-competent rock with indications of glacial sculpting or striae. This is true for the samples from Vardeklettane and Cottontoppen Jr. (Supplementary Fig. S1). Based on field evidence, it is likely that these bedrock outcrops were formed by episodic (subglacial or subaerial) removal of overlying regolith or weakened rock. Likewise, recent plucking is also possible for samples MAB01 and MAB05 ( Supplementary  Fig. S10) that have 10 Be and 26 Al inventories very similar to the Vardeklettane-samples, but for these samples, we lack material for 36 Cl measurements to test this hypothesis.
Our results indicate stability of the Heimefrontfjella escarpment around samples MH01 and 02 since 3.4-6.6 Ma (Figs. 4c and 5; Supplementary Table S5). These samples were collected at the edge of the escarpment, where exposed rock rises ca. 500 m above the coastal sector of the EAIS in DML on the western side, while being nearly flush with the polar plateau ice-sheet margin on the east. It is hard to imagine how these sites could be overridden by the EAIS in the current topographic configuration where increased accumulation upstream would more likely lead to increased ice flux through the escarpment outlet glaciers than ice flow across the escarpment at its highest parts. We propose that the ice-free conditions of the (generally <50 m wide) rim of the escarpment results from wind-scour and, perhaps, ablation of snow by radiation heating from surrounding rocks, locally creating a steady negative mass balance. Under these circumstances, exposure of these samples is controlled by escarpment retreat uncovering rock that may have been glacially sculpted and striated when the ice sheet was last warm-based in the Early-to-Mid Cenozoic 2,48,55 . Recent estimates indicate~20 km of escarpment retreat in DML since ice-sheet inception, of which the bulk occurred prior to the Miocene 50 . Our inversions show that nuclide inventories in samples on the escarpment crest are consistent with near-surface residence times since the Pliocene or Late Miocene, and <1-8% of ice burial since 1 Ma.

Conclusions
Inverse modeling of cosmogenic multi-nuclide datasets can be used to explore and quantify long-term ice-sheet history and erosion in glaciated regions characterized by low erosion rates, such as Antarctica. Long periods of exposure followed by episodic plucking of shallow bedrock layers, or stripping of regolith, complicates the interpretation of ice-burial histories when measuring only 26 Al and 10 Be. Our study demonstrates that both iceburial durations ( 36 Cl) and erosion-depth estimates ( 21 Ne) may be better resolved when additional nuclides can be measured in the same sample.
For the Heimefrontfjella escarpment setting, we find that nuclide inventories are consistent with a dominance of bedrock exposure (ice burial <20% of last 1 Myr for most samples) combined with generally shallow erosion depths (<0.02-0.9 m since 1 Ma). In contrast, a combination of a dominance of ice burial (up to 75-97%) and limited erosion (<0.1-1.0 m) since 1 Ma, provide the best fit to the data from Jutulstraumen, except at elevations >1700 m a.s.l. where the data are consistent with low erosion rates (<0.04-0.2 m) and little-to-no ice burial (<8-18%) since 1 Ma.
Based on these results we suggest that the escarpment in DML acts as a hinge-zone, where ice-dynamic changes resulting from migration of the grounding line is attenuated inland from the coastal portions of the EAIS, and where precipitation-controlled ice-thickness variations of the polar plateau ice taper off towards the coast. Finally, we suggest that long durations of bedrock exposure at the escarpment crest may reflect the local mass balance combined with slow escarpment retreat indicating escarpment stability since 3.4-6.6 Ma.

Methods
Heimefrontfjella samples. Haukelandnuten (Milorgfjella): Milorgfjella constitutes the easternmost extent of Heimefrontfjella (Fig. 1b), with the present-day ice sheet breaching the escarpment to the north-east. We collected two samples from striated (335-340°) quartz veins protruding from the gneissic bedrock at the ice-free rim on top of the escarpment at Haukelandnuten (MH01, 02). The samples were collected adjacent (<10 m) to a~500 m near-vertical north-facing bedrock wall.
Månesigden and Bowrakammen (Tottanfjella): The Månesigden-Bowrakammen ridge extends in a NW-SE direction, approximately perpendicular to the overall trend of the escarpment, and descends towards the coastal portion of the ice sheet ( Fig. 1). We collected four samples from Månesigden (MAB01, 02, 03, 05), and two from Bowrakammen (BRA01, 02) forming an elevation transect from 1580-1730 m a.s.l. The ice-sheet surface elevation differs by 160-180 m between the two sides of the ridge, yielding a wide range of elevations above the ice sheet (45-310 m; Supplementary Table S1). The local orthogneiss bedrock is weathered, often regolith covered, and displays up to 20 cm deep tafoni near MAB02. Samples were collected from the least weathered bedrock outcrops or from quartz veins (Supplementary  Table S1). Some sample sites preserve faint striations trending 220-255°(BRA01; Fig. 2b) and 318-328°(MAB02, 05), and bedrock around MAB05 was glacially sculpted. No erratic clasts were found above 1650 m a.s.l., below this elevation glacially transported clasts of apparently local lithology were observed in a few places.
Vardeklettane and 'Cottontoppen Junior' (Tottanfjella): We sampled two nunataks protruding the coastal portion of the ice sheet on a partially submerged ridge near the western-most end of Heimefrontfjella, also striking perpendicular to the escarpment. A sample was taken from a polished granite dyke protruding the weathered augengneiss bedrock on 'Cottontoppen Junior' (our informal name, CJ01; Fig. 2c) near the escarpment. Two bedrock samples were collected from ã 6 × 2 m quarzitic bedrock outcrop on Vardeklettane (VK01, 02; Fig. 2d) four kilometers from the escarpment. Although Vardeklettane is mostly regolith covered and weathered local bedrock comes apart in large slabs, the sampled sites appear glacially sculpted and largely intact. Faint striae were observed to trend 315°a nd a mafic erratic cobble was found in the vicinity, probably sourced from a nearby col (<200 m up-ice relative to established striae direction).
Study design. Cosmogenic nuclides accumulate in Earth's surface layers by cosmic-particle induced nuclear reactions and are removed through erosion and radioactive decay. The sensitivity to erosion is high due to the short attenuation length of the dominant spallation reaction (150-230 g cm −2 ) 56 . In intermittently glaciated regions, the interplay between periods of ice burial or subaerial exposure, and subglacial/subaerial erosion, regulates cosmogenic nuclide inventories. The multiple unknowns associated with complex exposure and erosion histories requires measuring multiple nuclides with different half-lives [15][16][17][18]57 . Traditionally, the interpretation of two-nuclide diagrams, casting the shorter-to-longer lived nuclide ratio (e.g., 26 Al/ 10 Be) against the longer-lived nuclide concentration (e.g., 10 Be) normalized to sea-level high-latitude production ( Supplementary Fig. S2), has been employed to constrain the duration of landscape burial beneath thick ice 58 . However, the underlying assumption of steady-state erosion is unlikely to be valid in polar regions where ice covers tend to be only sporadically erosive 58,59 . This is also true for Antarctica, where ice is inferred to be mostly cold-based and nonerosive 12,60,61 . To encompass the possibility of episodic plucking events, we employ a MCMC inversion approach 18,21,22 . This MCMC inversion approach enables us to explore the range of exposure and erosion scenarios that comply with the measured nuclide concentrations through iterative comparison between simulated and measured nuclide abundances.
Analytical methods. Laboratory procedures and 10 Be, 26 Al, and 36 Cl measurements on Jutulstraumen samples are described in ref. 13 . For the Heimefrontfjella samples, quartz separation and initial 10 Be- 26 Al chemistry was done at the Scottish Universities Environmental Research Center (SUERC). Rock samples were crushed, and quartz was isolated from the 250-500 μm sieved fractions by aqua regia, froth flotation, and magnetic separation before etching with dilute HF/HNO 3 in an ultrasonic tank for at least 3 days 62,63 . We assessed the purity of the quartz by induced coupled plasma optical emission spectrometry (ICP-OES) analysis. Samples were considered clean if aluminium concentrations were less than 250 μg/g quartz. In preparation for Al chemistry, processing blanks and samples with low Al were spiked with a Fischer Scientific 27 Al carrier to reach a total of~1000-1500 μg Al. For Be chemistry, all samples and blanks were spiked with (~200-220) μg 9 Be from an in-house carrier (CIAF PH9). A processing blank accompanied each batch of 8-14 unknown samples. After digestion in concentrated HF, an aliquot was removed from each sample to allow for a final total-Al determination on ICP-OES. Samples were then dried down and converted to chloride form preceding Be and Al isolation by ion chromatography. At this stage, samples were transferred to PRIME Lab at Purdue University, where samples were oxidized and mixed with Nbpowder and pressed into stainless steel cathodes for accelerator mass spectrometry (AMS) analysis.
All 10 Be/ 9 Be ratios were normalized to 07KNSTD 64 while 26 Al/ 27 Al ratios were normalized using the KNSTD standard 65 . All samples were corrected for contamination during the laboratory procedure using the processing blank associated with each batch of samples. Corrections are <5% for all samples and <0.5% for most samples. Uncertainties, including Be carrier concentration uncertainty, were propagated through to the final result. During the 26 Al/ 27 Al AMS run of our samples, one of the processing blank samples had an unusually high number of 26 Al counts, which was suggested to be caused by Mg interference in the AMS. As a result, all our 26 Al/ 27 Al ratios were corrected for Mg interference. For samples in this study, this correction was minor (<4%). Post-dissolution aliquots were dried down and re-dissolved in 5% HNO 3 , prior to ICP-OES determination of the total Al content (native + carrier). We prepared and measured one aliquot of the Co-Qtz-N laboratory inter-comparison material together with each batch of unknown samples 66 . The obtained results overlap within 1σ for both 10 Be (average = 2.45 ± 0.04 × 10 6 at g −1 (n = 4)) and 26 Al (average = 15.6 ± 0.4 × 10 6 at g −1 (n = 4)), and both values overlap the consensus values ( 10 Be: 2.53 ± 0.09 × 10 6 at g −1 and 26 Al: 15.6 ± 1.6 × 10 6 at g −1 ) from a round-robin exercise 66 (Supplementary Table S2).
Cosmogenic 21 Ne analysis was performed at the SUERC Noble Gas Isotope Laboratory on 50-100 mg splits of clean quartz following the procedure outlined in ref. 67 . Each sample was packed into Pt foil tubes and Ne was released by heating to 1300°C using a diode laser. Neon isotope composition was determined using a modified ThermoFisher ARGUS VI magnetic sector mass spectrometer. Repeat measurements of CREU-1 quartz standard material 68 run with our samples yielded a reproducibility of ±4%, which we apply as uncertainty for samples with only single 21 Ne measurements or where repeat measurements yielded standard errors <4%. For samples where repeat measurements overlap within 1σ we use the weighted mean, otherwise, we use the higher mass aliquot.
Feldspar-rich mineral separates for 36 Cl-analysis were obtained from the floatfraction following froth-flotation, magnetic components were removed, and the residue leached three times in 5% nitric acid in an ultrasonic bath. Sample dissolutions, Cl separations, and AMS measurements were performed at PRIME Lab using the same procedure as outlined in ref. 13 . The computed nucleogenic 36 Cl was subtracted from the total measured 36 Cl and uncertainties propagated through to final results (Supplementary Table S1). Sample compositions are shown in Supplementary Tables S3-S4.
Production rate calculations. We compute site-specific spallation production rates for 10 Be and 26 Al in quartz scaled for latitude and elevation according to the 'St' scaling scheme 69 , with sea-level high-latitude reference production rates of 4.01 atoms (g quartz) −1 yr −1 for 10 Be and 27.93 atoms (g quartz) −1 yr −1 for 26 Al 70 . We derive 21 Ne spallation production from the 21 Ne/ 10 Be ratio of 4.08 from ref. 9 . We use a spallation attenuation length of 155 g cm −2 (ref. 56 ) consistent with the atmospheric depths (~720-900 g/cm 2 ; ERA-40) and low cut-off rigidity at high latitudes (75 to 71°S). We compute parameters for 10 Be, 26 Al and 21 Ne production by muon interactions by fitting two exponentials to the 'model 1 A' nuclide-specific, site-atmospheric-pressure dependent muon production profile 71 . We calculate the parameters for 36 Cl production in feldspar-rich mineral separates and bulk-rock (sample VID03) using MATLAB codes from CRONUScalc v. 2.1 56 , again using the 'St' production rate scaling method 69 and fitting exponential functions to production depth profiles to increase computational speed. The 36 Cl production pathways include spallation reactions (Ca, K, Fe, and Ti), muon production (Ca and K), and neutron capture production ( 35 Cl). The computed 36 Cl production rates for all 36 Cl samples (n = 10) are shown in Supplementary Fig. S5. We note that the 'LSDn' scaling scheme generally yields better results, especially at high latitudes 72 , but there is currently no open-source framework available for calculating the 'LSDn' production rates for all four nuclides in a computationally consistent way. Furthermore, we note that any erosional flexural isostatic uplift due to trough incision would lead to an underestimation of the true near-surface residence times of our samples.
Inverse modeling of multi-nuclide inventories. Exposure and erosion histories are defined by a set of parameters, which are constrained through the MCMC inversion (Fig. 2). Two parameters define the exposure history: a global marine benthic δ 18 O ref. 23 threshold value, δ 18 O th , and the timing of last deglaciation, T dg , following the LGM. These two parameters are allowed to vary freely between models, but the latter in a relatively narrow time interval as it is often well determined by 36 Cl/ 26 Al/ 10 Be ages on nearby erratic boulders 13 (Table S1). For sample sites that may have remained ice-free during the LGM, the T dg parameter boundaries are set wide enough to allow full exposure during the LGM. To characterize the sample trajectory from depth within bedrock to the surface, and to allow for (glacial and subaerial) erosion to vary through time, an erosion history is parameterized using four timedepth tie-points (T dg , Z dg ; T 2 , Z 2 ; T 3 , Z 3 ; T 4 , Z 4 ) and the initial steady-state erosion rate E 0 (Fig. 2). Here, Z dg represents the depth of the sample at the time of the last deglaciation (T dg ), while T 2-4, and Z 2-4 represent three time and depth tie-points younger than 15 Ma and within 20 m depth. The multi-nuclide inventories can be computed from any combination of these ten parameters and the resulting cosmogenic-nuclide inventories compared to measured values. We use a Metropolis-Hastings type MCMC method to approximate the distribution of models that can explain the measured nuclide concentrations 21,73,74 . To ensure that the results are not biased by the starting position, we use 10 random walkers initiated at different positions in the model space. To ensure a wide initial search of the model space, each Markov chain starts with a burn-in phase of 10.000 iterations with a target acceptance ratio (proportion of accepted to rejected models) of 0.15. This is followed by (up to) 500.000 iterations with a target acceptance ratio of 0.4, until 100.000 models are accepted. The accepted models are used to construct posterior parameter distributions. The results of the individual walkers were checked for convergence, and to ensure that the model results are not unduly influenced by the prior choice of parameter space boundaries. Finally, we note that this inversion approach encompasses two simple scenarios that are often used for interpreting cosmogenic nuclides in bedrock surfaces: (i) steady-state subaerial erosion, and (ii) episodic erosion of a sufficient depth to remove previously formed cosmogenic nuclides followed by continuous subaerial exposure without erosion. However, the inversion approach described here also allows for exploring a much broader set of potential exposure and erosion scenarios.

Data availability
All cosmogenic nuclide data presented in this manuscript is available at https://doi.org/ 10.5281/zenodo.7422400.

Code availability
MATLAB code used for generating the MCMC inversions can be downloaded from https://github.com/jaluan/cosmo-inversion