Geothermal heat flux and basal melt rate in the Dome C region inferred from radar reflectivity and thermal modelling

Basal melt rate is the most important physical quantity to be evaluated when looking for an old-ice drilling site, and it strongly depends on the geothermal heat flux, which is poorly known under the East Antarctic ice sheet. The wetness of the ice-bed interface can be assessed from radar echoes on the bedrock, considering that a wet bedrock has a stronger reflectivity than a dry one. But, as the basal conditions depends on the climatic forcing lagged by the thick ice, the basal ice may 5 be cold today whereas it was in average temperate in the past. Accordingly, the risk of mismatch between present and past conditions must be evaluated, and more generally the geothermal heat flux and basal melt rate in the Dome C region, which is a promising old-ice target. Here, we run a 1D heat model over the last 800 ka in inverse mode to constrain the value of geothermal heat flux by assessing a critical ice thickness, i.e. the minimum ice thickness that would allow the local melting distribution 10 at present. A regional map of the geothermal heat flux is then inferred on a 80 km×130 km area, and shows a N-S-oriented gradient, with a value range of 48− 60 mWm−2. The forward model is then emulated by a polynomial function, to compute a time-averaged value of the basal melt rate, spatially variable over the region. Two main subregions appear to be free of basal melting because of the thin overlying ice, and a third one, north of Dome C, because of a low geothermal heat flux. 15


