Introduction

The East Antarctic Ice Sheet (EAIS) is the largest continental ice mass on Earth. With a freshwater storage equivalent to ~52.2 m of global sea level1, even minor changes in ice-sheet volume can strongly influence global sea level and climate dynamics2. It is therefore critical to understand its sensitivity to various forcing agents. Recent advances in satellite gravimetry and ice-sheet modeling have yielded refined estimates of EAIS mass balance3 and allowed for a re-examination of its response to global climatic changes4,5. A problem with satellite gravimetry is that the observed changes inevitably include the effect of glacial isostatic adjustment (GIA) of the Solid Earth driven by past ice sheet changes. This is the largest source of error in estimating net mass balance3,6. However, because space-geodetic observations of ice-sheet change only extend over the last few decades, EAIS behavior before the 1970s remains largely unknown. Modeling studies of ice-sheet sensitivity to atmospheric and oceanic warming and circulation, and sea-level rise also require well-constrained records of past ice-sheet changes for model validation and refinement7,8. Together, these factors highlight the need for long term (millennial-scale) glacio-geological records of ice surface elevation changes obtained from field-based studies. This is especially pertinent for understanding the large-scale processes impacting the EAIS9,10, which is crucial for direct comparison with present-day observations11.

Ice thickness changes over the EAIS since the Last Glacial Maximum (LGM) are of particular importance because large-amplitude climatic changes during this period allow for an exploration of the forcing mechanisms responsible for ice-sheet retreat over large scales. Inflow of warmer Circumpolar Deep Water (CDW) has been proposed as a trigger for retreat of the ice sheet across the continental shelf after the LGM12,13, and ice-shelf melt driven by CDW intrusion is considered the primary cause of the current negative mass balance of the Antarctic ice sheets14,15,16. Eventual ice-shelf collapse is expected to induce inland ice-sheet thinning because of accelerated ice outflow13,17,18. Another potential trigger for perimeter-wide ice sheet retreat is sea-level rise19,20. Ice loss due to the demise of most Northern Hemisphere ice sheets during the last deglaciation is estimated to have resulted in a sea-level rise of 80–130 m around Antarctica, and the inclusion of sea level forcing in ice-sheet models causes Antarctic ice sheets to retreat further, quicker, and earlier19. However, the relative importance of see-level rise and influx of warmer CDW for driving EAIS retreat remains uncertain.

In this study, we present a new, well-constrained ice-sheet thinning history for the period since the LGM for Dronning Maud Land (DML), one of the major sectors of the EAIS (Fig. 1a, b). We use the ice-thinning chronology obtained by surface exposure dating of glacially transported rocks from Gjelsvikfjella (Gjelsvik Mountains) to refine a reconstruction of the ice-load history and calculate regional sea-level changes in DML considering GIA effects. Based on these field and modeling results, we infer a regional, relative sea-level highstand at 9–8 thousand years ago (ka), with concomitant influx of warmer CDW triggering a breakup of ice shelves via the marine ice sheet instability mechanism, and subsequent rapid thinning of coastal ice sheet sectors in DML.

Fig. 1: Study area.
figure 1

a Antarctica and the Southern Ocean. Antarctic shading indicates ice sheet elevation, and Southern Ocean contours show the vertically integrated transport stream function78. Red (blue) colors indicate a cyclonic (anticyclonic) circulation, and the Weddell Gyre is highlighted by a white arrow. Contour intervals are 10 Sv (1 Sv = 1.0 × 106 m3 s−1). EAIS East Antarctic Ice Sheet, SR Shackleton Range, WS Weddell Sea, WG Weddell Gyre. The black box marks the location shown in b. b Overview of Dronning Maud Land and the location of the Gjelsvikfjella study area (black box, c) within the Jutulstraumen catchment (blue dashed line). The yellow circles and green star are locations where seawater temperatures are extracted from the atmosphere–ocean general circulation model40 and regional sea-level change is calculated, respectively. The blue dashed line indicates the area of the Jutulstraumen catchment used for ice sheet model experiments. The red arrow shows the approximate flow direction of warmer Circumpolar Deep Water (CDW). c Overview of Gjelsvikfjella, including Rabben, Grjotfjellet, and Terningskarvet. Arrows indicate the direction of ice flow. Small white boxes correspond to Fig. 2a–c. Blue eyes indicate approximate locations for the photographs of Fig. 3. The base map is from the Quantarctica GIS data package79 for QGIS. Ice velocities80 and ice surface contours81 (every 100 m) are overlain on the Landsat Image of Antarctica (LIMA)82.

Full size image

Results

Gjelsvikfjella is a group of mountains (nunataks) in the western part of Fimbulheimen, piercing the EAIS 100–150 km inland from the grounding line (Fig. 1b, c). The ice-sheet surface elevation on the coastal side of Gjelsvikfjella is ~1000 meters above sea level (m a.s.l.), rising as high as ~2500 m a.s.l. on the inland side, forming a mountain escarpment in the Jutulstraumen catchment. Depositional and erosional glacial landforms on the mountains Grjotfjellet, Rabben, and Terningskarvet (Fig. 1c) are used here to investigate past elevation changes of the EAIS.

