Geothermal heat flux of Ridge B region in Antarctica inferred from basal dry–wet distribution

Abstract

Ridge B is one of the least studied areas in Antarctica but has been considered to be a potential location for the oldest ice on Earth. Among important parameters for calculating where very old ice may exist, geothermal heat flux (GHF) is critical but poorly understood. Here, GHF is determined by quantifying the transitions between dry and wet basal conditions using a radioglaciological method applied to airborne radio-echo sounding data. GHF is then constrained by a thermodynamic model matched to the transitions. The results show that GHF in Ridge B varies locally and ranges from 48.5 to 65.1 mW m⁻², with an average value of 58.0 mW m⁻², which is consistent with the current known GHF constrained by subglacial lakes and derived from Vostok ice core temperature measurements. Our work highlights the value of considering local GHF when locating the oldest ice in this potential region or other regions.

1. Introduction

Finding a 1.5 million-year-old ice core is the key to resolving the mechanisms behind the major climate reorganization during the Mid-Pleistocene Transition (Van Liefferinge and others, 2018). Investigations of ‘old ice’ have been carried out in Dome Fuji, Vostok Station, Dome C and Titan Dome (Petit and others, 1999; Watanabe and others, 2003; EPICA community members, 2004; Karlsson and others, 2018; Beem and others, 2021). Ridge B is part of the main ice divide in the East Antarctic Ice Sheet (EAIS). Ice thickness varies from 2000 to 4000 m, which coupled with the very low accumulation of ice (Siegert, 2003; Leysinger Vieli and others, 2011) marks the region as having potential for containing the oldest ice in Antarctica (Van Liefferinge and Pattyn, 2013; Lipenkov and others, 2019; Cui and others, 2020a; Ekaykin and others, 2021). However, Ridge B is also one of the most underexplored areas in Antarctica, which has led to few glaciological assessments of this potential.

Geothermal heat flux (GHF) is key to predicting where the oldest ice may exist, as it is an important boundary condition of ice flow models (Larour and others, 2012; Golledge and others, 2014; Pittard and others, 2016; Seroussi and others, 2017; Reading and others, 2022). GHF can affect the ice-sheet behavior by controlling the freezing and melting of the ice-sheet bed; a process dominant in areas with low ice velocity such as Ridge B (Fahnestock and others, 2001; Joughin and others, 2009; Larour and others, 2012; Pittard and others, 2016). Ice flow models that predict basal temperatures can offer insights into likely locations of old ice, but these lack precision at Ridge B (Wolff, 2005; Brook and others, 2006; Van Liefferinge and Pattyn, 2013; Burton-Johnson and others, 2020) partly because GHF is poorly constrained. Direct and accurate measurements of GHF require deep boreholes which are rare in Antarctica (Carson and others, 2014; Burton-Johnson and others, 2020; Reading and others, 2022), hence an alternative means to evaluate GHF is needed.

At present, several GHF models covering Ridge B have been proposed based on seismological data and/or satellite/airborne magnetic data, most of which provide low spatial resolution estimations and show great differences (Shapiro and Ritzwoller, 2004; Maule and others, 2005; Purucker, 2012; An and others, 2015; Martos and others, 2017; Shen and others, 2020; Li and others, 2021; Lösing and Ebbing, 2021; Stål and others, 2021; Haeger and others, 2022). Martos and others (2017) used airborne magnetic data to obtain a high-resolution GHF model, but the data used have gaps across Ridge B (Golynsky, 2001; Martos and others, 2017; Golynsky and others, 2018). Li and others (2021) used new airborne magnetic data to infer the GHF at Ridge B, obtaining results that are significantly higher than all previous datasets. Siegert and Dowdeswell (1996) estimated the minimum GHF implied by the then-known subglacial lakes in Ridge B by assuming the ice base in the lakes is at the pressure melting temperature. Given known values of ice accumulation and ice thickness, GHF can be calculated through a simple thermodynamic model. Here, we estimate the GHF of Ridge B based on the airborne radio-echo sounding (RES) data collected by the International Collaborative Exploration of the Cryosphere by Airborne Profiling in Prince Elizabeth Land (ICECAP/PEL) project (Cui and others, 2018, 2020b). We report the new limit of GHF under the ice sheet in the Ridge B area and analyze the GHF anomaly. Based on the latest airborne RES data, we use an improved radioglaciological method (Lang and others, 2022) to detect the locations of pressure melting point (PMP) at the ice-sheet bed, and diagnose the distribution of subglacial dry and wet zones. Using this knowledge we then use a thermodynamic model to extract GHF.