Introduction
The Dome C region (East Antarctica) has been a region of strong interest for the paleoclimatic science during the last decades, with two deep ice cores having been retrieved, including the oldest one ever dated (Lorius et al., 1979;Jouzel et al., 2007).Furthermore, the Dome C region has been identified as one possible candidate to host even older ice in its vicinity (Fischer et al., 2013; Based on radar equation, Zirizzotti et al. (2012) proposed a method to recognize the wet areas using the Italian RES dataset collected around Dome C (http://labtel2.rm.ingv.it/antarctica/).The authors adopted a linear model of electromagnetic wave adsorption in the ice column, based on the EDC ice core analysis, to account for the amplitude differences (in dB) between the ice surface and the bottom echoes.In this paper, the thresholds considered to ascribe an echo to a wet or dry 75 basal contact were ≥ 7.7 dB and < 1 dB respectively.Figure 1 shows the distribution of dry and wet areas under Dome C, revealing interesting patterns: (i) the Concordia trench is characterized by the presence of wet points only, which is a consequence of the very thick ice; (ii) the tops of the two main bedrock reliefs (candidate A and candidates B − C − D) appear to be dry; (iii) almost all the northern points are dry, despite the very thick ice at this place; (iv) in between, and in particular 80 under the drilling site of EPICA, dry and wet points are scattered, showing no clear trend.
We emphasize that these observations show the bed conditions at present, whereas the historical conditions all along the glacial/interglacial periods should be investigated to evaluate the quality of a future old-ice drilling site.The absence of basal water today means a cold basal ice, but this may simply be a consequence of the very strong temperature signal of the Last Glacial Maximum 85 (LGM) lagged by the thick ice.The present cold state may not be representative of the whole glacialinterglacial periods, and thawing could have occured in the past under different climatic conditions.
Hence, this paper primarily aims to assess the risk of past temperate conditions at places known to be cold today.To this end, we present a 1D heat model forced by reconstructed climatic conditions, and run it in two ways.First, in inverse mode to infer the value of the GHF in the Dome C region, 90 using radar echoes as observational constraints.Second, in forward mode to compute the average past basal melt rate under the GHF inferred at the first step.We present an easy-to-use emulator of the past basal melt rate, which is a convenient thermal boundary condition for future modelling works, and key-criterium in the location of old-ice drilling sites.

Heat model 95
This section presents the heat model accounting for the relations between GHF, ice thickness and vertical advection of the ice.We first explain the main assumptions upon which the model is based.
The areas on which these assumptions are valid are presented in section 3.

1D-assumption 100
In a dome region, the horizontal velocities are very small (a few centimeters per year, (Vittuari et al., 2004)), so the horizontal advection terms can be neglected.Similarly, the deformational heat is several orders of magnitude smaller than the vertical advection term.Finally, due to the small aspect ratio (thickness/typical horizontal length), the temperature field is mainly vertically stratified, so that the horizontal diffusion can be neglected as well.Thus, the heat balance of ice has dependences on 105 the vertical direction only.

Water circulation at the base of the ice sheet
Since the wetness indicated by the basal echoes will be used as constraints on the inverse model, we first have to pinpoint the origin of the observed basal water.Whether the local ice is temperate and the water comes directly from the thawed ice, or the local (cold) ice energy budget was imbalanced 110 by the latent heat of a water flow coming from a temperate place, following the gradient of hydraulic potential.In the Dome C area, the ice surface is very flat and the hydraulic potential almost follows the bedrock heights, thus corresponding to the local maxima of hydraulic potential.Water cannot go upwards on these topographic features (Fig. 2).Hence, the water observed on the lee of a relief has a local origin.

115
Could the dry areas observed on Fig. 1 correspond to a thawing ice whose meltwater would have flown away?Meltwater can be driven out by different types of hydrological networks: connected cavities (Lliboutry, 1968;Kamb, 1987) or efficient structured channels network (Röthlisberger, 1972).Weertman (1972) showed that these channels cannot form upstream of the hydrological network, which is the case here, and would rather take the form of a film of water.Whatever the exact type 120 of structure (film or cavities), basal melting is a continuous process that permanently feeds the hydrological network.Upstream of this network, the thawed water cannot be driven out faster than the melt rate, so that some water would always remain at the base of the ice.As a consequence, we suggest that the areas seen as dry in Fig. 2 effectively correspond to cold ice.

GHF spatial variability 125
According to the little that we know about the geology, we may assume that the GHF is uniform on a ∼ 10 km-lengthscale.Under such an assumption, the presence of basal water is only the consequence The Cryosphere Discuss., doi:10.5194/tc-2017Discuss., doi:10.5194/tc- -23, 2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License. of a thicker ice, and not of GHF spatial variations.As such, thickness and reflectivity have to be locally correlated, and the model will be run for small areas where this correlation is reasonably ascertained.

Geometry and coordinate system
We consider here a one-dimensional vertical domain, oriented upwards along the z-axis.Instead of being set in the z-coordinates, the equations are expressed in reduced depth ζ = (s − z)/(s − b), s being the surface height, and b the bed height.In the ζ-coordinate system, the domain size remains the same whatever the ice thickness H, and no remeshing is needed in the resolution of 135 the finite difference scheme.Therefore the ice thickness is a simple parameter evolving with time depending on the accumulation forcing at the surface.The thickness evolution is accounted for by the conceptual model of Parrenin et al. (2007), wich emulates for the Dome C the 3D large scale simulations of Ritz et al. (2001).At each timestep, explicit expressions for the thickness and the bedrock elevation are solved, that depend on the accumulation rate and six tabulated parameters.

140
The typical difference between ice thicknesses of glacial and interglacial periods is 150 m.

Heat equation
The heat balance of ice only depends on the vertical, and is written as follows in the ζ-coordinate system (Ritz et al., 1997): where K is the thermal conductivity of the ice, c its heat capacity, and D the relative density w.r.t.
ice.The vertical velocity u ζ accounts for the true ice velocity, but also for the evolution of the ice thickness, which changes the relationship between z and ζ.

Velocity model
The vertical advection of ice acts in the heat balance by transporting cold towards the bed.Instead of 150 solving the equations of motion, it will be basically accounted for by a 1D shape function (Parrenin et al., 2007;Ritz, 1987): (2) The Cryosphere Discuss., doi:10.5194/tc-2017Discuss., doi:10.5194/tc- -23, 2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.
the temperature profile, since it controls the vertical advection of ice from the surface to the base.The lower p, the more non linear vertical velocity profile.In practice, above a typical value of p = 10, the profile is very close to being linear (Fig. 3), and higher values of p do not change the profile 160 significantly.The vertical velocity u ζ in the ζ-coordinate system is then expressed as follows in the case of a shape function in a dome region (Parrenin et al., 2007): which is equivalent to (5) 165

Ice thermal properties
The specific heat capacity and heat conductivity depend on the absolute temperature as follows (Cuffey and Paterson, 2010, p. 400): The presence of firn at the surface of the ice lowers the conductivity at the ice sheet surface, and this must be accounted for (Cuffey and Paterson, 2010, p. 401): where K i is the conductivity of the ice.The density profile of the firn layer is taken from Parrenin et al. (2007).the air content is about 0.1 cm 3 /g (Martinerie et al., 1992), whereas the nitrogen saturates in ice at ∼ 2.6 cm 3 g −1 under a pressure of 27 MPa (Wiebe et al., 1932).The air is far from saturation for ice sheets.Only the partial pressure of air in the ice P should be accounted for, so that the dependence 185 of the melting temperature T m on the pressure P (in MPa) and P is expressed as (Ritz, 1992): In ice sheets, P is of the order of 1 MPa.This is an unusual choice for such an important parameter, but we argue that Eq.( 9) is consistent with the temperature profile of the EPICA Dome C ice core, where the deepest measured temperature was 270.05K at 3223 m, 50 m above the bedrock. 190 The temperature profile can be extrapolated to the bedrock (similar to Dahl-Jensen et al. (2003) at North GRIP) to 271.04 K.The melting temperature computed with B = 0.098 K.Pa −1 would be 0.8 K too low, whereas, it is found to be 270.96K with Eq.( 9).As the ice moves very slowly in the dome region, we assume that the pressure is only isostatic.Once the temperature field is computed, the melt rate at the bottom can be explicitely known by where Φ g is the GHF and L the latent heat of ice (J kg −1 ).For cold basal ice, Eq. ( 10) is used as a Neumann boundary condition, equating m to zero.

Boundary condition at the ice surface
The atmospheric temperature forcing is continuously transferred through the whole column of ice, 200 so that the present thermal conditions at the bed is the result of the whole climatic history.The atmospheric paleotemperatures are known over the last 800 ka from the δD measurements on the EPICA Dome C ice core (Jouzel et al., 2007).Here we use the model and notations of Parrenin et al. (2007), linking the isotope content of the ice to the temperature difference w.r.t. a reference temperature T 0 : The influence of the chosen value for α must be examined, since uncertainties affect our knowledge of the isotopic thermometer on long timescales in Antarctica, and α may vary by -10% to +20% (Jouzel et al., 2003).On top of the nominal value of 1/6.04 K, α will be given two additional values of 0.13 and 0.20, chosen so that the maximum temperature difference over climatic periods with the  2007), who considered an exponential accumulation model: For Dome C, the value of the β coefficient was evaluated at 0.0156±0.0012by an inverse method constrained by known age markers.Hence, we will check the sensitivity of our model for three 215 values of β (0.0138, 0.0156 and 0.0175).The two extreme values being chosen so that the maximum accumulation difference over climatic periods with the nominal run is ±0.25 cm a −1 , which is a more intuitive way to express the sensitivity.
Regarding the choice of initial state, we considered that the duration needed for a step climatic signal to reach the bedrock is ∼ 10 ka, and ∼ 100 ka to stabilize.As we do not know any true initial 220 state for the temperature profile in the past, we decided to run the model over the whole reconstructed period (800 ka), so that the computation is independent of the initial state and the final vertical temperature profile is as realistic as possible.

Numerical method
Equations are solved with a finite difference implicit scheme on a 50-element mesh, and the temporal 225 dependence of the model is solved on a 1000-year timestep.These values were selected as trade-offs between result accuracy and computation speed.
The physical coupling between flow and heat content of the ice is accounted for by non linear iterations, in which updated values of m, K and c are computed, till the temperatures do not show a discrepancy larger than 10 −5 K between two iterations.At each timestep, the heat equation is solved 230 for a boundary condition corresponding to the thermal state of the basal ice (temperate or cold).The computed solution may be inconsistent with the imposed boundary condition, i.e. that the pressure melting point may be exceeded (cold → temperate), or the melt rate may become negative (temperate → cold).In these cases, the equation is solved again with the new consistent boundary condition.

Measurement spots 235
We define the critical ice thickness H c as the minimum ice thickness for which present basal melting becomes possible for a given GHF and a given p.The model is run for increasing values of the GHF, till the melting point is reached at present, so that we determine a unique tuple (H c , p, Φ g ).As the unknown is the GHF, an estimation of this critical ice thickness is needed, and this is done in two steps from the wet/dry map displayed on Fig. 1.

240
First, we determine spots for which ice thickness and reflectivity are correlated, i.e. for which the top of bed reliefs are dry and their lees are wet (see Fig. 2).The melting starts to be physically possible in between.Ten corresponding spots are selected (black rectangles in Fig. 1), that are hereafter Second, we evaluate the critical ice thickness for each spot.Along the radar line, going upwards along a bed relief, we select the wet point for which the ice is the thinnest, and the dry point for which the ice is the thickest.The same is done going downwards, so that four points are selected along the radar line (black arrows on Fig. 2).The critical ice thickness H c is estimated by the average of the 250 ice thickness measured at these points.The uncertainty on H c is taken as the half of the standard deviation of the four thicknesses, so that almost all the possible values (H c ± 2σ) range within the extreme thicknesses measured at the considered spot.For the large spots C3 and C6, this operation is done twice, on the two perpendicular radar lines going through them (radar lines C and 3, and C and 6 respectively), so that eight dry/wet points are selected. 255

GHF inversion
The two main input parameters that influence the presence or absence of water on the bed are the ice thickness H and the shape parameter p.Both of them are affected by uncertainties, and the corresponding 2D-parameter space is explored with the heat model using a Monte-Carlo method based on 200 (H, p) values spread along a Gaussian distribution.The critical thickness H c measured 260 on the map is used as a prior for H.The shape parameter p cannot be less than −1, and the inversion is therefore done on the derived parameter p : The prior value for p is taken as a gaussian distribution of mean p = 1.5 and standard deviation A kriging interpolation was carried out to evaluate the GHF between the evaluated spots.Unfortunately, they are scarce (ten points only), and unevenly distributed within the study area.As such, the computed experimental variogram was poorly described, and our confidence in the computed kriging standard deviation is low.Moreover, the presence of local variability at a scale of a few tens of kilometers cannot be excluded, so that the validity of such an interpolation is limited.However, we highlight the importance of such a map to show the regional trend that can be expected around Dome C.
Figure 5 shows the interpolation between the ten spots, in which the GHF field smoothly evolves along a N-S gradient, with a typical norm of 0.1 mW m −2 km −1 .At the two sites where the correlation between reflectivity and ice thickness is weak (E4 and E6), the inferred values of the GHF 285 do not mismatch with the N-S gradient, even if the value at E4 seems significantly lower than its neighbours.Given that the kriging standard deviation does not account for the uncertainty on the GHF estimation at the spots, the GHF at the EPICA drilling site is estimated at 54.5 ± 3.5 mW m −2 .
5.2 Forward mode: emulator of the basal melt rate averaged over the last 400 000 years Computing an average value of the basal melt rate is needed to assess the risk of losing the oldest ice 290 layers during the glacial/interglacial periods.To do so, the GHF field inferred at the previous step is now used as an input for the heat model.Running it transient allows the computation of the basal melt rate at each time step, which is done for a large range of H and p values, while α and β are kept at 1/6.04 and 0.0156 respectively.Because of the uncertainty on the initial state, the accuracy of our temperature profile is unclear during the first glacial cycles, and we average the computed melt rate 295 over the last 400 ka only.
The empirical relationship that links the average melt rate to the model parameters appears to be very regular (Fig. 6).The slope of the isomelt contour lines shows a close equivalence between GHF and ice thickness (1 mW m −2 ⇔ 60 m).The sensitivity of the melt rate on p seems to be equivalent respectively, are basically chosen as their expected values under the condition chosen for the other 315 variable.They are computed on the joint distribution given by inversion step as follows, where Φ i g is the GHF at the spot i, and Φi g its central value inferred by the inversion.
The figure 8 shows that the basal melt rate goes up to 1.0 mm a −1 where ice is very thick, and 320 vanishes over several spots.The candidates A, B, C and D appear to be melt-free, which was expected, since they benefit from a relatively lower ice thickness.Less expected is the presence of a potentially melt-free area sixty kilometers north from Dome C (point N8).The ice is quite thick there (∼3300 m), but the computed GHF is low enough to prevent ice from melting.For high GHF values (case a)), the basal melt rate is generally increased by a typical amount of ∼ 0.2 mm a −1 , in compar-

