A site for deep ice coring at West Hercules Dome: results from ground-based geophysics and modeling

Abstract

Hercules Dome, Antarctica, has long been identified as a prospective deep ice core site due to the undisturbed internal layering, climatic setting and potential to obtain proxy records from the Last Interglacial (LIG) period when the West Antarctic ice sheet may have collapsed. We performed a geophysical survey using multiple ice-penetrating radar systems to identify potential locations for a deep ice core at Hercules Dome. The surface topography, as revealed with recent satellite observations, is more complex than previously recognized. The most prominent dome, which we term ‘West Dome’, is the most promising region for a deep ice core for the following reasons: (1) bed-conformal radar reflections indicate minimal layer disturbance and extend to within tens of meters of the ice bottom; (2) the bed is likely frozen, as evidenced by both the shape of the measured vertical ice velocity profiles beneath the divide and modeled ice temperature using three remotely sensed estimates of geothermal flux and (3) models of layer thinning have 132 ka old ice at 45–90 m above the bed with an annual layer thickness of ~1 mm, satisfying the resolution and preservation needed for detailed analysis of the LIG period.

1. Introduction

The West Antarctic ice sheet (WAIS) is vulnerable to rapid retreat because it is grounded below sea level (e.g. Mercer, 1978; Bamber and others, 2009), and has collapsed at least once in the Pleistocene (Scherer and others, 1998). During the Last Interglacial (LIG, 116–132 ka), global eustatic sea level was higher than today, though the timing and magnitude of the sea level highstand remain uncertain. Estimates of the upper limit of sea level during the LIG range from 4 m (Dyer and others, 2021) to 9 m higher than today (Dutton and others, 2015). A collapse of the WAIS would help explain the LIG sea level, but there is no unequivocal evidence that this occurred.

Ice core records provide one possible means to obtain information about past changes in the WAIS. Existing ice core records from West Antarctica either do not contain ice dating to the LIG (e.g. WAIS Divide, Buizert and others, 2015) or, if they do, the ice is very near the bed and is difficult to interpret (e.g. Siple Dome, Brook and others, 2005; Roosevelt Island, Lee and others, 2020; Skytrain Ice Rise, Mulvaney and others, 2020). An ice core site from a location that would be sensitive to the climatic changes induced by WAIS collapse, but where the ice flow would not be directly affected, has the potential to elucidate the size and rate of change of the WAIS during the LIG.

Hercules Dome (86° S, 105° W) is located in East Antarctica ~400 km from the South Pole towards the Horlick and Thiel ranges in the Transantarctic mountains. The primary moisture source region for snowfall at Hercules Dome is the Pacific sector, such that the moisture pathway passes over the WAIS; in contrast, the moisture pathways for existing East Antarctic ice cores that preserve LIG records do not pass over the WAIS (Sodeman and Stohl, 2009; Nicolas and Bromwich, 2011). Hercules Dome is not likely to have changed elevation significantly in elevation from the LIG to present even in scenarios of WAIS collapse (Sutter and others, 2016), with changes of <75 m in some models (e.g. Pollard and DeConto, 2016; Garbe and others, 2020). On the other hand, the composition of snowfall at Hercules Dome should be sensitive to changes in the elevation of the WAIS (Steig and others, 2015). Simulations with multiple atmospheric models show that a lowered WAIS topography would cause a cyclonic anomaly which changes regional atmospheric circulation, bringing a greater flow of warm, moist air across the WAIS toward the South Pole (Steig and others, 2015; Holloway and others, 2016). Based on these simulations, an ice core from Hercules Dome would provide a climate record that is influenced by the size and extent of the WAIS.

Climate records from a deep ice core at Hercules Dome could be used in conjunction with other paleoclimate records to constrain the size and rate of change of the WAIS during the LIG. Indicators of past ice-sheet extent, such as subglacial bedrock exposure (e.g. Hein and others, 2016; Spector and others, 2018), lack information on the rate of change which is likely preserved in ice core records. A record from Hercules Dome would complement other ice core records both directly affected by WAIS topography changes, such as Skytrain Ice Rise (Mulvaney and others, 2020), and indirectly affected, such at Mount Moulton (Korotkikh and others, 2011). Indeed, Steig and others (2015) showed if the WAIS collapsed, surface temperatures at Mount Moulton would have warmed less, and Hercules Dome would have warmed more, than sites farther inland on the East Antarctic plateau. Therefore, existing ice core records from the East Antarctic plateau provide independent information against which the Hercules Dome record could be compared.

Hercules Dome was first identified as a promising location for a deep ice core during the 2001/02 US-ITASE traverse from West Antarctica to South Pole (Jacobel and others, 2005). Ice-penetrating radar was collected along the main traverse route as well as two perpendicular 50 km long tracks. The radar revealed clear internal layering with ice thicknesses that ranged from ~1500 to 2800 m. A 72 m firn core was also recovered. Analyses of the firn core show that annual layering and volcanic signals are well preserved and the 400 year average accumulation rate was 0.12 m a⁻¹ ice equivalent (Dixon and others, 2013).

We conducted a geophysical survey across the Hercules Dome region during the 2018/19 and 2019/20 austral summers. Recent satellite observations (Cryosat-2, ICESat2 and Maxar Imagery) revealed a complex surface topography (Fig. 1). There are three local surface divides, which we refer to here by their relative positions in polar-stereographic space with respect to grid north (Fig. 1). ‘West Dome’ (85.8° S, 102.9° W) lies furthest to the true north, and is connected to ‘East Dome’ by a ~70 km ridge. A secondary ridge extends approximately grid south from East Dome that we term ‘South Ridge’. East Dome, where the US-ITASE traverse visited, has large variations in bedrock topography and lacks a clear dome summit; East Dome is the current source region for South Ridge. We therefore focus on West Dome where our geophysical data suggest a simple flow regime and small bedrock topography variations (Fig. 2), and assess the suitability of West Dome as a site for recovering a deep ice core that preserves LIG ice. We focus on a 10 km long radar line (named x135 within our grid) which nearly transects the dome summit, and describe the results of ice-flow modeling constrained by near-surface, deep-profiling and phase-sensitive radar.

