Using ice-flow models to evaluate potential sites of million year-old ice in Antarctica

Finding suitable potential sites for an undisturbed record of million-year old ice in Antarctica requires slowmoving ice (preferably an ice divide) and basal conditions that are not disturbed by large topographic variations. Furthermore, ice should be thick and cold basal conditions should prevail, since basal melting would destroy the bottom layers. However, thick ice (needed to resolve the signal at sufficient high resolution) increases basal temperatures, which is a conflicting condition for finding a suitable drill site. In addition, slow moving areas in the center of ice sheets are also low-accumulation areas, and low accumulation reduces potential cooling of the ice through vertical advection. While boundary conditions such as ice thickness and accumulation rates are relatively well constrained, the major uncertainty in determining basal thermal conditions resides in the geothermal heat flow (GHF) underneath the ice sheet. We explore uncertainties in existing GHF data sets and their effect on basal temperatures of the Antarctic Ice Sheet, and propose an updated method based on Pattyn(2010) to improve existing GHF data sets in agreement with known basal temperatures and their gradients to reduce this uncertainty. Both complementary methods lead to a better comprehension of basal temperature sensitivity and a characterization of potential ice coring sites within these uncertainties. The combination of both modeling approaches show that the most likely oldest ice sites are situated near the divide areas (close to existing deep drilling sites, but in areas of smaller ice thickness) and across the Gamburtsev Subglacial Mountains.


Introduction
One of the major future challenges in the ice coring community is the search for a continuous and undisturbed icecore record dating back to 1.5 million years BP (Jouzel and Masson-Delmotte, 2010).The reason for such a quest is that the oldest part of the EPICA Dome C ice core has revealed low values of CO 2 from 650 000 to 800 000 yr ago (Lüthi et al., 2008), and is therefore out of phase with atmospheric temperature change.This questions whether such partial decoupling between the CO 2 record and climate had precursors over longer timescales (Jouzel and Masson-Delmotte, 2010).Marine records show evidence of a reorganization of the pattern of climate variability around 1 Myr ago, shifting from the "obliquity"-dominated signal, characterized by 40 000 yr of weak glacial-interglacial cycles, to the "eccentricity"-dominated signal with longer glacialinterglacial cycles (Lisiecki and Raymo, 2005).The origin of this major climate reorganization (the so-called mid-Pleistocene transition, MPT) remains unknown and may be intrinsic to a series of feedback mechanisms between climate, cryosphere and the carbon cycle (Jouzel and Masson-Delmotte, 2010).Alternatively, a recent study has demonstrated that climate oscillations over the past four million years can be explained by a single mechanism: the synchronization of nonlinear internal climate oscillations and the 413 000 yr eccentricity cycle (Rial et al., 2013).According to model calculations in conjunction with spectral analysis, Rial et al. (2013) find that the climate system first synchronized to this 413 000 yr eccentricity cycle about 1.2 million years ago, coinciding roughly with this MPT.A deep ice core covering a time span of more than one million years would shed a light on the mechanisms involved.ap of major drill sites and locations cited in the paper.Blue triangles depict the 379 l lakes (Wright and Siegert, 2012), major lakes are highlighted by black lines.