325
ison with case b).Hence, the extend of the melt-free spots is drastically decreased, however they still remain under this high GHF pessimistic assumption.The temporal evolution of the basal melt rate is also shown for two (H c , Φ g ) temperate configurations, leading to the same average value (Fig. 7).
The difference between minimum and maximum value of the basal melt rate is 0.6 mm a −1 for a 2770 m-thick ice, but additionnal 500 m of ice dampens this amplitude by a half, and smoothes out 330 high frequencies.The difference between present melt rates and maximum past values suggest that, more generally, it is possible that present cold basal ice was melting during the warmest periods.

Sensitivity tests
The influence of the climatic forcing is now investigated, by running the model for different slopes of the isotopic thermometer.The results for the GHF obtained by the inverse method are shifted 335 positively by 1.4 mW m −2 for α = 0.20 and negatively for α = 0.13.However, the average basal melt rate is changed at the opposite by 0.1 mm a −1 , so that the final difference in the inferred melt rate is only 0.04 mm a −1 (Table 2).Because the parameter α affects both the derivation of the GHF on the spots and the model melt rate, the two steps partly compensate each other when producing the melt rate map.The accumulations reconstructed at the Dome itself were equally used for the whole region, 345 whereas spatial variations surely affect the surface accumulation rate around Dome C (Frezzotti et al., 2005).The sensitivity of our model to a ±10 % variation of the surface accumulation results in a ±0.3 mW m −2 difference for the GHF, and to an opposite ±0.08 mm a −1 for basal melt rate.
The final sensitivity of the basal melt rate is thus ±0.05 mm a −1 .A higher accumulation rate results in a lower basal melt rate, which is expected.