Fig. 1. Overview map of the Hercules Dome geophysical survey in polar stereographic (PS) coordinates. The directions in the PS projection are used in the remainder of the paper. Arrow indicates direction of geographic north. The inset shows the two radar lines, x135 and x133, at West Dome. Surface elevation contours from the Reference Elevation Model of Antarctica (REMA, Howat and others, 2019) are shown at 10 m spacing for the overview and 2 m spacing for the inset.

Fig. 2. (a) GNSS derived surface elevation along x135 with ApRES sites shown by colored dots: n0 (yellow), n1000 (orange), n2000 (pink), n5000 (purple). (b) VHF radar of line x135 with arrows pointing to bright layer identifiable throughout survey. Data quality is lower to the north. (c) The depth of the bright layer on x135 agrees well with x133, which does not have the data gap or quality issues. VHF data along x133 is used for the accumulation forcing in the kinematic model. (d) HF radar of x135 showing internal layer and bed reflections. (e) Traced internal layers of x135 radargram shown in (d).

2. Geophysical data and model setup

The three radar systems used at Hercules Dome in the 2018/19 and 2019/20 field seasons as well as the kinematic model setup used to help interpret the observations are described in this section; Section 3 will present the results. The extent of each survey is indicated in Figure 1. We use stereographic grid directions, following the orientation of Figure 1, in the remainder of the paper. This paper focuses on the radar data collected near the West Dome summit. Two 10 km long lines across the divide were collected (Fig. 1, inset); x135 is closest to the dome summit and x133 is 2 km to the east.

2.1 Deep (HF) radio-echo sounding surveys

Radio-echo sounding (RES) has been used to map ice thickness and bed elevation, basal reflectivity and ice-sheet internal stratigraphy along two 10 km lines at West Dome (x135 in Fig. 2). We used a custom high-frequency (HF, 3 MHz) monopulse radar system (Welch and others, 2003, 2009; Christianson and others, 2012, 2016). The RES data were processed using standard techniques, including horizontal stacking, bandpass filtering, correction for antenna separation, geolocation from contemporaneous kinematic Global Navigation Satellite System (GNSS) data, interpolation to standard trace spacing and time-wavenumber migration using the open-source radar-processing and interpretation toolbox ImpDAR (Lilien and others, 2020). Englacial reflections were picked through the semi-automated digitization routine in ImpDAR. The two-way travel times (TWTT) were then converted to depth using a spatially uniform firn-density profile (ITASE 02-4, Dixon and others, 2013) and a Looyenga and others (1965) mixing model (Fig. 2d). There is a ~100 m data gap at 4.9 km on x135. The noise on the picked layers in Figure 2e is due to sample resolution and the picking algorithm in ImpDAR.

2.2 Shallow (VHF) RES surveys

To image the upper 20–200 m of the ice sheet, we used a Pulse-EKKO Pro radar system with very high-frequency (VHF, 100 MHz) antennas. Our VHF data-processing flow is similar to that described for the HF radar and includes bandpass filtering, correction for antenna separation, geolocation from contemporaneous kinematic GNSS data, interpolation to standard trace spacing and exponential range gain to increase the relative amplitude of deeper reflections. Englacial reflectors were traced using the semi-automated reflection digitization routine in ImpDAR. A distinct high-amplitude reflector at ~75 m depth (Fig. 2b), observed in all the shallow RES profiles, is used to infer the spatial pattern of accumulation at West Dome.

We convert the TWTT to depth with a depth–density profile based on the firn core density measurements; we fit a Herron and Langway (1980) model to the firn core measurements to produce a smooth depth–density profile. The bright layer is at 73.5 m depth at the closest intersection (within 1 km) of our VHF profiles and the US-ITASE 02-4 core (Dixon and others, 2013; Fig. 1). The core extends to 72 m depth and was collected 17 years prior to our VHF survey. The age at 73.5 m depth is 1600 CE, or 420 years before the VHF radar data were collected, where we have assumed a steady-state firn column and determined the additional age in the bottom 1.5 m from the same model fit used for the TWTT conversion.

For the x133 and x135 radar lines, the conversion of TWTT for the bright layer to depth (Fig. 2c) and the inference of accumulation rate from the layer depth both depend on the density profile. Because there is no firn core at West Dome, we use a Herron and Langway (1980) firn model to estimate the depth–density profile. We use fixed inputs of 370 kg m⁻³ for the surface density and −40°C for the surface temperature; we start with an initial value of 0.13 m a⁻¹ for the accumulation rate (all accumulation rates are given in ice equivalent and be converted to a kg m⁻² a⁻² using a density of 917 kg m⁻³), which is greater than the value inferred from the ITASE core because the average bright layer depth at West Dome is deeper than at the ITASE core site. We then infer the accumulation rate from the layer depth and the modeled density profile. To find a self-consistent solution, we iterate by calculating a new density profile along the radar line with the inferred accumulation rate at each trace, convert TWTT to depth and infer the updated accumulation rate with the new density profile. The solution becomes self-consistent within three iterations.

2.3 Automated phase-sensitive RES (ApRES)

Acquisitions of autonomous ApRES (Nicholls and others, 2015) were performed at four locations along the x135 line in the 2018/19 season, and repeated in 2019/20. ApRES measurements are used to determine the relative phase offset of englacial reflectors between the two acquisitions. ApRES processing is implemented in ImpDAR (Lilien and others, 2020) and follows the processing flow described by Brennan and others (2014). The wave speed is not corrected for the density variations of the firn and we do not interpret the data in the upper 150 m. The phase offset is converted to a relative range offset which determines a depth–velocity averaged over the period between the two measurements. Since phase is relative, the phase offset profile is unwrapped against a point of known velocity such as the bed or surface which then can shift the velocity profile by a constant.

2.4 Kinematic model setup

We employ a kinematic model to help interpret the internal structure of West Dome. In this section we describe the model setup and the application of the model follows in Section 4. The kinematic model uses inputs of the accumulation rate, bed topography and prescribed shape functions for the horizontal and vertical velocities (Nereson and others, 2002; Waddington and others, 2007; Koutnik and others, 2012). These model inputs are based on the geophysical observations described in Section 3. The ability to fit the measured layers provides information about the stability of the current flow regime near the divide and about the appropriate forcing for depth–age modeling.