Grjotfjellet (1690 m a.s.l.) is located west of the Jutulsessen nunatak and north of the Fjellimellom Glacier, a major artery through Gjelsvikfjella. Grjotfjellet has a series of well-defined moraine ridges on its western slope, here numbered as moraines I (highest elevation) to IV (lowest elevation). Moraine I occurs up to an altitude of ~1500 m a.s.l. half-way up the slope, and moraine IV spans altitudes of ~1370–1400 m a.s.l. and extends along and immediately east (30–50 m) of the contemporary glacial margin (Figs. 2a and 3a). The mountain slope above moraine I lacks distinct glacial landforms but many erratic cobbles and boulders (up to several meters high) are distributed across the weathered bedrock. The density of erratics generally decreases with increasing elevation.

Fig. 2: Detailed topographic maps of Grjotfjellet.
figure 2

Locations of samples and their 10Be exposure ages from Grjotfjellet (a), Rabben (b), and Terningskarvet (c) are shown along with corresponding elevation reference points on the contemporary ice sheet. Well-defined moraine ridges at Grjotfjellet are outlined by yellow lines with Roman numerals indicating moraine numbers (I–IV). 10Be exposure ages of erratic samples collected from the moraine ridges are labeled by their corresponding moraine numbers. Blue dashed lines indicate contemporary (bold) and inferred past (thin) locations of the ice sheet margin. Topographic contours are in meters above sea level (m a.s.l.), based on the 2 m-resolution Reference Elevation Model of Antarctica (REMA)53 and Earth Gravitational Model 200854. The background image is a morphometric protection index red relief image map83 produced from the same DEM with glacial-ice areas shaded in half-transparent blue. See Supplementary Data 1 and 2 and Supplementary Note 1 for more details on the samples, their exposure ages, and age interpretation.

Full size image
Fig. 3: Photographs of Grjotfjellet, Terningskarvet, and Rabben.
figure 3

Approximate sample site locations are shown by yellow circles along with 10Be exposure ages. a Overview of Grjotfjellet, with Terningskarvet in the background. 10Be exposure ages of erratic samples collected from three moraine ridges are labeled by their corresponding moraine numbers. b Sample locations on Grjotfjellet immediately adjacent to the contemporary ice sheet surface. Yellow dashed lines in a, b delineate well-defined moraine ridges at Grjotfjellet, with Roman numerals indicating moraine numbers that correspond to those in Fig. 2a. c Overview of Rabben. The surface elevation of the ice sheet on the western side of Rabben (where the photograph was taken) is ~150–200 m lower than that on the opposite side. d The 2690 m a.s.l. summit of Terningskarvet.

Full size image

Rabben is an isolated nunatak at 1440 m a.s.l. located in the northernmost sector of Gjelsvikfjella (Fig. 1c). Slithallet Slope, another glacial artery connecting the ice sheet plateau with its coastal sector, exhibits northward ice flow on the eastern side of the nunatak, whereas ice on its western side appears quasi stagnant, owing to its location on the lee side of Jutulsessen. Contemporary ice-surface elevations on the eastern side are 150–200 m higher than the western side (Figs. 2b and 3c). No moraines are recognized on Rabben, but there is an abundance of erratic cobbles and boulders sitting on the north-south trending ridge as well as on the western slope.

Terningskarvet, a series of flat-topped summits at similar elevations, comprises the highest part of Gjelsvikfjella, and sits more than 200 m above the present ice stream surface to the southwest (Fig. 2c). One of its flat-topped summits, at 2690 m a.s.l., consists of highly weathered bedrock with scattered relatively unweathered erratic cobbles and boulders of local gneiss type lithology (Fig. 3d). The current ice flow is blocked by Terningskarvet, except for a small local ice artery that just barely overpasses the mountain immediately south of the flat-topped summit on which we collected samples (Figs. 1c and 2c).