2. Data

2.1 Airborne RES data

Since 2015, ICECAP/PEL has surveyed the largest ‘data gap’ in Antarctica with the Chinese fixed-wing airborne platform ‘Snow Eagle 601’ (Cui and others, 2018, 2020b). Snow Eagle is equipped with a phase-coherent RES system, operates at a central frequency of 60 MHz and a peak power of 8 kW, making it capable of penetrating deep ice (>4 km) in Antarctica (Cui and others, 2020a). This study uses data from ICECAP/PEL collected during the 32nd, 35th and 36th Chinese National Antarctic Research Expedition (CHINARE) (2015/16, 2018/19 and 2019/20) (Fig. 1), processed by 2-D focused synthetic aperture radar processing algorithm (Peters and others, 2007). The study area covers the central region of Ridge B where ice divides converge.

Figure 1. Distribution of the survey lines used in this study. (a) Survey lines' distribution with ice flow velocity map (Rignot and others, 2011) as background, the aerial view in the upper left corner shows the location of the study area in Antarctica, using the MODIS Mosaic of Antarctica image map as the background. (b) Survey lines' distribution with surface elevation map (Fretwell and others, 2013) as background. (c) Survey lines' distribution with bed elevation map (Morlighem, 2020; Morlighem and others, 2020) as background. The study region is marked by purple box. The gray lines represent the ice divides (Creyts and others, 2014), and the blue triangle represents the known subglacial lakes (Livingstone and others, 2022). The survey lines are marked with black lines, which are part of the airborne RES data collected by the ICECAP/PEL project in seasons of 2015/16, 2018/19 and 2019/20. The white dashed line segment in (a) is the location of the survey line shown in Section 3.1 for the example of the subglacial dry–wet distribution.

2.2 Surface temperature and accumulation rate

The data used in this study include the annual average surface temperature and accumulation rate, which are from the latest version of a regional atmospheric climate model that is specifically adapted for using over Antarctica (RACMO2.3p2/ANT) (van Wessem and others, 2018). RACMO2.3 is a regional climate model developed by the Institute for Marine and Atmospheric Research Utrecht at Utrecht University, which combines the dynamical core of the High Resolution Limited Area Model version 6.3.7 and cycle CY33r1 with the physics package of the European Centre for Medium-Range Weather Forecasts Integrated Forecast System (van Wessem and others, 2014). RACMO2.3 can provide products at a horizontal resolution of 27 km and a vertical resolution of 40 levels (van Wessem and others, 2014). The accumulation rate and surface temperature used in this study were obtained by taking the average of monthly accumulation rate data and surface temperature data every 3 h over 40 years (1979–2019). A bilinear interpolation algorithm is used to interpolate the values needed across the region.

3. Inversion method of geothermal heat flux

3.1 Diagnosis of basal conditions

RES can be used to infer basal conditions and identify subglacial lakes on a regional scale, since the presence of water at the ice-bed interface is responsible for a remarkable increase in the amplitude of the reflected echoes (Siegert and others, 2005; Zirizzotti and others, 2010; Fujita and others, 2012). Lang and others (2022) produced a method to automatically identify the dry–wet transition zone (DWTZ), from which the dry and wet zones can be determined. In ice divide regions, by assuming that the rate of liquid water generation by basal melting in the upstream zone is greater than the discharge rate, and that the ice sheet is in thermal equilibrium, the DWTZ can be used to detect the locations where the bed is at PMP (Passalacqua and others, 2017; Lang and others, 2022). The detailed introduction of the method of Lang and others (2022) is as follows.

