Main

The West Antarctic Ice Sheet (WAIS) has long been considered susceptible to collapse. The capacity for the WAIS to rapidly change volume is linked to the marine ice-sheet instability—a nonlinear relationship between ice flux and ice thickness that drives reduced models of marine-terminating glaciers to either advance or retreat from grounding zones positioned on topography that deepens inland. The WAIS may be vulnerable to this instability in the marine basin confined by the Transantarctic Mountains, coastal mountains of Marie Byrd Land, the Whitmore block and the subglacial highlands of the Ellsworth and Thiel mountains where ice flows into Ross and Filchner–Ronne ice streams.

Despite the well-understood capacity for the WAIS to retreat into the interior, very little is known about the confining subglacial highlands that probably serve as ice-sheet refugia during deglaciation and as nucleation centres during readvance1,2. Three-dimensional model simulations of the WAIS, palaeo-sea-level change observations and reconstructions of ice volume during recent interglacials suggest that the ice sheet’s instability response may permit the ice sheet to deglaciate faster than it can reglaciate3,4. This asymmetry has been linked to large interior surface mass-balance fluxes required to promote WAIS regrounding across ocean-separated West Antarctic highland regions that include the Ellsworth Mountains and the coastal mountains of Marie Byrd Land2,5. Ice transport from the larger East Antarctic margin was instrumental for the early formation of the WAIS in the Miocene epoch. The ice-covered extensions of the Horlick and Thiel mountains near Hercules Dome have also long been recognized as key sites for ice-sheet regrowth1,6, probably facilitating ice transport to the WAIS when it was smaller than it is today.

Following previous studies that used subglacial morphology to understand ice-marginal behaviour when the Antarctic Ice Sheet was smaller5,7,8,9, we use new geophysical data collected near Hercules Dome to better constrain the past behaviour in the upstream reaches of the Ross and Filchner–Ronne ice streams (Fig. 1). Using high-frequency impulse and multi-element ice-penetrating radar, we produce three-dimensional (swath) bed topographies of the slow-flowing Antarctic interior (Figs. 1b and 2). The features revealed by the radar swaths combined with regional airborne gravity data6,10 provide new insight into tectonic controls on glaciation at the junction of East and West Antarctica and the initiation of past ice ages when ice near Hercules Dome flowed more rapidly than it does today.

Fig. 1: The subglacial topography and geology of West Antarctica.
figure 1

a, Overview of the basins and provinces with formation ages of outcrops in West Antarctica overlying the ice-sheet bed topography. b, An overview map of Hercules Dome with interpolated aerogravity surface anomaly shown with composite MODIS (moderate-resolution imaging spectroradiometer) imagery54. Also plotted are the radar profiles (light orange) from this study and flight lines of aerogravity data (light blue) collected with airborne radar observations as part of the PolarGAP survey6,10. Inset shows the locations of panels a and b in Antarctica. Black arrows denote flow axis inferred from bedforms and the valley axis imaged in subglacial topography. Dark orange lines indicate the approximate position of the basin troughs probably associated with past tectonism. Ma, million years ago. Panel a adapted with permission from ref. 40, Springer Nature Limited.

Survey overview and high-resolution subglacial topography mapping

Hercules Dome, located in East Antarctica at ~86° S, 105° W (Fig. 1a,b), is a compelling candidate location for a deep ice core11,12 that would provide information on the configuration of the WAIS during the last interglacial13, when global mean temperatures were as warm as they are today14,15. Ice-flow modelling studies show that the Hercules Dome region remained glaciated during the last interglacial period3, and climate model results show that significant lowering of the adjacent WAIS surface would cause atmospheric circulation changes that may be reflected in the isotopic composition of precipitation at Hercules Dome due to its proximity to the WAIS13,16. Ice-penetrating radar data collected across a transect of Hercules Dome during the US portion of the International Trans-Antarctic Scientific Expedition (ITASE) traverses revealed continuous layers to depths just hundreds of metres above the bed11,17. This supported conjectures that flow through parts of the dome has been stable throughout much of the last glacial period11. Satellite observations of the dome’s complex multi-axis structure18 support this possibility and suggest that the position of the dome is anchored by local subglacial topography (Figs. 1b and 2a).