We measured in situ 10Be and 26Al concentrations on 39 quartz-bearing erratics collected from Gjelsvikfjella; 26 from Grjotfjellet, 9 from Rabben, and 4 from Terningskarvet (see Methods and Supplementary Data 1 and 2). At Grjotfjellet, the erratic samples span a vertical transect of >300 m from the lowest-elevation moraine (IV) up to near the peak, and their 10Be exposure ages fall in the range from 1.13 ± 0.13 ka (sample G25, 1371 m a.s.l.) to 738 ± 54 ka (sample G01, 1674 m a.s.l.), with generally older ages from samples at higher positions (Figs. 2a, 3a, b and Supplementary Data 2). Two ages of 14.7 ± 1.0 ka and 20.4 ± 1.3 ka at 110–120 m above the ice sheet reference surface (samples G10 and G11) are both from the highest-elevation moraine (I). Given that the erratics collected at even higher elevations all yield ages older than 50 ka (Fig. 4a), it is likely that moraine I represents the LGM glacial highstand. A vertical profile of the ages shows a clear trend of younger ages towards the contemporary ice sheet surface for samples younger than LGM (Fig. 4a). The post-LGM ages are distinctly separated from older outliers, which are probably derived from recycled erratics with multiple exposure histories as reported elsewhere in Antarctica21. The ages indicate that the ice sheet surface was ~100 m higher than at present at ~8.4 ka, but rapidly thinned down to ~15 m above present at ~5.7 ka (Fig. 4b). The ages from moraine IV, located immediately adjacent to the contemporary ice sheet surface (G24–26), indicate that the ice sheet margin has been near its present position for at least the last 1000–2000 years. Although we did not sample between the ~5.7 ka position and Moraine IV, the two endmember interpretations of this distribution are that the ice sheet has either remained relatively stable since ~5.7 ka, or that it has been thinner than today prior to thickening to its present configuration.

Fig. 4: Surface exposure ages plotted against relative heights above the contemporary ice sheet surface.
figure 4

a Grjotfjellet for all samples. b Grjotfjellet for samples showing ice sheet lowering since LGM. Roman numerals indicate moraine numbers that correspond to those in Figs. 2a and 3a. c All samples from Terningskarvet. d Rabben samples showing ice sheet lowering since LGM, with the inset showing all samples. Markers surrounded by dashed circles are used for the linear regression analysis of the rapid deglaciation phase (Fig. 6). All uncertainties are reported as 1σ external errors.

Full size image

Aside from two pre-LGM outliers, 10Be ages from Rabben ridge, covering the range 6.47 ± 0.48–8.21 ± 0.56 ka, demonstrate that the nunatak was fully covered by ice immediately before ~8.2 ka. An apparent age of 11.1 ± 0.8 ka at ~70 m above the reference surface (sample R09) is from an erratic on the western slope rather than on the ridge (Figs. 2b and 3c), and we regard this sample as a recycled clast. This dataset shows a strikingly similar trend to that from Grjotfjellet: the ice sheet surface was ~70 m higher than the ice sheet reference surface at ~8.2 ka, but it rapidly thinned down to ~20 m above the reference surface at ~6.5 ka (Figs. 2b–4d). The relative height of ~70 m thus needs to be taken as a minimum value for the ice thinning since the LGM at Rabben.

Unlike the Holocene ages prevalent at Grjotfjellet and Rabben, all erratics from Terningskarvet yield much older 10Be ages between 774 ± 56 ka and 2460 ± 290 ka (Figs. 2c–4c). These exposure ages indicate that Terningskarvet has not been fully covered by the EAIS for the last several hundreds of thousands of years or has only been covered for relatively brief periods by non-erosive ice. Because we do not observe clear evidence of ice cover during the LGM, it is likely, but not certain, that the LGM glacial highstand remained below the summits in Terningskarvet.

Although there is an excellent match between the Holocene deglaciation history in Grjotfjellet and Rabben, we note that local ice sheet history should not be directly interpreted as reflecting regional changes. The ice sheet response at a specific site could be exaggerated, reduced, or even reversed by effects of local bedrock topography above and below the ice sheet22. To examine these effects at our sites, we performed high-resolution (with grid-size down to ~500 m) ice sheet model experiments for the entire Jutulstraumen catchment across western and central DML regions (Fig. 1b) using 1 km-grid bedrock topography (see Methods and Supplementary Note 2). A total of 75 experiments spanning a range of plausible basal sliding, rheological, and climate conditions, resulted in a wide variety of steady-state configurations of the ice sheet at the LGM. Of these, 53 experiments reproduced ice sheet configurations consistent with previous reconstructions of thicker margins and thinner interior geometries9,23,24 and with our inferences from Grjotfjellet (Fig. 4b), where model results yield ice surface elevations 80–205 m higher than at present. In those preferred models, the ice sheet is 120–250 m higher than at present at Rabben, which is also consistent with the observed ice-surface drop of >70 m since the LGM (Fig. 4d).

One important implication of our modeling is a regionally-differentiated response of the ice sheet to glacial conditions. Although a higher-than-present ice sheet surface prevails over the modeling domain downstream of the escarpment, regional differences are pronounced, varying from ~100 m higher than at present at Grjotfjellet and Rabben to as high as ~800 m above present along ice streams such as nearby Jutulstraumen (Fig. 5b). In contrast, the modeled elevation differences at Terningskarvet are ~30–50 m lower than at present, which supports the interpretation from the 10Be ages that this peak remained ice-free during the LGM. The model constrains the escarpment as a hinge zone for ice-sheet surface change below which the coastal sector of the EAIS has thinned since the LGM and above which the polar plateau has thickened to present elevations (Fig. 5b). Conversely, basal or rheological properties cannot be properly constrained from our results. Within the climate and parameter space explored with our high-resolution model setup, the largest spread in simulated ice sheet geometries was due to climate differences. However, we note that the combination of ice sheet model parameters in our best-performing ensemble member (Fig. 5; see “Methods”) is consistent with estimates from previous work (see Section 2.1.5 in the Supplementary Note 2)25,26. In summary, our model experiments show that Grjotfjellet and Rabben, by virtue of their location downstream of the escarpment, are representative of the regional coastal sector ice-sheet changes (i.e., thinning since the LGM), but not of the magnitude of change because this varies markedly over the region.

