Geothermal flux beneath the Antarctic Ice Sheet derived from measured temperature profiles in deep boreholes

The temperature at the Antarctic ice sheet bed and the temperature gradient in subglacial rocks have been directly measured only a few times, although extensive thermodynamic modelling has been used to estimate geothermal heat flux under ice sheet. During the last five decades, deep ice-core drilling projects at six sites – Byrd, WAIS Divide, Dome C, Kohnen, Dome F, and Vostok – have succeeded in reaching to, or nearly to, the bed in inland locations in Antarctica. When 15 temperature profiles in these boreholes and heat flow model are combined with estimations of vertical velocity, the heat flow at ice sheet base is translated to a geothermal heat flux of 117.8±3.3 mW m-2 at Byrd, 67.3±8.6 mW m-2 at Dome C, 79.0±5.0 mW m-2 at Dome F, and -3.3±5.6 mW m-2 at Vostok, close to predicted values. However, estimations at Kohnen and WAIS Divide gave flux of 161.5±10.2 mW m-2 and 251.3±24.1 mW m-2, respectively, far higher than that predicted by existing heat flow models. The question arises as to whether this high heat flow represents regional values, or if the Kohnen and 20 WAIS Divide boreholes were drilled over local hot spots.

. For a central point in the WAIS (78S, 110W), the average GHF is expected to be 110 mW m -2 . The GHF can also be estimated on the basis of geological information, where uniform values are attributed to large geologically homogeneous areas (An et al., 2015;Goodge, 2018;Llubes et al., 2006;Martos et al., 2015;Pollard et al., 2005).
Some studies suggest using remote methods to estimate GHF underneath the Antarctic ice sheet. For example, satellite 35 magnetic data showed that the GHF underneath the ice sheet varies from 40 to 185 mW m -2 and that areas of high GHF coincide with known current volcanism and some areas known to have ice streams (Fox Maule et al., 2005). In the central part of the EAIS, the average GHF was estimated to be in the range of 50 to 60 mW m -2 ; however, elevated GHFs were found along the WAIS-EAIS boundary and around the Siple Coast. Similarly, high GHFs were found around Victoria Land, Oates Land, and George V Land. Observations of crustal heat production within the continental crust underneath the 40 Lambert-Amery glacial system in East Antarctica also show high heat flux of at least 120 mW m −2 (Pittard et al., 2016).
Direct temperature measurement obviously produces the most reliable GHF estimations and can be used to verify results of preliminary thermal modeling and geological-geophysical studies. While over 10,000 heat flow measurements have been made globally, 90% are from Europe, North America, and southern Africa. South America, Asia and Australia have far fewer measurements, and Antarctica has virtually none (Neumann et al., 2000). Because drilling through the thick ice is 45 extremely complicated, time-consuming, and expensive, direct temperature measurements in Antarctic subglacial till/bedrock environments have only been conducted twice so far, both under the WAIS: at the subglacial Lake Whillans (285±80 mW m -2 ) (Fisher et al., 2015) and near the grounding zone of the Whillans Ice Stream (88±7 mW m -2 ) (Begeman et al., 2017), ~100 km apart (Fig. 1). The tremendous difference between these two adjacent sites suggests a high spatial variability of GHF. 50 More reliable GHF estimations under the Antarctic ice sheet can be made from available temperature profiles in ice boreholes. During the last five decades, deep ice-core drilling projects at six sites -Byrd (Ueda and Garfield, 1970), WAIS Divide (Slawny et al., 2014), Dome C (Augustin et al., 2007, Kohnen (Wilhelms et al., 2014), Dome F (Motoyama, 2007, and Vostok (Lukin and Vasiliev, 2014;Vasiliev et al., 2011) -have succeeded in reaching, or nearly reaching, the ice sheet bed at inland locations in Antarctica. Reported drill site conditions -snow accumulation rate, mean annual surface air 55 temperature, ice thickness, and drilling depth -are summarized in Table 1.
The Byrd and Kohnen holes encountered water at the base, which welled up into the holes. The borehole at Vostok penetrated the subglacial Lake Vostok at 3769.3 m, and here, water rose from the lake to a height of more than 340 m.
Drilling of other holes was stopped within 10-50 m from the bed. All these holes were temperature logged and provide a good opportunity to fill the gap in our knowledge of the GHF under the Antarctic ice.

Temperature and temperature gradient at the base of Antarctic Ice Sheet
Temperature in Byrd, WAIS Divide, Vostok, Dome C, Kohnen, and Dome F boreholes were measured using different devices and different methods. All boreholes were mechanically drilled and filled with kerosene-based drilling fluid.
Temperature profiles were then obtained by logging with custom-made borehole temperature loggers (Dome C, Kohnen, 65 Dome F, WAIS Divide, and Vostok) or a thermistor embedded into the drill (Byrd).
Temperature measurements in mechanically drilled boreholes filled with drilling fluids are difficult because the increasing temperature with depth can cause convective cell mixing and because of the long time required to get the sensor in equilibrium with the surrounding ice temperature (Clow, 2008). Despite this, we believe that measurement precision in all holes was better than 0.05 °C, and in many instances, did not exceed 0.01 °C. All depth data were recalculated to true 70 vertical depth.
Measured temperature profiles in four of the boreholes (Vostok, Dome C, Kohnen, and Dome F) increase nearly linearly with depth, as would be expected in locations with minimal accumulation and hence small vertical velocities (Fig. 2).
Vertical advection is much greater at the Byrd and WAIS Divide boreholes in West Antarctica; in these locations the upper part of the ice sheet is nearly isothermal, but at depth the temperature gradient is nearly the same as at the other sites. 75 Temperature gradients at the bed are 2.02-3.12 °C/100 m at Dome C, Kohnen, Dome F and Vostok and slightly higher in West Antarctica, 3.70-3.88 °C/100 m at Byrd and WAIS Divide (Table 2).
Ice thicknesses generated from extrapolation of temperature profile to the depth of pressure melting point in assumption of Clausius-Clapeyron slope of 0.0742 K/MPa (Cuffey and Paterson, 2010) are in good agreement with radar data, except for the WAIS Divide, where the difference is ~30 m. This could be caused by scintillations on the melted ice-bedrock interface. 80

GHF estimation model
A one-dimensional time-dependent energy-balance equation (Dahl-Jensen et al., 2003;Johnsen et al., 1995) is usually used to model the temperature distribution through the ice as a function of the climate conditions on the surface and the GHF from the bedrock: (1) 85 where T is the temperature as a function of z, x and t; z represents the vertical coordinate (at the ice sheet base, z=0, while at ice sheet surface, z=H); H is the ice thickness; x is the horizontal coordinate; t is the time; k is the conductivity of ice dependent on T; r is the density of ice; c is the specific heat capacity of ice dependent on T; w and are respectively the vertical velocity and horizontal velocity of ice sheet dependent on z and t.
We assume that ice sheet is always in thermal steady state and horizontal advection can be neglected (Robin, 1955). This 90 assumption reduces the non-steady-state heat-transfer equation to which can be rewritten as The vertical velocity in the ice is given by (Fischer et al., 2013) where ;<=% is the basal melt rate; Acc is the surface accumulation rate dependent on t; m is the form factor that accounts variation of horizontal velocity. 100 Classically, vertical velocity linearly depends on / (Cuffey and Paterson, 2010) and = 0. However, at an ice divide, the downward flow of ice is slower, for the same depth, than at locations away from the divide (Raymond, 1983). This reduces the cooling influence of vertical advection and increases the basal temperature. Within this near-divide zone, the form factor could be from 0.5 (Fischer et al., 2013) to 1.0 (Raymond, 1983). All discussed sites are located at, or near, ice divide, thus, we assume m = 1. 105 Substitution of Eq. (5) into Eq. (4) and integrating on the assumption of k=const gives the following temperature distribution in ice sheet at steady state: where F is the surface temperature; G #$ #' H I is the temperature gradient at ice sheet base; $ = / is the thermal diffusivity of ice. 110 Least square method is used to fit measured borehole temperatures with this equation. In fitting, the initial values of and ;<=% can only be guessed and this results in unavoidable uncertainty of fitting. To over this large uncertainty, a common genetic algorithm (GA) was designed to find optimal global solution of temperature fitting by limiting these unknown parameters changing in a predetermined range (Reeves and Rowe, 2002). In the GA, the crossover fraction is set to be 0.9 while the migration fraction is 0.2. To obtain accurate solution and save calculation time, 115 we set the population size to be 8000 and the generations to be 20. Usually, after 15 generations of iteration, the optimal solution can be found. All the calculations were performed with MATLAB R2014b software. For each deep borehole, the fitting experiments were trialed five times to avoid random error of GA. Then, the average value in the five fitting experiments was used as the GHF from bedrock into ice at selected site and the uncertainty ranges came from the difference between the maximum/minimum and the average GHF values.
Note that the vertical velocity is markedly affected by energy which is used to melt/freeze ice = ;<=% . Thus, the GHF will be: 125 where L is the specific latent heat for melting of ice.
In our model, we assume that k and c are constant and equal to values at the temperature of pressure melting point that can give better estimation of basal melting rate at the base of ice sheet (Fischer et al., 2013). In this case, 1 2 o2 o$ =-5.7e -3 K -1 . 135 Estimated vertical velocity profiles are shown on Fig.3. Fig. 4 shows the fitted temperature profiles compared with measured temperatures. Results of our estimates for basal melt/freezing rate and GHF are summarized in Table 2.
The temperature profile shows that the heat flows through the ice at six deep drilling sites in Antarctica must be >42.6-81.4 mW m -2 in order to match the observed temperature in the boreholes. The basal ice at all sites is at the pressure-melting point, and the amount of melt cannot be constrained by the energy balance equation alone. When the heat flow model is 140 combined with estimations of vertical velocity, the estimated heat flow can be translated to a GHF of 67.3-251.3 mW m -2 , far higher than that predicted by existing heat flow models.
Byrd. Our modeled GHF value at this location (117.8±3.3 mW m -2 ) is much higher than the first estimation made immediately after temperature logging (75.4 mW m -2 referenced by Ueda, 2007) simply because the latter one did not account for the basalt ice melt. The latest modeling (Martos et al., 2017) revealed high GHF at the location of Byrd Station 145 (132 mW m -2 with an error of ±5 mW m -2 ), compared with reference to previous models (An et al., 2015;Fox Maule et al., 2005;Van Liefferinge Pattyn, 2013). This is sufficiently close to our estimations.

Dome C.
The inverse approach to retrieving GHF from radar inferred distribution of wet and dry beds at the EPICA drilling site (Passalacqua et al., 2017) gave 54.5±3.5 mW m -2 , lower than estimates derived from borehole temperature profiling (67.3±8.6 mW m -2 ). Modeled GHF range (43-56 mW m -2 obtained by An et al., 2015;Fox Maule et al., 2005;Martos et al., 2017;Shapiro and Ritzwoller, 2004;Van Liefferinge Pattyn, 2013) is also 1.2-1.6 times lower than our estimates. However, a high spatial variation of GHF at Dome C area was found from radar-sounding data (Carter et al., 2009). The values of nearly 100 mW m -2 inferred for the southern shore of Concordia Subglacial Lake, approximately 50 km to the south of the drilling site, are also well outside modeled estimates.
Kohnen. The model with a standard GHF of 54.6 mW m −2 predicted a basal temperature 0.3 °C below the pressure melting 155 point at Kohnen (Huybrechtset al., 2007). However, subglacial water entering the borehole indicated that the actual GHF should be much higher than that assumed by the standard model. Under these circumstances, our estimates (161.5±10.2 mW m -2 ) approach to a certain degree the real heat flux.

Dome F.
A previously estimated GHF of 59 mW m -2 neglected the bottom ice melt rate (Hondoh et al., 2002) and thus is lower than our estimates (79.0±5 mW m -2 ). As the drill approached the base (approx. 10 m above), subglacial meltwater 160 leaked into the borehole and froze onto the drill, directly indicating that ice reaches the pressure melting point under high

GHF.
Vostok. The surface temperature-time curve for the upper bound of the present-day accumulation rate at Vostok corresponds to a GHF of 53 mW m -2 (Salamatin et al., 1998). We calculated that at the base of the ice sheet, the conductive flux is 42.6±0.4 mW m -2 while the latent heat flux from refreezing of the lake water is 46.3±5.6 mW m -2 . Thus, the GHF heat flux 165 at the base of the ice sheet has a negative value of -3.3±5.3 mW m -2 . This is in good agreement with the isotope studies that showed that the Vostok ice core consists of ice refrozen from Lake Vostok water, from 3539 m below the surface of the Antarctic ice sheet to its bottom (Jouzel et al., 1999).
At this stage we are not yet able to predict GHF at the bed of 600-m thick subglacial Lake Vostok because temperature profile in the lake is still indefinite. However, the DNA detection of thermophile bacterium in the near-base accretion ice 170 suggests the existence of near-bottom warm waters with temperatures as high as 50 °С (Bulat et al., 2012). If so, the GHF in the lake sediments can reach 200-240 mW m -2 . These values can be considered as paleo-GHF because microorganisms were picked up thousands of years ago but still actual accounting for long duration of geological processes.

WAIS Divide.
Preliminary estimations of GHF showed quite high values of 215±13 mW m -2 , depending on the actual ice thickness (Clow et al., 2012). Our estimates are even higher at 251.3±24.1 mW m -2 . The question arises as to whether this 175 high heat flow represents a regional value, or if the WAIS Divide borehole was drilled over a local hot spot. There is no depression in the local surface topography or drawdown in the subsurface layers detected by ice-penetrating radar, as would be expected over a local hot spot. Thus, the high heat flow at this location appears to represent a regional value. How far this heat flow high extends into the interior of the West Antarctic Rift System will require further work.
Prediction of the future behavior of the Antarctic Ice Sheet requires accurate ice-sheet models. However, GHF models based 180 on seismic tomography, radar data, magnetic field observations, the tectonic age and geological structure of the bedrock yielded mixed results at sites of deep ice-core drilling in Antarctica. Processing of available temperature profiles in ice boreholes showed that Antarctica is warmer than predicted. The GHF found at the sites of the EAIS lying on young oceanic crust are 1.1-2.5 times higher than the mean continental GHF. Currently, this phenomenon is also confirmed by tectonic https://doi.org/10.5194/tc-2020-32 Preprint. Discussion started: 4 March 2020 c Author(s) 2020. CC BY 4.0 License.
reconstructions (Carson et al., 2014). However, this warm base could be caused by radiogenic or tectonically induced heat 185 effects, and a coherent explanation for this phenomenon is still required. Most numerical models of the EAIS basal conditions assume the GHF to be 42-65 mW m -2 . However, the presence of basal meltwater beneath most of the Antarctic ice sheet requires GHF ≥ 80 mW m -2 (Budd et al., 1984), which is in good correlation with estimates from the deep drilling sites.
Variability of crustal thickness, hydrothermal circulation (Seroussi et al., 2017), magmatic intrusion (Van Wyk De Vries et 190 al., 2017, and thermal conductivity variability are the main contributors to the elevated and highly variable GHF in West Antarctica (Begeman et al., 2017). One of the first pieces of evidence for an "unreasonably high" GHF (>100 mW m -2 ) under WAIS came from temperature-depth profiles in a 480-m-deep borehole drilled at Crary Ice Rise, Antarctica (Bindschadler et al., 1990). The GHF at the vents of subglacial volcanoes in West Antarctica can be as high as 25 W m -2 and one such 23-km-wide caldera was revealed ~100 km to the south of the WAIS Divide drill site (Blankenship et al., 1993). 195 However, more precise GHF estimates and explanations for an elevated GHF would be possible after temperature logging and subglacial rocks studies from deep boreholes that are required to drill in Antarctica in the distant future.  1966-1968a 2006-2011d 1990-1998, 2005-2014f,g 1999-2004i 2002-2006k 2003-2007n Surface