To better characterize the subglacial landscape, we deployed multi-element swath radar19,20,21 and conventional high-frequency impulse ice-penetrating radar22,23 at Hercules Dome. Both radar systems were used to map englacial layering and along-track subglacial topography at ~5 m posting; where available, swath radar mapped cross-track topography at ~25 m posting in swaths ~2 km wide (see Fig. 1 for an overview of the impulse radar profiles and swath-mapped topography). These data cover previously unmeasured regions of the ice-sheet interior and can be connected to existing ice-penetrating radar data6,11,17 to extend our knowledge of the ice-sheet stratigraphy into the bottleneck between East and West Antarctica. We use polar stereographic (EPSG:3031) referenced directions (that is, grid east, grid west) to describe relative positioning within the survey. First, we characterize the large-scale features in the radar data before discussing smaller-scale landforms and evidence of subglacial water mapped in radar swaths. We then contextualize these new observations with existing aerogravity data collected by the PolarGAP survey6,10,24, past studies of the regional tectonic history and theoretical models of ice-sheet response time.

Subglacial valleys and bedforms

Both radars imaged deep U-shaped valleys at Hercules Dome, similar to those typically found in glaciated alpine environments (Figs. 24). The deepest and most pronounced valleys appear to be part of one contiguous feature stretching from the grid southwest region of the survey to the grid east. The dramatic relief (1,000–1,500 m over a few kilometres) continues to affect modern ice flow, promoting deflection greater than 5° in the englacial layer slopes in profiles crossing the valley (Fig. 4). Grid north of this feature, we image additional kilometre-scale features consistent with alpine glaciation, including hanging valleys and arêtes (Figs. 24).

Fig. 2: Subglacial alpine landscapes revealed in swath topography.
figure 2

a, Swath bed elevation overview created by geolocating off-nadir energy using the multi-element swath radar overlain over spline interpolated bed elevation55 in the white region shown in Fig. 1b. ae, Revealed in these radar swaths is evidence of U-shaped valleys, host to subglacial water (a); landscapes characteristic of glacier highlands (b) and alpine valleys (c); glacial sills (d, e) and lineations (e) within the Hercules Dome valley (d).

Within the large U-shaped valley system, Fourier analysis of swath radar topographies revealed younger sub-kilometre-scale landforms that appear to be associated with streaming flow and two candidate subglacial lakes that suggest the bed in the valleys is thawed (Figs. 3 and 4). The landforms include strongly tapered, elongated crags and tails and weakly tapered longitudinally symmetric bedforms consistent with drumlins, glacial lineations and glacial sills (Figs. 2 and 3 and Extended Data Figs. 13). Many of these features are thought to align with the flow direction during formation and can therefore be used to understand past ice-flow direction25. Across all radar swaths, windowed Fourier analysis of landform orientation suggests ice flowed rapidly along the subglacial valley axes (Fig. 3 and Extended Data Figs. 13).

Fig. 3: Swath image of the Hercules Dome subglacial valley with spectral roughness orientation.
figure 3

a, An overview of radar data collected at Hercules Dome. be, Spectral analysis of local subglacial topography in the glacier highlands (c,d) and valley floors (b,e). Each Fourier-analysed panel (b, c, d, e) corresponds to different areas shown and labeled in the map view (a). Black solid lines in Fourier-analysed panels of each outlined domain denote wavelength (in metres), and black dashed lines indicate orientation (azimuth with 45° marking northwest–southeast features and 135° marking northeast–southwest features). Colours capture the power (decibels) of spectral energy relative to the direction of the survey.

Fig. 4: Comparison of HD subglacial valleys and exposed postglacial alpine valley networks.
figure 4

