Where is the 1-million-year-old ice at Dome A?

Ice fabric influences the rheology of ice, and hence the age–depth profile at ice core drilling sites. To investigate the age–depth profile to be expected of the ongoing deep ice coring at Kunlun station, Dome A, we use the depth-varying anisotropic fabric suggested by the recent polarimetric measurements around Dome A along with prescribed fabrics ranging from isotropic through girdle to single maximum in a three-dimensional, thermo-mechanically coupled full-Stokes model of a 70 × 70 km² domain around Kunlun station. This model allows for the simulation of the near basal ice temperature and age, and ice flow around the location of the Chinese deep ice coring site. Ice fabrics and geothermal heat flux strongly affect the vertical advection and basal temperature which consequently control the age profile. Constraining modeled age–depth profiles with dated radar isochrones to 2∕3 ice depth, the surface vertical velocity, and also the spatial variability of a radar isochrones dated to 153.3 ka BP, limits the age of the deep ice at Kunlun to between 649 and 831 ka, a much smaller range than previously inferred. The simple interpretation of the polarimetric radar fabric data that we use produces best fits with a geothermal heat flux of 55 mW m⁻². A heat flux of 50 mW m⁻² is too low to fit the deeper radar layers, and 60 mW m⁻² leads to unrealistic surface velocities. The modeled basal temperature at Kunlun reaches the pressure melting point with a basal melting rate of 2.2–2.7 mm a⁻¹. Using the spatial distribution of basal temperatures and the best fit fabric suggests that within 400 m of Kunlun station, 1-million-year-old ice may be found 200 m above the bed, and that there are large regions where even older ice is well above the bedrock within 5–6 km of the Kunlun station.

1 Introduction