The reflectivity variation profile of the basal interface can be generated based on layer information of the surface and bed. By ignoring the transmission loss caused by multiple reflections between internal layers, the reflectivity variation of the basal interface can be expressed as follows:

(1)$$\Delta R = P_{{\rm ib}}-P_{{\rm ai}} + L_{{\rm Ga}}-L_{{\rm Gb}} + R_{{\rm ai}}-2C_{{\rm ai}}-R_{{\rm ir}}-2L_{\rm i}$$

where the subscripts representing different interface materials as follows: b represents the bed with an unknown condition, a represents air, i represents ice and r represents bedrock. In addition, P represents the reflected power, L _G represents the geometric spreading loss, C represents the one-way transmission loss at an interface, R represents the reflection loss at an interface, L _i is the ice absorption loss, and ΔR represents the variation in reflectivity of basal interface relative to frozen bedrock. The specific calculation method of P is given in Lang and others (2022). According to Eqn (1), ΔR corresponding to each azimuthal sample of the bed can be calculated. For the generated ΔR profile, the theoretical values of the threshold used to identify dry and wet locations in ΔR profile at this time are ΔR _wet = R _iw − R _ir and ΔR _dry = 0, respectively, where w represents water.

However, the theoretical value of the calculated echo power loss terms including transmission loss C and reflection loss R may have a regional error in the actual environment of the Antarctic ice sheet, and losses in the process of echo transmission may not be fully estimated; both of these conditions will lead to errors of ΔR obtained by taking these terms as inputs. The subglacial water bodies in the region can be used as reference to correct the identification threshold. The average value of ΔR profile of subglacial water body is calculated by:

(2)$$\Delta R_{{\rm wet}\hbox{-}{\rm c}} = \displaystyle{1 \over {n_2-n_1 + 1}}\sum\limits_{i = n_1}^{n_2} {\Delta R( i ) } $$

where ΔR _wet-c (dB) is the corrected wet threshold, n ₁ and n ₂ determine the range of subglacial water body. For the case of multiple subglacial water bodies in the region:

(3)$$\Delta R_{{\rm wet}\hbox{-}{\rm c}} = \displaystyle{1 \over m}\sum\limits_{\,j = 1}^m {\Delta R_{{\rm wet}\hbox{-}{\rm c}}( j ) } $$

where m represents the number of subglacial water bodies. Therefore, the newly corrected threshold can be used to identify wet locations in ΔR profile and can be specified as ΔR _wet‐c. At this time, the corrected threshold can be used to identify dry locations in the ΔR profile as ΔR _dry-c = R _ir − R _iw + ΔR _wet-c. So far, the dry and wet locations can be identified only by thresholds at a regional scale, but the dry–wet distribution in DWTZs is still not effectively estimated.

By taking the ΔR profile and the terrain profile as inputs, Lang and others (2022) proposed three groups of features extracted by a feature calculation window to describe the specificity of DWTZ relative to other areas, in order to drive the SVM classification model with an Radial Basis Function kernel to automatically detect DWTZ: (1) features for ΔR profile, (2) feature for terrain profile and (3) feature for both ΔR and terrain profiles. Therefore, DWTZs in each transect can be identified.

The final subglacial dry–wet distribution of the region was generated based on the identified DWTZs and the corrected dry and wet thresholds. First, the lowest point of ΔR profile in the window of a DWTZ was taken as the reference dry location. The location where the difference between the ΔR of the location and the ΔR of the reference dry location is greater than |ΔR _wet-c − ΔR _dry-c| was determined as the reference wet location, and the midpoint of the ΔR of the reference dry and wet locations was determined as the locations where the bed is at PMP to determine the dry–wet distribution within the DWTZ. Second, for other areas in ΔR profiles, dry and wet locations were identified using the dry and wet thresholds respectively, where ΔR ≤ ΔR _dry-c represents a dry location, and ΔR ≥ ΔR _wet-c represents a wet location. In this way, the complete distribution of the subglacial dry–wet distribution could be generated, the predicted dry locations imply frozen bedrock with a high probability, and the wet locations imply the presence of subglacial water with a high probability, and the location where the bed is at PMP represents the critical transition point from cold to temperate.