Fig. 5: Surface elevation difference between the best-fit LGM experiment and the contemporary ice surface.
figure 5

Positive values denote thicker ice during the LGM, while ice-free areas are in white (see Supplementary Note 2 for more details). a Dashed white line marks the modern grounding-line. Rectangle outlines the location shown in b. Study sites (Fig. 1c) are shown as white triangles. b Terningskarvet (TS) Grjotfjellet (GF), and Rabben (RB).

Full size image

Discussion

New 10Be exposure ages from Grjotfjellet and Rabben in Gjelsvikfjella provide clear evidence for ~100 m of rapid EAIS thinning between 8.4 and 5.7 ka while new high-resolution ice sheet modeling results show that these sites reflect the overall EAIS behavior in the coastal sector downstream of the mountain escarpment in the Jutulstraumen catchment.

Evidence for a thicker EAIS during the LGM has been reported elsewhere in and around DML. Along the northern side of the Shackleton Range (~30°W), immediately west of DML, the EAIS drains into the Ronne-Filchner Ice Shelf and in-situ cosmogenic 14C dates indicate that this portion of the ice sheet was 350–650 m thicker during the LGM27. In western DML, west of the Jutulstraumen ice stream, the EAIS is thought to be at least 100–200 m thicker than at present within the last 25 ka (Fig. 1b)23. In central and eastern DML, less ice thickening during the LGM is inferred: <50 m in the Wohlthat Massif28,29 and <100 m in the Sør Rondane Mountains30 (Fig. 1b). Our new data reinforces the hypothesis that the EAIS had a gradient of increased thickening towards the west (the Weddell Sea) along the coast of DML during the LGM23.

A study of the post-LGM deglaciation in the easternmost area of DML depicts a strikingly similar story to ours31. A total of 32 10Be dates from the southern Sôya Coast, Lützow–Holm Bay (Fig. 1b and Supplementary Fig. 1.1b) indicate that the EAIS during the LGM extended at least 10 km seaward from the contemporary ice margin and thinned by ~400 m from ~9 to ~5 ka, concurrent with ice-margin retreat. The higher magnitude of ice thinning at the Sôya Coast can be related to its coast-proximal setting adjacent to the Shirase Glacier (Fig. 5). A large-scale ice-sheet model simulation of the entire EAIS shows that the inland-coast ice sheet thinning gradient observed in our study area is common among all catchments in DML, including the Sôya Coast (see Supplementary Fig. 3.1). From our new chronological constraints and modeling, combined with previous observations, we argue that a rapid ice-sheet thinning occurred throughout the coastal regions of DML between the Early (9–8 ka) and Middle Holocene (6–5 ka).

To further examine the thinning phase, we apply a Monte Carlo linear regression analysis using the iceTEA toolkit32 to all the data from Grjotfjellet and Rabben that comprise a clear Holocene thinning phase. The linear regression analysis reveals that the composite best-fit thinning rate for the sites in Gjelsvikfjella is 0.05 m yr−1 (0.03–0.10 m yr−1 in a 95% confidence interval) with an onset age of 8.3 ka (7.6–9.1 ka in a 95% confidence interval; Fig. 6a and Supplementary Fig. 1.4a). We apply the same analysis to the selected age-elevation dataset for the Sôya Coast31 (Supplementary Data 3) and obtain a higher best-fit thinning rate of 0.15 m yr−1 (0.09–0.51 m yr−1 in a 95% confidence interval) with an onset age of 8.6 ka (7.4–10.2 ka in a 95% confidence interval; Fig. 6b and Supplementary Fig. 1.4b). These thinning rates are smaller than those reported from the Amundsen Sea sector of West Antarctica (>1 m yr−1), but comparable to those from the interior parts of the Weddell Sea and Ross Sea sectors33. The reported timings of the onset of post-LGM thinning vary vastly across Antarctica, from ~12 to ~3 ka33. Our estimated onset of thinning at Gjelsvikfjella and the Sôya Coast fall within this range. Although the timings of thinning onset overlap within errors, the slightly earlier median timing at the Sôya Coast may reflect gradual dynamic thinning from the coast to the ice sheet interior. In any case, the consistent timing of thinning onset at two sites that are ~1400 km apart could demonstrate synchronized ice-sheet thinning across the coastal sector of DML, which then implies a common external forcing mechanism.

Fig. 6: Comparison of ice-sheet thinning in Droning Maud Land using surface exposure ages and ice load models since the LGM.
figure 6