350
Given that the expression of the pressure melting point is an unusual choice in glaciology (Eq.9), the inverse method was reiterated with a more common value (B = 0.098 K.Pa −1 ) as a test.The inferred GHF values were found to be shifted by −0.6 mW m −2 compared to our expression.The order of magnitude of the results remains the same whatever the expression for the pressure melting point, as well as the regional pattern of GHF. 355

Forward mode: emulator of the critical ice thickness
The method performance is finally assessed by comparing modelled present melting areas to the ones observed by Zirizzotti et al. (2012).The relation between critical ice thickness and GHF is emulated by a one-variable polynom for p = 2: This allows to compute the difference between critical ice thickness and the present ice thickness (Bedmap2), so that negative values correspond to melting areas, and positive values to melt-free areas.We intend to broadly mimick the wet/dry pattern in areas where no critical thickness was measured (north of Dome C), so we tune on purpose the imposed GHF field, within the uncertainty range.The map presented on Fig. 9 was build for Φm g − 1 mWm −2 .Superimposed with the observa-365 tion data, much of the map is quite well assigned, even far from the ten measurement spots.Despite the lack of accuracy of the thickness description, many small-scale structures are well described (e.g.E6-7, G4, G-H8, H6, I9, L7-8).The steeper the bed, the sharper the limit between melting and no-melting areas; this limit appears very smooth in the central part.If this threshold between temperate and cold ice is somewhat offset locally, the spatial pattern of dry/wet areas is however Around Dome C, ice may flow over a hilly bedrock, and the velocity profile is not necessarily as smooth as a 1D synthetical shape function.However, the energy budget at the base of the ice actually depends on the total advection of cold from the top, not on the particular shape of the vertical velocity profile.Using a synthetic profile is just a convenient way to account for a realistic transport of cold towards the bed.
390 Some of the given confidence intervals are quite low (E4, E6, L7 and N8), and this is a consequence of the tiny altitude difference measured on Fig. 1 between the highest wet points and the lowest dry ones for a given spot.Since the correlation between ice thickness and reflectivity was weak, the confidence intervals at E4 and E6 are probably underestimated, and some local effects may not have been accounted for in this study (e.g.small relief and GHF variabilities).