An example of generating a subglacial dry–wet distribution containing the locations of where the bed is at PMP through the method proposed by Lang and others (2022) is shown in Figure 2. The survey line named TSH-GCX0g-R40a (hereinafter referred to as R40a) is taken as an example to illustrate how to generate the distribution of dry and wet zones, and locations of where the bed is at PMP. Figure 2a shows the transect of R40a, which shows that the bed layer is mainly composed of undulating bedrock and the subglacial lake $90^\circ {\rm E}$. Figure 2b shows the layer information of the ice surface and bed extracted from the transect of R40a. The ΔR profile calculated through layer information is shown in Figure 2c, and it can be seen that the ΔR values are higher in the subglacial lake and several other low-lying areas. Figure 2d shows the distribution of the subglacial dry and wet zones and locations of where the bed is at PMP. The cyan dots represent the freezing zone, the red dots represent the melting zone, the black dots represent the uncertain state, and the yellow dots represent the locations of PMP.

Figure 2. A subglacial dry–wet distribution of transect TSH-GCX0g-R40a containing the locations of where the bed is at PMP through the method proposed by Lang and others (2022). (a) Radargram of transect TSH-GCX0g-R40a. (b) The extracted information of transect TSH-GCX0g-R40a, the air-ice interface is marked in blue, the ice-bed interface is marked by purple. (c) The ΔR profile of transect TSH-GCX0g-R40a, the ΔR profile is marked by blue line, the dry threshold is marked by cyan line, the wet threshold is marked by red line. (d) The diagnosis results of dry and wet zones of transect TSH-GCX0g-R40a, the dry zone is marked in cyan, the wet zone is marked in red, unknown condition area is marked in black, and the locations where the bed is at PMP is marked in yellow.

In addition, a partial enlarged view of a DWTZ is shown in Figure 3, which is the right end of the lake $90^\circ {\rm E}$ in transect of R40a. Figure 3a shows the original transect of the DWTZ, Figure 3b shows the corresponding ΔR profile, and Figure 3c shows the distribution of the subglacial dry and wet zones. The purple dashed window represents the recognition window corresponding to the DWTZ. It can be seen that the DWTZ inside the purple window has undergone a transition from wet to dry from left to right, and the reflectivity has also changed from high to low. The midpoint of the local reflectivity inside the DWTZ window has been determined as the location where the bed is at PMP. Based on these locations the minimum ice thickness required for basal melting can be obtained by layer information of surface and bed, hereafter referred to as the Critical Ice Thickness (CIT) which can be determined by using the surface and bed layer elevations. Therefore, the distribution of dry and wet zones, and the CIT corresponding to each location where the bed is at PMP can be obtained.

Figure 3. A subglacial dry–wet distribution of a DWTZ in transect TSH-GCX0g-R40a containing the locations of where the bed is at PMP through the method proposed by Lang and others (2022). (a) The image of the DWTZ, the purple dashed window represents the recognition window corresponding to the DWTZ. (b) The ΔR profile of the DWTZ, the ΔR profile is marked by blue line, the dry threshold is marked by cyan line, the wet threshold is marked by red line. (c) The diagnosis results of dry and wet zones of transect TSH-GCX0g-R40a, the dry zone is marked in cyan, the wet zone is marked in red, unknown condition area is marked in black, and the locations where the bed is at PMP is marked in yellow.

3.2 Thermodynamic model