a, Profiles of U-shaped valleys we image with nadir-focused multichannel radar and high-frequency impulse radar at Hercules Dome (note candidate subglacial lake water in trough bottom of the first radar profile). b, Swath radar digital elevation model of the trough at Hercules Dome. c, Offset elevation model captured by lidar56 of the Yosemite Valley. (Panels b and c have the same vertical scale.) d, Profiles correspond to the bed elevation profiles from Hercules Dome radio echograms shown in panel a (blue, black and grey) with a profile of El Capitan (dashed black) in the Yosemite Valley (c,e). Many large-scale (>1 km) landforms observed in these profiles are characteristic of deglaciated landscapes. Credit: photograph in e, Carol M. Highsmith Archive, Library of Congress, Prints and Photographs Division.

Gravity-field anomalies

The ice-flow direction imprinted in the swath-imaged valleys implies ice flow into larger basins that underlie the upper reaches of West Antarctic ice streams. Using radar and gravity data collected by the PolarGAP survey6,24, we evaluate the hypothesis that the large basins host low-density sediment columns by calculating gravity anomalies associated with density contrasts consistent with the modern ice thickness, basin geometry and geology using methods described by refs. 26,27,28. Local geology is unknown in the region of Hercules Dome, and thus substrate density and basement depth were treated as free parameters informed by outcrop observations from the Ellsworth and Thiel mountains29 and density data available for other Antarctic sediment basins30,31. Free-surface gravity anomalies range between –90 and 90 mGal, with well-defined lows coincident with topographic troughs imaged with high-frequency radar (Figs. 1b and 4 and Extended Data Fig. 4). Forward gravity calculations indicate that these topographic basins are probably underlain by thick (1,000–2,000 m) low-density (2,400 kg m–3) sediment columns (Fig. 5 and Extended Data Fig. 5).

Fig. 5: Tectonic controls on past and present ice flow.
figure 5

a, High-frequency impulse and swath radar data with flow paths through the Hercules Dome valley into the Patuxent Trough. The regional surface velocities57 are plotted in the background. b,c, Gravity modelling results of the PolarGAP transect (b) across the Patuxent Trough (c). In b, blue curves show the gravity effect of ice thickness, purple curves indicate the sum of the Moho gravity effect and orange curves represent the sum of these contributions, with sediment density hypotheses that can be compared with observations (black). Because gravity solutions associated with substrate density and depth are non-unique, we held the sediment thickness constant and used different hypotheses for sediment density to evaluate the likely presence of a sedimentary basin hosted within the Patuxent Trough.

The pattern of satellite and aerogravity data suggests that Hercules Dome is associated with a positive gravity anomaly, with ice probably anchored by remnants of an igneous province or potentially old cratonic bedrock (Fig. 5). The upstream extensions of Mercer, Whillans, Support Force and Foundation ice streams continue into basins (topographic troughs and negative gravity anomalies) that isolate Hercules Dome from East and West Antarctica (Figs. 1 and 5 and Extended Data Figs. 4 and 5). The valleys imaged with swath radar appear to be tributaries of these larger basins and, together with the density contrasts and the morphologic context of the Neoproterozoic rifted margin, suggest that alpine glaciers fed into palaeo ice streams that flowed from the East Antarctic margin to the Ross and Filchner–Ronne seas.

Discussion

New constraints on the geologic evolution of Antarctica

Our new swath radar observations reveal an alpine subglacial landscape at the interior origins of Whillans, Mercer, Foundation and Support Force ice streams in the central Transantarctic Mountains in unprecedented detail. Most of the subglacial features we image are inconsistent with the current slow outflow from Hercules Dome and resemble deglaciated alpine landscapes, such as the Yosemite Valley and marine-proximal glaciated valleys found in the Antarctic Peninsula, Alaska and Norway (Fig. 4). Ice-core evidence13,15,32, ocean sediment cores15 and far-field sea-level data33,34 support a smaller WAIS during the last interglacial (130–115 kyr before present) and marine isotope stage 11 (420–360 kyr before present). Still, the durations of these interglacials were probably too short to promote the incision of >1-km-deep valleys at Hercules Dome. This suggests that they were probably carved in the Pliocene epoch or earlier when the ice sheet was much thinner and when flow did not override locally high topography, driving localized erosion and valley formation35.