Published by Copernicus
24 Fig. 1.Map of major drill sites and locations cited in the paper.Blue triangles depict the 379 subglacial lakes (Wright and Siegert, 2012), major lakes are highlighted by black lines.
Deep ice-core drillings have been carried out in the past in Antarctica, reaching back in time over several hundred thousands of years.Amongst the longest records are Vostok (Petit et al., 1999), EPICA Dome Concordia (EPICA community members, 2004), Dome Fuji (Watanabe et al., 2003), and EPICA Dronning Maud Land (EPICA community members, 2006); all of these sites are depicted in Fig. 1.The longest record is from EPICA Dome Concordia, which extends over 800 ka into the past.What those records have all in common is that they are all recovered in the center of the ice sheet, and given the fact that the Antarctic Ice Sheet has been relatively constant in size over the last 13 million years (DeConto and Pollard, 2003), they are consequently undisturbed by dramatic changes in ice flow, contrary to the longest records from the Greenland Ice Sheet (Johnson, 2001;NEEM community members, 2013).
In theory, and in the absence of basal melting, these deep Antarctic records could reach several million years back in time with layers getting infinitesimally thin near the bottom.In reality, however, all deep records lack the bottom sequence as they are all found to be at the pressure-melting point, and lower layers are melted away or heavily disturbed due to complex basal processes.Furthermore, resolving deep records not only requires that the bottom sequence is unaltered, but that the ice is sufficiently thick so that the gas signal can still be retrieved and analyzed with sufficiently high resolution in the bottom layers.
In this paper we use two thermodynamic ice-sheet models to infer suitable areas for retrieving long ice-core records.We first investigate the most influential parameters having an effect on ice-core record length.Secondly, we apply this simple concept to evaluate uncertainties in geothermal heat flow (GHF) and use this uncertainty to guide the search for suitable drilling locations.Thirdly, we carry out a sensitivity analysis with a three-dimensional thermodynamical model (Pattyn, 2010) to determine the sensitivity of basal conditions to uncertainties in GHF, guided by a priori knowledge of basal conditions through the geographical distribution of subglacial lakes.

Why obvious drill sites are unsuitable
Obvious places to look for oldest ice are the deepest parts of the ice sheet, where ice is thick, and accumulation rates are low.However, a thick ice cover insulates very well and keeps the geothermal heat from escaping to the surface.Furthermore, we know that at least 379 subglacial lakes exist under the Antarctic Ice Sheet (Fig. 1), which implies that large portions of bottom ice should be at the pressure-melting point (Smith et al., 2009;Pattyn, 2010;Wright and Siegert, 2012), and therefore destroying bottom layers.Most subglacial lakes occur in the so-called Lakes District (stretching between subglacial Lake Vostok and Wilkes Land in East Antarctica), characterized by a thick ice cover and also low geothermal heat flow (Shapiro and Ritzwoller, 2004;Pollard et al., 2005).Therefore, GHF is not the main culprit in causing subglacial melt.
The interplay between GHF and accumulation rates is very subtle, as high GHF increases basal temperatures, while high accumulation rates cool down the ice mass.To illustrate this we calculate the minimum GHF needed to reach the pressure-melting point at the bottom of any ice mass as a function of environmental parameters.This can easily be determined analytically (Hindmarsh, 1999;Siegert, 2000).Using the simplified model of Hindmarsh (1999), valid in the absence of horizontal ice advection due to motion, the minimum heat flow G min (mW m −2 ) needed to reach the pressure-melting point is where and where In Eqs.
The result is illustrated in Fig. 2, displaying the minimum GHF needed to reach the pressure-melting point at the base of an ice sheet as a function of ice thickness H and surface accumulation rate ȧ, based on Eqs.(1)-(3) for a mean surface temperature T s = −50 • C. Despite these low surface temperatures, pressure-melting point is reached for relatively low values of GHF, as long as the ice is thick and accumulation rates are small, which is rather typical for the interior parts of the East Antarctic Ice Sheet.For high accumulation rates, one needs a significantly higher GHF to reach melting point at the base for a given ice thickness.
Even the low GHF values in Fig. 2, in conjunction with low accumulation rates, are quite common over the central part of the East Antarctic Ice Sheet (van de Berg et al., 2006), which hampers the retrieval of a long time sequence under thick ice conditions.Despite the simplicity of the model, it can be applied to central parts of the Antarctic Ice Sheet, where horizontal advection is absent or negligible, to explore suitable drill sites as a function of known (or estimated) geothermal heat fluxes.
3 Uncertainties in Antarctic GHF and the location of oldest ice