According to the CIT corresponding to each location where the bed is at PMP, and subglacial dry–wet distribution, we are able to use traditional glaciological methods to constrain the GHF. Since Ridge B is the ice divide in the middle of EAIS, the horizontal velocity is nearly <2 m a⁻¹ (Fig. 1a), therefore the horizontal advection can be ignored (Rignot and others, 2011; Van Lieffering and Pattyn, 2013). Similarly, the low rates of horizontal ice flow imply that strain heating from vertical shear should be small, and we neglect that term as well, along with horizontal diffusion, leaving us with a 1-D steady-state thermodynamic model (Passalacqua and others, 2017):

(4)$$\displaystyle{{k_i} \over {\rho _{\rm i}c_{\rm i}H^2}}\displaystyle{{\partial ^2T} \over {\partial \xi ^2}}-\displaystyle{{u( \xi ) } \over H}\displaystyle{{\partial T} \over {\partial \xi }} = 0$$

where T is temperature, t is time, H is ice thickness corresponding to the ice-sheet bed, ξ is the normalized vertical coordinate, and ξ = 0 on the bed, ρ _i = 910 kg m⁻³ is the density of ice, c _i = 2009 J kg⁻¹ K⁻¹ is the specific heat capacity, k _i = 2.1014 W m⁻¹ K⁻¹ is the thermal conductivity (Cuffey and Paterson, 2010; Van Lieffering and Pattyn, 2013). Following Passalacqua and others (2017) a temperature-dependent shape function is used to determine the vertical distribution of u(ξ), the change in vertical velocity relative to the bed. The boundary conditions of the 1-D steady-state thermodynamic model are: T = T _s at ξ = 1 (where T _s is mean surface temperature), u = 0 at ξ = 0 and u = −a at ξ = 1 (where −a is the accumulation rate), and dT/dξ = −GH/k _i at ξ = 0 (where G is GHF). GHF is obtained from the bed temperature gradient obtained by solving the thermodynamic equation (Hindmarsh, 1999; Van Lieffering and Pattyn, 2013; Passalacqua and others, 2017). According to the CIT of the location of PMP, the temperature of PMP, T _pmp, can be calculated (Pattyn, 2010):

(5)$$T_{{\rm pmp}} = T_0-\gamma H\xi $$

where T ₀ = 273.15 K, γ = 8.7 × 10⁻⁴ km⁻¹. For points at which we assume the basal temperature just reaches PMP, we treat GHF as a free parameter and iteratively solve Eqn (4) to find the value of GHF that brings the basal temperature up to the PMP, with the other boundary conditions (ice thickness, surface temperature, and surface accumulation rate) held constant.

3.3 Methods to calculate GHF and evaluate uncertainty

As mentioned above, GHF corresponding to the locations where the bed is at PMP is calculated by running the thermal model. These positions are distributed discretely in the region (marked by black cross in Fig. 4), thus a Kriging interpolation (Oliver and Webster, 1990) is used for spatial interpolation to obtain a preliminary GHF model.

Figure 4. (a) The final GHF model in Ridge B area. The white circles marked by label ① and label ② represent the area with relatively low GHF values and high GHF values, respectively. (b) Diagnosis results of dry and wet zones in Ridge B area. The bedrock digital elevation model in the background is from BedMachine v2 (Morlighem, 2020; Morlighem and others, 2020). (c) The uncertainties of the final GHF model. The gray lines represent the ice divides (Creyts and others, 2014), the cyan point represents dry zone, and the red point represents wet zone. The known subglacial lakes (Livingstone and others, 2022) in the region are marked by pink triangles, the locations where the bed is at PMP are marked by yellow crosses, and the supplementary correction points are marked by purple crosses.

In order to constrain the GHF more accurately, according to the distribution of the dry and wet zones, additional restrictions are imposed on the GHF model. The subglacial dry–wet distribution includes information on the distribution of subglacial water bodies and local frozen zones in the Ridge B region. The subglacial water bodies in the region can be used to limit the minimum value of the local GHF, and the lowest bed elevation of the local frozen zone can be used to limit the maximum value of the local GHF; an approach that has been widely used in previous works (Siegert, 2000; Van Liefferinge and Pattyn, 2013; Fudge and others, 2019). These end member estimates are used to adaptively modify the preliminary GHF model by finding locations where the GHF values limited by local freezing zones and subglacial water bodies contradict the preliminary GHF model. Thus, the correction points are locations where the GHF value in the preliminary GHF model is less than the minimum GHF value limited by subglacial water bodies, and locations where the GHF value in the preliminary GHF model exceeds the maximum GHF value limited by local freezing zones. Having supplemented the maximum/minimum GHF value of GHF correction points of the initial assessment, we re-perform Kriging interpolation to generate the final GHF model.