Within these valleys, we identify bedforms associated with streaming flow that are thought to require faster basal sliding speeds for formation than the ice surface speeds we observe at Hercules Dome today25,36,37. Most of our observations of glacial bedforms come from deglaciated environments and include features formed or modified during glacier retreat. Features associated with deglaciation (kettles, kames, eskers and moraines) represent the ‘death mask’ of a retreating ice sheet and have been used to understand past changes in ice-sheet extent. Similarly, observed bed features consistent with fast flow in glaciated environments that are not currently flowing quickly allow us to identify past regions of fast flow and candidate nucleation centres that promoted past ice-sheet advance. The landforms and large-scale topography all indicate that ice flowed rapidly along the trough axis at Hercules Dome. This flow orientation is at times orthogonal to the current outflow from the local ice divide. The inferred past flow direction suggests that ice from Hercules Dome fed into larger basins beneath palaeo ice streams associated with negative gravity anomalies. Even today, interior ice flow is influenced by these basins as the surface is depressed and surface velocities increase in these regions6 (Extended Data Fig. 4).

Evidence of past streaming flow in the interior of the continent and thick sediment deposits within the downstream basins are also consistent with the region’s tectonic history38,39,40,41. The host rock that anchors Hercules Dome and the nearby southern Transantarctic Mountains formed in a Neoproterozoic–Cambrian metasedimentary basin24,42. The outcrop ages of Thiel Mountain bedrock43, sedimentological, structural and geochemical studies42, and gravity and magnetic data10,24 confirm that the bedrock was deformed during the Cambro–Ordovician Ross orogeny. The major basins are imaged only in cross section by the PolarGAP aerogeophysical survey (Fig. 5 and Extended Data Fig. 4), but they resemble graben or wrench basins. The tributary valleys at Hercules Dome feed into these basins, which geographically align with the West Antarctic ice streams44. Ice flow is clearly topographically steered along the basins as it flows from the East Antarctic continental margin to West Antarctica6.

The potential fault-bounded nature of the basins is difficult to connect to a single stage of tectonism as any established crustal fault could reactivate in successive tectonic events. Large-magnitude extension/transtension occurred in the Late Cretaceous when the West Antarctic rift system and the Marie Byrd Land block separated from the East Antarctic Ross orogen as an active subducting margin and magmatic arc40. Over a span of just 20 Myr, the movement of these blocks contributed to over 600 km of extension in the West Antarctic interior, triggering melt and emplacement of the upper crust that promoted thermal subsidence to form the WAIS’s quintessential over-deepened basin38,40 and references therein. The faults formed during the development of the West Antarctic Rift System have been linked to rifting across the Siple Coast45, a region that also bears geophysical evidence of Neogene tectonism46. Such basin-bounding structures influence the position of marginal ice streams farther from the interior44,47. The position of the upper extension of the Mercer and Whillans ice streams and the palaeo-alpine glacier networks we image at Hercules Dome appear to be controlled by these larger tectonic basins.

Implications of the geology for ice-sheet evolution

Marine outlet glaciers advance or retreat in response to flux imbalances between discharge across glacier grounding zones and the integrated accumulation upstream48,49. Spreading ice shelves confined to fjords or pinning points can stabilize a growing ice sheet by buttressing grounded ice, which, in concert with sufficiently positive surface mass balance, can lead to additional ice flow from the glacier interior, resulting in further glacier advance. The pace of this change depends on how quickly interior surface mass-balance changes can be transported to the grounding zone. This ‘glaciological filter’, associated with properties of conveying ice streams, plays a vital role in determining the pace of early ice-sheet advance.

Guided by observations of swath-mapped lineations and topographic ridges in currently fast-flowing glaciated environments under Thwaites Glacier21,50, we infer a contrast in basal resistance fields associated with lineated and highland features of the Hercules Dome topography. We then estimate basal resistance for these two substrate types and incorporate them into an idealized outlet-glacier model51,52 to characterize how the fast-flowing valley glacier networks in our observations (Figs. 2 and 4) would have affected the interior ice-sheet response to perturbations in surface mass-balance forcing (Methods and Fig. 6).

