Articles | Volume 12, issue 8
The Cryosphere, 12, 2773–2787, 2018

Special issue: Oldest Ice: finding and interpreting climate proxies in ice...

The Cryosphere, 12, 2773–2787, 2018

Research article 30 Aug 2018

Research article | 30 Aug 2018

Promising Oldest Ice sites in East Antarctica based on thermodynamical modelling

Promising Oldest Ice sites in East Antarctica based on thermodynamical modelling
Brice Van Liefferinge1, Frank Pattyn1, Marie G. P. Cavitte2,3, Nanna B. Karlsson4,a, Duncan A. Young2, Johannes Sutter4,5, and Olaf Eisen4,6 Brice Van Liefferinge et al.
  • 1Laboratoire de Glaciologie, Université libre de Bruxelles, CP 160/03, Avenue F.D. Roosevelt 50, 1050 Brussels, Belgium
  • 2Institute for Geophysics, Jackson School of Geosciences, University of Texas at Austin, Austin, Texas, USA
  • 3Department of Geological Sciences, Jackson School of Geosciences, University of Texas at Austin, Austin, Texas, USA
  • 4Alfred Wegener Institute Helmholtz-Centre for Polar and Marine Research, 27568 Bremerhaven, Germany
  • 5Climate and Environmental Physics, Physics Institute, and Oeschger Centre for Climate Change Research, University of Bern, Bern, Switzerland
  • 6Department of Geosciences, University of Bremen, Bremen, Germany
  • anow at: Geological Survey of Denmark and Greenland (GEUS), Oster Voldgade 10, 1350 Copenhagen, Denmark