We choose a kinematic model for two primary reasons. First, the existing radar coverage is limited. The primary radar line is 10 km in length, which results in a distance of only ~3 ice thicknesses between the divide and downstream boundary. Using dynamic models requires extending the lower boundary to ~10 ice thicknesses from the divide, so that the area of interest is not affected by the edge boundary (Martín and others, 2006). Second, dynamic models have difficulty reproducing divide flow, particularly in areas with variations in bed topography (Goel and others, 2018). Therefore, we use a kinematic model where we can specify the shape of velocities and the position of the divide.

The horizontal velocities, u, are calculated as:

(1)$$u( {x, \;\hat{z}, \;t} ) = \bar{u}( {x, \;t} ) \phi ( {x, \;\hat{z}, \;t} ) $$

where $\bar{u}$ is the depth-averaged velocity, x is the distance along the flowline, $\hat{z}$ is the relative height above the bed, t is the time and ϕ is the horizontal velocity shape function. We use the Lliboutry (1979) approximation:

(2)$$\phi ( {\hat{z}} ) = \displaystyle{{\,p + 2} \over {\,p + 1}} \Big ( {1-{( 1-\hat{z}) }^{\,p + 1}} \Big) $$

where p is the shape factor. The corresponding vertical velocity shape function is:

(3)$$\psi ( {\hat{z}} ) = \left({1-\displaystyle{{\,p + 2} \over {\,p + 1}}( {1-\hat{z}} ) + \displaystyle{1 \over {\,p + 1}}{( 1-\hat{z}) }^{\,p + 2}} \right).$$

The depth-averaged horizontal velocity is calculated from the ice flux from upstream, based on the accumulation rate and ice thickness. We assume no divergence in the ice flow (i.e. the flowband width is constant) based on comparisons with the measured surface velocities along x135 (Fig. S1). The vertical velocity is calculated as

(4)$$w = ( {\dot{b}-\dot{S}} ) \psi -\bar{u}\phi \left({\displaystyle{{\partial B} \over {\partial x}} + \hat{z}\displaystyle{{\partial H} \over {\partial x}}} \right)-\bar{u}H\mathop \int \limits_0^{\hat{z}} \displaystyle{{\partial \phi } \over {\partial x}}\partial \hat{\zeta }$$

where w, ψ, ϕ vary in x,$\;\hat{z}$ and t; $\dot{b}, \;\dot{S}, \;\;\;\bar{u}$ vary in x and t; B and H vary in x and $\dot{b}$ is the accumulation rate, $\dot{S}$ is the change in surface elevation, B is the bed elevation and H is the ice thickness.

An example of the modeled velocity field is shown in Figure 3. The details of our model results are presented in Section 4, but we present the individual velocity components, the three terms on the right-hand side of Eqn (4), to aid in visualizing the model construction. By representing the velocity field with these three terms, we can separate the influence of the bed from dynamics within the ice column to better understand the cause of variability in ice-flow along profile.

Fig. 3. (a) Modeled vertical velocity field using Eqn (4) with the divide at 4.7 km (black triangle). (b) First term of Eqn (4), which has the 1-D effects and is set by the accumulation rate and shape function. (c) Second term of Eqn (4), which expresses the vertical velocity for bed-parallel horizontal flow. (d) Third term of Eqn (4), which expresses the impact of the varying horizontal shape function. All panels use the same color scale to show the relative importance of the terms but have contours of different values to emphasize the magnitude of the terms.

We refer to the first term on the right-hand of Eqn (4) side as the 1-D term, which is governed by the surface vertical velocity. The surface vertical velocity (w _s) is given by the difference between the accumulation rate ($\dot{b}$) and the change in the surface elevation ($\dot{S}$): $w_{\rm s} = \dot{b}-\dot{S}$. We do not consider basal melt (see Sections 3.4 and 4.4).

We refer to the second term on the right-hand side of Eqn (4) as the 2-D term since it ensures bed and surface parallel flow.

The third term on the right-hand side in Eqn (4) accounts for variations in the shape of the horizontal vertical velocity, which in our case is the transition from divide to flank flow.

The bed elevation and ice thickness (Section 3.2) as well as the accumulation rate (Section 3.3) used in the kinematic model are derived from the radar data. The prescribed shape functions are based on vertical velocity profiles fit to the ApRES data (see Section 3.4). The ApRES measurements are at only four locations and we must define the vertical velocity continuously along the flowline. To do so, we define one shape function for the divide and one for the flank, and prescribe a smooth transition between them. We follow Nereson and Waddington (2002) and assume the width of divide flow influence can be approximated as:

(5)$$\beta ( x ) = \exp \left({-\displaystyle{{x^2} \over {2\sigma^2}}} \right)$$

where σ is the half-width of the divide zone for the vertical velocity and the transition for the horizontal velocity is then:

(6)$$\alpha ( x ) = \displaystyle{{\sqrt {2\pi } \sigma } \over {2x}}{\rm erf}\left({\displaystyle{x \over {\sqrt 2 \sigma }}} \right)$$

The transition in the vertical velocity occurs in 1–2 ice thicknesses from the divide, based on modeling studies (Raymond, 1983) and observations of the width of Raymond arches (e.g. Matsuoka and others, 2015). The transition in horizontal velocity occurs over a longer distance. The value of p for the shape functions (Eqn (3)) are 0.3 for the divide and 7 for the flank. The divide-zone half-width is 2800 m, ~1.5 ice thicknesses.

3. Geophysical observations

The geophysical observations reveal the internal structure and bed of West Dome and provide inputs and constraints for the kinematic model. We first determine the internal stratigraphy from the HF radar reflections. Second, we use the HF radar to determine the bed elevation and assess ambiguities in the bed picks using the kinematic model and measured internal layers. Third, we find the accumulation rate across the dome using a bright layer in the VHF radar which can be dated to 420 years before 2020, when the data were collected, near the US-ITASE shallow core. Finally, we assess the englacial vertical velocity profiles and identify a change in shape with distance from the present ice divide.

3.1 Internal layer geometry