In order to evaluate the uncertainty in the inferred GHF, we employed the method proposed by Fudge and others (2019) by systematically changing the input parameters and running the thermal model. Considering the following factors, in order to calculate the final GHF uncertainty, we calculate the average value of a set of model uncertainties corresponding to surface temperature, accumulation rate and CIT, and then combine them in quadrature to give the total uncertainty. The thermal equilibration time of the ice sheet should be on the order of 100 ka or longer. We use Vostok ice core data (Petit and others, 1999) to calculate a 100 ka averaged surface temperature and accumulation rate. Then, we identify the difference between a 100 ka averaged surface temperature and present-day surface temperature, and a ratio between the long-term average accumulation rate and the present-day accumulation rate at the Vostok ice core location. By doing this we find a temperature uncertainty of 4.74 K and an accumulation rate uncertainty of 26.6%. We considered the possibility of other potential uncertainties and ultimately varied the surface temperature by 5 k and the accumulation rate by 28%. In addition, we have more accurately determined the locations of where the bed is at PMP through the method proposed by Lang and others (2022), which reduces the uncertainty in CIT for each PMP site when compared to the method used by Passalacqua and others (2017). However, in order to achieve a more realistic assessment of the uncertainty of GHF, we still vary the CIT by 3%.

4. Results

We have identified the locations where the bed is at PMP corresponding to the DWTZs (which is at the transition from cold to temperate, marked by yellow crosses in Fig. 4), and the locations of the correction points (marked by purple crosses in Figs 4a, b), and display the subglacial dry–wet distribution (marked by red and cyan points in Fig. 4b, representing the wet and dry zones, respectively). We consider that there is some consistency between the diagnosis results of dry and wet zones and the current list of subglacial lakes in Ridge B (marked by pink triangles) (Livingstone and others, 2022), especially the subglacial lake $90^\circ {\rm E}$ (in left-top corner of Fig. 4b marked by label $90^\circ {\rm E}$, the bedrock is lower with darker color) (Bell and others, 2006). It can be seen that the distribution of the wet zone along two radar profiles crossing the lake $90^\circ {\rm E}$ corresponds to the very low bedrock elevation (Fig. 4b). Most of the PMP locations exist in valleys or the transition zone from mountain peak to valley, and a few exist in the transition zone from gentle slopes to flat areas. Some black crosses occur in the middle of cyan or red lines, because there is a small line segment of the opposite color that cannot be seen at the region scale. In addition, some regions experience the boundary between dry and wet zones, but they have not been marked by yellow crosses because the pattern of subglacial dry–wet distribution in these regions does not correspond to the characterization of the DWTZ.

The final GHF model is shown in Figure 4a, which ranges from 48.5 to 65.1 mW m⁻², with an average value of 58.0 mW m⁻². In the region, there is a main ice divide extending from the south side to the north side of the region. We note that the final GHF model shows a relatively low value on the west side of the main ice divide, a relatively high value on the north side of the main ice divide, and a relatively low value on the northwest side of the region. Combining uncertainty estimation and GHF distribution, we consider that in the Ridge B region, the GHF value is minimal around label ①, ~50–52 mW m⁻², and the GHF value is maximal around label ②, ~63–65 mW m⁻².

The uncertainty of GHF estimation as shown in Figure 4(c) includes the uncertainty generated by Kriging interpolation and model input parameters. We found that the uncertainty of GHF in the study area is <8 mW m⁻², with an average value of 7 mW m⁻². We found that the uncertainty of parameters has a greater impact on the overall level of uncertainty in the region, and the uncertainty of Kriging interpolation has a greater impact on the spatial distribution of uncertainty in the region.