395
We cannot exclude any assignment error of the dry points, if ever a small water film is present but not detected.If so, we would only be able to assess a lower boundary for the GHF.Water remaining at the base enables basal sliding, which can be suspected by observing the slope of the ice surface.
The surface slope is the source of motion, so that no sliding enhances a steeper slope, and sliding, a smaller slope.Most of our measurement spots are standing on slopes locally steeper than the 400 regional slope, so that it is likely that sliding does not occur there, giving an additional hint that the base is cold (Table 1).If the basal ice is anyway temperate despite this favourable clue, the absence of detection of the meltwater likely means that the water film is very thin and the melt rate is very small.The average inferred GHF would be only offset by a small amount, and the regional GHF gradient would be unchanged.

Sensitivity to parameters
The sensitivity of our results to the climatic forcing shows that the confidence interval on the GHF at the 10 spots could be larger than presented so far, possibly by ∼ 2 mW m −2 .However, for the archiving process, the truly important parameter is the basal melt rate, which is much less sensitive.
For example, a lower α correspond to a lower inferred GHF (inverse run), but the average melt 410 rate is higher for the same GHF (direct run).The latter compensates two thirds of the effect of the former, so that the melt rates computed for the ten measurement spots are quite robust to our lack of knowledge on the climatic forcing.
As the present surface accumulation pattern shows a N-S gradient, we wonder if our GHF gradient could come from an artefact of our method, which does not account for the spatial variations of 415 accumulation.As the sensitivity of the GHF on the surface accumulation is rather limited, accounting for its variations would not radically change our results.Furthermore, we do not know how stable is the pattern of accumulation over the glacial/interglacial periods, so that assuming its stability would results in unnecessary additionnal uncertainties.where they find a positive heat flux anomaly and a high melt rate on its flanks (∼ 2.5 mm a −1 ).