Fig. 6: Implications of troughs and valleys on glacier response time.
figure 6

a, Illustration of an early WAIS with an open Ross Sea and alpine glaciers flowing from Hercules Dome into the neighbouring rift system as viewed from the Filchner–Ronne sector looking towards the Ross sector. b, Geologic cross section of the evolution of the West Antarctic rift system and the WAIS in the Late Cretaceous and the onset of glaciation. c, The equilibrium interior ice thickness estimated by an idealized model, associated with different parameterizations of glacier sliding. Stars indicate the resistances and associated thickness corresponding to rough highland-like features. Diamonds indicate the resistances and thickness inferred for lineated regions. The dashed black line indicates the approximate depth of the glacially incised trough, which provides a bound for the ice thickness when the trough was formed according to theory predicted by ref. 35. d, Interior transport timescale associated with different parameterizations of glacier sliding. Panel b adapted with permission from ref. 40, Springer Nature Limited.

Sampling a range of different sliding parameters and equilibrium estimates of interior surface mass balance that produce the same glacier length, we calculate characteristic ice thickness and glacier response time using the idealized model (Methods). We find that interior ice transport over a basal environment consistent with landforms we image at Hercules Dome reduces the equilibrium thickness and the characteristic timescale for glacier advance compared with the response time inferred for the modern configuration of the ice sheet (Fig. 6). The response times of the glaciers and ice streams that characterize committed changes in glacier length and interior thickness in response to climate perturbations44,45 also describe the efficiency of the glacier’s interior flux response to climate variability.

We see the filtering properties of the ice sheet expressed directly in the englacial layering in the trough that bisected Hercules Dome. Layers recorded in the high-frequency radar data collected at Hercules Dome and dated using legacy high-frequency ice-penetrating radar data traced to the South Pole ice core reveal patterns that change coherently with the glacially incised trough. Outside the glacially incised trough, the layers remain bed conformal while reflectors that cross the trough dip steeply below a dated package of reflectors that expand below surface layers that conform to a 13 kyr layer (Extended Data Figs. 6 and 7). This change in the shape of dated layers marks a transition when the ice sheet retreated from its last maximum position grounded across much of the Ross and Filchner–Ronne to the late Holocene configuration we are familiar with today. The relic glacially incised troughs that were probably formed during early WAIS formation continue to affect the flow behaviour of the ice-sheet interior and the response of the interior flux of the ice sheet to climate variability.

Conclusions

Our observations highlight the important and often underappreciated influence of Antarctica’s geologic history on the ice sheet’s capacity to change volume. Radar-imaged valleys at Hercules Dome probably hosted palaeo ice streams that conveyed ice to West Antarctica, potentially playing an important role in both seeding the early growth of the WAIS and facilitating subsequent readvances. The positions of sedimentary basins, probably controlled by the tectonic history of the continent, appear to precondition the past flow patterns we see imprinted on the floors of glacially incised valleys. Connections between past rifting and ice-stream flow have been observed near the present-day ice-sheet margins44,53 but had not been made in the interior, where streaming flow has implications for ice-sheet cycles. We conclude that past rifting that preceded or was activated during the extension event responsible for the over-deepened WAIS basin also facilitated early ice-sheet advance through interior basin-guided ice streams from mountain highlands adjacent to Hercules Dome. These features of the subglacial landscape affect the englacial layering and indicate that the glacially incised trough and tectonically controlled basins have continued to affect the flow behaviour of Hercules Dome over the most recent glacial cycle.

Methods

Impulse radar topography