We traced 22 internal layers from ~400 m below the surface to 200 m above the bed (Figs 2d, e). Layers are consistent and traceable along the full profile to within 400 m of the bed. Below 400 m, there is a difference in the character of the layers on the north and south sides of the divide. To the north, there is a bed-conformal and traceable layer between 100 and 200 m above the bed that extends to the northern end of the profile. Between that layer and the ice bottom, there appears to be coherent layering just north of the divide, but that signal is lost as the ice thickens along the profile. The lack of continuous layering in this interval is not surprising given the challenges of resolving deep layers (e.g. Drews and others, 2009): returns from the deep layers are weaker due to greater attenuation; the layers vary more in depth as they drape the irregular bedrock and the relatively long wavelength (~55 m) of the HF system integrates a larger number of small-scale electrical contrasts which can cause more complex reflections. To the south, there is a package of layers that diverge from the traced layer at 400 m, indicating a local thickening of the layer package that sits between 200 and 400 m above the bed at the divide as you move south along the profile. These layers lose coherence as they approach the ice bottom.

The internal layers cannot be dated with existing information from ice cores because there are no continuous radar profiles between the two West Dome radar profiles and the ITASE ground-based or airborne radar datasets. The only layer that can be dated is the bright reflector in the VHF data (Section 2.2, Fig. 2b) because it can be unambiguously matched by correlation to a radar line within 1 km of the dated ITASE core despite the lack of a continuous link. The internal reflectors of the HF radar profiles at West Dome are not distinct enough to be matched to the East Dome and South Ridge HF radar profiles.

3.2 Bed topography

We imaged the bed of West Dome along the x135 and x133 radar lines using the HF radar. The current surface divide is located above a bedrock high (Fig. 2). The ice thickness along x135 ranges from 1510 m below the divide summit to 1750 m at ~8 km along the radar line. The bed topography is smoother than at East Dome (Jacobel and others, 2005). The bed reflection is visible along the full profile (Fig. 4), although there are two locations where the bed reflection is ambiguous. The most important of these is at 6 km along the radar line. The brightest reflection is at 1640 m depth; however, there is also a weaker return from 1580 m depth. This smaller amplitude reflection is likely an off-nadir return and may indicate a decameter-scale bump at the bed close to the profile.

Fig. 4. Interpreted bed location is shown (red). There is an ambiguity in the bed reflection at 6 km (arrow).

To evaluate the significance of such a feature on the internal layering, we run the kinematic model with bed topographies based on both the original pick and by picking the decameter-scale bump at ~6 km. We also use two different low-pass filters of 500 and 1000 m, resulting in four different potential boundary conditions. The observed layers (and the associated model-predicted layers) deeper than 1000 m are shown in Figure 5 because they are the most affected by the bed boundary condition. The age of the internal layers is unknown so we cannot compare the modeled and measured layers based on age; instead, we compare the observed layers to a modeled layer with the same average depth. The modeled layers produced using the 500 m filtered beds show greater variability than the observed layers while the 1000 m filtering produces internal layering with smoothness comparable to the measured layers. Using the bed with 1000 m filtering and the decameter-scale bump, the model is able to reproduce the slight dip in the deepest layer at 5.5 km along the radar line. Therefore, we use the 1000 m filtered bed with the bump in all of the model runs used to help interpret the internal layering (Section 4).

Fig. 5. Black lines are observed internal layers. The bottom layers are the four different beds used in the kinematic model. The two brown lines have no bump at 6 km; the two blue lines have the bump added. Low-pass filtering of 500 m is shown as solid lines and 1000 m is shown as dash dot. The modeled internal layers have line color and type that match the bed.

3.3 Accumulation rate

The accumulation rates inferred from the VHF radar profiles along x135 and x133 are shown in Figure 6. There is a ~2 km gap in the x135 line because of poor radar quality. The inferred accumulation rates are similar between the two lines. The highest accumulation rates of 0.14 m a⁻¹ are in the north and are relatively constant from 7 to 10 km along the radar lines. The accumulation rates decrease by ~25% to <0.115 m a⁻¹ at the south end of the radar line. Because the data from x133 is of higher quality and agrees well with the high-quality portions of x135, we use the accumulation rate inferred from x133 as the forcing for the model runs used to help interpret the internal layering (Section 4).

Fig. 6. Accumulation rate inference from the bright layer dated to an age of 420 years (1600 CE). Profile x135 has poor data quality in the center so x133 is used for model forcing given the close agreement between the two. Note that the x133 profile is slightly offset at 5 km (Fig. 1, inset) such that we interpolate linearly from 4.8 to 5.1 km.

3.4 Englacial vertical velocities

Divides with frozen beds have vertical strain rates which are greater near the surface and smaller near the bed compared with locations of larger surface slopes and shear strain rates. We refer to this characteristic shape of the vertical velocity near divides as divide-flow, and the more linear vertical velocity away from divides as flank-flow. The difference in the vertical velocity patterns yields Raymond arches (Raymond, 1983) when the divide position is stable. Englacial vertical velocities were derived using ApRES acquisitions in the 2018/19 and 2019/20 field seasons at four sites shown in Figure 2a with the goal to (1) identify whether the characteristic shape of the vertical velocity is characteristic of those with a frozen bed and (2) provide constraints for kinematic modeling of internal layers.

The vertical velocity profiles are plotted by relative height above the bed in Figure 7a (where 0 represents the bed and 1 the ice surface) and by absolute height above the bed in Figure 7b. A useful way to visualize whether divide flow is occurring is to perform a linear fit of the upper portion of the vertical velocities and to compare the extension of the linear fit to the bed with the measured vertical velocities (Fig. 7a). The fit is calculated by linear regression between a relative height of 0.5–0.9. The near-surface portion is excluded because of the influence of firn compaction on the velocity profile above 0.9. For n5000, ~5 km north of the divide, the linear fit matches well for all velocities except for those below 0.05. This near-linear increase of vertical velocity with height above bed is typical of ice away from a divide (Kingslake and others, 2014). At sites n0 and n1000, there is a noticeable divergence from the linear fit, as the vertical velocities show little change below 0.3H.