5. Discussion

The GHF value range in our results is within the uncertainty range of other large-scale model estimates of GHF (Fig. 5). However, the spatial distribution of GHF is different. Comparison of GHF map in this study with other GHF datasets is shown in Figure 5 and Table 1. The comparison shows that our GHF estimation in Ridge B area is higher than that of Shapiro and Ritzwoller (2004) (Fig. 5e, the regional mean value is ~47.3 mW m⁻² and the std dev. is ~1.3 mW m⁻²), Shen and others (2020) (Fig. 5f, the regional mean value is ~47.5 mW m⁻² and the std dev. is ~1.2 mW m⁻²) and An and others (2015) (Fig. 5g, the regional mean value is ~53.4 mW m⁻² and the std dev. is ~1.4 mW m⁻²). These estimations show lower overall GHF values and have lower spatial resolution. Kang and others (2022) found that in the Antarctic Lambert–Amery Glacial system, the higher GHF estimated by Li and others (2021) and Martos and others (2017) is more consistent with the distribution of subglacial lakes in the whole region, and our average value in Ridge B area is just between them. In addition, Stål and others (2021) (Fig. 5d, the regional mean value is ~56.5 mW m⁻² and the std dev. is ~2.8 mW m⁻²) predict high GHF values (>63 mW m⁻²) in southeast part of Ridge B region. There is a known GHF (50–56 mW m⁻²) in the lake Vostok (Salamatin and others, 1998; Dmitriev and others, 2016), which is more consistent with the overall trend of GHF of Martos and others (2017) and our estimation. Our results are closer to the GHF estimation of Martos and others (2017) (Fig. 5a, our regional mean value and std dev. are ~58.0 and ~3.1 mW m⁻² respectively, and Fig. 5b, their values are ~56.4 and ~3.1 mW m⁻² respectively), and also show similar spatial distribution of GHF. However, the lack of airborne magnetic data in the region resulted in lower spatial resolution of the result of Martos and others (2017).

Figure 5. The comparison of the GHF results of this study with other existing datasets in Ridge B region. The six GHF models use the same color bar to represent the range of results. The known subglacial lakes (Livingstone and others, 2022) in the region are marked by pink triangles, and the ice divides (Creyts and others, 2014) are marked by gray lines. (a) GHF results of this study. (b) GHF results from Martos and others (2017). (c) GHF results from Li and others (2021). (d) GHF results from Stål and others (2021). (e) GHF results from Shapiro and Ritzwoller (2004). (f) GHF results from Shen and others (2020). (g) GHF results from An and others (2015).

Table 1. Comparison of GHF map in this study with other GHF datasets

Li and others (2021) estimated the GHF based on airborne high-resolution magnetic data, and the overall GHF results are higher than those from other existing datasets in Ridge B region (Fig. 5c, the GHF estimation of Li and others (2021) ranges from 55 to 82 mW m⁻², the regional mean value is ~69.7 mW m⁻², and the std dev. is ~6.1 mW m⁻²). Especially in the intersection zone of ice divides, the highest GHF result of 78–82 mW m⁻² is obtained. Such GHF values are enough to make a large area of melting on the ice-sheet bed at the ice thickness of ~3 km under the condition of ignoring horizontal diffusion, horizontal thermal friction and deformation heat. But the diagnosis results of our dry and wet zone do not show the corresponding phenomenon (Fig. 4b), although there may be unknown subtle drainage networks that do not form obvious water layers.

Using a 1-D vertical heat-transfer equation, a GHF of ~54 mW m⁻² is sufficient to keep most of the subglacial lakes near Ridge B maintaining their thermal state under pressure (Siegert and Dowdeswell, 1996; Wright and Siegert, 2012). The value of the locations of the subglacial lakes in the GHF model we reported matches this value. In particular, for a typical lake $90^\circ {\rm E}$ with an ice thickness of ~4000 m in the region, we have obtained an average GHF of ~58 ± 7 mW m⁻², which can maintain it.