a, b Results of linear regression analysis to determine a best-fit profile of ice surface lowering in Gjelsvikfjella (this study), western-central DML and the Sôya Coast, Lützow-Holm Bay, eastern DML31 (Fig. 1b) by using the iceTEA toolkit32. The exposure ages used for the analysis in Gjelsvikfjella are circled in Fig. 4. All exposure ages from the Sôya Coast31 except for bedrock and recent moraine samples are used for the analysis (see Supplementary Data 3). Exposure ages are shown with 1σ external errors and are plotted against relative heights above the contemporary ice sheet surface for Gjelsvikfjella and elevation for the Sôya Coast. Gray lines represent all models that fit the data, the solid black line is the best-fit to the data, and dotted black lines represent the 95% confidence interval. Rates of lowering were generated from 5000 iterations through randomly sampled points using 2σ internal uncertainties. c, d Ice load models for each area based on ICE6G42 and its refined version based on the field-based ice thinning histories for Gjelsvikfjella (c) and the Sôya Coast (d).

Full size image

The timing of EAIS thinning across DML postdates major climatic events since the LGM. Indeed, the Meltwater Pulse 1 A (MWP-1A) sea-level rise of 15.7–20.2 m in 500 years at ∼14.650 ka34,35 and global and regional warming of the atmosphere and the ocean36,37 clearly preceded the initiation of the rapid ice-sheet thinning (Fig. 7). Atmospheric temperatures estimated from ice cores at Dome Fuji38 and EDML39 peak at 12–11 ka, and then exhibit a gentle decline towards the present (Fig. 7d). A coupled atmosphere-ocean general circulation model shows that the ocean temperature at a water depth of 450 m, a typical depth of CDW37, has a similar trend40 (Fig. 7c and Supplementary Note 5). Although the increased ocean temperature is expected to have carried more heat to the coast of DML via intrusion of warmer CDW as part of the Weddell Gyre16,41, our documented thinning of the EAIS in DML was initiated 2000–4000 years later, indicating the presence of another controlling process.

Fig. 7: Post-LGM ice-sheet thinning in Droning Maud Land compared with paleoclimatic proxies from 20 ka to present.
figure 7

a The 95% confidence intervals of ice surface lowering in Gjelsvikfjella (red, left axis) and the Sôya Coast (blue, right axis) from Fig. 6, and original ICE6G ice-load models for Gjelsvikfjella and the Sôya Coast. All uncertainties are reported as 1σ external errors. b Relative sea levels at 0˚E, 71.5˚S calculated by an original ICE6G ice-load model42 and its refined version, compared with the eustatic sea level42. The thin lines are relative sea levels based on different earth models (applicable range for upper and lower mantle viscosity), and the bold line corresponds to the best model76 (see “Methods”). c Modeled Ocean temperatures at 450 m water depth40, extracted for several locations offshore of DML (see Fig. 1b). A depth of 450 m is selected because this is the typical depth where CDW is located37. d δ18O record from Dome Fuji38 (blue, right axis) and EDML39 (purple, left axis) ice cores using the Antarctic Ice Core Chronology 201284. The gray and red vertical bars represent the deglaciation phase used for refining the GIA ice load model, and the timing of Meltwater Pulse 1A (MWP-1A)34, respectively.

Full size image

We propose regional sea-level rise, determined by global sea level and regional GIA, as a trigger of rapid EAIS thinning in DML. The solid Earth’s delayed viscoelastic response (uplift) due to regional ice-load changes leads to a sea-level highstand, which we hypothesize culminated in rapid ice-sheet thinning. We calculate the regional sea-level change at the current grounding line of the Jutulstraumen catchment (0˚E, 71.5˚S, Fig. 1b) using the ice-load history model ICE6G42 with varying mantle rheology, and find that sea level exhibits a broad peak at 8–9 ka (Fig. 7b). However, deglaciation histories for Grjotfjellet, Rabben, and the Sôya Coast derived from the original ICE6G model do not match with those derived from the field-based data (Fig. 6c, d). Although the magnitudes of thinning are largely consistent, the initiation of thinning at ~11 ka in the ICE6G model clearly predates the onset of thinning constrained by field evidence by ~2000 years. We therefore recalculate the regional sea-level change based on a refined ICE6G ice-load model, where the magnitude of thinning remains unchanged, but the timing of ice-sheet thinning is delayed to fit the field evidence from the ice-sheet catchments that encompass Gjelsvikfjella and the Sôya Coast (Antarctic drainage systems 5–7; Supplementary Fig. 4.1a; see “Methods” and Supplementary Note 3). The refined regional sea-level history exhibits a clear highstand of at most 30 m above current sea level at 9–7 ka (Fig. 7b). The delayed timing of deglaciation leads to a delayed viscoelastic response (uplift) of the Solid Earth, thereby causing a more prominent overshoot of the regional sea-level rise relative to the present. This overshoot originates from the competing effects of the global sea-level rise of the melting Northern Hemisphere ice sheets and regional crustal uplift driven by the thinning Antarctic ice. A more pronounced regional sea-level highstand that coincides with the initiation of rapid ice-sheet thinning indicates that this process played a key role in the pacing of this important and swift response of the EAIS.