Fig. 7. Modeled and measured vertical velocity profiles for the four ApRES locations labeled by their approximate distance in m from the divide. The horizontal spacing of the vertical lines in both panels is 0.05 m a⁻¹. The positions are shown in Figures 1 and 2, with n0 in yellow, n1000 in orange, n2000 in pink and n5000 in purple. (a) ApRES vertical velocities; black lines are a linear fit to the interval of 0.5H to 0.9H. (b) Lliboutry fits to the vertical velocities where the surface vertical velocity (excluding firn compaction) is required to be within 0.02 m a⁻¹ of the 420 year average accumulation rate, the vertical component of ice flow over uneven bedrock is included and a constant shift for the vertical velocity profile is a free parameter. The transition from divide-flow with low p values to flank-flow with high p values occurs with increasing distance from the divide (see Section 2.4).

The four sites were revisited in the 2021/22 field season. Sites n0, n1000 and n5000 showed nearly identical vertical velocities for the 19/20–21/22 period as the 18/19–19/20 period; however, n2000 did not (Figs S2–S5). Therefore, we focus our discussion on the differences between n0 and n1000 near the present divide and n5000, which is outside the divide region. A discussion of the vertical velocities using the 2021/22 field data are in given in the Supplementary material.

To quantitatively describe the vertical velocity profiles, we fit the vertical velocities using two simplified versions of Eqn (4). First, we follow Kingslake and others (2014) and fit with the following 1-D term,

(7)$$w( {\hat{z}} ) = w_{\rm s}\psi{\hat{z}} )$$

The value of p in the calculation of ψ (Eqn (4)) is a useful measure of the shape of the vertical velocity and can be used to characterize the influence of the divide. Small values of p indicate divide-like flow, with increased vertical strain near the surface and reduced vertical strain near the bed. Higher p values indicate more uniform vertical strain with depth.

Since the ApRES measures phase, we can only infer a relative velocity, and we must specify one absolute velocity. We assume no motion at the bed and set the average of the velocities within 100 m of the bed to zero. This allows an estimate of the surface vertical velocity, excluding firn compaction, which should be similar to the ice-equivalent accumulation rate. The best-fit values of p and w _s are given in Table 1; the velocities within 100 m of the bed and within 150 m of the surface are not used in the fit. The inferred surface vertical velocity exceeds the inferred 420 year average accumulation rates from the VHF radar data, except at n2000. The inferred surface vertical velocity at n5000 is 0.17 m a⁻¹, nearly 0.04 m a⁻¹ (30%) greater than the 420 year average accumulation rate of 0.13 m a⁻¹. We suspect that this is because n5000 has the greatest horizontal velocity (~0.5 m a⁻¹) and is expected to have the largest component of the vertical velocity due to horizontal flow over the bedrock.

Table 1. Fits to ApRES vertical velocities

Bold values are the parameters inferred.

To account for ice flow over an irregular bed, we also fit the vertical velocities with

(8)$$w( {\hat{z}} ) = w_{\rm s}\psi + \bar{u}\phi \displaystyle{{\partial B} \over {\partial x}}$$

We neglect the component of the vertical velocity due to lateral variation in the shape of the horizontal velocity profile because it is the smallest term (Fig. 3). We estimate the depth-averaged horizontal velocity from the balance velocity. The bed slopes are calculated on a 200 m length scale using the unfiltered bed (Fig. S6).

In the fit using Eqn (8), the surface vertical velocity is constrained to within 0.02 m a⁻¹ of the measured accumulation rate. In addition, the fit finds a uniform shift of the vertical velocity profile since the velocities are relative; this relaxes the constraint we applied when using Eqn (7), where the average velocity in the bottom 100 m was set to 0, and allows us to determine if the vertical velocities are consistent with the glaciological setting. The fit values are also shown in Table 1.

The largest difference in the fit parameters between Eqns (7) and (8) is at n5000, where a reasonable fit with a surface vertical velocity of 0.151 m a⁻¹ is possible because of the negative (upward) vertical velocities near the bed caused by flow up the bedrock slope (e.g. Hills and others, 2022). The two measurements closest to the divide, n0 and n1000 are similar to each other with the value of p being indicative of divide-flow; n5000 has a large p value indicative of flank-flow. Because of the uncertainty in the measurements at n2000 and the lack of measurements in between n2000 and n5000, it is difficult to determine the width of the divide zone.

4. Kinematic model results

Our kinematic model and data comparison focuses on five objectives: (1) fit the shallowest layer using the kinematic model to determine how well the modern observations reflect late Holocene dynamics and forcing; (2) assess differences in layer depths in the divide region to evaluate the presence of a Raymond arch; (3) fit all layers across the full radar profile to evaluate the likelihood of divide migration within our domain; (4) determine the likelihood of basal melt and (5) use our best fit model to produce a depth–age scale that captures the annual-layer thickness near the bed.

4.1 Fit to shallowest layer

Before examining the model fit for all of the internal layers, we focus on the shallowest HF radar layer which has a depth of ~430 m and a modeled age of ~4 ka. We assess whether the modern observations of accumulation and vertical velocity are likely to have been constant for the past ~4 ka. The accumulation rate pattern is a 420 year average while the measured vertical velocities are a snapshot over 1 year.

Using the modern patterns of accumulation and vertical velocity for the past 4 ka produces a layer which is too shallow beneath the divide and too deep north of the divide (Fig. 8). This suggests that the modern accumulation and vertical velocity patterns have not been constant for the past 4 ka, which may indicate that the divide has not been stable at its current position. To evaluate this possibility, we test scenarios of divide onset occurring in 1 ka increments. Prior to the divide onset, we assume spatially uniform flank-flow vertical velocities and uniform accumulation. The modeled layers and the misfit are shown in Figure 8. The modeled layers fit the measured layer best with a 2 ka onset of modern conditions although there is still a significant discrepancy on the south edge of the profile.

Fig. 8. Model fit to the shallowest HF (deep) radar layer. The legend gives the spatially averaged misfit to the measured layer (in m) and the age at which the modern vertical velocity and accumulation rate patterns are applied (in ka). The divide zone has a width of 1.5H. The modeled age of the layer is ~4 ka.