Correspondence: Brice Van Liefferinge (


To resolve the mechanisms behind the major climate reorganisation, which occurred between 0.9 and 1.2 Ma, the recovery of a suitable 1.5 million-year-old ice core is fundamental. The quest for an Oldest Ice core requires a number of key boundary conditions, of which the poorly known basal geothermal heat flux (GHF) is lacking. We use a transient thermodynamical 1-D vertical model that solves for the rate of change of temperature in the vertical, with surface temperature and modelled GHF as boundary conditions. For each point on the ice sheet, the model is forced with variations in atmospheric conditions over the last 2 Ma and modelled ice-thickness variations. The process is repeated for a range of GHF values to determine the value of GHF that marks the limit between frozen and melting conditions over the whole ice sheet, taking into account 2 Ma of climate history. These threshold values of GHF are statistically compared to existing GHF data sets. The new probabilistic GHF fields obtained for the ice sheet thus provide the missing boundary conditions in the search for Oldest Ice. High spatial resolution radar data are examined locally in the Dome Fuji and Dome C regions, as these represent the ice core community's primary drilling sites. GHF, bedrock variability, ice thickness and other essential criteria combined highlight a dozen major potential Oldest Ice sites in the vicinity of Dome Fuji and Dome C, where GHF could allow for Oldest Ice.

1 Introduction

The relationship between the variations in atmospheric CO2 and atmospheric temperatures, determined from oxygen isotope records, is increasingly better understood through a wealth of marine and lacustrine records recently recovered (Kawamura et al.2017). However, characterising this relationship on short timescales, with direct sampling of the palaeo-atmosphere, requires a temporal resolution that can only be obtained from ice core records, which currently only go back as far as 800 ka (EPICA community members2004; Parrenin et al.2007). In particular, there is a strong interest in constraining greenhouse gas forcings between 0.9 and 1.2 Ma, a period during which glacial–interglacial periodicity changed from 40 to 100 ka cycles (Lisiecki and Raymo2005; Snyder2016) but without explained natural forcings (Clark et al.2006; Elderfield et al.2012; Imbrie1993). To resolve the mechanisms behind the major climate reorganisation during the MPT, the recovery of suitable 1.5 million-year-old ice core samples is fundamental. This old ice would provide us with unique and crucial insights into air composition as well as the isotopic composition and dust content of the ice throughout the MPT.

In order to retrieve a 1.5 million-year-old ice core in the centre of the Antarctic Ice Sheet (Wolff et al.2005), the base of the ice sheet should not have experienced melting or refreezing processes during this period (Fischer et al.2013; Wolff et al.2005). Furthermore, even in regions where basal melting can be considered to be insignificant, complex processes of mixing or folding due to rough bedrock topography can cause perturbations in ice flow over the bedrock and make accurate dating of the ice difficult or even impossible (Bell et al.2011). These processes have impacted the NEEM ice core analysis in Greenland (NEEM community members2013) as well as the signal of the deeper part of the EPICA Dome C ice core (Tison et al.2015). Moreover, in order to recover an interpretable climate signal, present-day ice surface velocities should remain below a certain threshold (less than 1 to 2 m a−1 for the horizontal surface velocities), so that ice travels as little as possible horizontally. Finally, ice should be as thick as possible in order to preserve a resolvable and thus an interpretable record within the deeper layers. In 2013, Fischer et al. (2013) and Van Liefferinge and Pattyn (2013) evaluated the conditions necessary for retrieving an old ice core record and highlighted candidate sites with potential 1.5 million-year-old ice in Antarctica. They stressed the importance of collecting more and spatially denser ice thickness data over Oldest Ice candidate sites to reduce uncertainties in model results, both in terms of basal ages and temperature conditions.

The major uncertainty in determining basal melting and basal temperature gradients stems from our limited knowledge of the spatial distribution of the geothermal heat flux (GHF) underneath the Antarctic Ice Sheet. As direct measurements are challenging, due to the presence of the thick ice cover, several approaches exist to derive GHF distributions based on limited data (An et al.2015; Fox-Maule et al.2005; Martos et al.2017; Purucker2013; Shapiro and Ritzwoller2004). All methods infer GHF from properties of the crust and the upper mantle and therefore provide average GHF values without accounting for crustal gradients in GHF. Furthermore, from an ice-sheet modelling perspective, it is crucial to know basal temperature gradients at the ice–bedrock interface and not GHF within the crust. Shapiro and Ritzwoller (2004) extrapolated the GHF from a global seismic model of the crust and the upper mantle. Fox-Maule et al. (2005) derived the GHF from satellite magnetic measurements, and Purucker (2013) provided a GHF data set as an update of the latter. More recently, An et al. (2015) analysed the Earth's mantle properties from a new 3-D crustal shear velocity model to calculate crustal temperature and surface GHF. Their GHF values for East Antarctica deviate by ±10 mW m−2 compared to Shapiro and Ritzwoller (2004), which used a similar method. They found very low GHF values, ∼40 mW m−2, in areas close to Dome C, Dome Fuji and Dome Argus, as well as across the Gamburtsev Subglacial Mountains. Their model, however, is invalid for GHF exceeding 90 mW m−2, but these high values concern the young areas of the crust, mainly in West Antarctica and the Transantarctic Mountains. Finally, Martos et al. (2017) provided the first high-resolution heat flux map on a 15 km by 15 km grid derived from the spectral analysis of a continental compilation of airborne magnetic data. Generally low values of GHF are found in East Antarctica with respect to West Antarctica. This data set estimates the GHF across all candidate sites with variations of up to 20 % from other data sets (Dome F, 65±12 mW m−2; Dome C, 58±12 mW m−2; and Dome A, 55±11 mW m−2). The five data sets differ both in absolute values as well as in their spatial distribution of GHF.

On smaller spatial scales, those particularly relevant to the search for Oldest Ice, GHF is constrained by using models based on ice-penetrating radar, on scales ranging from 1 to 100 km2 (Parrenin et al.2017; Passalacqua et al.2017). Spatially localised features such as lakes and deep ice-core drillings have to be taken into account when attempting to constrain GHF. Subglacial lakes have been documented under the ice of the Antarctic Ice Sheet through the collection of radar and seismic data. An ever-increasing number of lakes have been identified and the current count is close to 415 (Smith et al.2009; Wright and Siegert2012; Young et al.2016). However, more lakes are suspected to exist in currently unsurveyed areas. With respect to deep ice core drill sites, only a few drillings have reached the actual ice–bedrock interface, enabling a direct measurement of the GHF (or at least the basal temperature gradient), i.e. Vostok (Parrenin et al.2004; Petit et al.1999), EPICA Dome C (EPICA community members2004; Parrenin et al.2007), Dome Fuji (Fujii et al.2002; Hondoh et al.2002; Watanabe et al.2003), and EPICA Dronning Maud Land (EPICA community members2006; Ruth et al.2007). All drillings revealed a basal temperature close to or at pressure melting point (pmp).

Since the initial efforts to identify areas of 1.5 million-year-old-ice sites (Fischer et al.2013), a lot of progress has been made in predicting such candidate sites through the collection of detailed ice-penetrating radar data (Cavitte et al.2016; Steinhage et al.2013). Models focussing on divide-adjacent areas and using these radar data also add confidence to the probability of detecting Oldest Ice (Parrenin et al.2017). Furthermore, the mechanisms that control the geometry and the ice volume as well as Antarctic Ice Sheet stability are also increasingly better understood (Pollard et al.2015; Shakun et al.2015). Shakun et al. (2015) put forward the strong coupling between ice volume and temperature over climatic cycles from planktonic δ18O records. Pollard et al. (2015) put forward new mechanisms of hydrofracturing and ice cliff failure, producing a rapid retreat of the ice sheet during past warm periods.

Dense ice-penetrating radar data recently collected over Dome Fuji and Dome C have been instrumental, not only to better constrain the most promising candidate Oldest Ice sites, but also to eliminate some of the modelled candidate sites. Processes active at the base of the ice sheet visible from radar data reduce chances of recovering Oldest Ice. This is the case for areas where significant subglacial water networks have been observed, such as seen in the vicinity of Dome Argus (Wolovick et al.2013), or where subglacial lakes or subglacial trenches have been detected (Wright and Siegert2012; Young et al.2016). Since the aim is to avoid melting at the base while preserving a sufficient ice-core resolution close to the basal-ice layers, this poses an additional problem: ice acts as an insulator, and therefore the greater the ice thickness, the warmer the ice at the base. However, thick ice is needed in order to sufficiently resolve the climatic signal at depth. Conversely, where freezing mechanisms (ice flow divergence, ridge-line freezing) or a reduced ice thickness prevent basal ice from melting, it has been shown that the probability for recovering Oldest Ice is greater, such as in the Gamburtsev Mountains (Creyts et al.2014) or around Dome Fuji and Dome C (Passalacqua et al.2017; Young et al.2017). On a more detailed scale, Cavitte et al. (2018) have shown that near Dome C the snow accumulation pattern is rather stable in time, leading to limited variations in surface topography over the last glacial cycles.

Obviously, the selection of candidate sites will be made, building on radar data criteria. However, since our current radar coverage of the ice sheet interior is currently limited to small, localised areas, it is essential to use thermodynamic models to complement these radar data to characterise basal conditions. In addition, models have the advantage of highlighting areas of interest on small and large scales (Passalacqua et al.2017; Pattyn2010; Van Liefferinge and Pattyn2013). In the context of the Oldest Ice initiative, recently collected radar data and modelling advances have highlighted three candidate areas in particular: Dome Fuji, Dome C and the Dome Argus area, even though radar data still need to be refined for the latter (Sun et al.2014; Wolovick et al.2013). Furthermore, logistical issues cannot be ignored when deciding the next deep ice core drilling site.

So far, thermomechanical modelling has been based on steady-state temperature fields (Pattyn2010; Van Liefferinge and Pattyn2013). However, previous interglacials were demonstrated to have had higher surface temperatures than today (Lisiecki and Raymo2005; Snyder2016), which, in combination with thicker ice (Pollard and DeConto2009; Pollard and Deconto2016) could impact basal temperatures and therefore the basal ice record. Given the size of the Antarctic Ice Sheet and the low vertical advection rates in the interior during prolonged glacial periods, steady-state conditions probably overestimate the probability of melting bed conditions. Here, we use a transient 1-D thermodynamical model to determine whether inland basal conditions over the last 1.5 million years remained frozen and to determine in particular the threshold value of GHF (Gpmp) to satisfy these conditions. Our calculated threshold values are then statistically evaluated through a comparison with existing GHF data sets and their uncertainties (An et al.2015; Fox-Maule et al.2005; Martos et al.2017; Purucker2013; Shapiro and Ritzwoller2004). The obtained probability distribution of ice that has remained frozen over the last 1.5 Ma is used to refine areas of potential Oldest Ice, both on a global scale, to re-examine previous distributions of potential Oldest Ice, and on a local scale, to focus specifically on the Dome Fuji and Dome C districts using the new radar coverage (Karlsson et al.2018; Young et al.2017) and previously suggested constraints (Fischer et al.2013).

Figure 1 Top: GHF from (a) to (d), Martos et al. (2017), Purucker (2013), Shapiro and Ritzwoller (2004) and An et al. (2015) data sets. GHF values are centred on the value of 51 mW m−2. This value corresponds to the average of the GHF data sets within our area of interest. The green line outlines areas with surface velocities <2 m a−1 (Pattyn2010). GHF anomalies are limited to the 2700 m a.s.l. surface elevation contour. The blue triangles locate Dome Fuji and Dome C ice cores (from top to bottom). Refer to Fig. 4 for the northing–easting polar stereographic grid and latitude–longitude coordinates. Bottom: (e) surface-temperature (C) reconstruction adapted from Snyder (2016) near Dome C; (f) ice thickness (m) reconstruction from Pollard and DeConto (2009) at Dome C.


2 Thermodynamical modelling

2.1 Steady-state model

Van Liefferinge and Pattyn (2013) analytically determined the minimum geothermal heat flux necessary to reach the pressure melting point at the base of the ice (Hindmarsh1999; Siegert2000). While a positive GHF increases the temperature at the base of the ice, the surface accumulation cools down the ice from the top. It follows that thick ice combined with low accumulation requires a low GHF to avoid melting from occurring at the base. Although accumulation is relatively well constrained, this is not the case for GHF. In addition, available data sets (An et al.2015; Fox-Maule et al.2005; Martos et al.2017; Purucker2013; Shapiro and Ritzwoller2004) have relatively large errors. In the vicinity of divide areas, GHF uncertainty is 55 %–70 % and 40 % for the Purucker (2013) and Fox-Maule et al. (2005) data sets and for the Shapiro and Ritzwoller (2004) data set respectively in our regions of interest. Van Liefferinge and Pattyn (2013) also used the GHF values from available data sets to calculate the basal temperature and highlight areas with basal melting by running a thermomechanical steady-state ice-flow model. The result was a map of mean basal temperatures on an ensemble model of 15 runs.

2.2 Transient model description

In this paper we solve the vertical temperature profile over time, taking into account vertical diffusion and advection,

(1) T t = k ρ i c p 2 T z 2 - w T z ,

where k is the thermal conductivity, ρi is ice density, cp is the heat capacity of ice, T is the ice temperature, and w the vertical ice velocity. Values for these parameters are listed in Table 1. In divide-adjacent areas, horizontal advection and horizontal heat conduction may be safely neglected. In areas with a relatively smooth bed, as in that case, horizontal conduction is much lower than vertical conduction (Hindmarsh2018, 1999) and horizontal advection and horizontal heat conduction can also be safely neglected. Therefore, here the model is limited to the interior slow-moving areas of the Antarctic Ice Sheet. In this perspective, the 2 m a−1 surface velocity contour is highlighted on all the figures (Pattyn2010) and used as the cut-off surface velocity in the search for Oldest Ice (Van Liefferinge and Pattyn2013). This avoids the need for a correction of the climate signal due to upstream ice advection. The basal boundary condition for a cold base bed is given by

(2) T b z = - G k ,

where G is the geothermal heat flux. At the surface, the temperature is defined by a Dirichlet condition, i.e. T=Ts. The vertical velocity considers a profile based on simple shear using Glen's flow law with a flow exponent n=3 (Hindmarsh1999; Pattyn2010),

(3) w ( ζ ) = - a ˙ ζ n + 2 - 1 + ( n + 2 ) ( 1 - ζ ) n + 1 ,

where ζ=(s-z)/h, h is the ice thickness, s is the surface elevation, and a˙ is the surface accumulation rate.

2.3 Model forcing

Atmospheric forcing is applied in a parameterised way, based on the observed fields of surface mass balance (accumulation rate) obtained from the RACMO regional atmospheric climate model over the period 1980–2004, calibrated with observed surface mass balance rates (van de Berg et al.2006; van den Broeke et al.2006) and surface temperature (van den Broeke2008). For a change in background (forcing) temperature ΔT, corresponding fields of precipitation a˙ and atmospheric temperature Ts are defined by Huybrechts et al. (1998) and Pollard and DeConto (2012).


where γ is the atmospheric lapse rate and δT is 10 C (Pollard and DeConto2012). The subscript “obs” refers to the present-day observed value. Any forcing (increase) in background then leads to an overall increase in surface temperature corrected for elevation changes according to the environmental lapse rate γ. Surface elevation changes with time are obtained from changes in ice thickness with time obtained from a model that takes into account isostatic adjustment, given s(t)=b+H(t), where b is the bed elevation and H(t) the time-varying ice thickness, defined by

(6) H ( t ) = H 0 + ( H p ( t ) - H 0 p ) ,

where H0 is the present-day ice thickness from Bedmap2 (Fretwell et al.2013), updated with the local high-resolution available data for Dome Fuji and Dome C (Karlsson et al.2018; Young et al.2017), and Hp is the ice thickness variation in time obtained from ice sheet modelling over the last 2 Ma (Pollard and DeConto2009). Finally, background temperature changes ΔT(t) are taken from the reconstruction of Snyder (2016), discussed in Sect. 4.1 and scaled to Dome C ice-core temperature reconstruction (Parrenin et al.2007) (Fig. 1).

Table 1Model parameters and constants.

Download Print Version | Download XLSX

2.4 Limit values of GHF

The model is applied on a 5 km by 5 km grid size for the whole of the interior Antarctic Ice Sheet and on a 500 m by 500 m grid size for the two detailed analyses of Dome Fuji and Dome C, with 40 layers in the vertical. For each grid point within our Antarctic domain, the temperature profile is forced with changes in ice thickness, surface temperature and surface mass balance for a given GHF value. This is then repeated for a series of GHF values (Fig. 2), varying around a standard value of 51 mW m−2.

We define Gpmp as the threshold value of GHF for which basal melting may occur during the last 1.5 Ma. Gpmp is determined using a second-degree polynomial fit function between GHF and the maximum basal temperature over the period of the last 1.5 Ma of each run as illustrated at the Dome C site in Fig. 3. GHF values that generate basal temperatures at the pressure melting point are not used for the fit function. To constrain the contribution of the ice cover to the basal temperature variation, the model is further run with uncertainties in the ice thickness. The chosen uncertainty corresponds to a 10 % thicker and thinner ice sheet, which in our areas of interest is equivalent to a variation of 450 to 250 m in elevation. These differences are larger than the variations in ice thickness of the Pollard and DeConto (2009) reconstruction (between 50 and 250 m). A thicker ice cover insulates more than a thinner one and prevents heat flow from escaping as quickly to the surface. Our Gpmp calculation indicates a variation of 6 % to 8 % for the threshold GHF due to the variation in ice thickness. For example, our value of the threshold GHF calculated at Dome C is 51.6 mW m−2. With a higher and a thinner ice cover, these values reach 48.1 and 55.9 mW m−2 respectively, representing a variation of 6.6 % in threshold GHF. This calculation highlights the non-negligible role of the ice thickness on Gpmp variations and therefore also shows the reduced impact of uncertainties in the GHF data sets on the calculation of the basal temperature.

2.4.1 Constraints on GHF

The presence of subglacial water, in the form of lakes or even wet sediments, can inform the basal temperatures, as it implies the pressure melting point has been reached. This allows for local constraints on models (Pattyn2010; Van Liefferinge and Pattyn2013). Subglacial lakes are used as in Van Liefferinge and Pattyn (2013) to constrain the GHF data sets. Lakes are considered to be at the pressure melting point, which implies a local GHF value equal to or larger than the Gpmp. In order to calculate the probability of a frozen bed, a Gaussian function is applied to match the GHF data set with the Gpmp at lake positions. The value of GHF corresponding to 0.95 (2σ) of the probability of the cumulative distribution function (CDF) is used at lake locations (Gpmpc). On the margin of the influence area, which is 20 km or, if known, the size of the subglacial lake (Young et al.2017), a threshold value corresponding to the Gpmp is applied. Two cases are possible: the GHF from the data set is higher than the Gpmp or the GHF is lower. For both cases, a correction (Gc) is made as follows (Pattyn2010):

(7) G c ( x , y ) = G + G pmp c - G exp - x 2 + y 2 20 2 ,

where (x,y) is the horizontal distance in kilometres from the respective lake positions. Without this correction, subglacial lake areas would have GHF values corresponding to a frozen bed.

Figure 2 Basal temperature evolution at Dome C over a period of 1.5 Ma calculated for an ensemble of GHF values illustrated by colours. Red colours indicate high GHF values, which induce temperatures close to the pressure melting point. Blue colours show low GHF values which lead to colder basal temperatures. The GHF values on the illustration are restricted between 22 and 52 mW m−2.


2.5 Constraints on Oldest Ice candidate sites

Until now, we have described how to calculate the probability of frozen conditions at the bed. However, the presence of Oldest Ice at depth only allows for a limited range of key ice parameters. Furthermore, we argue that the flatness of the bed should also be taken into consideration as it can affect ice flow and compromise stratigraphic integrity (NEEM community members2013; Tison et al.2015). We introduce this bed topography constraint in the form of the standard deviation of the spatial bedrock topography variability (σb) (Pattyn2017; Pollard et al.2015; Young et al.2017), which we can assimilate to the roughness of the bed, calculated over an area of 5 km by 5 km. We also introduce a criterion on the maximum distance from radar data, as the density of the radar coverage strongly influences the calculated of the roughness of the bed (σb).

A previous Antarctic-wide analysis (Van Liefferinge and Pattyn2013) used a limit of 2 m a−1 for the horizontal surface velocity, an ice thickness larger than 2000 m and cold basal conditions as acceptable ranges for the occurrence of Oldest Ice. We modify this approach by (i) restricting the parameters' range of values, (ii) taking into account the Gpmp instead of the basal temperature, (iii) adding a σb basal roughness threshold value of 20 m, which implies a relatively smooth bed topography for Dome Fuji and Dome C areas over a radial distance of 2500 m, and (iv) including thresholds of 4 and 2 km on the maximum distance from radar data (for Dome Fuji and Dome C). In addition, (v) we use a minimum H value of 2500 m, as we consider that a minimum H value of 2000 m could be inadequate to obtain a sufficient age resolution at the base of the ice column. We also (vi) use a 1 m a−1 threshold for horizontal surface velocities to limit the influence of ice flow. Finally, the choice of a drill site for an Oldest Ice core will have to be within reasonable distance from radar data in order to provide the necessary upstream constraints when reconstructing the ice core's age-depth chronology. This is already taken into account in the constraints listed above.

Figure 3Example model result for a location near Dome C. The polynomial fit (black line) indicates the value of GHF needed to keep the bed frozen (corrected for pressure melting). In this example, 13 values were used for the fit calculation (no H or surface-temperature uncertainty). The blue lines represent calculated GHF values, applying a H uncertainty from the top to the bottom of −10 % and +10 % and the red lines applying surface-temperature variation of −5 and +5C from the top to the bottom.


Figure 4Map of Gpmp, i.e. the maximum GHF needed to keep a frozen base over 1.5 Ma. Colours represent the magnitude of the GHF (mW m−2). The colour scale's central GHF value, in yellow, is 51 mW m−2. The small black triangles locate the subglacial lakes (Smith et al.2009; Wright and Siegert2012). The green line outlines areas with surface velocities <2 m a−1 (Pattyn2010).


Figure 5The red curve is the cumulative distribution function (CDF) based on the mean and standard deviation of Shapiro and Ritzwoller (2004) GHF data set at Dome C. The blue line is the obtained threshold Gpmp of 51.6 mW m−2 and the indicated probability of having a GHF less than that value.


3 Results

3.1 Large-scale GHF probability distributions

Threshold Gpmp values for the interior of the East Antarctic Ice Sheet for any ice slower than 2 m a−1 is displayed in Fig. 4. According to Fig. 4, Gpmp varies between 20 and 100 mW m−2. Two regions can clearly be distinguished on the map, one with lower values (in blue), located between South Pole and Dome C, and one with higher values located between the Gamburtsev Mountains and Dome Fuji (in red). The difference between both regions is ∼10 mW m−2. This means that the Dome Fuji area allows for higher values of GHF to keep the bed frozen, compared to, for instance, the Dome C area. The Gamburtsev Mountains area differs markedly from other regions, with a high Gpmp between 70 and 100 mW m−2 (due to thinner ice cover and lower surface temperatures than at Dome C and Dome Fuji), while the Vostok region presents the lowest Gpmp values.

The resulting map of Gpmp (Fig. 4) is compared to the published data sets of GHF (An et al.2015; Martos et al.2017; Purucker2013; Shapiro and Ritzwoller2004). Given that each data set has uncertainties associated to the GHF values, a normal probability density function (PDF) and a normal cumulative distribution function CDF (Fig. 5) can be constructed based on the mean and standard deviation of those values. In our case, the CDF can be interpreted as the probability that Gpmp equals or exceeds the GHF of data sets. If the Gpmp is lower than the GHF, the probability of having temperate basal conditions is lower.

Our obtained Gpmp values are then matched against the CDF (see blue line on Fig. 5) to calculate the probability that ice remained frozen over the last 1.5 Ma. The process is repeated for each of the data sets and the resulting probability is shown in Fig. 6. The An et al. (2015) data set appears to exhibit very low GHF values in comparison with the other data sets, especially in the dome regions (Fig. 1), which lead to very low probabilities of reaching the pmp at the domes. The probability map is therefore not shown as it is not a major constraint. On a global scale, GHF values are generally higher in the Martos et al. (2017) data set, which results in a overall lower probability of having a frozen bed, which is more coherent with the basal temperature map proposed by Pattyn (2010) and Van Liefferinge and Pattyn (2013). Although the regions with very high or very low Gpmp values highlighted on the Gpmp distribution map stand out on the three maps, they are most pronounced in the Shapiro and Ritzwoller data set. The Dome C region is interesting since its values are close to the 0.5 threshold between temperate and freezing conditions. The Dome Fuji region has a higher probability of being below the pressure melting point with probabilities between 0.3 and 0.4 on both Shapiro and Ritzwoller (2004) and Purucker (2013) data sets. Regarding the Martos et al. (2017) data set, our analysis is more contrasted. The probability of having a non-frozen bed is much higher in the northern part of Dome Fuji region than in the south.

Figure 6Probability that ice reached the pressure melting point over the last 1.5 Ma according to the GHF data sets from Martos et al. (2017), Purucker (2013) and Shapiro and Ritzwoller (2004) from left to right.


3.2 Small-scale GHF probabilities and Oldest Ice: Dome Fuji and Dome C

The Dome Fuji and Dome C regions are analysed with the same model but applied at a significantly higher spatial resolution (Figs. 7 and 8). As expected, results are in line with the previous continental-scale analysis and Dome Fuji generally exhibits higher values for Gpmp compared to the Dome C region.

The subglacial highlands 40 km south-west of Dome C, informally named “Little Dome C”, and to the north of the subglacial Concordia Trench show the lowest probability of being temperate at the base (∼0.2), regardless of the presence of subglacial lakes, with a Gpmp around 56 and 61 mW m−2 respectively.

Two sites also emerge in the vicinity of Dome Fuji with Gpmp values around 66 mW m−2 and a probability to be non-frozen of 0.1. The first site is located to the northwest of the Dome Fuji region (centred on 76 S30 E, northing–easting 1230∕665 km). The other site is located along a topographic feature characterised by a relatively low ice thickness to the southeast of Dome Fuji. As we can see in Figs. 7 and 8, the Dome Fuji and Dome C sites are well constrained by the recently collected new radar surveys (Karlsson et al.2018; Young et al.2017), thereby avoiding interpolation errors in ice thickness measurements and in the Gpmp calculation.

The basal topography (Figs. 7 and 8) shows a smoother bed (lower bed roughness) at Dome Fuji than at Dome C, enhanced by the differences in radar data cover density. In some cases, the steepest slopes are found near bedrock highs (highlighted by a high σb), and ease with distance from the summit. This is most visible in the vicinity of Dome Fuji (Fig. 7). In the Little Dome C region, the lowest slopes (with a σb around 15 m) are located towards the edges of the subglacial highlands and in the troughs, also shown by Young et al. (2017). In both regions, basal topography also displays flat-topped mountains. In the Dome Fuji region, these plateaus also correspond to high and constant Gpmp values as well as lower ice thickness. However, our results are strongly dependent on radar data coverage density, as will be explained in Sect. 5.

Regarding the Oldest Ice candidate sites at Dome Fuji and Dome C, the red and blue areas on the Figs. 7d and 8d locate the most promising sites, with more conservative parameter values comprising thicker ice and slightly lower ice velocities in red than those in blue. As the ice thickness is mostly larger than 2500 m around Dome C, only red areas are shown. We do not show the same base maps for both data sets because ice thickness has the strongest influence on our model results at Dome Fuji, while the bed roughness is more relevant for Dome C due to the difference in spatial radar data resolution (Figs. 7d and 8d). In the vicinity of Dome Fuji, we note three types of promising sites: (1) extended areas such as those centred on northing–easting 1210∕665 km, 1015∕814 km and 1030∕875 km; (2) several smaller sites scattered in the vicinity; and (3) areas enclosing domes such as those centred on northing–easting 1090∕885 km, 990∕848 km or 1150∕830. All three types fulfil all conditions. In the vicinity of Dome C, we also highlight a number of promising sites. These are scattered over the Little Dome C subglacial highlands and along a transect from the Dome C ice core and the northing–easting -860/1315 km (see Fig. 11). In comparison to Dome Fuji, the Dome C sites are less extensive in area, probably due to the denser radar survey coverage.

4 Discussion

Figure 7 (a) Map of Gpmp results for the Dome Fuji area with the 1-D model calculated on a 500 m × 500 m grid, given in a WGS 84 northing–easting coordinate system (km). (b) Probability that ice reached the pressure melting point over the last 1.5 Ma according to Shapiro and Ritzwoller (2004). The small black triangles locate subglacial lakes and the circle locates the Dome Fuji ice-core site. (c) Standard deviation of bedrock variability; (d) ice thickness from Karlsson et al. (2018). In blue are potential locations of Oldest Ice with H>2000 m, σb<20 m, a probability that ice reached the pressure melting point <0.4, surface velocity <2 m a−1, distance to radar lines <4 km. In red are potential locations with H>2500 m and a surface velocity <1 m a−1. The green dashed lines outline the new radar survey (Karlsson et al.2018).


Figure 8 (a) Map of Gpmp results for Dome C from the 1-D model calculated on a 500 m × 500 m grid, given in a WGS 84 northing–easting coordinate system (km). (b) Probability that ice reached the pressure melting point over the last 1.5 Ma according to Shapiro and Ritzwoller (2004). The small black triangles locate subglacial lakes (Young et al.2017) and the circle locates the Dome C ice-core site. (c) Standard deviation of bedrock variability. (d) Basal topography from Young et al. (2017) with in red, potential locations of Oldest Ice for H>2500 m, σb<20 m, a probability that ice reached the pressure melting point at less than 0.4 and from a distance to radar lines of less than 2 km (right). The green dashed lines locate the 2015–2016 radar survey. LDC locates the Little Dome C area as defined by Parrenin et al. (2017) and Cavitte et al. (2018). Figure 11 is a detailed view of the Little Dome C region.


Knowledge of GHF values at the ice–bedrock interface is a crucial boundary condition for ice flow modelling, yet it remains the most difficult parameter to measure in situ. Constraining this parameter is therefore essential. From an ice-sheet modelling perspective, it is more realistic to know the temperature gradient at the ice-bed interface than a specific GHF value at the interface. The thermal gradient inside the bedrock has an impact on heat availability to the ice (Lowrie2007) as does the thermal inertia of the bedrock. Ritz (1987) shows that bedrock temperature will reach equilibrium after thousands of years, on the scale of several climatic cycles, after a change in ice surface temperature. However, knowing the composition, the thickness and the thermal conductivity of the bedrock is also a challenge. At first approximation, we can use a GHF value without taking into account crustal thickness. This simplification is frequently made in glacier and ice sheet models (Huybrechts and de Wolde1999; Huybrechts et al.1996; Pattyn2003; Pollard et al.2005; Ritz et al.1997), particularly in steady state.

At present, the only constraints on basal GHF are provided by remote measurements and modelling approaches. In this work, we quantify the GHF needed to reach the pmp (Gpmp) and therefore do not calculate an absolute value of the GHF. To do so, we provide constraints on Gpmp, the threshold value of GHF leading to basal melting, by taking into account the glacial–interglacial history of the ice sheet over 1.5 Ma. Our results generally agree with those of Parrenin et al. (2017) and Passalacqua et al. (2017). However, because of the difference in spatial scales, a more detailed comparison is beyond the scope of this paper. We will now discuss the influence of the key parameters (surface-temperature variations, δa˙ and δH) on determining Gpmp on locating Oldest Ice candidate sites. For ease of analysis, Fig. 9 summarises their variations. We will demonstrate that the spatial variability of the distribution and the probabilities of Gpmp (Fig. 4) are directly related to these parameters.

4.1 Surface-temperature forcing

Since no high-resolution reconstructions of surface temperature currently exist over glacial–interglacial timescales, we have chosen to use present-day surface temperatures (van den Broeke2008) generally used in models (Pattyn2010; Van Liefferinge and Pattyn2013) scaled by the Snyder (2016) surface-temperature reconstruction. The Snyder (2016) data set is based on a multi-proxy database and modelling, predicting warmer surface temperatures previous to 800 ka than Lisiecki and Raymo (2005). This global surface-temperature data set is controversial as it may overestimate the Earth System Sensitivity to greenhouse gases and hence the global-mean surface temperature (Schmidt et al.2017). In our case, taking into account warmer surface temperatures between 2 Ma and 800 ka represents a conservative boundary condition and therefore increases our confidence in our predictions of Oldest Ice candidate sites. A higher surface temperature will result in a decrease in the advection of cold temperatures towards the base of the ice, therefore decreasing the Gpmp and so reducing the probability of finding Oldest Ice. However, as explained in Sect. 3, the higher the value of the GHF, the higher the attenuation of the surface-temperature variations with depth. We have shown that an error in surface temperature has a lower or equal effect to an error in ice thickness. The use of Snyder (2016) as a surface-temperature boundary condition should not affect our predictions of Oldest Ice candidate sites as the values lie in the error range (Fig. 3). In addition, absolute differences in temperature from one reconstruction to the next (of the order of a few C) are dwarfed by differences between a glacial and an interglacial period (of the order of 14 C, Fig. 9). And finally, taking into account all past climate variations (accumulation and ice thickness variation) reduces the GHF required to reach the pmp and so contributes to a conservative estimate of Oldest Ice candidate sites. The probability maps of frozen bed conditions obtained (Figs. 7b and 8b) refine the Oldest Ice candidate sites first described in Van Liefferinge and Pattyn (2013).

4.2 Limits on Gpmp calculation

Surface temperatures and accumulation rates are spatially relatively homogeneous in our regions of interest (Fig. 9b and d). This is not the case for ice thickness (Fig. 9c). We can clearly see areas where the mean ice thickness over the last 1.5 Myr is relatively high, more than 3500 m (Dome C, Vostok) and other areas where the mean ice thickness is lower (Gamburtsev Mountains), with differences over time of the order of 1000 m. This is also the case for the Lakes District and Dome C areas, where ice thickness variations are large, around 200 m. In our model this thickness variation corresponds to a Gpmp decrease or an increase of 10 mW m−2.

The variations between the higher and the lower accumulation are around 0.03 m a−1 for the whole period. In a very extreme scenario, we can also consider this fluctuation to be the maximum error of our simulation. This gives us Gpmp variation of the orders of +5.5 and −6.4 mW m−2, for an increase and a decrease in accumulation respectively. For changes in surface temperature, the differences in Gpmp are then +7.6 and −7.2 mW m−2. The combination of the two errors indicates variations of +13.1 and −12.7 mW m−2. Whereas H dominates our result for the Gpmp calculation, errors in accumulation and surface temperature can also have a major impact on the Gpmp, of the order of 25 % of the Gpmp value.

Dome Fuji and Dome C are interesting locations to look at in detail as they provide direct measurements of the temperature profile from ice core measurements. It is therefore possible to deduce the GHF at the base under present conditions. Our Gpmp value at Dome Fuji is 57.3 mW m−2, which agrees with values previously calculated by Seddik et al. (2011) and Hondoh et al. (2002) of 60 and 59 mW m−2, by taking into account a small amount of basal melting. The GHF calculated from the temperature gradient from the Dome Fuji deep ice core is 51.5 mW m−2. In comparison, the four GHF data sets show relatively low values (between 40 and 59 mW m−2), except the value of 65 mW m−2 provided by Martos et al. (2017). In the Dome Fuji region, the bed has a high probability of being close to the pmp, according to the four data sets, but the calculated value of the basal temperature gradient at Dome Fuji points to conditions that are likely non-frozen. This comparison increases the confidence in our approach, as does the comparison with the steady-state approach, described in the following section.

Figure 9 Palaeo-reconstructions for the Antarctic Ice Sheet over the last 1.5 Ma. (a) Ice thickness variation from Pollard and DeConto (2009). The colour scale is truncated at 300 m. (b) Surface mass balance changes (m a−1) related to surface-temperature variations (Pollard and DeConto2012). (c) Mean ice thickness (m) from Pollard and DeConto (2009). (d) Reconstructed variation in the amplitude of surface temperature (C) forced by the multi-proxy database and modelling of Snyder (2016).


4.3 Steady-state model comparison

The spatial variability of GHF noted in the Gpmp maps derived above is also visible in the results of the Van Liefferinge and Pattyn (2013) simple model and the new recalculation of Karlsson et al. (2018). In both models, we observe a clear spatial variability of the GHF, which mainly reflects the ice thickness of the ice sheet. In the GHF histogram calculated for both models (transient and steady state), we clearly note a difference in the mode of the distributions. The difference is 5 mW m−2, with lower values for the steady-state model. This is also what emerges from Fig. 10, which shows differences in GHF between the Van Liefferinge and Pattyn (2013) simple model and the calculated Gpmp. The major part of the highlighted region (Fig. 10) shows GHF values corresponding to the pressure melting point, which are ∼5 mW m−2 lower for the steady-state model except in the area between Dome C and South Pole, where the difference is in some places slightly positive or close to zero, explained by a lower a˙ and δH. We attribute the lower values in our previous steady-state model (simple model) to failure in taking into account variations in surface temperature coupled with changes in thickness. A steady-state model will produce an amplitude of basal temperature variations similar to the surface-temperature variations, but a transient model leads to a much smaller amplitude, of the order of 3 C at the base for a surface variation of 14 C. If ice thickness is reduced, advection of cold temperatures at the surface is increased, which in the transient model implies a decreased basal temperature and a higher value of Gpmp.

Figure 10 Comparison of the GHF needed to reach the pressure melting point calculated using the Van Liefferinge and Pattyn (2013) simple model and the Gpmp calculated in this paper (mW m−2). Blue colours (negative values) indicate that we need a higher GHF for the Gpmp to reach the pressure melting point and vice versa for the positive values (in red).


5 Conclusion and implications for Oldest Ice candidate sites

Figure 11 Detail of promising sites around Little Dome C. Potential locations of Oldest Ice are shown in red for H>2500 m, σb<20 m, a probability that ice has reached the pressure melting point at less than 0.4 km and from a distance to radar lines of less than 2 km. In pink are potential locations of Oldest Ice with the same parameter values but a value of σb<30 m. The background displays basal topography from Young et al. (2017). The blue dashed line locates the elongated bedrock feature discussed in the paper. The green dashed lines locate the 2015–2016 radar survey. LDC and the dashed rectangle locate the Little Dome C area as defined by Parrenin et al. (2017) and Cavitte et al. (2018).


Although there is a large number of parameters that can influence the presence or absence of Oldest Ice at depth, our modelling approach identifies and constrains the key parameters. The obtained Gpmp probability maps have a strong dependency on the spatial resolution of the input data sets, namely the horizontal resolution of the GHF data sets and the horizontal spacing of the radar surveys used. Additionally, a direct comparison of Dome Fuji and Dome C is precluded by the difference in spatial resolution of their respective radar data sets. We note that the bed roughness (σb) is lowest in regions where radar line spacing is the widest, clearly visible on the marginal regions of the σb maps (Figs. 7c and 8c). To take into account the influence of the resolution of the radar surveys, we restrict ourselves to a maximum radial distance from any radar line of 4 km for the Dome Fuji region and 2 km for the Dome C region, so as to take into account the horizontal resolution of their respective radar surveys and reduce the uncertainty in bed topography roughness.

5.1 Dome Fuji

Dome Fuji shows lower σb values on average due to the low density of the radar coverage. The region shows high Gpmp values combined with a thin ice cover. Therefore, the spatial variations in ice thickness, H, dominate the distribution of Oldest Ice for this region. The most promising Oldest Ice candidate sites in the vicinity of Dome Fuji are located on the edges of subglacial mountains, which have the advantage of offering a thicker H and a lower σb while keeping cold conditions at the base. Plateau areas also show potential Oldest Ice sites, but these are less promising due to their lower age resolution as a result of their thinner ice cover.

5.2 Dome C

Dome C is characterised by higher values of σb on average due to the radar coverage's higher spatial density. We note that our σb distribution is similar to that calculated by Young et al. (2017), which adds confidence in our results.

The bed roughness and the Gpmp probability distributions have the strongest influence on the location of candidate sites for this region. In the vicinity of Dome C, potential candidate sites are found in two areas in particular (Fig. 11): near Little Dome C, where σb values are low and the probability of a frozen bed is high, as well as along a transect from the Dome C ice core and the northing–easting -860/1315. The geometry of this elongated bedrock feature is evocative of a raised fault block (also referred to as a horst in geology), which, if confirmed, implies an uplifted but relatively flat bedrock surface. This is promising because it offers a wider area with appropriate ice thicknesses for the recovery of Oldest Ice.

We conclude that, following the analysis of the recent radar data surveys and our modelling efforts, both regions remain interesting as Oldest Ice drilling sites. This work highlights a number of candidate locations that will benefit from the collection of additional geophysical data and modelling.

Data availability

The data that support the findings in this study are available from the corresponding author. Gpmp data and the location of promising Oldest Ice sites are available in Matlab © and GIS formats.

Author contributions

BVL and FP developed the model and ran the analysis. MGPC and DAY provided the processed radar data for the Dome C region. NBK, JS and OE provided processed radar data for the Dome Fuji region. OE coordinated the BE-OI project. The manuscript was written by BVL, FP and MGPC with relevant input from all co-authors.

Competing interests

The authors declare that they have no conflict of interest.

Special issue statement

This article is part of the special issue “Oldest Ice: finding and interpreting climate proxies in ice older than 700 000 years (TC/CP inter-journal SI)”. It is not associated with a conference.


We would like to first thank the editor Eric Larour and both reviewers for their constructive reviews on the paper. This publication was generated in the frame of Beyond EPICA-Oldest Ice (BE-OI). The project has received funding from the European Union's Horizon 2020 research and innovation programme under grant agreement no. 730258 (BE-OI CSA). It has received funding from the Swiss State Secretariate for Education, Research and Innovation (SERI) under contract number 16.0144. It is furthermore supported by national partners and funding agencies in Belgium, Denmark, France, Germany, Italy, Norway, Sweden, Switzerland, The Netherlands and the United Kingdom. Logistic support is mainly provided by AWI, BAS, ENEA and IPEV. The radar survey at Dome C was supported by NSF grant PLR-1443690, and ADD for logistic. Brice Van Liefferinge received a “Fonds Van Buuren” grant to complete this research. Computational resources have been provided by the Shared ICT Services Centre, Université libre de Bruxelles. We would like to thank David Pollard for the data and the fruitful discussions.

The opinions expressed and arguments employed herein do not necessarily reflect the official views of the European Union funding agency, the Swiss Government or other national funding bodies.

This is UTIG contribution #3261.

This is BE-OI publication number #6 .

Edited by: Eric Larour
Reviewed by: Ralf Greve and one anonymous referee


An, M., Wiens, D. A., Zhao, Y., Feng, M., Nyblade, A., Kanao, M., Li, Y., Maggi, A. and Lév`̂que, J. J.: Temperature, lithosphere-asthenosphere boundary, and heat flux beneath the Antarctic Plate inferred from seismic velocities, J. Geophys. Res.-Sol Ea., 120, 8720–8742,, 2015. a, b, c, d, e, f, g

Bell, R. E., Ferraccioli, F., Creyts, T. T., Braaten, D., Corr, H., Das, I., Damaske, D., Frearson, N., Jordan, T., Rose, K., Studinger, M., and Wolovick, M.: Widespread Persistent Thickening of the East Antarctic Ice Sheet by Freezing from the Base, Science, 331, 1592–1595, 2011. a

Cavitte, M. G., Blankenship, D. D., Young, D. A., Schroeder, D. M., Parrenin, F., Lemeur, E., Macgregor, J. A., and Siegert, M. J.: Deep radiostratigraphy of the East Antarctic plateau: connecting the Dome C and Vostok ice core sites, J. Glaciol., 62, 323–334, 2016. a

Cavitte, M. G. P., Parrenin, F., Ritz, C., Young, D. A., Van Liefferinge, B., Blankenship, D. D., Frezzotti, M., and Roberts, J. L.: Accumulation patterns around Dome C, East Antarctica, in the last 73 kyr, The Cryosphere, 12, 1401–1414,, 2018. a, b, c

Clark, P. U., Archer, D., Pollard, D., Blum, J. D., Rial, J. A., Brovkin, V., Mix, A. C., Pisias, N. G., and Roy, M.: The middle Pleistocene transition: characteristics, mechanisms, and implications for long-term changes in atmospheric pCO2, Quaternary Sci. Rev., 25, 3150–3184,, 2006. a

Creyts, T. T., Ferraccioli, F., Bell, R. E., Wolovick, M., Corr, H., Rose, K. C., Frearson, N., Damaske, D., Jordan, T., Braaten, D., and Finn, C.: Freezing of ridges and water networks preserves the Gamburtsev Subglacial Mountains for millions of years, Geophys. Res. Lett., 41, 8114–8122, 2014. a

Elderfield, H., Ferretti, P., Crowhurst, S., Mccave, I. N., Hodell, D., and Piotrowski, A. M.: Evolution of Ocean Temperature and Ice Volume Through the, Science, 337, 704–710,, 2012. a

EPICA community members: Eight Glacial Cycles from an Antarctic Ice Core, Nature, 429, 623–628, 2004. a, b

EPICA community members: One-to-One Coupling of Glacial Climate Variability in Greenland and Antarctica, Nature, 444, 195–198, 2006. a

Fischer, H., Severinghaus, J., Brook, E., Wolff, E., Albert, M., Alemany, O., Arthern, R., Bentley, C., Blankenship, D., Chappellaz, J., Creyts, T., Dahl-Jensen, D., Dinn, M., Frezzotti, M., Fujita, S., Gallee, H., Hindmarsh, R., Hudspeth, D., Jugie, G., Kawamura, K., Lipenkov, V., Miller, H., Mulvaney, R., Parrenin, F., Pattyn, F., Ritz, C., Schwander, J., Steinhage, D., van Ommen, T., and Wilhelms, F.: Where to find 1.5 million yr old ice for the IPICS “Oldest-Ice” ice core, Clim. Past, 9, 2489–2505,, 2013. a, b, c, d

Fox-Maule, C., Purucker, M. E., Olsen, N., and Mosegaard, K.: Heat Flux Anomalies in Antarctica Revealed by Satellite Magnetic Data, Science, 309, 464–467, 2005. a, b, c, d, e

Fretwell, P., Pritchard, H. D., Vaughan, D. G., Bamber, J. L., Barrand, N. E., Bell, R., Bianchi, C., Bingham, R. G., Blankenship, D. D., Casassa, G., Catania, G., Callens, D., Conway, H., Cook, A. J., Corr, H. F. J., Damaske, D., Damm, V., Ferraccioli, F., Forsberg, R., Fujita, S., Gim, Y., Gogineni, P., Griggs, J. A., Hindmarsh, R. C. A., Holmlund, P., Holt, J. W., Jacobel, R. W., Jenkins, A., Jokat, W., Jordan, T., King, E. C., Kohler, J., Krabill, W., Riger-Kusk, M., Langley, K. A., Leitchenkov, G., Leuschen, C., Luyendyk, B. P., Matsuoka, K., Mouginot, J., Nitsche, F. O., Nogi, Y., Nost, O. A., Popov, S. V., Rignot, E., Rippin, D. M., Rivera, A., Roberts, J., Ross, N., Siegert, M. J., Smith, A. M., Steinhage, D., Studinger, M., Sun, B., Tinto, B. K., Welch, B. C., Wilson, D., Young, D. A., Xiangbin, C., and Zirizzotti, A.: Bedmap2: improved ice bed, surface and thickness datasets for Antarctica, The Cryosphere, 7, 375–393,, 2013. a

Fujii, Y., Azuma, N., Tanaka, Y., Nakayama, M., Kameda, T., Shinbori, K., Katagiri, K., Fujita, S., Takahashi, A., Kawada, K., Motoyama, H., Narita, H., Kamiyama, K., Furukawa, T., Takahashi, S., Shoji, H., Enomoto, H., Sitoh, T., Miyahara, T., Naruse, R., Hondoh, T., Shiraiwa, T., Yokoyama, K., Ageta, Y., Saito, T., and Watanabe, O.: Deep Ice Core Drilling to 2503 m Depth at Dome Fuji, Antarctica, Mem. Natl Inst. Polar Res., Spec. Issue, 56, 103–116, 2002. a

Hindmarsh, R.: Ice-sheet and glacier modelling, in: Past Glacial Environments, Elsevier, 605–661,, 2018. a

Hindmarsh, R. C. A.: On the Numerical Computation of Temperature in an Ice Sheet, J. Glaciol., 45, 568–574, 1999. a, b, c

Hondoh, T., Shoji, H., Watanabe, O., Salamatin, A. N., and Lipenkov, V.: Depth-Age and Temperature Prediction at Dome Fuji Station, East Antarctica, Ann. Glaciol., 35, 384–390, 2002. a, b

Huybrechts, P. and de Wolde, J.: The Dynamic Response of the Greenland and Antarctic Ice Sheets to Multiple-Century Climatic Warming, J. Climate, 12, 2169–2188,<2169:TDROTG>2.0.CO;2, 1999. a

Huybrechts, P., Payne, A., and The EISMINT Intercomparison Group: The EISMINT Benchmarks for Testing Ice–Sheet Models, Ann. Glaciol., 23, 1–12, 1996. a

Huybrechts, P., Abe-Ouchi, A., Marsiat, I., Pattyn, F., Payne, A., Ritz, C., and Rommelaere, V.: Report of the Third EISMINT Workshop on Model Intercomparison, European Science Foundation (Strasbourg), 1998. a

Imbrie, J.: On the Structure and Origin of Major Glaciation Cycles. 2. The 100,000-Year Cycle, Paleoceanography, 8, 699–735, 1993. a

Karlsson, N. B., Binder, T., Eagles, G., Helm, V., Pattyn, F., Van Liefferinge, B., and Eisen, O.: Glaciological characteristics in the Dome Fuji region and new assessment for “Oldest Ice”, The Cryosphere, 12, 2413–2424,, 2018. a, b, c, d, e, f

Kawamura, K., Abe-Ouchi, A., Motoyama, H., Ageta, Y., Aoki, S., Azuma, N., Fujii, Y., Fujita, K., Fujita, S., Fukui, K., Furukawa, T., Furusaki, A., Goto-Azuma, K., Greve, R., Hirabayashi, M., Hondoh, T., Hori, A., Horikawa, S., Horiuchi, K., Igarashi, M., Iizuka, Y., Kameda, T., Kanda, H., Kohno, M., Kuramoto, T., Matsushi, Y., Miyahara, M., Miyake, T., Miyamoto, A., Nagashima, Y., Nakayama, Y., Nakazawa, T., Nakazawa, F., Nishio, F., Obinata, I., Ohgaito, R., Oka, A., Okuno, J., Okuyama, J., Oyabu, I., Parrenin, F., Pattyn, F., Saito, F., Saito, T., Saito, T., Sakurai, T., Sasa, K., Seddik, H., Shibata, Y., Shinbori, K., Suzuki, K., Suzuki, T., Takahashi, A., Takahashi, K., Takahashi, S., Takata, M., Tanaka, Y., Uemura, R., Watanabe, G., Watanabe, O., Yamasaki, T., Yokoyama, K., Yoshimori, M., and Yoshimoto, T.: State dependence of climatic instability over the past 720,000 years from Antarctic ice cores and climate modeling, Science Advances, 3, e1600446,, 2017. a

Lisiecki, L. E. and Raymo, M. E.: A Pliocene-Pleistocene Stack of 57 Globally Distributed Benthic δ18O Records, Paleoceanography, 20, PA1003,, 2005. a, b, c

Lowrie, W.: Fundamentals of geophysics, 2nd Edn., Cambridge university press, 2007. a

Martos, Y. M., Catalan, M., Jordan, T. A., Golynsky, A., Golynsky, D., Eagles, G., and Vaughan, D. G.: Heat flux distribution of Antarctica unveiled, Geophys. Res. Lett., 44, 11417–11426,, 2017. a, b, c, d, e, f, g, h, i, j

NEEM community members: Eemian interglacial reconstructed from a Greenland folded ice core, Nature, 493, 489–494,, 2013. a, b

Parrenin, F., Remy, F., Ritz, C., Siegert, M., and Jouzel, J.: New Modelling of the Vostok Ice Flow Line and Implication for the Glaciological Chronology of the Vostok Ice Core, J. Geophys. Res., 109, D20102,, 2004. a

Parrenin, F., Dreyfus, G., Durand, G., Fujita, S., Gagliardini, O., Gillet, F., Jouzel, J., Kawamura, K., Lhomme, N., Masson-Delmotte, V., Ritz, C., Schwander, J., Shoji, H., Uemura, R., Watanabe, O., and Yoshida, N.: 1-D-ice flow modelling at EPICA Dome C and Dome Fuji, East Antarctica, Clim. Past, 3, 243–259,, 2007. a, b, c

Parrenin, F., Cavitte, M. G. P., Blankenship, D. D., Chappellaz, J., Fischer, H., Gagliardini, O., Masson-Delmotte, V., Passalacqua, O., Ritz, C., Roberts, J., Siegert, M. J., and Young, D. A.: Is there 1.5-million-year-old ice near Dome C, Antarctica?, The Cryosphere, 11, 2427–2437,, 2017. a, b, c, d, e

Passalacqua, O., Ritz, C., Parrenin, F., Urbini, S., and Frezzotti, M.: Geothermal flux and basal melt rate in the Dome C region inferred from radar reflectivity and heat modelling, The Cryosphere, 11, 2231–2246,, 2017. a, b, c, d

Pattyn, F.: A New 3D Higher-Order Thermomechanical Ice-Sheet Model: Basic Sensitivity, Ice-Stream Development and Ice Flow across Subglacial Lakes, J. Geophys. Res., 108, 2382,, 2003. a

Pattyn, F.: Antarctic Subglacial Conditions inferred from a Hybrid Ice Sheet/Ice Stream Model, Earth Planet. Sc. Lett., 295, 451–461, 2010. a, b, c, d, e, f, g, h, i, j

Pattyn, F.: Sea-level response to melting of Antarctic ice shelves on multi-centennial timescales with the fast Elementary Thermomechanical Ice Sheet model (f.ETISh v1.0), The Cryosphere, 11, 1851–1878,, 2017. a

Petit, J. R., Jouzel, J., Raynaud, D., Barkov, N. I., Barnola, J. M., Basile, I., Bender, M., Chappellaz, J., Davis, M., Delaygue, G., Delmotte, M., Kotlyakov, V., Legrand, M., Lipenkov, V. Y., Lorius, C., Pepin, L., Ritz, C., Saltzman, E., and Stievenard, M.: Climate and Atmospheric History of the Past 420,000 Years from the Vostok Ice Core, Antarctica, Nature, 399, 429–436, 1999. a

Pollard, D. and DeConto, R. M.: Modelling West Antarctic Ice Sheet Growth and Collapse Through the Past Five Million Years, Nature, 458, 329–332,, 2009. a, b, c, d, e, f

Pollard, D. and DeConto, R. M.: Description of a hybrid ice sheet-shelf model, and application to Antarctica, Geosci. Model Dev., 5, 1273–1295,, 2012. a, b, c

Pollard, D. and Deconto, R. M.: Contribution of Antarctica to past and future sea-level rise, Nature, 531, 591–597,, 2016. a

Pollard, D., DeConto, R. M., and Nyblade, A. A.: Sensitivity of Cenozoic Antarctic Ice Sheet Variations to Geothermal Heat Flux, Global Planet. Change, 49, 63–74, 2005. a

Pollard, D., Deconto, R. M., and Alley, R. B.: Potential Antarctic Ice Sheet retreat driven by hydrofracturing and ice cliff failure, Earth Planet. Sc. Lett., 412, 112–121,, 2015. a, b, c

Purucker, M.: Geothermal heat flux data set based on low resolution observations collected by the CHAMP satellite between 2000 and 2010, and produced from the MF-6 model following the technique described in Fox Maule et al. (2005), available at:, last access: 23 March 2013. a, b, c, d, e, f, g, h, i

Ritz, C.: Time Dependent Boundary Conditions for Calculation of Temperature Fields in Ice Sheets, IAHS Publ., 170, 207–216, 1987. 

Ritz, C., Fabre, A., and Letréguilly, A.: Sensitivity of a Greenland Ice Sheet Model to Ice Flow and Ablation Parameters: Consequences for the Evolution Through the Last Glacial Cycle, Clim. Dynam., 13, 11–24, 1997. a

Ruth, U., Barnola, J.-M., Beer, J., Bigler, M., Blunier, T., Castellano, E., Fischer, H., Fundel, F., Huybrechts, P., Kaufmann, P., Kipfstuhl, S., Lambrecht, A., Morganti, A., Oerter, H., Parrenin, F., Rybak, O., Severi, M., Udisti, R., Wilhelms, F., and Wolff, E.: “EDML1”: a chronology for the EPICA deep ice core from Dronning Maud Land, Antarctica, over the last 150 000 years, Clim. Past, 3, 475–484,, 2007. a

Schmidt, G. A., Severinghaus, J., Abe-ouchi, A., Alley, R. B., Broecker, W., Brook, E., Etheridge, D., Kawamura, K., Keeling, R. F., Leinen, M., Marvel, K., and Stocker, T. F.: Brief Communications Arising Overestimate of committed warming Snyder replies, Nature Publishing Group, 547, E16–E17,, 2017. a

Seddik, H., Greve, R., Zwinger, T., and Placidi, L.: A full Stokes ice flow model for the vicinity of Dome Fuji, Antarctica, with induced anisotropy and fabric evolution, The Cryosphere, 5, 495–508,, 2011. a

Shakun, J. D., Lea, D. W., Lisiecki, L. E., and Raymo, M. E.: An 800-kyr record of global surface ocean δ18O and implications for ice volume-temperature coupling, Earth Planet. Sc. Lett., 426, 58–68,, 2015. a, b

Shapiro, N. M. and Ritzwoller, M. H.: Inferring Surface Heat Flux Distributions Guided by a Global Seismic Model: Particular Application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224, 2004. a, b, c, d, e, f, g, h, i, j, k, l, m, n

Siegert, M. J.: Antarctic Subglacial Lakes, Earth-Sci. Rev., 50, 29–50, 2000. a

Smith, B. E., Fricker, H. A., Joughin, I. R., and Tulaczyk, S.: An Inventory of Active Subglacial Lakes in Antarctica Detected by ICESat (2003–2008), J. Glaciol., 55, 573–595, 2009. a, b

Snyder, C. W.: Evolution of global temperature over the past two million years, Nature, 538, 226–228,, 2016. a, b, c, d, e, f, g, h

Steinhage, D., Kipfstuhl, S., Nixdorf, U., and Miller, H.: Internal structure of the ice sheet between Kohnen station and Dome Fuji, Antarctica, revealed by airborne radio-echo sounding, Ann. Glaciol., 54, 163–167,, 2013. a

Sun, B., Moore, J. C., Zwinger, T., Zhao, L., Steinhage, D., Tang, X., Zhang, D., Cui, X., and Martín, C.: How old is the ice beneath Dome A, Antarctica?, The Cryosphere, 8, 1121–1128,, 2014. a

Tison, J.-L., de Angelis, M., Littot, G., Wolff, E., Fischer, H., Hansson, M., Bigler, M., Udisti, R., Wegner, A., Jouzel, J., Stenni, B., Johnsen, S., Masson-Delmotte, V., Landais, A., Lipenkov, V., Loulergue, L., Barnola, J.-M., Petit, J.-R., Delmonte, B., Dreyfus, G., Dahl-Jensen, D., Durand, G., Bereiter, B., Schilt, A., Spahni, R., Pol, K., Lorrain, R., Souchez, R., and Samyn, D.: Retrieving the paleoclimatic signal from the deeper part of the EPICA Dome C ice core, The Cryosphere, 9, 1633–1648,, 2015. a, b

van de Berg, W. J., van den Broeke, M. R., Reijmer, C. H., and van Meijgaard, E.: Reassessment of the Antarctic Surface Mass Balance using Calibrated Output of a Regional Atmospheric Climate Model, J. Geophys. Res, 111, D11104,, 2006.  a

van den Broeke, M. R.: Depth and Density of the Antarctic Firn Layer, Arct. Antarct. Alp. Res., 40, 432–438, 2008. a, b

van den Broeke, M. R., van de Berg, W. J., and van Meijgaard, E.: Snowfall in Coastal West Antarctica much greater than previously assumed, Geophys. Res. Lett., 33, L02505,, 2006. a

Van Liefferinge, B. and Pattyn, F.: Using ice-flow models to evaluate potential sites of million year-old ice in Antarctica, Clim. Past, 9, 2335–2345,, 2013. a, b, c, d, e, f, g, h, i, j, k, l, m, n, o

Watanabe, O., Jouzel, J., Johnsen, S., Parrenin, F., Shoji, H., and Yoshida, N.: Homogeneous Climate Variability across East Antarctica over the Past Three Glacial Cycles, Nature, 422, 509–512, 2003. a

Wolff, E., Brook, E., Dahl-Jensen, D., Fujii, Y., Jouzel, J., Lipenkov, V., and Severinghaus, J.: The oldest ice core: A 1.5 million year record of climate and greenhouse gases from Antarctica, available at: (last access: November 2017), 2005. a, b

Wolovick, M. J., Bell, R. E., Creyts, T. T., and Frearson, N.: Identification and control of subglacial water networks under Dome A, Antarctica, J. Geophys. Res.-Earth, 118, 140–154, 2013. a, b

Wright, A. P. and Siegert, M. J.: A Fourth Inventory of Antarctic Subglacial Lakes, Antarct. Sci., 24, 659–664, 2012. a, b, c

Young, D., Schroeder, D., Blankenship, D., Kempf, S. D., and Quartini, E.: The distribution of basal water between Antarctic subglacial lakes from radar sounding, Philos. T. R. Soc. A, 374,, 2016. a, b

Young, D. A., Roberts, J. L., Ritz, C., Frezzotti, M., Quartini, E., Cavitte, M. G. P., Tozer, C. R., Steinhage, D., Urbini, S., Corr, H. F. J., van Ommen, T., and Blankenship, D. D.: High-resolution boundary conditions of an old ice target near Dome C, Antarctica, The Cryosphere, 11, 1897–1911,, 2017. a, b, c, d, e, f, g, h, i, j, k

Short summary
Our paper provides an important review of the state of knowledge for oldest-ice prospection, but also adds new basal geothermal heat flux constraints from recently acquired high-definition radar data sets. This is the first paper to contrast the two primary target regions for oldest ice: Dome C and Dome Fuji. Moreover, we provide statistical comparisons of all available data sets and a summary of the community's criteria for the retrieval of interpretable oldest ice since the 2013 effort.