We propose that the incipient deglaciation of the DML margin of the EAIS was initiated by ~11 ka when the inflow of warmer CDW reached the ice front (Fig. 8a, b). However, in the DML region, the deglaciation was apparently limited to coastal areas and did not propagate into the ice-sheet interior because the sea-level highstand was not prominent enough to allow inflow of CDW beneath the ice shelves. While the ocean temperature remained relatively high (Fig. 7c), sea level continued to rise locally until ~9–8 ka, reaching probably ~15 m and possibly >20 m higher than at 12–11 ka (Fig. 7b). The sustained sea-level rise likely promoted grounding-line retreat by increased ice-shelf buoyancy and facilitated warm CDW inflow beneath the ice shelves, thereby enhancing basal melting and calving via the marine ice sheet instability mechanism, culminating in rapid EAIS thinning that propagated >100 km inland (Fig. 8c). This interpretation is supported by a marine sediment record recovered from the Lazarev Sea, offshore Schirmacheroase (Schirmacher Oasis), central DML (Fig. 1b), which shows the establishment of fully marine conditions during the Early Holocene43 after grounding line retreat.

Fig. 8: Schematic image of a thinning of the EAIS across DML.
figure 8

a The EAIS advanced due to sea-level fall during the LGM19. b Warming of the CDW initiated retreat of the ice shelf, as observed in the Weddell Sea sector20. c Regional sea-level highstand triggered enhanced basal melting and calving of ice shelves, facilitated by the inflow of warmer CDW, which peaked at 9–8 ka and continued to ~5 ka, resulting in grounding line retreat and ice-sheet thinning. d Sea-level fall due to the GIA of the Antarctic continent stabilized ice shelves and/or may have caused a readvance of the EAIS, supported by pinning on ice rises and ice rumples along the DML margin44.

Full size image

The subsequent sea-level fall from Mid-to-Late Holocene (Fig. 7b), due to the delayed Solid Earth’s response to surface unloading, likely stabilized ice shelves and reduced calving. Ice rises and ice rumples scattered along the coast of DML likely play an important role in this stabilization process44 (Fig. 8d). The above effects of regional GIA on the local ice sheet stabilization process are thought to be responsible for the termination of the deglaciation phase and may have even caused a readvance and thickening in the coastal sector of DML. The moraines adjacent to the contemporary ice surface found at several locations in Gjelsvikfjella, including moraine IV at Grjotfjellet (Fig. 3b), probably resulted from this stabilization and/or readvance/thickening. The notion of a readvance is supported by marine sediments recovered offshore of the Riiser-Larsen Ice Shelf, western DML (~16.5°W; west of Fig. 1b), which revealed a local readvance of the ice shelf in the Late Holocene45. Additional evidence for the EAIS readvance in DML was found in the southern part of Sôya Coast, where lacustrine sediment embedded with Late Holocene molluscan shell fossils directly underlies the contemporary ice sheet46. Although a late Holocene readvance of the EAIS needs to be further verified, several lines of geological and geomorphological evidence support stabilization or a slight readvance of the EAIS in DML during the sea-level fall phase following the Middle Holocene45,46.

Holocene readvance of the Antarctic Ice Sheet has also been suggested for the Weddell Sea47, Ross Sea47,48,49, Amundsen Sea50, and Denman Glacier51 sectors, although direct dating49 is limited due to the subglacial position of the relevant materials. As for DML, Mid-to-Late Holocene relative sea-level fall is likewise invoked to explain these apparent readvances47,50, although for the Ross Sea sector, regional climate change has been advocated as an alternative driving mechanism52.

Aided by high-resolution ice-sheet modeling, our new ice thinning chronologies from Gjelsvikfjella along with previously published data from the Sôya Coast31, demonstrate that abrupt regional-scale deglaciation of the EAIS occurred across coastal DML from the Early to Middle Holocene. This event was probably due to a regional sea-level highstand that facilitated the inflow of warmer CDW. Although higher regional sea level is thought to have been a direct final trigger, we emphasize that the preceding warming of CDW was an important prerequisite, and the mixed effect of these two factors pushed the ice sheet towards a regional tipping point. Our hypothesized mechanism needs to be tested by future studies, including field-based reconstructions of ice-sheet history at other sites, analysis of marine sediments off the coast of DML, and more detailed coupled atmosphere, ocean, and ice-sheet modeling. The location and geological setting of DML provides an excellent study area for advancing the understanding of the EAIS response to global climate change, and for refining the future trajectory of the Antarctic cryosphere.

Methods

Field campaign and rock sampling for surface exposure dating