Data sets and model setup
The above simple model is applied to central areas of the Antarctic Ice Sheet that are characterized by slow ice motion (i.e., the vicinity of ice divides).To do so, ice thickness is taken from the recent BEDMAP2 compilation (Fretwell et al., 2013) and resampled on a 5 km grid.Surface-mass balance is obtained from van de Berg et al. (2006) and van den Broeke et al. (2006), based on the output of a regional atmospheric climate model for the period 1980 to 2004, and calibrated using observed mass balance rates.Surface temperatures derived from van den Broeke ( 2008), based on a combined regional climate model, calibrated with observed 10 m ice temperatures.Using these data sets enables us to calculate the minimum required GHF to reach the pressure-melting point at the bed.Since this is a vertical-column model with no horizontal advection, it is only valid for divide areas.The simulations are therefore only carried out for regions with horizontal velocities smaller than 2 m yr −1 .Ice sheet velocities in divide areas are determined based on balance velocities, stating that the mass of ice flowing out of any area within the horizontal domain ∇ H (x, y) exactly equals the sum of the inflow and the ice accumulated over the area (Budd and Warner, 1996;Fricker et al., 2000;Le Brocq et al., 2006): with and where 1 and 0 are the bottom and the surface of the ice sheet (m a.s.l.), respectively.Integrating Eq. ( 4) over the whole surface of the ice sheet, starting at the ice divides, one obtains the vertically averaged horizontal balance velocities v H = (v x , v y ).Details of this procedure are given in Pattyn (2010).
Using the above data sets, the minimum geothermal heat flow G min from Eq. ( 1) needed to reach the pressure-melting point at the bed is calculated for the areas of the Antarctic Ice Sheet where horizontal flow velocities are < 2 m yr −1 and where ice thickness H > 2000 m.This ice thickness is considered to be the lower limit for possible recovery of a million-year old climate signal (Fischer et al., 2013) .The limit of 2 m yr −1 for the horizontal surface velocities was chosen from a simple age calculation using a Lagrangian algorithm.This showed that for a surface velocity of 2 m yr −1 , the maximum span of the origin of the ice can be several hundreds of kilometers (up to 1000 km), which will definitely complicate the interpretation of the climatic signal.Therefore, larger values should definitely be avoided.Smaller values of surface velocities lead to a similar pattern of potential oldest ice sites, but generally smaller in extent, which hamper a good visualization on a continental scale.The calculated values of G min are compared below to other GHF databases.
Several data sets of derived GHF underneath the Antarctic Ice Sheet exist.The first one (G 1 ) uses a global seismic model of the crust and the upper mantle to guide the extrapolation of existing heat-flow measurements to regions where such measurements are rare or absent (Shapiro and Ritzwoller, 2004).The second GHF database (G 2 ) stems from satellite magnetic measurements (Fox Maule et al., 2005).Their values of GHF are in the same range as Shapiro and Ritzwoller (2004), but the spatial patterns are markedly different, and the (G 2 ) values are considerably higher in many regions.The third data set (G 3 ) represents a recent update of (G 2 ) derived by Purucker (2013).This uses low-resolution magnetic observations acquired by the CHAMP satellite between 2000 and 2010, and produced from the MF-6 model following the same technique as described in Fox Maule et al. (2005).While the technique is similar, GHF values are considerably lower than the latter, and even lower than those derived from the seismic model (Shapiro and Ritzwoller, 2004).