For these tests, we have assumed that the divide in the surface topography causes both the vertical velocity pattern and the accumulation pattern. The vertical velocity pattern is primarily caused by the low deviatoric stress beneath the divide although the development of ice fabric may also play a role (Pettit and others, 2007; Martín and others, 2009); however, it is less clear how closely the surface topography of the divide and the accumulation pattern are coupled. More complicated temporal and spatial variations in the accumulation pattern and divide position can improve the fit; however, without evidence to guide a priori expectations for how the patterns should develop in time, the improvements would not necessarily be meaningful. An improved fit with more complicated variations in the forcings would reinforce the primary conclusion that the modern conditions are unlikely to have been stable for the past ~4 ka. In the following sections, we explore the fit to all of the internal layers to extend our understanding of the divide position and accumulation pattern to older ages.

4.2 Evaluating the presence of a Raymond arch

Although the ApRES vertical velocities indicate that divide-flow is occurring, there is not an obvious Raymond arch beneath the divide. However the measured internal layers are shallower in the ice column at ~6 km than might be expected based on the bed topography (Fig. 2d). In areas of flat bedrock, the amplitude of a Raymond arch can be determined by comparing the depth of an internal layer beneath the divide to the depth on the flank outside of the divide region (e.g. Roosevelt Island, Conway and others, 1999), but the identification of Raymond arches is challenging in areas with greater ice thickness, lower accumulation rates and without smooth bedrock (e.g. Goel and others, 2018). Since the irregular bedrock at our site does not allow a simple comparison between the depth of a layer beneath the divide and on the flank, we use the kinematic model.

We compare the height of the measured layer to one from the kinematic model forced with a spatially constant vertical-velocity shape function (p = 4, Fig. 9). The difference between the measured and modeled layers at 5.5 km along the radar line has both similarities and differences to the Raymond arch amplitude at Roosevelt Island (Fig. 9b). The Raymond arch at Roosevelt Island began forming ~3 ka (Conway and others, 1999; Martín and others, 2006). The onset of a layer mismatch begins deeper in the ice column, at a normalized height of ~0.7 at West Dome compared to 0.9 at Roosevelt Island. There are also large misfits in the layer depths outside of the divide zone at West Dome, particularly to the south (Fig. 9a). Thus, the existence of a Raymond arch is uncertain. In Section 4.3, we model the full set of internal layers to better understand the observed layering.

Fig. 9. (a) Modeled internal layers (blue) assuming uniform flank flow (p = 4) compared to measured layers (black). The bed is shown in brown. (b) The difference between modeled and measured layers (blue circles) at 5.75 km (black vertical line) compared to the Raymond arch amplitude of Roosevelt Island (red squares) (Conway and others, 1999).

4.3 Model fit to all layers

The kinematic model fit to the shallowest layer (Section 4.1) suggested that the modern divide and accumulation pattern have not been stable for more than a couple of thousands of years while the possibility of a Raymond arch exists (Section 4.2) offset to north from the present divide. Here we evaluate the fit to all the layers for test scenarios of divide migration and variations in the accumulation pattern. The primary goals are to evaluate the likelihood of divide migration within our domain and identify the appropriate forcings for the thermal and depth–age modeling.

The kinematic model fit to the full depth of measured layers is shown in Figure 10 for two scenarios: (A) using the measured accumulation and vertical velocity forcings for all ages and (B) using the measured forcings for only the past 2 ka (i.e. divide onset at 2 ka, see Section 4.1) and spatially uniform patterns of accumulation and vertical velocity for ages older than 2 ka. The RMSE, calculated as fraction of layer depth offset, is improved by a factor of ~2 when the modern patterns are used for only the past the 2 ka. The resulting misfit shows a pattern of the modeled layers being too deep to the north of the divide, and too shallow to south.

Fig. 10. (a, b) Bed (brown), measured internal layers (black) and modeled internal layering (blue) using the modern spatial patterns of vertical velocity and accumulation rate for all ages (a) and only the past 2 ka (b, uniform vertical velocity and accumulation rate before 2 ka). (c, d) Relative misfit (measured depth – modeled depth)/measured depth. Panels (a) and (b) show only 13 of the 22 measured layers for clarity while panels (c) and (d) show all 22 layers.

There are two primary changes in the model forcing that will yield deeper layers to the south, shallower layers to the north and thereby improve the fit: (1) move the region of divide flow to the north and (2) reverse the accumulation gradient such that there is more accumulation to the south and less to the north.

We performed model runs with the divide shifted north to 6 km along the flowline for different ages prior to 2 ka. The divide then migrates from 6 to 4.7 km (the modern divide position) between 2.5 and 2 ka. The best fit occurs with a divide onset at 6 km at 6 ka. In Figure 11, we show the misfit for divide onset at 6 ka as well as for 9 ka; the misfit for the base scenario (Fig. 10b) of divide onset at 2 ka at 4.7 km along the flowline is also shown. In the scenarios with divide onset prior to 2 ka, the modeled layers are shallower to the north because of the increased layer thinning in the upper portion of the ice sheet. The modeled layers are also deeper to the south, better matching the measured layers, because the average depth of the modeled layer is constrained to match the average depth of the measured layer. Divide onset at 6 ka results in better fits from 5 to 7 km and from 0 to 4 km along the flowline throughout most of the depth; however, the misfit increases from 4 to 5 km, the modern divide region. The scenario with divide onset at 9 ka further increases the misfit beneath the modern divide, which drives an overall increase in the misfit for the full flowline. The timing of the divide onset should be interpreted with caution because the layers are not dated and the improvement in fit is modest. These tests with the kinematic model indicate that the divide is unlikely to have been established prior to the mid-Holocene.

Fig. 11. Relative misfit as calculated in Figures 10e, f for kinematic model scenarios with no divide onset prior to 2 ka (a), divide onset at 6 ka at 6 km (b) and divide onset at 9 ka at 6 km (c).