425
The heat anomaly cannot be induced with our method, since the GHF is interpolated at this place.
This illustrates that our method is reasonably reliable where we do have an estimation of the critical thickness (corresponding to interesting dry, low-reflective cold spots), and at least describes a GHF regional background elsewhere.
At the EPICA drilling site, we compute an average melt rate of 0.32 ± 0.25 mm a −1 , and this 430 value seems to be slightly lower than, but however consistent with the value of 0.56 ± 0.19 mm a −1 already inferred by Parrenin et al. (2007), which was computed over longer timescales.The range of possible basal melt rates inferred from our method seems realistic enough to contain its effective local value, at a place where there is no observation of the critical thickness yet.The basal melt rate inferred for the different old-ice candidates, where we have observations of H c , is thus reliable.The In our point-by-point method, there is no horizontal regularization constraints to make the GHF being realistic at a regional scale.Yet, the inferred heat flow at the ten spots shows a certain spatial correlation, which is a good omen of plausibility.More fundamentally, we can explain at first order a relatively complex pattern of wet/dry areas on a whole region with a single physical key (a smoothly-445 varying GHF field) without any 2D water routing model, or horizontal heat transport depending on both ice flow and bedrock topography.This application of the parcimony principle of Occam's razor supports the validity of our 1D approach, and means that the main physical mechanisms have been accounted for.Of course, a higher-dimension ice flow and hydrological model, over a refined bedrock, would be necessary to go one step beyond to describe the water variability more closely, 450 but more assumptions regarding the model parameters would also lead to higher model uncertainties.
The GHF confidence intervals given for the surroundings of Dome C are now about five to ten times lower than the ones previously available.In the frame of the oldest-ice project, we suggest that Φm g + 2σ Φ is a realistic higher boundary for the local GHF value.

455
The basal melt rate map (averaged over the last 400 000 years) shows that the old-ice candidates appear to be melt-free all along glacial/interglacial periods.On the contrary, the northern part of Dome C is broadly a thawing area, whereas Zirizzotti et al. (2012) observed no presence of water, in spite of the local thick ice.In this region, the climatic signal of the last deglaciation has reached the bedrock later than in regions of thinner ice (Fig 9, contour lines).If the basal ice was cold at the LGM, 460 the pressure melting point has not been reached yet since the deglaciation signal has just reached the bed for the last ∼ 2 000 years.For thinner ice, the deglaciation signal is close to be already past; if the basal ice is still cold now, it should be the case over the whole glacial/interglacial periods as well.
The time lag for the climatic signal to reach the bed then appears to be a self-sufficient parameter to detect cold-based areas on long time scales.However, this only criteria would not be sufficient to 465 detect all of these areas, or to discuss their spatial extent.

Consequences for the old-ice targets
As expected, the regions where the basal ice is supposed to be cold over the glacial/interglacial periods are the candidate A and the candidates B, C and D. These three last places may undergo a lower GHF than candidate A and could be better places to exclude the possibility of basal melting. is not significant enough to make candidates B, C and D first choices, but they remain spots of great interest.Three-dimensional modelling will be performed to study the regional ice flow towards these sites in more details.
Our study also suggests that a certain spot in the northern part of Dome C may undergo a low enough GHF to prohibit basal melt rate over long time scales.As the ice is significantly thicker 480 than at other previously considered candidates, a very old ice could be retrieved with a much better resolution than elsewhere.Considering that the clue leading to this conclusion is not cross-checked (one observation only), we emphasize the importance of carrying out additional survey to assess the validity of this statement (e.g.ground/airborne radar).
Finally, the evolution of the benthic δ 18 O past sea level proxy (Lisiecki and Raymo, 2005), and 485 its close similarity with the one of the EDC reconstructed temperatures, show that the mean air temperature at Dome C was probably higher by ∼ 2 K before −800 ka than after.As a consequence, the mean basal melt rate was probably higher by 0.1 mm a −1 as well.Even if we conlude this study on favorable statements, we should bear in mind that basal ice has maybe undergone little melting on tops and flanks of bed reliefs during the Mid-Pleistocene.