Results
In view of the large uncertainty in GHF estimates, we combined all three data sets into two databases (i.e., a mean GHF, G, and a standard deviation, σ G ).The latter is calculated based on the inter-data set variability and the standard deviation given for the Shapiro and Ritzwoller data set in the following way: The calculated values of G min are directly compared to the map of mean GHF.For G min < G, the observed GHF is too elevated to prevent the bottom ice from reaching pressure melting and most likely (within error bounds) the ice is temperate.For G min > G the minimum GHF needed to reach pressure melting at the base is higher than the value reported.Of course, this information needs to be further evaluated against the dispersion between the GHF data sets, represented by σ G .The result is shown in Fig. 4, where the rectangular area points to the potentially most suitable conditions in terms of basal temperature (i.e., the largest excess of minimum GHF above actual GHF in combination with the lowest variability between the three GHF data sets).Although the limits of the rectangle are arbitrarily chosen, they assure that the probability of reaching cold ice at the bed is sufficiently high.The furthest to the right in Fig. 4, the colder the bed because a significantly higher GHF than observed is needed to make the bed temperate; the lower the value of σ G , the more  likely there is a small spread (hence reduced uncertainty) in GHF, so that the observed value is likely.The thickest ice, as expected, corresponds to zones that are temperate (negative values of G), while for large positive G and small σ G , ice is also the thinnest.
These restrictions (combined with the ice-flow speed limit and minimum ice thickness) mean that only very few areas in the central part of the Antarctic Ice Sheet can be considered likely to host cold-bed conditions.The largest zone is situated near Dome Argus on top of the Gamburtsev Subglacial Mountains (Fig. 5).However, subglacial mountain ranges are likely characterized by an uneven bed topography (Bell et al., 2011), which may also hamper the interpretation of the paleo-climatic signal (Grootes et al., 1993).Other potential areas are situated around Dome Fuji as well as on Ridge B, between subglacial Lake Vostok and Dome Argus.The Dome Concordia area seems less prone to cold basal conditions, due to the large uncertainty in GHF and the thick ice, which makes temperate conditions more likely.This is corroborated, in reality, by the abundance of subglacial lakes around Dome Concordia.

Thermomechanical ice-flow modelling
The simple thermodynamic model used in the previous section neglects horizontal advection, which, even in the interior of the Antarctic Ice Sheet, plays a significant role in determining the thermal properties of the ice-bed interface.We therefore extend the model that does not include horizontal flow and present a more advanced thermomechical icesheet model to calculate basal temperatures for a set of given boundary conditions and applied to the whole Antarctic Ice Sheet.Moreover, we try to reduce uncertainties in GHF by incorporating actual information on bed properties, such as the geographical distribution of subglacial lakes.

Model description
The thermodynamical model used for this purpose is the same as described in detail in Pattyn (2010).The major differences are related to the way the horizontal flow field is calculated.Moreover, we use a new series of data sets on ice thickness (see Sect. 3.1) and geothermal heat flow.
The thermodynamic equation for the temperature distribution in an ice mass is given by where T is the ice temperature (K) and v = (v x , v y , v z ) is the three-dimensional ice velocity vector (m yr −1 ).The last term on the right-hand side represents internal heating rate per unit volume (Pattyn, 2003), where ε and σ are effective strain rate and effective shear stress, respectively.Horizontal diffusion is neglected, and the temperature field is considered to be in steady state ( ∂T ∂t = 0).Boundary conditions for Eq. ( 7) are the surface temperature T s and a basal temperature gradient, based on the geothermal heat flux: where G is the geothermal heat flux entering the base of the ice sheet and the second term on the right-hand side of Eq. ( 8) is heat produced due to basal sliding.τ b is the basal shear stress, and can be defined as τ b = −ρgH ∇ H s, where ∇ H s is the surface slope.Whenever the pressure-melting point is reached, the temperature in the ice is kept at this value T pmp = T 0 − γ (s − z).