An impulse ice-penetrating radar system operating at a centre frequency of 3 MHz was used to map bed topography and ice-sheet internal layering at nadir. See ref. 23 for comprehensive system details. This radar system is an evolutionary advance from earlier impulse radars22,58 with increased bandwidth, faster digitization, higher sampling frequency and increased dynamic range23. Impulse radar data processing steps include bandpass filtering (eighth-order Butterworth bandpass from 0.5 to 5.0 MHz), time correction for antenna separation, geolocation from dual-frequency Global Navigation Satellite System data, interpolation to a precise trace spacing of 5 m and along-track time–wavenumber migration (assumed radar wave speed in the ice of 169 m μs–1). An adaptive horizontal filter that subtracts the average low-frequency component in a moving 100-trace boxcar window with a vertical exponential taper was also applied at shallow depths to improve imaging of shallow stratigraphy. The basal reflector was digitized using a semi-automated routine that identified the Ricker wavelet corresponding to the ice–bed interface23,59,60.

Swath radar topography

Conventional ice-penetrating radars have limited or no ability to localize energy in the cross-track direction, so cross-track resolution has historically been dictated by the choice of cross-track line spacing when designing surveys. Swath processing of a multi-element antenna array allows energy localization in the along- and across-track directions. In the along-track direction, phase differences from sequential acquisitions are used to position subsurface returns, which results in a conventional nadir-focused image. However, by combining target position from phase differences in sequential along-track acquisitions with phase differences between receiver elements in the cross-track antenna array, energy can be positioned in both directions simultaneously, allowing three-dimensional mapping of the subsurface19,21. We applied the Multiple Signal Classification algorithm implemented in the Center for Remote Sensing and Integrated Systems toolbox19 to estimate the direction of arrival signals. Following ref. 61, we digitized the reflector in the cross-track direction to geolocate off-nadir reflection information from multichannel frequency-modulated pulsed radar data. This algorithm produced fine-resolution (~25 m posting) bed topography maps over large swaths (~2 km wide).

Sediment density and layer thickness gravity modelling

To identify spatial patterns in the aerogravity data (Fig. 1b), the gravity disturbance line data (gravity anomaly using ellipsoid heights) distributed by the PolarGAP survey were gridded using a continuous curvature spline with a tension factor of 0.3 at a horizontal resolution of 1 km. After large-scale features of interest were identified, we conducted forward modelling to identify whether sediment was present beneath topographic basins. The forward modelling of the gravity anomalies for bathymetry and sediment-layer thickness and density is based on the gravity-anomaly calculations of ref. 26. In this method, the domain of interest is discretized into rectangular prisms composed of ice, sediment and crystalline bedrock. The prisms used in the discretization have fixed horizontal dimensions (5 km × 5 km), but their depth and density can vary according to the different sediment distribution hypotheses we test. We assume a sediment density consistent with seismic and borehole data (2,400 kg m–3) and vary the basement depth according to the position of the basins identified by ref. 6. The total gravity anomaly at a given point is computed as the sum of the contributions from all prisms. We account for regional gravity-anomaly fields by calculating the long-wavelength Bouguer gravity anomaly from satellite gravimetry and assume that this gravity anomaly is caused by variations in crustal thickness. We compare solutions for the gravity anomaly with observations and the null hypothesis that the basins contain no sediment. In all cases, we improve the fit to the gravity observations when we include sediment in the basins compared with the simulation where we assume there is no sediment.

The effect of fault-bounded ice streams on timescales of ice-sheet advance

We use a reduced marine outlet-glacier model51 to assess what past streaming flow near Hercules Dome implies for past regional ice dynamics. The model describes an idealized outlet-glacier system with two dynamical stages. In the first stage, accumulation (S) falls in the interior catchment of length, L. This surface mass-balance flux L × S is balanced by ice flow towards the outlet-glacier margins (Q). The interior flux, Q, has the general form

$$Q={\left(\frac{{\rho }_{{\mathrm{i}}}g}{C}\right)}^{n}\frac{{H}^{\alpha }}{{L}^{\gamma }}.$$
(1)

where ρi is ice density, g is acceleration due to gravity and C is a friction coefficient assuming Weertman-type sliding; n is the flow-law exponent, H is the average interior ice thickness, and α and γ depend on the partitioning of stress accommodated by internal deformation and basal traction.

The interior reservoir of the outlet glacier drains according to this interior flux until it reaches a second conceptual reservoir near the grounding line. Here, stresses change as the glacier goes afloat, and the characteristic thickness hg corresponds to the flotation thickness. The grounding-zone flux, Qg, is given by