7 Conclusions
The geothermal heat flux is a poorly constrained geophysical parameter in Antarctica, despite its crucial influence on the ice flow properties and old-ice archiving.Around Dome C, the available continental estimations are currently of limited benefit, and we lack a more precise local estimation.
We have presented a simple inverse method based on a 1D heat model, constrained by amplitude 495 analysis of RES echoes recorded at ice/bedrock interface trying to discriminate wet and dry areas on the bed.Assuming that the GHF is locally uniform, the presence of basal water is only a function of ice thickness.The critical ice thickness, for which the pressure melting point is reached today, is inferred from wet/dry thresholds used in the analysis of RES amplitude data (Zirizzotti et al., 2012).The heat model, accounting for the whole history of the ice (temperature, accumulation and 500 thickness evolution), is inverted for this critical thickness, to retrieve the value of the GHF and of the time-averaged basal melt rate for ten spots around Dome C.
Our method is valid in dome areas, where horizontal advection and diffusion can be neglected, but its principle could be adapted for other regions with a more complex physical scheme.However, it assumes that the origin of the observed basal water is local, so that it is better suited for flat regions,

505
where there is no upwards water transport on the bed reliefs.Moreover, in places where horizontal ice flow is significant, deformational heat should be accounted for in the energy balance of ice.
Furthermore, we show that the role of the ice thickness appears to be dual.Of course, in average, it limits the heat diffusion towards the surface and enhances basal melting.But under an evolving The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.
climate, it lags the temperature forcing, so that the base of a thin ice is more sensitive to climate 510 evolution, while the base of a thick ice has just begun to be under deglaciation and may still be cold today.Hence, as the LGM was one of the coldest climatic conditions ever recorded, the lag effect of the ice thickness must absolutely be accounted for to correctly interpret present basal conditions.
The interpolated map of the GHF shows a regional gradient, oriented North-South.Where no critical thicknesses have been measured, the GHF values are consistent with the pattern of dry and 130 conditions at the bottom At the ice/bed interface, the thermal conditions depend on whether the pressure melting point is reached or not.For thawing ice, the temperature simply equals the melting temperature.The melting temperature of pure ice linearly depends on the pressure following a Clapeyron law, for which the corresponding coefficient is B = 0.074 K Pa −1 .For glacier ice, a coefficient of 0.098 K Pa −1 The CryosphereDiscuss., doi:10.5194/tc-2017Discuss., doi:10.5194/tc--23, 2017     Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.

The
Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.denoted by the indexes of their central point.The correlation between ice thickness and reflectivity is shown on Fig. 4. For two spots (E4 and E6), the correlation is weak, but we retain them as they 245 are the only ones available in the central part of our study area.
σ p = 0.3 so that the corresponding values for p are mainly distributed between common values of 1 265 and 10.For a given present-day critical ice thickness H c and a given p, the 1D heat model is run for increasing values of Φ g (by steps of 0.25 mW m −2 ), and the first Φ g value that allows melting for the present time is selected as the local GHF.The thicker the ice, the lower the GHF needed to allow for melting.The explored value range for Φ g is 40-70 mW m −2 .mode: geothermal heat fluxThe derived values of the mean GHF for the ten measurement spots are within a range of 48-60 mW m −2 (Table1), and the highest values of the GHF are found south of Dome C. Given the model assumptions, the inferred value for the GHF at two potential drilling sites (C6 and H1) are respectively 59.3 ± 2.2 and 53.9 ± 3.3 mW m −2 .275 9 The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.
to a range of 2 mW m −2 of GHF.An additional flux of 1 mW m −2 correspond to an increased melt 300 rate of 0.09 mm a −1 .To account for this relationship by a multivariable polynomial, it seems natural to make m depend linearly on Φ g , and quadratically on H and p.Over the positive-melt-value domain, we computed by a least-square minimization method the following relation:m = − 5.148 × 10 −7 H 2 + 4.688 × 10 −3 H + 89.519Φ g + 3.08 × 10 −3 p 2 − 5.887 × 10 −2 p − 14.335.(14)305The performance of the emulator is assessed by comparing the melt rate output from the thermal model with the one computed with the emulator.The average absolute error is 0.014 mm a −1 , and the corresponding standard deviation is 0.016 mm a −1 , so that the errors due to the emulator are The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.significantely lower than the corresponding uncertainties due to the GHF and p.This polynomial function is thus considered precise enough for our purpose.310 We now call Φg the GHF field and Φm g the one based on the average values derived at the spots.The averaged basal melt rate in the past was inferred for two pessimistic cases: a) a high value of the GHF ( Φm g + 2σ Φ ), and b) a low value of p (p − 2σ p ), leading to little advection of cold ice.As Φ g and p are not independent within the inversion process, the values of p and Φg , for a) and b) ., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.For extreme values of the β coefficient (accumulation model), the results are shifted positively by 0.3 mW m −2 for β = 0.0138 and negatively for β = 0.175.Similarly to α, the average basal melt rate is changed at the opposite by 0.01 mm a −1 , so that the final difference in the inferred melt rate is only 0.02 mm a −1 .