Finding a continuous and undisturbed million-year-old ice core record in the Antarctic has been identified by the International Partnership for Ice Core Sciences (IPICS) as one of the most important scientific challenges in ice core research in the near future (http://www.pages.unibe.ch/ini/end-aff/ipics/white-papers, last access: 4 May 20018). This is because the last eight glacial cycles are characterized by irregular cycles of roughly 100 ka in length. Both climate and greenhouse gases co-vary closely. Between 900 ka and 1.2 Ma BP glacial cycles are more regular and are paced at significantly higher frequencies (Lisiecki and Raymo, 2005). The relationship between greenhouse gases and ice sheet growth and decay during these times is presently unknown since it can only be derived from the atmospheric record archived in an Antarctic ice core covering this time interval.

The search for a continuous and undisturbed stratigraphic record containing 1 million-year old ice has also interested and challenged the ice modeling community for several decades (e.g., Van Liefferinge and Pattyn, 2013). Potential sites require thick ice, low accumulation rate, and cold (that is frozen) basal conditions. However, thick ice increases basal temperatures and may lead to basal melting. Geothermal heat flux is largely unknown in Antarctica and estimates have relatively large uncertainty (Van Liefferinge and Pattyn, 2013), which in turn is the major uncertainty of determining the basal thermal state. Basal temperature is very sensitive to geothermal heat flux, and potentially locally variable in mountainous terrain (e.g., Van Liefferinge and Pattyn, 2013); localized basal melting and freezing consequently strongly affects vertical velocity and the age profile (e.g., Sun et al., 2014; Parrenin et al., 2017). Ice fabric is also an important factor in determining the speed of vertical advection in the ice sheet which consequently controls both the basal temperature and the age profile. Depth-varying ice fabric will especially influence the age profile of the deeper ice layers where the base is frozen, although the fabric will not strongly change the temperature profile in the ice.

Van Liefferinge and Pattyn (2013) suggested that the oldest ice sites are most likely situated near the divide areas (in some cases, close to existing deep drilling sites, but in areas of smaller ice thickness) and across the Gamburtsev Subglacial Mountains.

Dome A is the top of the East Antarctic ice sheet and above the underlying Gamburtsev Mountains. Being near the center of the East Antarctic, at an altitude of about 4092 m a.s.l. the mean annual temperature (that is the measured temperature 10 m below the surface) at Dome A is −58.5 ^∘C, the lowest mean annual surface temperature on Earth (Hou et al., 2007). Ice flow in this region is very slow measuring less than 0.3 m a⁻¹ (Yang et al., 2014). The average snow accumulation rate is also small, at about 25 mm ice equivalent a⁻¹ over the past several centuries (AD 1260–2004) (Jiang et al., 2012). Therefore, the Dome A region has good potential regarding the recovery of the oldest ice in an ice core (e.g., Xiao et al., 2008).

Kunlun station (80^∘25^′01${}^{\prime \prime }$ S, 77^∘06^′58${}^{\prime \prime }$ E, 4087 m a.s.l.) was located where the ice thickness is maximal in the vicinity of Dome A specifically for deep ice core drilling to acquire high-resolution records approaching 1 million years in length (Cui et al., 2010). However the mountainous terrain of the Gamburtsev Mountains causes basal melting and refreezing in some places (Bell et al., 2011), which may lead to the loss of the oldest ice, and furthermore complicates the stratigraphic record.

Sun et al. (2014) modeled Dome A ice flow, temperature, and age by applying a full-Stokes model to the summit region where detailed surface radar profiles are available – we use the same domain here. As ice fabric information was not available, Sun et al. (2014) used some simple formulations to define an envelope of possible fabric effects: isotropic and prescribed anisotropic ice fabrics that vary the evolution from isotropic to single maximum at 1∕3 or 2∕3 depths. Using these fabrics resulted in basal ages varying by 500 000 years despite age–depth profiles being constrained by dated radar isochrones in the upper one third of the ice sheet. However, Wang et al. (2017) recently presented spatial variations in ice fabric across Dome A obtained from polarimetric radar data in a 30 × 30 km² grid around Kunlun station. Four distinct ice fabric layers were identified and their ages at Kunlun station found by tracing dated internal ice-sheet layering from the Vostok ice core drilling site.

In this study, we utilize the observed ice fabric determined at Kunlun station along with several prescribed alternative anisotropic ice fabrics in a three-dimensional, thermo-mechanically coupled full-Stokes model to simulate age–depth profiles, and improve the results of Sun et al. (2014). We also use the more plentiful recent measurements including dated radar isochrones at Kunlun station to elucidate the stability of the region on glacial timescales, and the localized variability in geothermal heat flux. Our approach contrasts with that recently used to explore possible ancient ice around Dome C (Parrenin et al., 2017) where a 1-D flow model was used in conjunction with extensive radar profiles.

2 Domain, data, and mesh

The surface and topography ice thickness in the Dome A region come from both airborne and ground-based measurements. The Antarctic Gamburtsev Province Project (AGAP) surveyed the region with flight lines 5 km apart orientated in the north–south direction (Bell et al., 2011) with perpendicular lines every 33 km. Ground-based surveys were done by the 21st and 24th Chinese National Antarctic Research Expedition (CHINARE) in a 30 × 30 km² square along lines typically a few km apart, the along-track radar resolution is 125 m (Sun et al., 2014; Cui et al., 2010). A stake network was also established for ice motion using differential GPS receivers, and data collected in 2008 and 2013 (Yang et al., 2014). More recently polarimetric radar observations were also collected on 5 km-spaced ground-based survey grid (Wang et al., 2017) using a 179 MHz radar and orthogonally orientated antennae with 17 m along-track spacing. These survey observations deduced the existence of four to six layers of different ice fabric, and we make use of these data for the ice dynamics simulation.

The Vostok ice core provides absolute dates for radar internal reflections, and since these radar reflections are often continuous over hundreds of kilometers they can provide age–depth profiles over an extensive region of the ice sheet (e.g., Wang et al., 2017). Isochrones in two 150 MHz airborne radar transects collected by the Alfred Wegener Institute were tracked from Vostok to Dome A (Sun et al., 2014; Wang et al., 2017; Fig. 1a) providing 12 dated layers at Kunlun station. We select the 153.3 ka isochrones for detailed analysis making use of polarimetric radar data collected in a set of four triangles centered on Dome A and 160 km in length.

The modeled domain is a 70 × 70 km² square centered at Kunlun station (Fig. 1). The surface and bedrock topographic data in the 70 × 70 km² domain come from AGAP while in the 30 × 30 km² domain we combined the AGAP and CHINARE data. Crossover analysis of radar lines shows 96 % of differences in both surface and bed elevations were less than 150 m. The surface is flat but the bedrock has gradients in excess of 20 % (Fig. 1). The domain was divided into 21 vertical layers with the lower 6 having logarithmic spacing and the bottommost layer representing 0.3125 % of ice thickness. The mesh contains 48 811 elements and 51 940 nodes.

Figure 1 (a) The locations of Dome A (black circle), Kunlun (black +), Vostok (black x), and the 70 × 70 km² study region (black box). The background is surface elevation with 100 m contour intervals. The two radar lines that connect the Vostok drill site to Dome A are shown in red. (b) The 70 × 70 km² finite element mesh in the vicinity of Dome A is projected on a polar stereographic map with standard parallel at 71^∘ S and central meridian at 0^∘ E. The background is bedrock elevation. The boundaries of the inner region and the whole region are shown, with the inner 30 × 30 km² region centered on Kunlun station displaying a 300 m resolution, and the outer region a 3 km resolution. There are 21 terrain-following vertical layers with thinner layers near the base. The bar in the center denotes the drilling site at Kunlun station.

3 Model

3.1 Field equations

We used the open source finite element method package Elmer/Ice (Gagliardini et al., 2013; http://elmerice.elmerfem.org) to solve the complete three-dimensional, thermo-mechanically coupled full-Stokes model across the model domain (Fig. 1). The following equations define the momentum, mass conservation, and temperature of the ice:

$$\begin{array}{}\text{(1)}& {\displaystyle }& {\displaystyle }\mathit{\rho }\mathrm{\nabla }\cdot \mathit{\sigma }+\mathit{\rho }\mathit{g}=\mathrm{0},\text{(2)}& {\displaystyle }& {\displaystyle }\mathrm{\nabla }\cdot \mathit{u}=\mathrm{0},\text{(3)}& {\displaystyle }& {\displaystyle }\mathit{\rho }c\left({\displaystyle \frac{\partial T}{\partial t}}+(\mathit{u}\cdot \mathrm{\nabla })T\right)=\text{div}\left(\mathit{\kappa }\right(T\left)\mathrm{\nabla }T\right)\end{array}$$

Equation (1) is the Stokes equation denoting the balance for linear momentum: the acceleration (inertia force) is negligible; $\mathit{\sigma }=\mathit{\tau }-p\mathbf{I}$ is the Cauchy stress tensor, where τ is the deviatoric stress tensor that has a non-linear constitutive relationship with the strain rate tensor $\dot{\mathit{\epsilon }}=\frac{\mathrm{1}}{\mathrm{2}}(\mathrm{\nabla }\mathit{u}+\mathrm{\nabla }{\mathit{u}}^{\text{T}})$ which will be discussed in more detail below; p is the pressure; and I is the identity matrix. Equation (2) is the incompressibility condition which implies the conservation of mass. Equation (3) is the heat transfer equation which comes from the conservation of energy. u and T denote ice flow velocity and ice temperature, respectively; ρ, c and κ are the respective density, heat capacity, and heat conductivity of ice; and g is acceleration due to gravity. We neglect strain heating of the ice by internal deformation.

The age of the ice, A, at any point in the ice is governed by

$$\begin{array}{}\text{(4)}& {\displaystyle \frac{\partial A}{\partial t}}+\mathit{u}\cdot \mathrm{\nabla }A=\mathrm{1}.\end{array}$$

3.2 Rheology

We use a non-linear anisotropic constitutive relation between the deviatoric stress tensor τ and strain rate tensor $\dot{\mathit{\epsilon }}$ following Gillet-Chaulet et al. (2006) and Martín and Gudmundsson (2012):

$$\begin{array}{ll}\text{(5)}& {\displaystyle }\dot{\mathit{\epsilon }}=& {\displaystyle \frac{\mathrm{1}}{\mathrm{2}{\mathit{\eta }}_{\mathrm{0}}}}\left(\mathit{\beta }\mathit{\tau }+{\mathit{\lambda }}_{\mathrm{1}}{\mathit{a}}^{\left(\mathrm{4}\right)}:\mathit{\tau }\right.{\displaystyle }& {\displaystyle }\left.+{\mathit{\lambda }}_{\mathrm{2}}(\mathit{\tau }\cdot {\mathit{a}}^{\left(\mathrm{2}\right)}+{\mathit{a}}^{\left(\mathrm{2}\right)}\cdot \mathit{\tau })+{\mathit{\lambda }}_{\mathrm{3}}({\mathit{a}}^{\left(\mathrm{2}\right)}:\mathit{\tau })\mathbf{I}\right),\end{array}$$

where a⁽²⁾ and a⁽⁴⁾ are the second- and fourth-order orientation tensors of ice fabric, respectively; I is the identity matrix; the symbols ⋅ and : are the contracted product and the double contracted product, respectively; and the three λ are expressed as

$$\begin{array}{ll}\text{(6)}& {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{1}}=\mathrm{2}\left(\mathit{\beta }{\displaystyle \frac{\mathit{\gamma }+\mathrm{2}}{\mathrm{4}\mathit{\gamma }-\mathrm{1}}}-\mathrm{1}\right),{\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{2}}=\mathrm{1}-\mathit{\beta },\\ {\displaystyle }& {\displaystyle }{\mathit{\lambda }}_{\mathrm{3}}=-{\displaystyle \frac{\mathrm{1}}{\mathrm{3}}}\left({\mathit{\lambda }}_{\mathrm{1}}+\mathrm{2}{\mathit{\lambda }}_{\mathrm{2}}\right).\end{array}$$

The parameter β is the ratio of the shear viscosity parallel to the basal plane to that in the basal plane, and should be significantly smaller than 1 since ice crystals deform primarily by shear in the basal plane. The parameter γ is the ratio of the viscosity in compression or traction along the c axis to that in the basal plane, and is close to 1 (Gillet-Chaulet et al., 2006). η₀ denotes the basal shear viscosity,

$$\begin{array}{}\text{(7)}& {\mathit{\eta }}_{\mathrm{0}}={\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}A(T{)}^{-\frac{\mathrm{1}}{n}}{\left({\displaystyle \frac{\mathrm{1}}{\mathrm{2}}}\text{tr}\left({\dot{\mathit{\epsilon }}}^{\mathrm{2}}\right)\right)}^{\frac{\mathrm{1}-n}{\mathrm{2}n}},\end{array}$$

where n is the power-law exponent and taken as 3, “tr” denotes trace, and A(T) is the rate factor described by the Arrhenius law (Cuffey and Paterson, 2010),

$$\begin{array}{}\text{(8)}& A\left(T\right)=A_{\mathrm{0}}\exp \left(-{\displaystyle \frac{Q}{R{T}_{\text{h}}}}\right).\end{array}$$

Here the coefficient A₀ is the prefactor, which takes 3.985 × 10⁻¹³ Pa⁻³ s⁻¹ at temperatures below −10 ^∘C and 1.916 × 10³ Pa⁻³ s⁻¹ at temperatures between −10 and 0 ^∘C; T_h denotes Kelvin temperature adjusted for melting point depression: $T_{\text{h}}=T+\mathit{\beta }p$ where $\mathit{\beta }=\mathrm{9.8}\times {\mathrm{10}}^{-\mathrm{8}}$ KP${}_{\text{a}}^{-\mathrm{1}}$; Q denotes the activation energy for creep, which takes 60 kJ mol⁻¹ at temperatures below −10 ^∘C, and 139 kJ mol⁻¹ at temperatures between −10 and 0 ^∘C; R= 8.314 J mol⁻¹ K⁻¹ is the gas constant. In Eq. (5), a⁽²⁾ and a⁽⁴⁾ are defined as

$$\begin{array}{ll}\text{(9)}& {\displaystyle }& {\displaystyle }{\mathit{a}}^{\left(\mathrm{2}\right)}=\oint f\left(\mathit{c}\right)\mathit{c}\otimes \mathit{c}d\mathit{c}=<\mathit{c}\otimes \mathit{c}>,{\displaystyle }& {\displaystyle }{\mathit{a}}^{\left(\mathrm{4}\right)}=\oint f\left(\mathit{c}\right)\mathit{c}\otimes \mathit{c}\otimes \mathit{c}\otimes \mathit{c}d\mathit{c}=<\mathit{c}\otimes \mathit{c}\otimes \mathit{c}\otimes \mathit{c}>,\end{array}$$

where f(c) is the normalized orientation distribution function (ODF) of the c axes, c, with $\oint f\left(\mathit{c}\right)d\mathit{c}=\mathrm{1}$; therefore, the sum of the diagonal components of a⁽²⁾ equals 1. In order to reduce the number of variables, we use the invariant-based optimal fitting closure approximation (IBOF) proposed by Chung and Kwon (2002). The components of a⁽⁴⁾ are approximated as functions of those of a⁽²⁾,

$$\begin{array}{ll}\text{(10)}& {\displaystyle }{a}_{ijkl}^{\left(\mathrm{4}\right)}=& {\displaystyle }{\mathit{\beta }}_{\mathrm{1}}\text{Sym}\left({\mathit{\delta }}_{ij}{\mathit{\delta }}_{kl}\right)+{\mathit{\beta }}_{\mathrm{2}}\text{Sym}\left({\mathit{\delta }}_{ij}{a}_{kl}^{\left(\mathrm{2}\right)}\right){\displaystyle }& {\displaystyle }+{\mathit{\beta }}_{\mathrm{3}}\text{Sym}\left(a_{ij}^{\left(\mathrm{2}\right)}{a}_{kl}^{\left(\mathrm{2}\right)}\right)+{\mathit{\beta }}_{\mathrm{4}}\text{Sym}\left({\mathit{\delta }}_{ij}{a}_{km}^{\left(\mathrm{2}\right)}{a}_{ml}^{\left(\mathrm{2}\right)}\right)\\ {\displaystyle }& {\displaystyle }+{\mathit{\beta }}_{\mathrm{5}}\text{Sym}\left(a_{ij}^{\left(\mathrm{2}\right)}{a}_{km}^{\left(\mathrm{2}\right)}{a}_{ml}^{\left(\mathrm{2}\right)}\right)+{\mathit{\beta }}_{\mathrm{6}}\text{Sym}\left(a_{im}^{\left(\mathrm{2}\right)}{a}_{mj}^{\left(\mathrm{2}\right)}{a}_{kn}^{\left(\mathrm{2}\right)}{a}_{nl}^{\left(\mathrm{2}\right)}\right),\end{array}$$

where “Sym” denotes the symmetrical part of its argument and β_i are six functions of the second and third invariants of a⁽²⁾. Further following Chung and Kwon (2002), we assume β_i are polynomials of degree 5 in the second and third invariants of a⁽²⁾ and use the coefficients computed by Gillet-Chaulet et al. (2006) so that the computed a⁽⁴⁾ from Eq. (9) fits the fourth-order orientation tensor given by the ODF from Gagliardini and Meyssonnier (1999).

3.3 Ice fabric

There are several typical types of fabric in the ice sheet: random ice-crystal fabric, perfect single pole (or single maximum) fabric, and vertical girdle fabric. The evolution of the fabric depends on the specific history of stress conditions experienced by the ice as it travels through the ice sheet. The fabric is represented by the three eigenvalues of the orientation tensor (e.g., Martín and Gudmundsson, 2012), a₁₁, a₂₂, and a₃₃. Sun et al. (2014) used three simple fabric distributions, but here we include radar observations of fabric to produce the following four archetypes of fabric in the central 30 × 30 km² domain:

1. Isotropic fabric (random ice-crystal fabric): $a_{\mathrm{11}}=a_{\mathrm{22}}=a_{\mathrm{33}}=\frac{\mathrm{1}}{\mathrm{3}}$.

2. Single maximum fabric (perfect single pole): $a_{\mathrm{11}}=a_{\mathrm{22}}=\mathrm{0}$, a₃₃=1.

3. “Girdle fabric” meaning a smooth linear transition from isotropic at the surface to single maximum at some transition depth, z_s. Sun et al. (2014) used $z_{\text{s}}=\mathrm{1}/\mathrm{3}$ and 2∕3 depth, and thence to the ice base;

4. “Kunlun fabric” meaning using measured ice fabric layer depths at Kunlun station. Wang et al. (2017) defined six layers for the Kunlun fabric. Here we experimented with subsets of layers.

The Wang et al. (2017) lowermost T5 and T6 layers are rather weak and indistinct in most of the survey grid. The layer T1 is relatively flat, T2 is relatively flat over half the region but has large spatial variation in depth over the other half, while the T3 and T4 layers have large spatial variation of depth and are even missing in some locations around their survey grid. Experiments with four layers T1–T4 show a slightly larger vertical velocity at Kunlun station and deeper depths for the 153.3 ka isochrone (Sect. 4.2), in addition to almost the same horizontal flow velocity as with just the top two layers T1 and T2 in our simulations. So here we present simulation results based on a fabric model using just the two upper layers. At Kunlun station T1 is present from the surface to a depth of 807.3 m, corresponding to ages of 0 to 57 ka, and T2 is present from 807.3 to 1226.2 m, corresponding to ages of 57 to 106 ka. The ice is isotropic in T1, and we assume a linear transition from isotropic at the T1 interface with T2 to single maximum at the base of the T2 layer. We then use single maximum for all ice below T2. Wang et al. (2017) do not present unique solutions for the fabric variation in their layers, nor define how the transition from isotropic to single maximum occurs with depth, so the assumptions we make here are perhaps the simplest, but not the only possible interpretations of the fabric data. The three eigenvalues of the orientation tensor for the fabric archetypes we examine are shown in Fig. 2.

Figure 2 Fabric as a function of depth (at Kunlun station) for girdle fabric with $z_{\text{s}}=\mathrm{1}/\mathrm{3}$ (dashed curve) and $z_{\text{s}}=\mathrm{2}/\mathrm{3}$ (dotted curve), and Kunlun fabric (solid curve). The values of a₁₁ (equals a₂₂) are in blue in (a), while a₃₃ is in black in (b). For single maximum $a_{\mathrm{11}}=a_{\mathrm{22}}=\mathrm{0}$, a₃₃=1; and for isotopic ice $a_{\mathrm{11}}=a_{\mathrm{22}}=a_{\mathrm{33}}=\frac{\mathrm{1}}{\mathrm{3}}$.

3.4 Boundary conditions

The ice surface is assumed to be stress-free and changes in atmospheric pressure and wind stress are neglected as follows:

$$\begin{array}{}\text{(11)}& \mathit{\sigma }\cdot \mathit{n}{\mathrm{|}}_{\mathrm{surface}}=\mathrm{0},\end{array}$$

where σ is the Cauchy stress tensor and n the unit normal vector pointing outwards.

The present-day surface temperature is −58.5 ^∘C, while it is likely about 10 ^∘C warmer than that during the Last Glacial Maximum (LGM) over the East Antarctic plateau (Ritz et al., 2001). The viscosity of the ice would change over time and is not a linear function of temperature. Sun et al. (2014) found that none of the simulations using a surface temperature of −68.5 ^∘C matched well with the dated radar isochrones at Kunlun station, and we confirm this finding with the extended set of dated isochrones extending to 2∕3 ice depth. While glacial period temperatures were likely warmer on average than −68.5 ^∘C they were certainly colder than present day. Sun et al. (2014) explain the poor fits for cold surface temperature simulations as being due to the key role of warm interglacials in determining the vertical velocity profile of the ice because of the exponential Arrhenius dependence on temperature of the ice viscosity (Eq. 8), along with much higher accumulation rates during interglacials. We therefore prescribe surface temperature to be the present value of −58.5 ^∘C in this study.

We run the model with a no-slip condition at the bed. We would however, expect that sliding might occur where there is melting at the bottom. However, surface speeds in our study region are very small (a mean speed of ∼ 11 ± 2.5 cm a⁻¹, Yang et al., 2014) and well matched to the model results we show later from ice deformation without basal sliding (Sect. 4.4); hence basal sliding must only be a small fraction of the total velocity, and should not affect the results we show.

For a cold base (temperature below the pressure melting point), a Neumann-type boundary condition is applied for the basal temperature:

$$\begin{array}{}\text{(12)}& \mathit{\kappa }\left(T\right)\mathrm{\nabla }T\cdot \mathit{n}{\mathrm{|}}_{\mathrm{bed}}=G,\end{array}$$

where G denotes the geothermal heat flux. For a warm base, (temperature reaching the pressure melting point), the basal melting rate (i.e., the vertical velocity w) is calculated by

$$\begin{array}{}\text{(13)}& w={\displaystyle \frac{G-\mathit{\kappa }\left(T\right)\mathrm{\nabla }T\cdot \mathit{n}{\mathrm{|}}_{\mathrm{bed}}}{\mathit{\rho }L}},\end{array}$$

where L denotes the latent heat of ice.

Geothermal heat flux is the most significant unknown boundary condition. Van Liefferinge and Pattyn (2013) produce a map of the broad-scale heat flux and its uncertainty based on three different estimates, which gives about 50 ± 25 mW m⁻² in the Dome A region. Experiments by Sun et al. (2014) suggest a reasonable spatial pattern of basal melting can be obtained using geothermal heat fluxes in the range of 50–60 mW m⁻², with values less than about 45 mW m⁻² producing little or no basal melt; this is in apparent conflict with the radar observations of Bell et al. (2011). Here, we make our simulations with either constant 50, 55, or 60 mW m⁻² heat fluxes across the domain.

The age of ice at the surface is set to zero. This is not necessarily trivial given the low accumulation rates and low temperatures at Dome A, but there is no evidence from radar that the area was an ablation region (Siegert et al., 2003) with negative accumulation at any time in the past.

At the model domain sidewalls we use an adiabatic boundary (i.e., vanishing normal component) for heat flux and a hydrostatic pressure condition from the surrounding ice.

Figure 3 The best-fit simulations (solid lines) at Kunlun station using geothermal heat fluxes of 50, 55, and 60 mW m⁻². Modeled age–depth profile (a) and vertical velocity – depth profile (b). The black points denote the dated radar internal reflection horizons tracked from the Vostok ice core site, using a 37 m firn correction (based on the EDML ice core density profile, Ruth et al., 2007) subtracted from the radar depths to convert to the ice-equivalent model scale. The shaded cyan band shows an envelope of acceptable fits to the radar isochrones and age profile with depth, however the dashed (60 mW m⁻²) line likely has surface velocities that are too high.

4 Simulations and results

We ran steady-state simulations with present day climate forcing and fixed geometry. We used three values of geothermal heat flux 50, 55, and 60 mW m⁻², and the four different types of fabrics described in Sect. 3.3. The model equations detailed in Sect. 3 were solved numerically with the Elmer/Ice model. In detail, we first computed an isotropic steady-state solution of the velocity and temperature fields for a linear rheology (power-law exponent n=1 in Eq. 7). Secondly, we used these results as initial conditions for an isotropic fabric steady-state run with n=3. Thirdly, the isotropic results were used as initial conditions for each of the anisotropic fabric steady-state runs. Finally, the age equation was solved and integrated for 1.5 million years by a semi-Lagrangian method (Martín and Gudmundsson, 2012), using the previously obtained steady-state velocity profile.

We initially ran simulations with a geothermal heat flux of 50 mW m⁻², then we used a restart from that thermal condition, for the second set of simulations with a geothermal heat flux of 55 mW m⁻², and finally did the same again using 60 mW m⁻².

4.1 Modeled age at Kunlun station

We define a best fit in age profile by the rms age error of the simulations from the dated radar isochrones. In Fig. 3 we plot these best fit fabrics for each of the three geothermal heat fluxes. In addition to the age error we can also usefully estimate model performance by the surface vertical velocity. At steady state, surface vertical velocity equals surface accumulation rate. The average accumulation during the past 800 ka is 17.7 mm i.e., a⁻¹ using the EPICA Dome C record (Bazin et al., 2013), which is very close to what the three best fit simulations achieve (Fig. 3b; Table 1).

With a geothermal heat flux of 50 mW m⁻², the best fit is a girdle fabric with $z_{\text{s}}=\mathrm{2}/\mathrm{3}$. The modeled age–depth profile is a noticeably poor fit with the deeper radar isochrones although it matches well in the shallow part. With a geothermal heat flux of 55 mW m⁻², the simulation using Kunlun fabric is the best fit; and with 60 mW m⁻², the simulation using a girdle fabric with $z_{\text{s}}=\mathrm{1}/\mathrm{3}$ is best. Furthermore, this 60 mW m⁻² girdle fabric $z_{\text{s}}=\mathrm{1}/\mathrm{3}$ is the best match overall to the measured data and gives a basal age of 687 ka (Table 1).

We want to bracket the possible age–depth profile, and make best use of the polarimetric radar observations of fabric. Therefore we use the simulation with Kunlun fabric and a geothermal heat flux of 55 mW m⁻² as an upper bound of basal age (831 ka). For the lower bound we choose the measured Kunlun fabric with a geothermal heat flux of 60 mW m⁻² because the lower geothermal heat fluxes seem to produce poor fits while this simulation nicely brackets the best fit overall, although the simulated surface vertical velocity is higher than expected (Table 1). Using the abovementioned gives a lower bound on basal age of 649 ka.

Table 1 Modeled and observed isochrone ages and surface velocities at Kunlun.

4.2 Spatial variability of fabric

We examine how the spatial variation in depth of the 153.3 ka radar isochrone along a track centered at Dome A and passing Kunlun station (Fig. 4a) can be simulated with the fixed fabrics that define the best fits in Fig. 3. We define misfit using a robust measure, which involves the median of the absolute difference between the modeled and measured depths. Figure 4b shows that among the three best fit simulations, the 50 mW m⁻² simulation has the largest misfit of 360 m, while the misfit of the other simulations are all less than 180 m, with the best overall fit (93 m) using the lower bound basal age simulation of Kunlun fabric with G=60 mW m⁻².

Wang et al. (2017) observed large spatial variability in the depth of the ice fabric layers T2–T4. As noted in Sect. 3.3, we applied the depth of the top two fabric layers T1–T2 at Kunlun station to the whole region. Wang et al. (2017) show that the depth of the T2 layer is relatively constant in the region to the north of Kunlun station (which includes triangle 4 in Fig. 4a), but is much deeper in the southern region (including triangles 2, and 3). Smaller T2 layer depth would result in slower simulated ice vertical velocity, hence older ice. Therefore, the modeled depths of the 153.3 ka isochrones to the south of Kunlun are underestimated. The underestimation is larger where the ice is thicker (triangle 3) and smaller in triangle 2 where the ice is thinner, see Fig. 4b.

Figure 4 (a) Measurement tracks (black curve) consisting of four triangles, for the depth of the 153.3 ka isochrone layer; the common point of the four triangles is Dome A; the white cross is Kunlun station; the background is bedrock elevation in a 70 × 70 km² region. (b) The measured (black) and modeled (colored) depths of the 153.3 ka isochrone using the simulations shown in Fig. 3. The distance coordinate in (b) starts from Dome A and follows the tracks of triangles 1–4, with Dome A passes marked by vertical black lines, and Kunlun station marked as a vertical magenta line in triangle 3.

4.3 Modeled age at depth in the central region

Using measured Kunlun fabric and geothermal heat fluxes of 55 and 60 mW m⁻², the modeled age at 95 % depth in the central 30 × 30 km² region is shown in Fig. 5. The age dependence on ice depth is such that deep ice that melts has relatively young ages at 95 % depth, as does thin ice. Melting removes old ice at the base, while thin regions have all their very old ice very close to the bed. There are many more locations where the age simulation reaches the 1.5 Ma limit under the 55 mW m⁻² (Fig. 5c) than under the 60 mW m⁻² (Fig. 5d) heat flux reflecting the more widespread basal melting. The maximum age is reached at depths as shallow as 2000 m under both heat fluxes (Fig. 5b), showing that a shrewd (or lucky) choice of location may recover very ancient ice even under the higher heat flux. But there are no locations with the oldest ice at depths above 2600 m with the 60 mW m⁻² heat flux, and above about 2800 m with 55 mW m⁻² heat flux (Fig. 5b). As we discussed earlier in Sect. 4.2, the age in triangle 3 are probably overestimated (Fig. 5c and d). The modeled age inside triangle 4 has the most confidence.

Figure 5 95 % depth in the central 30 × 30 km² model domain (a) and modeled age of the ice at this depth using Kunlun fabric and a geothermal heat flux of 55 (c, blue dots in b) and 60 (d, red dots in b) mW m⁻² and surface temperature of −58.5 ^∘C. The areas with no basal melt are arbitrarily limited to an age of 1.5 Ma. The black cross is Kunlun station. The black line marks the polarimetric radar survey route triangles marked in Fig. 4.

At the Greenland summit drill site, the GRIP ice core contains small (cm-scale) overturned folds 200 m above bedrock (Taylor et al., 1993), at Dome C stratigraphic continuity was lost only 60 m above the bed (Tison et al., 2015). Although the bedrock topography is smoother in central Greenland than around Dome A, ice sheet temperatures are warmer, vertical velocities higher, and the potential of summit migration over glacial cycles probably greater than the Dome A region. The GRIP ice core is in a similar dynamical pure stress (vertical compression-only) regime as Dome A, but it is not a perfect analogy. Dome C may be a better analogy but using a conservative approach we map the age of the ice 200 m above bedrock in Fig. 6. There is ice at least 1-million-year-old ice simulated on the side slopes of the valley below Kunlun station; the closest to Kunlun station being found directly below a point about 380 m away under 55 mW m⁻², and 1 km away under 60 mW m⁻² heat fluxes. However these positions are less reliable parts of the domain than the area to the north of Kunlun in triangle 4, where old ice is about 5 km away.

Figure 6 Modeled age at a height of 200 m above the bedrock using Kunlun fabric and a geothermal heat flux of 55 mW m⁻² (a) and 60 mW m⁻² (b) in the standard grid coordinate system (unit: km, see Fig. 1), and WGS 1984 latitude and longitudes (inclined grid with the South Pole to the lower left). Kunlun station is marked by a black plus sign. The black curve is the 1 Ma age contour. The thick black line marks the polarimetric radar survey route triangles marked in Fig. 4.

Figure 7 (a) Surface topography with contours, and the measured (black arrows; Yang et al., 2014) and modeled surface horizontal velocity (see legend for details) near Kunlun station. Kunlun station is marked by a black circle. The coordinate system is WGS 1984 plotted using Antarctic polar stereographic with standard parallel at 71^∘ S and central meridian at 0^∘ E. The inset box is a close-up of one velocity datum showing the differences between ice fabrics. (b) Modeled surface vertical velocities (unit: cm a⁻¹) using Kunlun fabric, a geothermal heat flux of 55 mW m⁻², and a surface temperature of −58.5 ^∘C. Note the region plotted in (b) is the central 30 × 30 km² area while in (a) it is the larger 70 × 70 km² region.

4.4 Modeled surface velocity comparison with observation

Yang et al. (2014) calculated the surface horizontal velocity field at 12 survey stakes around Dome A using repeated GPS measurements, and found a mean speed of ∼ 11 ± 2.5 cm a⁻¹, with a maximum velocity of 29 ± 1 cm a⁻¹ and a minimum surface velocity of 3.1 ± 2.6 cm a⁻¹. The modeled surface horizontal velocities from the four best-fit simulations are very similar to each other (Fig. 7a), and are very close to the observed in both magnitude and direction. There is less variability between the four different simulated velocities than with the observed velocities. Thus the fabric cannot be usefully determined by the horizontal surface velocity components.

The modeled surface vertical velocity distribution is shown in Fig. 7b and, as discussed earlier (Sect. 4.1), may be compared with local accumulation. Within the central 30 × 30 km² domain almost all surface velocities are within ±50 % of the value at Kunlun station. There are some larger differences near the border of the larger 70 × 70 km² domain, with small parts even displaying upward velocities. This is likely an indication of the model transition zone flow to the surrounding ice sheet rather than a real effect.

Local accumulation is associated with precipitation, small-scale surface topography over the flat interior of the ice sheet, and wind-driven post-depositional processes (e.g., Frezzotti et al., 2005; Ding et al., 2011). Recent and paleo-surface accumulation rates across Dome A have been measured (e.g., Hou et al., 2007; Ding et al., 2011) and show that the Dome A area has the lowest accumulation rate and smallest spatial variability along a transect from the coast to the summit. This is because it is the coldest and highest region, with smooth topography, it is also furthest from the coast, and has the lowest surface wind speeds. The variations in vertical velocity from the model are not prescribed by surface weather, but determined by mass conservation, and hence reflect advection processes in the ice sheet. Any differences from measured accumulation indicate that the ice sheet is out of steady-state balance. As shown in Fig. 3 there are only small differences in vertical velocity for the best fit fabrics for each of the three geothermal heat fluxes we use, although the lower age bound using a 60 mW m⁻² heat flux produces a value at Kunlun which is too large. Hence, although the vertical velocity does not, in practice, constrain the ice fabric, it can help eliminate a geothermal heat flux that is too high.

Figure 8 Basal temperature relative to pressure melting point using Kunlun fabric and a geothermal heat flux of 55  mW m⁻² (a) and 60  mW m⁻² (b), and the hydraulic potential (c). The bedrock areas at pressure-melting point in (a) and (b) are surrounded by a white contour. Kunlun station is marked as white plus sign. The black curve in (a) and (b) shows the outline of the valley defined as the 1500 m contour in bedrock elevation.

4.5 Modeled basal melt and temperature

Basal temperature depends on surface accumulation rate, ice thickness, and basal geothermal heat flux. Since we use fixed geometry, the surface accumulation rate equals surface vertical velocity. As shown in Fig. 7b, the spatial variation of surface vertical velocity is very small in the central 30 × 30 km² region. Therefore, the high temperature area is located along the valley where the ice is thick (Fig. 8). Using Kunlun fabric and a geothermal heat flux of 55 mW m⁻², the basal ice at the Kunlun station drill site is predicted to be at pressure-melting point (Fig. 8a), along with most of the large valley. There is, however, simulated to be cold basal ice within a kilometer from Kunlun station (Fig. 8a). The spatial extent of melting is considerably larger using a geothermal heat flux of 60 mW m⁻² (Fig. 8b), with several of the side valleys now simulated to melt.

Bell et al. (2011) show extensive melt and refreezing features in the Gamburtsev Mountains. Refreezing is driven by ice thickness gradients pushing water up slope to cooler regions where is can refreeze. This is most likely where a bedrock ridge occurs across the general direction of water flow driven by hydraulic potential. No refreezing features were observed within the domain we model here. Surface slopes in the summit region of Dome A are very low (Fig. 7a), so the hydraulic potential (Fig. 8c) of water at the bed is essentially governed by the bed slopes. Calculation of hydraulic potential shows that this is indeed the case and water flow should be along the valley in the vicinity of Kunlun drill site. The oldest ice closest to Kunlun (Fig. 6) is expected to be perpendicular to this flow direction, on the valley walls or in the regions without basal melt in Fig. 8.

5 Uncertainties

Our approach here is relatively sophisticated in terms of ice models presently in use, but there are several limitations that almost certainly mean that details of the simulation will be wrong. We make the key assumption that the ice sheet is in steady state, and the surface geometry is fixed, which means the surface accumulation rates balance the vertical velocity and are also fixed in time. However, the basal thermal condition is sensitive to the ice thickness although other simulations of the whole Antarctic ice sheet suggest that elevation changes at Dome A have been less than 50 m over glacial cycles (Ritz et al., 2001; Saito and Abe-Ouchi, 2010). Transient simulations with varying geometry and surface accumulation rate in the past 800 ka would improve the model result.

We used a spatially constant geothermal heat flux. Although geothermal flux may vary over kilometer scales, it seems unlikely in East Antarctica. For example, Carson et al. (2014) suggest heat flow may vary by a factor of > 150 % over 10–100 km length scales in East Antarctica. Passalacqua et al. (2017) explored variation in heat flux around Dome C using data from radar surveys, and prescribed a uniform geothermal heat flux over 10 km scales. Schroeder et al. (2014) similarly inferred geothermal heat flux variability from radar surveys over Thwaites Glacier in West Antarctica, which is proximal to the Mount Takahe volcano that was active during the Quaternary, finding heat fluxes could double over ranges of about 20 km. We do not expect any recent magmatic activity in the Gamburtsev Mountains, and the situation of Dome C is probably a reasonable analogue. However there is simply no data to constrain the heat flux around Dome A, and hence modeled thermal structure, ice viscosity, and age–depth profile. Liefferinge and Pattyn (2013) explored the uncertainty in existing geothermal heat flux datasets and their effect on basal temperature with a spatial resolution of 5 km. The basal temperature was calculated using the steady-state thermodynamic equation in which ice flow velocity is calculated from the shallow-ice approximation. The mean geothermal heat flux of the three existing datasets at Dome A is about 45 mW m⁻², with root mean square error of about 20 mW m⁻². Their modeled basal temperature at Dome A is about −10 ^∘C corrected for the dependence on pressure with a root mean square error of about 6 ^∘C. Due to the coarse resolution (5 km) used in the whole Antarctic simulations of Liefferinge and Pattyn (2013), the modeled basal temperature does not have obvious spatial variation across the Dome A region at scales of hundreds of kilometers.

The Gamburtsev Mountain is characterized by large spatial variability in bedrock topography, which means that a full-Stokes model which considers all the stress components is better able to capture the ice dynamics than the shallow-ice approximation is (e.g., Zhao et al., 2013). In our study, large variations in basal temperature are simulated using a full-Stokes model run at around 500 m resolution. The basal thermal state is then very sensitive to the geothermal heat flux (Sun et al., 2014), which we explored using 45, 50, 55, and 60 mW m⁻²; these values span the broad range suggested by Liefferinge and Pattyn (2013).

We also use a spatially constant fabric across our entire model domain, with transitions between fabrics at two fixed depths taken from those measured at Kunlun station by Wang et al. (2017). As discussed in Sect. 4.2, this leads to lower confidence in the age of the basal ice in the region south of Kunlun than to the north. This furthermore means that we have more confidence in finding very old ice in the slightly further away northern region of Fig. 6 than to the south of Kunlun.

Our results suggest spatial variability in basal melting, which may introduce basal accretion in places (Bell et al., 2011), although there is no radar evidence of any basal accretion features in the vicinity, the model could be improved by adding basal hydrology. Basal melting may also introduce sliding at the ice/bed interface, which we explicitly excluded in the model, however, comparison with observed horizontal velocities suggests that this is not an issue. Indeed extraction of sliding rates from inverse modeling using observed velocities would be extremely difficult at Dome A given the very low speeds which make satellite interferometry impossible, in addition to the sparse network of GPS locations.

Surface measurements of horizontal velocity do not constrain fabric information in the ice sheet. The influence of fabric is felt in the deeper ice not near the surface. Hence accurate estimates of fabric must rely on observations from the deeper layers, such as radar isochrones, or potentially vertical velocity profiles from phase sensitive radar. These observations together with a flow model allow geothermal heat flux and hence basal temperatures to be estimated over extended regions where assumptions of unchanging heat flux and fabric hold. Testing this hypothesis by tracking the depths of a 150 ka isochrone with the model suggests that fabric and heat flux variations are not very fast on 10 km horizontal scales, but that localized basal melt may complicate this diagnostic method.

The special ice flow conditions at ice divides often leads to the presence of Raymond arches (Raymond, 1983), where older ice is at shallower depth than it is several ice thicknesses away from the divide. These features are visible as uplifted radar internal reflections in profiles across the divide. The strongest Raymond arches are found in high-accumulation coastal domes where the bed is cold and flat and the ice column is closer to isothermal (e.g., Hindmarsh et al., 2011). However, bed topography is complex at Dome A and Raymond arches are not seen in the observed radar profiles. Furthermore, our ice dynamics package, Elmer/Ice, includes all the physics needed to produce the Raymond effect, but we also detect no such feature in transects across the flow divide. We believe this is because the Raymond arch is being obscured by a combination of rugged basal topography and thermal structure. The strong thermal gradient in the ice sheet tends to reduce the Raymond effect: the tendency of the non-Newtonian rheology to produce a stiff layer near the bed where strain rates are low is counteracted by the tendency of warm temperatures to produce softer ice at depth. The viscosity of the basal ice under the dome is softer than the viscosity of the super cold ice near the surface, but it is still much stiffer than the basal ice away from the dome, causing the old ice to be somewhat up-warped under the ridge. Moreover, the high basal melt rates of 2–3 mm a⁻¹ at Kunlun station draw down ice and obscure the Raymond effect.

Very old and deep ice near bedrock is likely to have experienced vertical mixing via various mechanisms: boudinage between layers with different rheology, small scale non-laminar flow, or regelation around any bed irregularities (Taylor et al., 1993). Although in central Greenland mixing was limited to areas closer than 200 m above the bed, mixing may scale with the vertical relief in the area, which would be very large in the case of the Kunlun site if the ice dome location has migrated by 10 km or more over time. However, the coherence of the radar isochrones to at least 2∕3 ice depth from Vostok through the Gamburtsev Mountains to Dome A suggests that vertical mixing to the topographic scale of the mountains has not occurred. Furthermore analysis of the EPICA Dome C ice core revealed continuous stratigraphy to within 60 m of bedrock (Tison et al., 2015), and Parrenin et al. (2017) use that as a basis for locating ice up to 1.5 Ma old in the Dome C region. Comparing our Fig. 6 with the analysis in Parrenin et al. (2017) shows far more locations as having ice at least 1.5 Ma further than 200 m from the bed in the vicinity of Dome A than at Dome C. The nearest such ice to the Concordia station is about 10 km away, compared with 0.5–5 km from Kunlun station.

6 Summary and conclusions

Using the constraints of observed ice fabric from polarimetric radar observations, depths of dated internal isochrones, along with reasonable estimates of surface vertical velocity allows us eliminate both geothermal heat fluxes lower than 50 and higher than 60 mW m⁻² at Kunlun station. The lower heat flux together with observed fabric produces poor fits to dated radar isochrones deeper than half ice depth. The higher heat flux produces a vertical velocity at the surface that is too fast and inconsistent with good fits to measured accumulation rate and to the dated isochrones.

The best fits to the isochrones and surface velocities rather closely constrain the range of basal ages at the Kunlun drilling site to about 650–830 ka, with the upper end more likely than the lower because the lower age bound comes from an unrealistic 60 mW m⁻² heat flux. The spatial variability of age at 95 % ice thickness illustrates the non-linear dependence on ice thickness. Ice that is too deep lacks old ice due to melting, ice that is too thin leads to old ice being too close to the bed to be useful for ice coring.

Reasonable ice core stratigraphy may be preserved to 200 m above bed, as is the case in central Greenland, or 60 m in the case of Dome C, so we determined locations having 1-million-year-old-ice at least 200 m above the bed. Using our favored values for geothermal heat flux and ice fabric, ice this ancient may be found by vertical drilling within 400 m of the present Kunlun drill site; indeed this location would contain much older ice since it seems to be frozen to the bed. However we have more confidence in our simulation of ancient ice about 5–6 km to the north of Kunlun station than the closer sites to the south. Near-basal ice this close to Kunlun may be accessible with a straight forward repositioning of the drilling site (Talalay et al., 2017) rather than the logistics base. Hydraulic potential suggests that the regions of old ice near Kunlun would not contain refrozen melt water from the deeper valleys. Multiple cores from the same borehole may be recovered, which would allow for the sampling of different climate periods in detail as basal melting effectively stretches the relative younger ice. Thus the Kunlun station is well suited to provide the longest continuous stratigraphic record from Antarctica.