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, Reference Welch and Jacobel2003, Reference Welch, Jacobel and Arcone2009; Christianson and others, Reference Christianson, Jacobel, Horgan, Anandakrishnan and Alley2012, Reference Christianson2016). 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, Reference Lilien, Hills, Driscol, Jacobel and Christianson2020). 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, Reference Dixon2013) and a Looyenga and others (Reference Looyenga1965) 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 (Reference Herron and Langway1980) 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, Reference Dixon2013; 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 (Reference Herron and Langway1980) 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, Reference Nicholls2015) 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, Reference Lilien, Hills, Driscol, Jacobel and Christianson2020) and follows the processing flow described by Brennan and others (Reference Brennan, Lok, Nicholls and Corr2014). 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, Reference Nereson and Waddington2002; Waddington and others, Reference Waddington, Neumann, Koutnik, Marshall and Morse2007; Koutnik and others, Reference Koutnik and Waddington2012). 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, Reference Martín, Hindmarsh and Navarro2006). Second, dynamic models have difficulty reproducing divide flow, particularly in areas with variations in bed topography (Goel and others, Reference Goel, Martin and Matsuoka2018). 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 (Reference Lliboutry1979) 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 (Reference Nereson and Waddington2002) 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, Reference Raymond1983) and observations of the width of Raymond arches (e.g. Matsuoka and others, Reference Matsuoka2015). 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.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, Reference Drews2009): 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, Reference Jacobel and Welch2005). 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, Reference Raymond1983) 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, Reference Kingslake2014). 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 (Reference Kingslake2014) 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, Reference Hills2022). 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.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, Reference Pettit, Thorsteinsson, Jacobson and Waddington2007; Martín and others, Reference Martín, Gudmundsson, Pritchard and Gagliardini2009); 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, Reference Conway, Hall, Denton, Gades and Waddington1999), 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, Reference Goel, Martin and Matsuoka2018). 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, Reference Conway, Hall, Denton, Gades and Waddington1999; Martín and others, Reference Martín, Hindmarsh and Navarro2006). 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, Reference Conway, Hall, Denton, Gades and Waddington1999).

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, Reference Koutnik2016), 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 (Reference Fudge, Biyani, Clemens-Sewall and Hawley2019) and Buizert and others (Reference Buizert2021). 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, Reference Shen, Wiens, Lloyd and Nyblade2020), 68 ± 25 mW m⁻² (Stal and others, Reference Stal, Reading, Halpin and Whittaker2021) and 66 ± 11 mW m⁻² (Martos and others, Reference Martos2017). 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, Reference Fudge, Biyani, Clemens-Sewall and Hawley2019).

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, Reference Veres2013). 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, Reference Cuffey2016) and −6°C for SPICEcore (Kahle and others, Reference Kahle2021). 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, Reference Fudge2016) and SPICEcore (Kahle and others, Reference Kahle2021) 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 (Reference Stal, Reading, Halpin and Whittaker2021), 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.