370
often respected in Fig.9(e.g.G-H8, G9, J6).Given the uncertainties, only one spot appears to be undoubtly inconsistent (F8), but this is explained by a difference of vertical accuracy of the data (local gap of ∼ 100 m between Bedmap2 and the radar data).The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.The heat model allow an easy detection at the base of the ice of the very strong climatic signal of the LGM, and the corresponding time lag for basal temperature appears to linearly depend on the ice 375 thickness for the range 2700-3500 m (Fig. 9) by the following relation ∆t = 10H − 16610.(18)Theisocontours of the time lag ∆t are simply parallel to the thickness contours, but allow a temporal interpretation of the wet/dry patterns.Around 10 ka are required for the signal to reach the tops of bed reliefs, but almost twice this value for a 3400 m-thick ice.Hence, considering the 380 duration of a deglaciation, the thermal state of the basal ice corresponds to very different climatic periods depending on the overburden ice thickness.

6. 1 . 3
Consistency with published data and measurements 420 Our GHF values are comparable to the ones inferred by Fox Maule et al. (2005) and Purucker (2013), but, as their uncertainties are quite large, they cannot be considered as a litmus test.For the northern part of our domain, the relatively low values of GHF (∼ 50 mW m −2 or less) is closely consistent with the estimation of Carter et al. (2009) at the same place, except on a bed relief (candidate D),

435
ice temperature gradient measured in the borehole is −0.022655K m −1 at the base, corresponding to a net heat flux of 48 mW m −2 going into the ice.The present corresponding amount of melting for our interpolated GHF is 0.7 mm a −1 , which seems consistent with the time-averaged value ofParrenin et al. (2007) regarding the fact that the deglaciation signal has already reached the bedrock for 4 ka.440 14 The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.6.1.4Structure of the GHF field

470
However, they are separated from the topographic dome by the Concordia trench.The deep ice that could be drilled at these places have likely crossed the trench, and the stratigraphy may have been affected by this topographically-disturbed region.Furthermore, the ice velocity is probably higher than over the candidate A because of the steeper slope starting from the dome, and this may The CryosphereDiscuss., doi:10.5194/tc-2017Discuss., doi:10.5194/tc--23, 2017     Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.enhance basal disturbance.Given the confidence intervals, the GHF difference with candidate A 475

515
wet points, in particular in the northern region, which appears to be dry today, and hosts a site of old ice potential.All the previously considered old-ice candidates are very likely cold-based, or to have undergone very little basal melting.We hope this work will provide the community with helpful and realistic thermal boundary conditions for any ice flow modelling in the Dome C region.The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.The Cryosphere Discuss., doi:10.5194/tc-2017-23,2017 Manuscript under review for journal The Cryosphere Discussion started: 14 March 2017 c Author(s) 2017.CC-BY 3.0 License.

Figure 1 .
Figure 1.Wet and dry points (resp.blue and red), adapted from Zirizzotti et al. (2012).The yellow star shows the EPICA drill site, distant from the topographic dome by 1.4 km.Side numbers and letters are used to locate any intersection point on the radar grid.Large-size letters stand for the main sites of interest (candidates A, B, C and D).Projection: WGS84/Antarctic Polar Stereographic -EPSG:3031 (m).The north direction is bottom right.

Figure 2 .Figure 3 .Figure 4 .
Figure 2.This figure illustrates how water and basal topography are linked in a region of flat surface for a uniform local GHF.The blue line shows the presence of water at the bed, and the direction of its flow.The brown area corresponds to a strip of unknown thermal conditions.

Figure 7 .
Figure 7. Evolution of the basal melt rate, depending on ice thickness and geothermal heat flux.

Table 2 .
Sensitivity of the GHF (inverse runs, [mW m −2 ]) and basal melt rate (forward runs, [mm a