$${Q}_{{\mathrm{g}}}={{\Omega }}{h}_{{\mathrm{g}}}^{\beta },$$
(2)

where Ω and β depend on assumptions about glacier sliding and ice-shelf buttressing (for example, refs. 45,58). Two coupled equations describe the evolution of H and L as the glacier adjusts towards a balance of accumulation, ice flow and discharge:

$$\frac{{\mathrm{d}}H}{{\mathrm{d}}t}=S-Q/L-H/({h}_{{\mathrm{g}}}\times L)(Q-{Q}_{{\mathrm{g}}})$$
(3)
$$\frac{{\mathrm{d}}L}{{\mathrm{d}}t}=\frac{1}{{h}_{{\mathrm{g}}}}(Q-{Q}_{{\mathrm{g}}}).$$
(4)

Descriptions of these kinematic fluxes (L × S, Q, Qg) with time series of forcing (for example, perturbations in S or Qg) form a closed system of ordinary differential equations that describe glacier response to climate forcing. These equations can be linearized about small perturbations to equilibrium length and interior thickness to understand the dependence of the glaciological filtering properties of marine outlet glaciers51. Solutions to this linear system of equations in the case of stable geometries reduce to the sum of two linearly independent exponential functions with analytic expressions for eigenvalues that can be ordered to characterize the slow, τs, and fast, τf, response times of the conveying ice stream51,52. We focus on interior perturbations in cumulative surface mass-balance variability because the early advance of marine-terminating outlet glaciers must be driven by accumulation51. The system of equations can be projected onto its eigenmodes and sorted according to the timescales of the glacier’s length and thickness response, where ms is the slower mode and mf is the faster mode:

$$\frac{{\mathrm{d}}{m}_{{\mathrm{s}}}}{{\mathrm{d}}t}=-1/{\tau }_{{\mathrm{s}}}{m}_{{{{\mathrm{s}}}}}+{\sigma }_{{\mathrm{s}}}\,f(t)$$
(5)
$$\frac{{\mathrm{d}}{m}_{{\mathrm{f}}}}{{\mathrm{d}}t}=-\frac{1}{{\tau }_{{\mathrm{f}}}}{m}_{{\mathrm{f}}}+{\sigma }_{{\mathrm{f}}}\,f(t).$$
(6)

We focus on the slow mode, associated with the adjustment of interior ice, in our parametric evaluation of glacier advance sensitivity, which ultimately paces the full response of stable advance or retreat on millennial and longer timescales. The slow response time can be estimated51 as

$${\tau }_{{\mathrm{s}}}=\frac{\bar{H}}{\bar{S}-{s}_{T}\alpha }(\alpha +1/m+1-{s}_{T})$$
(7)

where sT is a parameter that describes the stability of the glacier response based on its geometry and grounding-line dynamics (via β). Total flux in the model is dictated by assumed accumulation and catchment area; we can then solve equation (1) for the implied thickness and use this thickness to calculate the slow response time. The equilibrium thicknesses and response times were calculated for combinations of sliding parameters and accumulation that resulted in the same equilibrium length using an implementation of the Adam gradient-descent inverse method. Refer to Extended Data Table 1 for a comprehensive list of parameter descriptions and values.

With a model that appropriately complements the weak constraints on ice geometry at the time of glaciation and a metric for evaluating the glaciological filter associated with the early ice-sheet response, we can understand how the valley troughs and trough basins associated with fast outflow affect ice-sheet nucleation and regrowth. Full-Stokes model simulations of Thwaites Glacier using high-resolution swath-mapped grids indicate that lineated features are typically associated with low basal drag environments compared with rough ridge-like features50, such as the glacial highlands we image at Hercules Dome. Here we evaluated the sensitivity of characteristic ice thickness and response time to different parameterizations of basal sliding and specifically assessed how the new features that we image at Hercules Dome, typically associated with rapid ice flow, may reflect conditions that promoted ice-sheet advance between collapse cycles.