Velocity field
Ice sheet velocities are obtained from a combination of satellite-derived velocities from radar interferometry and modelled velocities.Satellite-derived velocities are available for almost the entire continent (Rignot et al., 2011), but are only relevant in the coastal areas and for fast ice flow.Generally speaking, the error associated with the slow flowing areas is substantially higher than 100 % (Rignot et al., 2011).Furthermore, in the vicinity of the South Pole, interferometric velocities are lacking due to the sun-synchronous orbit of satellites.To fill in the gaps and to guarantee a continuous flow field for our simulations, a heuristic method was implemented that uses interferometrically derived velocities for flow speeds above 100 m yr −1 and modelled velocities for flow speeds below 15 m yr −1 .Modelled velocities are derived from balance velocities, described in Sect.3. Between 15 and 100 m yr −1 , both modelled and interferometric velocities are combined as a fraction of flow speed, in order to keep the transition between both data sets as smooth as possible and to guarantee a correct flow direction.Similar to Pattyn (2010), a shelfy-stream model is used to correct for the ice flow over large subglacial lakes and basal sliding is only allowed when the base is temperate or within a range of 1 K of subfreezing temperatures.The three-dimensional horizontal velocities are then determined from the shallow-ice approximation (Hutter, 1983), by where basal sliding v b is represented by a Weertman sliding law (Pattyn, 2010).The vertical velocity field is derived from mass conservation combined with the incompressibility condition for ice.Given an ice sheet in steady state, a simple analytical expression can be obtained, based on the horizontal balance velocities (Hindmarsh, 1999;Hindmarsh et al., 2009).Expressed in local coordinates, and in the absence of subglacial melting, this leads to The numerical solution of the model is detailed in Pattyn (2010).For all experiments, n = 3 was used, which corresponds to the isothermal case.However, in the thermomechanically coupled case, the exponent is larger (Ritz, 1987), which results in a different shape of the vertical velocity profile.Therefore, advection being more concentrated to the surface, this leads to warmer basal conditions compared to the isothermal case.However, this effect is most pronounced in areas where horizontal velocity gradients ∂v x /∂x, ∂v y /∂y are more important.Since we concentrate on the central areas of the ice sheet, this bias (underestimation of basal temperatures) will have a limited effect.

Input data calibration
Major input data sets are already described in Sect.3. In this section, we will focus on the improvements made to the initial GHF data sets in order to reduce uncertainty in GHF.
Direct measurements of GHF are very rare, and are usually obtained from temperature measurements in boreholes of deep ice-core drillings.Basal temperature gradients in observed temperature profiles of deep boreholes, compared with values from the three GHF data sets, show rather large discrepancies (Pattyn, 2010).Therefore, the three GHF data sets were corrected using observed basal temperature gradients, surface temperature and accumulation rates, in such a way that modeled temperature profiles match as closely as possible with the observed ones (Pattyn, 2010).This type of correction is made for sites where temperature profiles are available: Byrd (Gow et al., 1968), Taylor Dome (G. Clow and E. Waddington, personal communication 2008), Siple Dome (MacGregor et al., 2007), Law Dome (Dahl-Jensen et al., 1999;van Ommen et al., 1999), Vostok (Salamatin et al., 1994;Parrenin et al., 2004), South Pole (Price et al., 2002), Dome Fuji (Fujii et al., 2002;Hondoh et al., 2002), EPICA Dome C (Parrenin et al., 2007, C. Ritz, personal communication 2008), and EPICA DML (Ruth et al., 2007).The applied method consists of determining the difference between observed (o) and corresponding database values and to adapt a Gaussian function for a sufficiently large area of influence.For a variable in the database X (either surface accumulation, surface temperature or geothermal heat flux), its corrected value X c based on an observation X o is obtained by where (x, y) is the horizontal distance from this observed position (0,0).The area of influence is dictated by σ and calculations were performed for σ = 0, 20, 50, 100, and 200 km.
A σ -value of 0 means that no correction is carried out.Larger spans describe potential influence areas, and give a wider range than those explored in Pattyn (2010).As such, by tuning GHF (constraining the vertical temperature gradient) and surface-mass balance (constraining vertical advection), the difference between modeled and observed temperature profiles is less than 2 K.The remaining difference is still due to horizontal advection, which is a model output, as well as past changes in surface temperature that were not taken into account in the model.