A detailed geomorphological and glacial geological survey in Gjelsvikfjella, central DML, was carried out during December 2015 - February 2016 as a part of the 57th Japanese Antarctic Research Expedition. Based from Troll Station, a Norwegian Antarctic research station located in the western part of Gjelsvikfjella (Fig. 1c), we used Ski-doos and a helicopter to access the study areas.

Thirty-nine samples for surface exposure dating with cosmogenic 10Be and 26Al were collected from erratics at Grjotfjellet, Rabben, and Terningskarvet (Figs. 2 and 3 and Supplementary Fig. 1.2). Sample lithologies consisted mostly of gneisses, some granites, and one pegmatite (Supplementary Data 1). To identify sites most suitable for exposure dating, we carefully considered the glacial geologic and geomorphologic settings, and erratic lithology. Size and stability were considered when selecting erratics, to minimize the likelihood of post-depositional overturning. We principally chose boulder-sized and flat-shaped erratics stably resting on bedrock or compact till. We also collected minimally weathered erratics to avoid recycled materials with near-surface residence periods spanning multiple glaciations. We calculated topographic shielding factors for each location with the CRONUS-Earth online calculators (https://hess.ess.washington.edu/) using horizon measurements in twelve directions with a handheld level.

Sample coordinates (latitude and longitude) were measured with a handheld GNSS device (GPSMAP 62s, GARMIN), and the elevations were taken from the 2-m-grid digital elevation model (DEM) produced from the Reference Elevation Model of Antarctica (REMA)53 and Earth Gravitational Model 200854 using the sample coordinates from the GNSS. We also quantified the relative height of each sample site from the contemporary ice sheet surface using the same DEM. We first manually delineated the contemporary ice sheet margin (bold blue dashed lines in Fig. 2), excluding areas around wind scoops and wind-blown snow ridges. We then chose the point on the ice margin that is closest to each sample site and took it as the corresponding reference point. In terms of the switch from westward to southwestward projections between G11 and G09 at Grjotfjellet, we also considered the glacial geomorphology. The southern part of moraine I, where G10 and G11 were collected, extends nearly north–south, suggesting that those sample sites should be referenced to the west rather than to the southwest. There are no moraine features on the slope above G11 that can inform us of former ice flow directions, and so reference is made to ice surface elevations to the southwest simply because it is substantially closer.

We have no quantitative information about snow cover at Gjelsvikfjella, but we assume that a very low mean surface snow mass balance of 111 mm yr−1 for the Asuka drainage basin55 (Fig. 1) is a reasonable estimate for our study area. Most areas of Gjelsvikfjella including our sample sites are snow-free in summer probably because the low albedo of the nunataks promotes the complete removal of small amounts of accumulated snow through sublimation. Therefore, time-averaged thicknesses of overlying snow/ice are inferred to be at least one order of magnitude less than the attenuation length of fast neutrons (Λ ≈ 150 g cm−2 = 150 cm in water/ice), and thus we assume that the snow shielding is negligible at our sample sites.

Sample details are listed in Supplementary Data 1 and photographs of typical sampling sites are shown in Supplementary Fig. 1.2.

Laboratory and analytical techniques

Sample preparation for 10Be and 26Al was performed at the National Institute of Polar Research (NIPR) and Purdue Rare Isotope Measurement Laboratory (PRIME Lab), Purdue University, following established methods56,57. The upper 3–5 cm of the collected erratic samples were crushed and sieved to 250–500 μm size to obtain clean quartz grains for surface exposure dating. Aliquots of the dissolved samples were analyzed for 27Al by inductively coupled plasma mass spectrometry (ICP-MS) at NIPR and inductively coupled plasma-optical emission spectrometry (ICP-OES) at PRIME Lab. The weights of the samples and carriers, and 10Be/9Be and 26Al/27Al ratios, all measured at PRIME Lab, are listed in Supplementary Data 2.

AMS measurements were normalized to the 07KNSTD Be standard material58 with nominal 10Be/9Be = 2.851 × 10−12 and to the KNSTD Al standard material59 with nominal 26Al/27Al = 1.818 × 10−12. Procedural blanks and samples were both spiked with 0.1–0.5 mg 9Be carrier. For Al analysis, blanks were spiked with 1.0–1.5 mg 27Al carrier (Al carrier, 1000 ppm), and the samples were spiked with ~1.7 mg 27Al carrier with the value depending on the native Al content of the sample. Blanks ranged from 1.7 × 10−15 to 35.4 × 10−15 [10Be/9Be] and 0.2 × 10−15 to 8.8 × 10−15 [26Al/27Al]. Concentrations reported in Supplementary Data 2 involved subtraction of background 10Be and 26Al atoms measured in procedural blanks. Uncertainties include propagated AMS sample/lab-blank uncertainties, and a 3% uncertainty for the stable 27Al measurement.

Exposure age calculations

Surface exposure ages were calculated using the CRONUS-Earth V3 (http://hess.ess.washington.edu/math/v3/v3_age_in.html) online calculator60. Exposure ages were calculated using a 10Be half-life of 1.387 Ma61,62, an 26Al half-life of 0.705 Ma59, CRONUS-Earth global production rates63, and the LSDn nuclide-specific scaling scheme64. Density of rock samples was assumed to be 2.7 g cm−3. In this study, we used 10Be exposure ages to construct deglaciation chronology because the production rate of 10Be is better constrained than that of 26Al65. Exposure ages from 26Al were only used to check their correspondence to 10Be-based exposure ages (within error ranges), examine the possibility of inherited nuclides from recycling of erratic clasts through multiple glacial cycles, and validate the experimental accuracy of nuclide quantification (Supplementary Figs. 1.3 and 1.4). All ages are apparent exposure ages, computed under the assumption of continuous exposure at Earth’s surface and zero erosion of the rock surface. Apparent exposure ages will represent the timing of last deglaciation only if nuclide concentrations, which potentially accumulated in the erratics during previous exposure periods, were below detection levels at the time of deposition and if erratics escaped self-shielding by post-depositional rolling or surface erosion (self-shielding by the mass of eroded rock). We made no attempt to account for minor variations in nuclide production rates caused by elevation changes due to GIA30. Uncertainties in apparent exposure ages in this study include scaling-, production rate-, and analytical uncertainties (based on laboratory/AMS measurements).

High-resolution ice sheet modeling

We used the ice-flow model Úa66 to produce steady-state LGM ice-sheet geometries for the Jutulstraumen catchment, extending its ice front to the continental shelf67 (blue dashed lines in Fig. 1b). To encompass the uncertainties in LGM climate and model parameters, we performed a total of 75 model experiments, forced by a range of CMIP5 and CMIP6 climate models, under a wide range of plausible basal sliding and ice rheology distributions (see Supplementary Note 2 for more details).

To estimate the change in ice surface elevation over the sampled sites, we computed the difference between the ice sheet surface elevation in the LGM simulations and the contemporary ice sheet surface elevation. This difference is computed as the difference between the LGM elevation at the averaged sample coordinates for each site, and the contemporary BedMachine-Antarctica1 ice sheet surface elevation at the averaged reference-point coordinates of each site. This yielded a total of 75 points for each site, which we used to compute a histogram of elevation differences (blue bars in Supplementary Fig. 2.9). We repeated the same process for a subset of experiments (n = 53; brown bars in Supplementary Fig. 2.9), where (i) at least 90% of all cosmogenic nuclide sample sites were covered by ice during the LGM (discarding too-thin geometries), (ii) the inland LGM elevation was lower than at present (discarding too-thick geometries), and (iii) the modeled grounding line advanced to the reconstructed LGM location for this region9. We reported the range of elevation changes as the 25th to 75th percentiles in the frequency distribution for this subset (Supplementary Fig. 2.9). From the 75 experiments, we chose the best-fitting one to show the pattern of elevation changes (Fig. 5), based on the available postglacial surface-exposure ages and the inland surface-elevation changes inferred from the Dome Fuji ice core (see Supplementary Note 2 for more details).

Regional sea-level modeling

We calculated the past regional sea-level change of DML through GIA modeling68,69, which allows calculation of the Solid Earth’s viscoelastic deformation due to rearranged ice and water masses based on theoretical equations70 incorporating the Earth’s rotation feedback, shoreline migration, and gravitational effects71,72,73. The modern topography from the ETOPO dataset74 was used for regions north of 60°S whereas we used 5-minute resolution Bedmap275 for the remaining areas to the south. Ice-load history was based on the ICE6G model42, which was resampled to 15 min spatial resolution for the calculation coupled with detailed bathymetry data. Because the original ICE6G spatial resolution is 1 degree, we interpolate the gridded data of ICE6G to calculate GIA-induced sea levels using the nearest-neighbor method. In the Earth model, we assumed a 100-km-thick elastic lithosphere, which is underlain by an upper mantle with a viscosity of (0.8–2) × 1021 Pa s, and a lower mantle with a viscosity of (3–30) × 1021 Pa s76.

We also calculated the sea-level change based on the refined ICE6G model that was fitted to our field constraints. For the Antarctic drainage systems 5–777, in which Gjelsvikfjella and the Sôya Coast are located (Supplementary Fig. 4.1), we delayed the timing of the main phase of ice-sheet thinning from ~11 to 9 ka, after which the ice-sheet surface is linearly decreased down to the present level at 5 ka (Fig. 6). The magnitude of thinning was kept unchanged. To examine the model sensitivity, we also calculated the sea-level changes based on three other versions of refinement in which (i) the timing of ice-sheet thinning is delayed to 10 ka rather than 9 ka, (ii) the refinement is introduced to all the Antarctic drainage systems rather than to the limited drainage systems of 5–7, and (iii) a combination of (i) and (ii). We found that the results are consistent among all refined models (see Supplementary Note 4 for more details).