To assess the impact of the spatial pattern of accumulation, we performed two tests to compare with the base scenario in Figure 10b. In both scenarios, the modern pattern of vertical velocities began at 2 ka with uniform flank-flow for older ages. The first scenario kept the modern accumulation gradient fixed for all ages; the second scenario flipped the modern accumulation gradient for ages older than 2 ka (Fig. S7). Using the modern accumulation gradient for all ages increases the relative misfit from 2.1 to 3.2%, while using a flipped accumulation gradient increases the relative misfit from 2.1 to 2.2%. As indicated above, a reversal in the accumulation gradient results in deeper layers to the south and shallower layers to the north, yielding a similar misfit as the base scenario. Different accumulation gradients and/or timing of changes in the past accumulation gradients can yield better fits. This issue of non-uniqueness is inherent in inferring accumulation rates from deep internal layers near ice divides (e.g. Koutnik and others, 2016), and is more challenging when the layers have not been dated. Scenarios with an onset of the flipped accumulation gradient in the Holocene, and uniform accumulation for older ages, produce a worse fit than when the flipped gradient is applied for all ages. We have no a priori expectation that a reversal in the gradient should have occurred during the early-to-mid-Holocene when West Antarctic climate was relatively stable. While past changes in the accumulation gradient cannot be well constrained, it is unlikely that the modern pattern has persisted for more than a few thousand years.

Overall, the kinematic modeling shows that the divide position was unlikely to be stable before the mid-Holocene and that the pattern of accumulation is difficult to constrain. A key conclusion for the depth–age relationship is that because the divide is unlikely to have been stable, a flank-flow vertical velocity is most appropriate for modeling the depth of the older ice.

4.4 Temperature modeling and basal thermal state

We use a 1-D ice-and-heat flow model to estimate the internal and basal temperature at West Dome. The model is described in Fudge and others (2019) and Buizert and others (2021). We take two approaches. First, we estimate the internal temperature based on the 65 mW m⁻² average of three recent remotely sensed estimates of geothermal flux which are tightly clustered: 62 ± 10 mW m⁻² (Shen and others, 2020), 68 ± 25 mW m⁻² (Stal and others, 2021) and 66 ± 11 mW m⁻² (Martos and others, 2017). This provides a prior estimate for the temperature profile for use in planning ice core drilling and analyses. Second, we find the minimum geothermal flux at which melting begins (e.g. Fudge and others, 2019).

The ice-and-heat flow model is forced by time-varying surface temperature, accumulation rate and vertical velocity shape factor histories. We use an ice thickness of 1600 m, which is held constant. The histories are based on the two closest deep ice cores, the WAIS Divide ice core (WDC) and the South Pole ice core (SPICEcore); these histories are extended through the LIG using the EPICA Dome C core (EPICA, 2004; Veres and others, 2013). The temperature histories are calculated using the below equation:

(9)$$T( t ) _{{\rm Model}} = T( 0 ) _{{\rm Herc}} + ( {T{( t ) }_{{\rm core}}-T{( 0 ) }_{{\rm core}}} ) $$

where T(t)_Model is the surface temperature forcing for the thermo-mechanical model, T(0)_Herc is the modern temperature at Hercules Dome and is −41°C based on a 10 m firn temperature measurement in 2019, T(t)_core is the inferred temperature history from the one of the ice cores, T(0)_core is the modern surface temperature of the ice core site and t is the model time. The Last Glacial Maximum (LGM)-modern change is −11°C for WDC (Cuffey and others, 2016) and −6°C for SPICEcore (Kahle and others, 2021). The accumulation rate histories are calculated as

(10)$$\dot{b}( t ) _{{\rm model}} = \dot{b}( 0 ) _{{\rm Herc}} \times \displaystyle{{\dot{b}{( t ) }_{{\rm core}}} \over {\dot{b}{( 0 ) }_{{\rm core}}}}$$

Using the same notations as for the surface temperature in Eqn (9). The modern accumulation rate at West Hercules Dome, $\dot{b}( 0 ) _{{\rm Herc}}$, is 0.13 m a⁻¹. Both WDC (Fudge and others, 2016) and SPICEcore (Kahle and others, 2021) have inferred LGM accumulation rates of ~40% of modern.

The englacial temperature profiles are shown in Figure 12. The geothermal flux necessary for melting at present is 85 ± 7 mW m⁻² using the SPICEcore scaling and 89 ± 7 mW m⁻² using the WDC scaling. The uncertainty was calculated assuming an end-member modern surface temperature of −39°C (2°C warmer) and modern accumulation rate of 0.11 m a⁻¹ (0.02 m a⁻¹ lower). The minimum geothermal flux to produce melting is considerably larger than the remotely sensed estimates inferred by three different studies. Our estimated range of the geothermal flux necessary for basal melting overlaps with the uncertainty of only one of these estimates, Stal and others (2021), and only because of its large (25 mW m⁻²) uncertainty.

Fig. 12. Temperature profiles for West Dome. Solid lines use the average of the remotely sensed geothermal flux estimates (65 mW m⁻²); dashed lines assume the bed is at the pressure melting point. The surface temperature and accumulation rate histories are scaled from values at SPICEcore (blue) and WAIS Divide (red) using modern values of −41°C and 0.13 m a⁻¹.

4.5 Depth–age relationship

Because the measured internal layers cannot be dated, preliminary depth–age relationships must be assessed from modeling. We use the onset of the divide at 2 ka (Sections 4.1 and 4.3) and three different choices for the vertical velocity shape factor (4, 7 and 10) for ages prior to 2 ka. We use the two accumulation histories described in Section 4.4 which are based on WDC and SPICEcore. For the kinematic model forcing, the fractional accumulation is multiplied by the modern spatial accumulation-rate pattern.

The resulting depth–age relationships at 5.5 km along the x135 flowline are shown in Figure 13. The choice of vertical-velocity shape factor has a larger influence than the choice of accumulation history. Ice of 132 ka age, the start of the LIG, ranges from 47 to 87 m above the bed. The average LIG annual-layer thickness is ~1 mm.

Fig. 13. Depth–age relationships at 5.5 km along the flowline. Accumulation forcing is scaled to SPC (SPICEcore, solid) or WDC (WAIS Divide, dashed). Three different vertical velocity profiles for flank flow (p = 4, 7 or 10) are shown for each accumulation forcing. The height above the bed for ice of 132 ka age ranges from 47 to 87 m.

5. Discussion

5.1 Potential for a deep ice core