Subglacial lake correction
Numerous subglacial lakes have been identified from radioecho sounding.An initial inventory contained 145 lakes (Siegert et al., 2005), and more than 230 have been added since (Bell et al., 2006(Bell et al., , 2007;;Carter et al., 2007;Popov and Masolov, 2007;Fricker et al., 2007;Fricker and Scambos, 2009;Smith et al., 2009), such that at least 379 subglacial lakes of varying size are now known to exist (Wright and Siegert, 2012).Subglacial lakes are usually identified from radio-echo sounding (RES) in which they are characterized by a strong basal reflector and a constant echo strength corresponding to a smooth surface or have been identified through surface elevation changes using satellite altimetry, theorized to be the surface expression of rapid drainage or filling of subglacial lake sites (Fricker et al., 2007;Pattyn, 2011).Subglacial lakes are used to constrain the GHF data sets, considering them to be at the pressure-melting point.As such, we calculate the minimum GHF needed to reach the pressure-melting point using (1) for any position of a subglacial lake.The value for G min thus obtained is a minimum value, which means that if at that location the database p: Mean basal temperature according to the ensemble of 15 experiments (see text etails), corrected for the dependence on pressure.The color scale is truncated at ttom: Root Mean Square Error (RMSE, • C) according to the same ensemble.The are generally small, and tend to correspond to a higher RMSE.contains a higher value, the latter is retained.Spatial corrections are subsequently applied using the Gaussian function defined in Eq. ( 11) for different areas of influence as defined above.

Ensemble model results
The temperature field in the ice sheet was calculated for 15 different sets of boundary conditions: the three data sets of GHF (Shapiro and Ritzwoller, 2004;Fox Maule et al., 2005;Purucker, 2013), and each of the data sets corrected for subglacial lakes and existing temperature profiles for influence area size σ = 0 (no correction), 20, 50, 100, and 200 km.The result is given in Fig. 6, representing the mean basal temperature of the 15 experiments, corrected for the dependence on pressure melting, and the root mean square error (RMSE) corresponding to the different experiments.
Low values of RMSE correspond to zones where the correction is effective and the difference between the experi-   ments is low, or areas that despite the variability in GHF are always at the pressure-melting point.This is the case for the central part of the West Antarctic Ice Sheet, as well as extensive zones in the Lakes District, where the dense network of subglacial lakes keeps the bed at melting point.
In the ensemble experiments, the relation between accumulation (vertical advection), ice thickness and basal temperature is less straightforward than with the simple model.The focus of the full model is to reduce uncertainties on GHF using proxy data, and therefore the RMSE guides us towards more suitable sites.Figure 7 summarizes the most suitable drilling areas based on the full model for flow speeds < 2 m yr −1 , ice thickness H > 2000 m, and a basal temperature < −5 • C. The color scale denotes the RMSE based on the ensemble experiments.We deliberately excluded basal temperatures higher than −5 • C, a value considered to be sufficiently far away from the melting point in view of our model approximations.Suitable areas characterized by low values of RMSE (hence smaller spread in basal temperatures according to the ensemble experiments) are found near existing ice-core sites where a temperature gradient is at hand (Dome Concordia, Dome Fuji and Vostok).Since all three sites are at or close to the pressure-melting point at the base, suitable cold-based sites do not coincide exactly with the icecore locations, but lie nearby in locations where ice is thin enough to reduce basal ice temperatures.
Similarly to the simple model, suitable sites (sufficiently low basal temperature) are found in the Gamburtsev Mountain region as well as along Ridge B. However, the ensemble analysis results in a larger range of basal temperature due to either the lack of basal temperature gradient constraints and/or the absence of subglacial lakes.Both regions are characterized by relatively low basal temperatures (Fig. 6) and are unlikely to reach the pressure-melting point, despite the large ential locations of cold basal conditions in areas with ice thickness H > 2000 m, and flow speeds < 2 m yr -1 according to the simple model (depicted in Fig. 5) and the model (depicted in Fig. 7).
31 according to the simple model (depicted in Fig. 5) and the ensemble model (depicted in Fig. 7).
RMSE due to -mainly -differences between the GHF data sets.

Discussion and conclusions
Since both the simple and ensemble model results are complementary (but not totally independent) in nature, they can be combined to form a joint data set in order to investigate common grounds.The analysis is limited to flow speeds < 2 m yr −1 and ice thickness H > 2000 m, which are considered as suitable conditions for retrieving and resolving ice older than one million years.Given the uncertainty in GHF originating from the large dispersion between the different data sets (both spatially and in terms of absolute values), we apply a set of constraints to select the suitable sites for preservation of million year-old ice: (i) the minimum GHF needed to reach the pressure-melting point should be at least 5 mW m −2 higher than the mean value from the combined GHF data sets; (ii) the variability between the GHF data sets for a given site expressed by the standard deviation σ G should be < 25 mW m −2 ; (iii) the mean basal temperature according to the ensemble model calculations should be < −5 • C (but lower values are favored).Results are displayed in Fig. 8.We explored different values for these constraints, but the general pattern remains the same.The main effect is the stronger the constraint, the smaller the areas, but the geographical distribution is not altered.
Due to the velocity and ice thickness constraints, all sites are situated near the ice divides.Not surprisingly, areas near the major drill sites and where temperature profiles are available (Dome Fuji, Dome Concordia, Vostok and South Pole) are also retained.These are not the sites themselves, but zones of smaller ice thickness in their vicinity.Finally, suitable areas are found across the Gamburtsev Mountains and Ridge B (between Dome Argus and Vostok).The former is characterized by a much larger spatial variability in bedrock topography, while the latter may suffer from sparse constraints on ice thickness (according to Fig. 2 in Fretwell et al., 2013).
Subglacial topography is a key factor in determining suitable sites for oldest ice.Given the strong relationship between basal temperatures and ice thickness, as depicted by Fig. 2, it is quite likely to find suitable cold-based spots in the vicinity of deep ice-core sites that have the bottom ice at or near the pressure-melting point.Areas that should be avoided are those in which a large number of subglacial lakes are found, such as the Lakes District, where even low values of GHF are sufficient to keep the ice at pressure melting.
Another factor that may influence basal conditions is the glacial-interglacial history of the ice sheet and the timescales needed for the ice sheet to adapt thermally to different climates.Moreover, the temperature calculations made in this study are based on present-day observed parameters of surface temperature, ice thickness and accumulation rate.To test this effect, we calculated the minimum geothermal heat flow G min needed to keep the base at the pressure-melting point for environmental conditions that are the mean for a time span covering a glacial-interglacial cycle.We reduced the surface temperature T s by 6 K, reduced the surface accumulation rate ȧ to 60 % of its current value, and reduced ice thickness H by 100 m, which is appropriate for the divide areas.The results are surprisingly similar to the previously calculated values, and are therefore not shown separately.The main reason is that for this spread of values the reduced accumulation rate (which reduces vertical advection, hence warms the bottom ice layers) is largely counteracted by the decrease in surface temperature.However, we note that both calculations (present-day and mean glacial-interglacial) relate to steady-state conditions, which in reality is not the case.For instance, Rogozhina et al. (2011) demonstrate that for the Greenland Ice Sheet, basal temperature differences between an ice sheet initialized by a steady simulation (as in this study) and those generated by a paleoclimatic simulation can be up to 4.5 • C.
We suggest that the results presented here should not be used as a sole guide in the process of detecting suitable coldbased areas for retrieving a long ice-core record, due to a number of factors that were not taken into account: 1. Areas characterized by subglacial mountains or other bedrock relief variability may be thermally conducive to the preservation of ancient ice, but the topographic variability may well hamper the deciphering of the climate signal due to complex processes, such as ice overturning (NEEM community members, 2013) or refreezing (Bell et al., 2011).
2. The upper limit on the flow velocity of 2 m yr −1 may be too high for reconstructing the climate signal without having to rely heavily on ice-flow models for corrections due to upstream advection.In theory, ice could have traveled over several hundreds of kilometers before reaching the ice-core site (Huybrechts et al., 2007).Taking into account shifts in ice divides over glacial-interglacial periods would also influence the flow direction over time.
3. The spatial variability of GHF may in reality be much higher than represented in the three GHF data sets.
4. Areas where bedrock elevation data are unavailable (or where interpolation is based on sparse data) may be wrongly classified in the above analysis, and some suitable areas thus overlooked.
In summary, this paper gives an overview of the factors that influence the basal thermal conditions of the Antarctic Ice Sheet, which are useful to guide the search for potential deep drilling sites for IPICS oldest ice (more than one million years old) records.The two complementary thermal models that were employed virtually lead to similar results: most suitable sites are situated in the vicinity of the ice divides and close to areas where deep drillings have been carried out in the past.Another suitable area is in the vicinity of the Gamburtsev Subglacial Mountains.Ice thickness is found to be a major limiting factor, since too thick ice may lead to temperate basal conditions.This is the main reason why most of the current deep drillings have been found at or close to pressure melting point at the base.While this paper gives an overview of continental-scale basal conditions of the Antarctic Ice Sheet, the processed data sets from both the simple (G, σ G ) and the full model (T , RMSE T ) are made available online together with simple MatLab scripts to allow for a more detailed search/zoom for potential sites, based on the figures presented here.

Fig. 2 .
Fig. 2. Minimum GHF (mW m −2 ) needed to keep the bed at the pressure-melting point as a function of surface accumulation rate (ice equivalent, IE) and ice thickness and in the absence of horizontal advection.Results are shown for a mean surface temperature of T s = −50 • C.
) Both are depicted in Fig. 3. High values of σ G indicate a large dispersion between the three data sets.These are essentially found in West Antarctica and along the Transantarctic Mountains.The lowest values are restricted to the central parts of the East Antarctic continent.