We found that GHF in Ridge B region has spatial variability on a small scale, like other research results on GHF in local areas (Carter and others, 2009; Schroeder and others, 2014; Passalacqua and others, 2017). Our estimates can identify local features of GHF that were previously undiscovered, the spatial variability of GHF can occur at a small scale, and a few locations in our GHF model can change ~7 mW m⁻² on a scale of ~50 km. This phenomenon mainly occurs at the intersection of ice divides in the middle of the region; the potential cause may be the subtle changes in crustal heat or geological materials, or differing geologic histories of magmatic emplacement or differences in past or ongoing hydrothermal circulation (Burton-Johnson and others, 2020). In addition, due to the significant impact of small changes in GHF on ice-sheet melting, obtaining spatial variability of GHF at a finer scale can play an important role in simulating subglacial ice melting and water distribution (Colgan and others, 2022; McCormack and others, 2022; Shackleton and others, 2023).

The new GHF model can provide effective help for addressing the search for another drilling location for an oldest ice core. The ice of more than 1 Ma found near Vostok indicates that there may be very old undisturbed ice near the ice divides area at the upstream of the lake Vostok (Ekaykin and others, 2021). Therefore, a more accurate and high spatial resolution GHF model covering Ridge B region can help us study history of basal melting by providing more accurate boundary conditions for complex 3-D ice flow models to locate the oldest ice. In addition, the new GHF model can also be used to provide more precise estimate of basal temperature, constrain the basal melting in the region and study the development of unconsolidated water-saturated sediments and subglacial hydrological network (Rémy and others, 2004; Llubes and others, 2006; Ashmore and Bingham, 2014; Burton-Johnson and others, 2020).

In the selection of interpolation schemes for obtaining regional GHF distribution, traditional statistical interpolation methods have limitations in dealing with geographical problems. Like Passalacqua and others (2017), we only considered the Kriging interpolation in geostatistics, while uncertainty estimation is influenced by interpolation methods. Our research highlights the importance of exploring more reliable spatial interpolation methods, especially for areas with sparse data. In addition, the uncertainty of our GHF model in the southeast of the region is relatively small, and the potential reason may be that the distribution of PMP locations is dense and uniform, while the uncertainty in the northwest of the region is relatively large, the potential reason may be that the radioglaciological method used to find the positions of PMPs and estimate the GHF (Lang and others, 2022) is applicable to the area with undulating terrain in the ice divide area (Passalacqua and others, 2017). It is difficult to determine whether there is a freezing and melting transition in relatively flat terrain, which means that the method may not be applicable in such places (Passalacqua and others, 2017). In addition, the uncertainty in the central and eastern parts of the region is also relatively large; the potential reason may be that there are no survey lines here. If the uncertainty of GHF model is due to insufficient measurements, our research also highlights the necessity of further RES survey in the region.

6. Conclusions

With a recently developed radioglaciological method (Lang and others, 2022), we obtained the distribution of PMP positions and the distribution of dry and wet zones at the ice-sheet bed in Ridge B region based on the RES data collected in ICECAP/PEL project, and then built a new high-resolution model of GHF in Ridge B region through a thermodynamic model. GHF in Ridge B region varies locally and ranges from 48.5 to 65.1 mW m⁻², with an average value of 58.0 mW m⁻². We introduced the method of Fudge and others (2019) to evaluate GHF uncertainty, resulting in an average uncertainty of ~5 mW m⁻². Our GHF values reveal the higher spatial variability than previous models in the region and are consistent with the current known GHF constraints for subglacial lakes in the region and the GHF derived from the Vostok ice core and fit best with respect to the mean values to the GHF models of Martos and others (2017), Stål and others (2021) and An and others (2015). This study highlights the need to take variability of local GHF on a smaller spatial scale into account when locating the oldest ice in Ridge B and other potential regions, as well as studying subglacial hydrology and geology.