One goal of drilling a deep ice core at Hercules Dome is to recover ice that includes the onset of the LIG (132 ka), and preferable Termination II (the transition into the LIG from the previous glacial maximum) as well. To preserve ice of this age, there must be limited or zero basal melt and no more than minor disturbances to the stratigraphy. The current pattern of vertical velocity across the divide at West Dome suggests the ice is frozen to the bed. This is supported by the ice-and-heat flow modeling (Section 4.4). With a geothermal flux of 65 mW m⁻², the ice at the bed is ~−10°C, which would allow subglacial bedrock drilling after completion of the ice core.

The presence of a frozen bed today also implies no significant impact from basal melt at older ages, even if the bed is near the pressure melting point. In the model runs where we found the maximum geothermal flux before basal melt occurred, the bed remains frozen at all ages with the WAIS Divide forcing. For the SPICEcore scenario, there is a ~3000 year period of basal melt that occurs at the end of the LGM (i.e. from ~25 to 22 ka). The maximum melt rate is 0.04 mm a⁻¹, and the total predicted melt is 80 mm. While these two scenarios do not encompass the full range of possibilities, they indicate that significant basal melt is unlikely to impact the age of the basal ice.

The impacts of changes in ice thickness are difficult to evaluate because it is not clear if Hercules Dome was thicker or thinner at the LGM. A thicker ice sheet at the LGM would cause a warmer basal temperature and likely increased basal melt, while a thinner ice sheet would have the opposite effect. The modest amount of basal melt would not remove ice from the LIG. In fact, basal melt rates up to ~0.8 mm a⁻¹ would preserve thicker packets of ice from the LIG, although the ice would be closer to the bed. While the evidence for a frozen bed and limited basal melt may not be as definitive as a measured basal temperature, our study suggests that basal melt is unlikely to have removed ice older than 132 ka from the West Dome site.

The preservation of stratigraphic order is more difficult to assess. The layering at West Dome is resolvable until just a few tens of meters above the bed, but the deepest layer that is traceable along the full radar profile is ~400 m above the bed. Ice from the start of the LIG (132 ka) is likely between 47 and 87 m above the bed, and thus continuity of the deep stratigraphy is important for recovering ice of this age. The deep internal reflectors are coherent all the way to the ice bottom just north of the divide, with more bed conformal layering on the north side of the divide generally, suggesting the north is the preferable location for a core.

The internal layers suggest that the current divide position has not been stable beyond the mid-Holocene; however, the internal layers alone are not enough to constrain the extent of divide migration. Modest variations in the divide position of a few kilometers are enough to prevent the development of a clear Raymond arch above an irregular bed. The lack of a clear Raymond arch at West Dome is not surprising because the characteristic time of ~13 ka, ($H/\dot{b}$), is substantially larger than for the thinner and higher accumulation coastal domes where Raymond arches are widely observed (Matsuoka and others, 2015).

Our best-fit models indicate that average annual layer thicknesses are ~1 mm during the LIG. With continuous-flow analysis methods that are able to make measurements with ~0.5 cm resolution (e.g. Röthlisberger and others, 2000; Osterberg and others, 2006; Gkinis and others, 2010; Sigl and others, 2016; Jones and others, 2017), measurements of the chemistry and isotopes during the interglacial will have approximately decadal resolution. For gas sampling with a typical sample length of ~10 cm, this translates to a maximum of centennial resolution. The largest challenge from thin annual layering will not be resolution but the limited ice availability for ice of older ages. One approach to overcome the challenge of a thinner packet of LIG ice is to drill one or multiple replicate cores with a deviation from the main borehole.

One limitation of the shallower ice at West Dome is that borehole thermometry (e.g. Cuffey and others, 2016) to constrain the magnitude of the glacial–interglacial temperature change becomes increasingly difficult. However, new techniques such as water isotope diffusion lengths, and firn thickness and gas-age ice-age difference (Δage) constraints, have been pioneered to provide alternate ways of addressing this question (e.g. Gkinis and others, 2014; Holme and others, 2018; Kahle and others, 2021; Buizert and others, 2021).

The geophysical survey at Hercules Dome indicates that a continuous climate record through the LIG could likely be obtained at West Dome. West Dome is the highest elevation portion of Hercules Dome, and is the most pronounced dome with the simplest ice-flow setting. The ice thickness varies relatively little at West Dome, with a mean depth of ~1600 m, which minimizes the risk of disrupted layering from complicated ice flow. The ice temperature at the bed is likely below freezing, preserving old ice near the bed.

5.2 Future plans

A full-field season focused on the West Dome region is planned for 22/23. Future measurements will include more extensive HF and VHF surveys. Additional ApRES measurements will be made to both more fully characterize the englacial vertical velocities and to constrain ice fabric variations. These future measurements will improve both the spatial resolution and coverage of the data to more precisely define the best ice core site within the West Dome region. Other locations at Hercules Dome are less promising for a deep ice core, and are not prioritized for future research. At East Dome, where the US-ITASE traverse visited (Jacobel and others, 2005), the basal topography and ice flow appears prohibitively complex. South Ridge has moderately thicker ice but the ice flow is from East Dome and over complex bedrock topography.

6. Conclusions

The surface elevation of Hercules Dome is not defined by a single topographic high, but rather three features with local divides. We have focused on West Dome, which is at the highest elevation in the Hercules Dome region and has the simplest ice-flow regime. Coherent internal layers are imaged to within tens of meters of the bed. The ice is of moderate thickness (~1600 m) and the bed is likely frozen; a geothermal flux of ~20 mW m⁻² greater than suggested by regional estimates would be required to produce melting. The internal layer pattern cannot be well fit with a kinematic ice-flow model if modern accumulation and vertical velocity patterns are held constant in the past; the lack of fit indicates that the modern divide position has likely been established during the Holocene. As with other inland divides in Greenland and Antarctica, modest variations in divide position would prevent the development of a prominent Raymond arch. Our best-fit modeled depth–age relationship suggests ice from the LIG is between 45 and 90 m above the bed, and that the average annual layer thickness for that time period is ~1 mm. The modeled depth–age relationship is insensitive to the Holocene flow history because the age of the deep ice depends primarily on the average shape of the vertical velocity profile for the Last Glacial–Interglacial cycle. The lack of modern basal melt implies no significant melt in the past and therefore ice from the LIG should not have been removed. In conclusion, West Dome is a promising location for a deep ice core that will preserve a climate record sensitive to the size of WAIS during the LIG.