Fig. 3 .
Fig. 3. Top: Mean GHF G (mW m -2 ) based on GHF estimates by Purucker (2013), F et al. (2005) and Shapiro and Ritzwoller (2004).Bottom: Standard deviation σ G o datasets.The magenta triangles are the major drill sites (from top to bottom): Dome F Argus, South Pole and Dome Concordia.

Fig. 3 .
Fig. 3. Top: Mean GHF G (mW m −2 ) based on GHF estimates by Purucker (2013), Fox Maule et al. (2005) and Shapiro and Ritzwoller (2004).Bottom: Standard deviation σ G of the GHF data sets.The magenta triangles are the major drill sites (from top to bottom): Dome Fuji, Dome Argus, South Pole and Dome Concordia.

Fig. 4 .
Fig. 4. Scatter plot of G = G min −G versus σ G for all points with ice thickness H > 2000 m and horizontal flow speed < 2 m yr −1 .The color scale depicts ice thickness for each of the grid points.Negative values of G show where the pressure-melting point is reached, hence basal melt occurs.Positive values mean that the minimum required heat flow to reach the pressure-melting point is higher than the mean of the three GHF data sets.Points lying within the rectangle are likely to be cold-based, taking into account the variability of GHF.

Fig. 5 .Fig. 5 .
Fig. 5. Potential locations of cold basal conditions in areas with ice thickness H (colorbar) and horizontal flow speeds < 2 m yr -1 , for ∆G > 5 mW m -2 and σ G < 25 mW as calculated with the simple model. 29

Fig. 6 .
Fig. 6.Top: Mean basal temperature according to the ensemble of 15 experiments (see text for more details), corrected for the dependence on pressure.The color scale is truncated at −10 • C. Bottom: root mean square error (RMSE, • C) according to the same ensemble.The cold areas are generally small, and tend to correspond to a higher RMSE.

Fig. 7 .
Fig. 7. Potential locations of cold basal conditions in areas with ice thickness H > 20 horizontal flow speeds < 2 m yr -1 and basal temperatures as calculated with the fu -5 • C. The colorbar denotes the RMSE ( • C) based on the ensemble calculations. 30

Fig. 7 .
Fig. 7. Potential locations of cold basal conditions in areas with ice thickness H > 2000 m, and horizontal flow speeds < 2 m yr −1 and basal temperatures as calculated with the full model < −5 • C. The colorbar denotes the RMSE ( • C) based on the ensemble calculations.

Fig. 8 .
Fig. 8. Potential locations of cold basal conditions in areas with ice thickness H > 2000 m, and horizontal flow speeds < 2 m yr −1according to the simple model (depicted in Fig.5) and the ensemble model (depicted in Fig.7).