Articles | Volume 16, issue 9
The Cryosphere, 16, 3619–3633, 2022
The Cryosphere, 16, 3619–3633, 2022
Research article
08 Sep 2022
Research article | 08 Sep 2022

Evaluation of six geothermal heat flux maps for the Antarctic Lambert–Amery glacial system

Evaluation of six geothermal heat flux maps for the Antarctic Lambert–Amery glacial system
Haoran Kang1, Liyun Zhao1,2, Michael Wolovick3, and John C. Moore1,4,5 Haoran Kang et al.
  • 1College of Global Change and Earth System Science, Beijing Normal University, Beijing 100875, China
  • 2State Key Laboratory of Earth Surface Processes and Resource Ecology, Beijing Normal University, Beijing 100875, China
  • 3Glaciology Section, Alfred Wegener Institute, Bremerhaven, Germany​​​​​​​
  • 4CAS Center for Excellence in Tibetan Plateau Earth Sciences, Beijing 100101, China
  • 5Arctic Centre, University of Lapland, Rovaniemi, Finland

Correspondence: Liyun Zhao ( and John C. Moore (


Basal thermal conditions play an important role in ice sheet dynamics, and they are sensitive to geothermal heat flux (GHF). Here we estimate the basal thermal conditions, including basal temperature, basal melt rate, and friction heat underneath the Lambert–Amery Glacier system in eastern Antarctica, using a combination of a forward model and an inversion from a 3D ice flow model. We assess the sensitivity and uncertainty of basal thermal conditions using six different GHF maps. We evaluate the modelled results using all observed subglacial lakes. The different GHF maps lead to large differences in simulated spatial patterns of temperate basal conditions. The two recent GHF fields inverted from aerial geomagnetic observations have the highest GHF, produce the largest warm-based area, and match the observed distribution of subglacial lakes better than the other GHFs. The modelled basal melt rate reaches 10 to hundreds of millimetres per year locally in the Lambert, Lepekhin, and Kronshtadtskiy glaciers feeding the Amery Ice Shelf and ranges from 0–5 mm yr−1 on the temperate base of the vast inland region.

1 Introduction

The Lambert–Amery system in eastern Antarctica is believed to be relatively stable against climate change and has changed little over several decades of observations (King et al., 2007). However, there is also evidence of extensive subglacial canyons and lakes (Fretwell et al., 2013; Jamieson et al., 2016; Cui et al., 2020a). Subglacial canyons and lakes are conduits for subglacial water, transporting subglacial meltwater to the coast through complex hydrologic routing, which may change on relatively fast timescales (Malczyk et al., 2020). Jamieson et al. (2016) report a large subglacial drainage network in Princess Elizabeth Land (PEL), which would transport water from central PEL to the coast passing the Lambert–Amery region. Subglacial water can affect the ice flow (Stearns et al., 2008; Diez et al., 2018), influence the dynamical stability and basal mass balance (Gudlaugsson et al., 2017), and may enhance basal melt of ice shelves (Le Brocq et al., 2013).

Ice temperature is an important factor in the rheology of ice (Budd et al., 2013) and ice flow. Whether the basal ice is at the melting point influences the movement of the ice to a great extent. Ice at the melting point can lead to water flowing along hydraulic gradients and accumulating in local depressions (Fricker et al., 2016). The meltwater lubricates the ice/bed interface or saturates any sediment till layer, allowing higher ice velocities via basal sliding. For instance, the rapid retreat of Thwaites and Pope glaciers in the Amundsen Sea sector of western Antarctica is being facilitated by high heat flow in the underlying lithosphere (Dziadek et al., 2021). This bed–ice linkage forms the basis for making inferences on basal conditions via surface observations (Pattyn, 2010) or relict landforms (e.g. Näslund et al., 2005).

The ice temperature is controlled by deformational heat generated from strain within the ice, advection of heat due to lateral ice motion and the descent rate of ice from the surface, conduction of heat through the ice, and frictional heating from basal sliding. Ice temperature is hard to evaluate because of the scarcity of in situ measurements, which are typically obtained from boreholes that are very rarely drilled through the Antarctic ice sheet. Geothermal heat flux (GHF) is an important boundary condition for ice temperature simulation and is generally the largest source of uncertainty. Hence, geophysical survey methods are used to indirectly map GHF. To date, GHF datasets have been estimated from seismic models (Shapiro and Ritzwoller, 2004; An et al., 2015; Shen et al., 2020) and derived from airborne magnetic surveys (Li et al., 2021; Martos et al., 2017) and satellite geomagnetic data (Maule et al., 2005; Purucker, 2012).

Extensive ice-penetrating radar data have been collected recently over Princess Elizabeth Land (PEL; Fig. 1d), including the eastern part of the Lambert–Amery system (Cui et al., 2020a). This fills in large data gaps from older surveys and provides the basis for our study. The radar surveys reveal  1100 km long canyons (Fig. 1c) that are incised hundreds of metres deep into the subglacial bed that extend from the Gamburtsev Subglacial Mountains (GSM) to the coast of the Western Ice Shelf (WIS). Li et al. (2021) collected airborne magnetic data that can be combined with radar ice thicknesses and estimated depths at which the bedrock reaches its Curie temperature to invert for the geothermal flux. The resulting higher-resolution dataset (Li et al., 2021) implies a larger heat flux than previous estimates in this region. Furthermore, recently discovered subglacial lakes, including potentially the second-largest subglacial lake in Antarctica, add evidence for more widespread basal melting in the region than was thought based on the much sparser earlier survey data (Cui et al., 2020b). The complex subglacial topography, relatively high geothermal heat flux, and subglacial lakes imply a complex distribution of basal thermal conditions and subglacial water networks. These heterogenous basal conditions will have shaped much of the ice flow and mass balance of the Lambert–Amery system. This motivates us to investigate how the basal thermal conditions inferred from the new high-resolution topography dataset (Cui et al., 2020a) can be reconciled with surface ice velocities and existing geothermal heat flow maps.

Ice sheet models can be used to simulate the dynamics and thermodynamics of the ice sheet. Glaciologists have combined ice sheet models with measurements of vertical temperature profiles or thawed basal states to constrain the GHF of the ice sheets (e.g. Pattyn, 2010; Rezvanbehbahani et al., 2019). In the Lambert–Amery glacial system, Pittard et al. (2016) suggest that ice flow is most sensitive to the spatial variation in the underlying GHF near the ice divides and along the edges of the ice streams.

In this study, we simulate ice basal temperatures and basal melt rates in the Lambert–Amery system using the new high-resolution digital elevation model, along with six different published GHF maps as forcing for an offline coupling between a basal energy and water flow model and a 3D full Stokes ice flow model. We evaluate the quality of the resulting basal temperature field incorporating the Stokes model estimates of ice advection, strain, and frictional heating under the different GHF maps using all available observed subglacial lakes and surface velocities. Hence, we make inferences on which GHF maps yield the best match with observations in the region.

2 Regional domain and datasets

Our modelled domain is part of the Lambert–Amery system. It consists of two drainage basins, namely the Lambert Glacier basin and the American Highland basin, along with about half of the Amery Ice Shelf (Fig. 1). The 2D domain boundary outlines are defined by the inland ice catchment basin boundary, the central streamline, and the ice front of Amery Ice Shelf. The central streamline was chosen by selecting a point at the confluence of Lambert Glacier and Lepekhin Glacier and then advecting that point downstream to the ice front using the observed velocity field. The margins of the inland subbasin and the central streamline of the Amery Ice Shelf were chosen as boundaries because the mass flux across them is assumed to be zero by definition.

Figure 1The domain topography and location with domain boundary overlain. (a) Surface elevation. (b) Ice thickness. (c) Bed elevation. (d) The location of our domain in Antarctica. The solid black curve is the outline of the study domain, including the central streamline of Amery Ice Shelf and the boundary of inland subbasins based on drainage basin boundaries defined from satellite ice sheet surface elevation and velocities (Mouginot et al., 2017; Rignot et al., 2019). The solid red and white curves in panels (a)(c) are the grounding line and margin of Antarctic, respectively (Morlighem et al., 2020). The dotted black curve is the dividing line between Lambert Glacier basin and the American Highland basin. The dotted red curves in panels (b) and (d) are the boundary of ice thickness data from Cui et al. (2020a), inside which we incorporate data from Cui et al. (2020a) and outside from MEaSUREs BedMachine Antarctica, version 2. The white stars in panel (c) denote the locations of observed subglacial lakes (Wright and Siegert, 2012; Cui et al., 2020b), and the region within the white line at (1800 E, 300 N) is potentially the second-largest subglacial lake in Antarctic. The red arrows in panel (c) indicate the routing through the deep subglacial canyon system from GSM to WIS. The subbasins names of Lambert–Amery system are labelled in panel (d), with ML for MacRobertson Land basin, FG for Fisher Glacier basin, MG for Mellor Glacier basin, LG for Lambert Glacier basin, AH for American Highland basin, and AIS for Amery Ice Shelf.

The surface elevation, bedrock elevation, and ice thickness from Cui et al. (2020a) are used in most of the domain (Fig. 1b; Table 1), with additional data are from Making Earth System Data Records for Use in Research Environments (MEaSUREs) BedMachine Antarctica, version 2, at a resolution of 500 m (Morlighem et al., 2020). The bed elevation is calculated by subtracting the ice thickness from the surface elevation.

The surface ice velocity data are obtained from MEaSUREs InSAR-based Antarctic ice velocity Map, version 2, with a resolution of 450 m (Rignot et al., 2017). Data were largely acquired during the International Polar Years, from 2007 to 2009 and between 2013 and 2016. Additional data acquired between 1996 and 2016 were used as needed to maximize coverage.

Ice sheet surface temperature data are prescribed by ALBMAP v1 (an improved Antarctic dataset for high-resolution numerical ice sheet models) with a resolution of 5 km (Le Brocq et al., 2010a) and come from monthly estimates inferred from advanced very high-resolution radiometer (AVHRR) data averaged over 1982–2004. Subglacial lake locations are from the fourth inventory of Antarctic subglacial lakes (Wright and Siegert, 2012), with the addition of the newly discovered lakes (Cui et al., 2020b).

A total of six GHF datasets (Fig. 2; Table 2) are used in this study. All the datasets are interpolated into the same 2.5 km resolution.

Table 1Datasets used in this study.

Download Print Version | Download XLSX

Table 2The six GHF datasets used in this study.

Download Print Version | Download XLSX

Figure 2The spatial distribution of GHF over our domain, as described in Fig. 1. See Table 2 for the GHF map details.


3 Model

Our goal is to infer the basal thermal conditions, including basal temperature and basal melt rate, in the domain. Geothermal heat flux, englacial heat conduction, and basal friction heat are the main heat sources that determine the basal thermal conditions. Therefore, we need to model both ice flow velocity and stress for basal friction heat and ice temperature for englacial heat conduction.

We solve an inverse problem by a full Stokes model, implemented in Elmer/Ice, to infer the basal friction coefficient such that the modelled velocity best fits observations (Gagliardini et al., 2013). Using the best-fit basal friction coefficient, we obtain the ice flow velocity, stress, and basal friction heat. A proper initial vertical ice temperature profile subject to thermal boundary conditions is needed to solve the inverse problem. To obtain it, we use a forward model that consists of an improved shallow ice approximation (SIA) thermomechanical model with a subglacial hydrology model (Wolovick et al., 2021). The forward model uses the modelled velocity direction and basal slip ratio from the full Stokes inverse model to constrain its solution. We do steady-state simulations by coupling the forward and inverse models. We will describe the forward model in Sect. 3.1 and the inverse model in Sect. 3.2 and then the coupling in Sect. 3.3.

3.1 Forward model

The forward model consists of a thermomechanical steady state model using an improved shallow ice approximation (SIA) in equilibrium with the subglacial hydrological system (Wolovick et al., 2021). It has internal consistency between the following three components: ice flow, ice temperature, and basal water flux. The numerical model requires the following three coupled components to be consistent with one another: (1) integration for balance flux and englacial temperature downhill along the ice surface, (2) integration for basal water flux and freezing rate downhill in the hydraulic potential, and (3) rheology and shape function computations to determine the distribution of ice flux and shear heating. The model performs a fixed point iteration for consistency between these three components. In addition. we improve on the model used in Wolovick et al. (2021) by combining the observed velocity field, the velocity field from the full Stokes model, and the surface gradient direction to compute a merged surface flow direction field. The observations are used where flow is fast, Elmer/Ice modelled velocity is used where flow is slow, and the surface gradient is only used near the margins of the domain where the Elmer/Ice velocity field is not reliable (Fig. 3). The simulation is done on a finite difference mesh with resolution of 2.5 km.

The surface accumulation rate we used in the forward thermal model is the mean of Arthern et al. (2006) and Van de Berg et al. (2005). Both were accessed through the ALBMAP_v1 dataset (Le Brocq et al., 2010b).

One key complexity is how to deal with basal thermal boundary condition. At the bottom of ice shelves, we set basal temperature equal to the pressure melting point. At the bed of grounded ice, the boundary condition can be either a Dirichlet or Neumann condition, depending on the basal melting and subglacial water conditions. The basal boundary conditions are given by the following:


where k(T) is the temperature-dependent thermal conductivity of ice, m is the basal melt rate, Tm is the pressure-dependent melting temperature, and G is GHF, taking the six GHF datasets listed in Table 2. The thermal condition will switch from Neumann (Eq. 1) to Dirichlet (Eq. 2) if the basal temperature exceeds the pressure-dependent melting point. The opposite switch from Dirichlet to Neumann is determined by the hydrology model if there is insufficient water input to supply a large freezing rate.

Figure 3Surface velocity direction fields, in degrees clockwise from the grid north. The first row shows the direction from the surface gradient (a), Elmer/Ice modelled velocity (b), and the observed velocity direction (c). The middle row (d–f) shows the three corresponding weighting fields (the sum of these weights is 1). The bottom row shows the difference between the direction of surface gradient and Elmer/Ice modelled velocity (g), the difference between the observed velocity direction and Elmer/Ice modelled velocity (h), and the merged velocity field used in the forward model (i).

One improvement on the method from Wolovick et al. (2021) is that a temperate basal ice layer with non-zero thickness is permitted in our model in the case that the modelled basal ice temperature reaches the pressure melting point. We do this using a weak-form solution in which the volumetric englacial melt rate rises steeply as temperature exceeds the melting point. The englacial melting absorbs latent heat and serves to limit temperature rise. We parameterize the increase in volumetric melt rate as an exponential function of temperature with a 1 K e-folding temperature and a prefactor given by the englacial strain heating and the latent heat of fusion. All englacial meltwater generated this way is assumed to immediately drain to the bed.

Another key component of the forward model is the shape function determining the distribution of horizontal velocity with depth. We also improve the shape function in Wolovick et al. (2021) by including the basal slip ratio, u^b=ub/u, where ub is the basal velocity magnitude, and u is the vertically averaged horizontal velocity magnitude. The slip ratio is taken from the full Stokes inverse model. Other than the addition of a spatially variable slip ratio, the shape function calculation is unchanged from Wolovick et al. (2021).

3.2 Inverse model with full Stokes model

The spatial distribution of basal friction in the domain is modelled by solving an inverse problem using the three-dimensional full Stokes model, Elmer/Ice, an open-source finite element method package (Gagliardini et al., 2013). The inverse model is based on adjusting the spatial distribution of the basal friction coefficient to minimize the misfit between simulated and observed surface velocities. The modelled velocity is obtained by solving the full Stokes equation, which includes conservation equations for both the momentum and mass of the ice, as follows:


where τ​​​​​​​ is the deviatoric stress tensor, p is the isotropic pressure, ρi is ice density, g is the acceleration due to gravity (0, 0, −9.81) m s−2, and v is ice velocity. According to Glen's flow relation, deviatoric stress is related to the strain rate tensor, ε˙, which can be described by τ=2ηε˙, where the effective viscosity of the ice, η, is sensitive to the temperature-dependent flow rate factor A(T) calculated using an Arrhenius equation (Cuffey and Paterson, 2010). The ice temperature distribution comes from the forward model in Sect. 3.1.

3.2.1 Mesh generation and refinement

First, we use Gmsh software (Geuzaine and Remacle, 2009) to generate an initial 2D horizontal footprint mesh with the boundary described in Sect. 2. Then we refine the mesh by an anisotropic mesh adaptation code called the Mmg library (, last access: 2 August 2022). The resulting mesh is shown in Fig. 4 and has minimum and maximum element sizes of approximately 1000 and 8000 m. The 2D mesh is then vertically extruded using 10 equally spaced, terrain-following layers.

Figure 4The refined 2D horizontal domain footprint mesh (a). Boxes outlined in panel (a) are shown in detail and overlain with surface ice velocity in panels (b) and (c) and with ice thickness in panel (d).


3.2.2 Boundary condition

The ice surface is assumed to be stress free. At the ice front, the normal stress under the sea surface is equal to the hydrostatic water pressure. On the lateral boundary, the normal stress is equal to the ice pressure applied by neighbouring glaciers and the normal velocity is assumed to be 0. The bed for grounded ice is assumed to be rigid, impenetrable, and fixed over time. Since we perform a stress balance snapshot in the full Stokes model, we do not need to prescribe surface mass balance or basal mass balance in the boundary conditions.

The normal basal velocity is set to 0 at the ice–bed interface. The linear sliding law is used to describes the relationship between the basal sliding velocity, ub, and the basal shear force, τb, on the bottom of grounded ice, as follows:

(5) τ b = C u b .

To avoid non-physical negative values, C=10β is used in the simulation. We call β the basal friction coefficient rather than C. C is initialized to a constant value of 10−4 MPa m−1 yr (Gillet-Chaulet et al., 2012) and then replaced with the inverted C in subsequent inversion steps.

3.2.3 Surface relaxation

We relax the free surface of the domain by a short transient run to reduce the non-physical spikes in initial surface geometry (Zhao et al., 2018). The transient simulation period here is 0.5 year, with a time step of 0.01 year.

3.2.4 Inversion and improvement for basal friction coefficient

Taking the results from the surface relaxation as our ice geometry, we use an inverse model to retrieve the basal friction coefficient, the deviatoric stress field, and ice velocity field. The inverse model adjusts the spatial distribution of the basal friction coefficient to minimize the value of the cost function (Morlighem et al., 2010), which is defined as the difference between the simulated surface velocity and the observed, as follows:

(6) J 0 = Γ s 1 2 | u | - | u obs | 2 d Γ ,

where Γs is the ice surface, and u and uobs are the simulated and observed surface velocities.

To avoid over-fitting of the inversion solution to non-physical noise in the observations, the following regularization term,

(7) J reg = 1 2 Γ s C x 2 + C y 2 d Γ ,

is added to the cost function, and then the total cost function is defined as follows:

(8) J tot = J 0 + λ J reg ,

where λ is a positive regularization weighting parameter. An L-curve analysis (Hansen, 2001) has been carried out for inversions to find the optimal λ by plotting the term Jreg as the function of J0. The optimal value of 1010 is chosen for λ to minimize J0.

Basal friction in reality depends on basal temperature, i.e. it is relatively large on cold beds, since the ice is frozen, and small on warm beds, where the basal temperature reaches the pressure melting point, allowing the ice to slide (Greve and Blatter, 2009). However, in the inverse model, the basal friction coefficient (Eq. 5) is adjusted to match velocity observations without regard to basal temperature, which leads to unrealistic noise manifested as local spikes in modelled basal friction heat (Fig. 5a).

We improve the parameterization of β via C in Eq. (5) (Sect. 3.2.2) by considering basal temperature Tbed, as follows:

(9) β new = β old + α T m - T bed ,

where βold is from the inverse model, α is a positive factor to be tuned, and Tm is the pressure-dependent melting temperature. βnew equals βold at a bed with temperate ice and is larger than βold at a bed with ice temperature lower than Tm. We tune α in the range of [0.1, 2], with an interval of 0.1, and find that the local spikes in modelled friction heat become fewer (Fig. 5) as α increases from 0.1 to 1 but stay almost constant with α from 1 to 2. Therefore, we take α to be 1 and use the parameterization of βnew in Eq. (5) in all the simulations. Using Eq. (9), the difference in the simulated and observed surface velocity is unchanged over the region, except for some parts of the inland boundary.

Figure 5Comparison of modelled basal friction heat with basal friction coefficient βold (a) and βnew with α=1 (b), as driven by Martos et al. (2017) GHF. The white square is enlarged.

3.2.5 Basal melt rate

Based on the inverted basal velocity and basal shear stress, we can calculate the basal friction heat. We then produce the basal melt rate using the thermal equilibrium, as follows (Greve and Blatter, 2009):

(10) M = G + u b τ b + k ( T ) d T d z ρ i L ,

where M is the basal melt rate, G is GHF, ubτb is the basal friction heat, -kTdTdz is the upward heat conduction, ρi is the ice density, and L is the latent heat of ice melt. Geothermal heat and frictional heating from basal slip warm the base, while the upward heat conduction to the interior cools the base. Note that the basal melt rate can be either positive (melting) or negative (freezing), depending on the heat balance.

3.3 Experimental design of coupled simulations

We design the coupled simulations in an eight-step scheme for coupling the forward model and inverse model, similar to Zhao et al. (2018):

  1. We run the forward model with the velocity direction taken from a mixture of the surface gradient and surface velocity observations and obtain an initial modelled englacial temperature (Fig. 3).

  2. We do surface relaxation in Elmer/Ice with the englacial temperature from step 1.

  3. Taking the results from step 2 as the initial state, we do an inversion in Elmer/Ice, using the modelled englacial temperature from step 1, to obtain a modelled surface velocity best fit to the observed surface velocity. The modelled surface velocity will remove some artefacts in the observed field.

  4. We run the forward model using the velocity directions derived by merging the Elmer/Ice modelled velocity, the surface gradient, and the surface velocity observations (Fig. 3). We use the modelled velocity from the full Stokes inverse model to constrain the basal slip ratio, then constrain the rheology and shape function in the forward model. Then we obtain an updated modelled englacial temperature.

  5. We run the inverse model in Elmer/Ice with the improved englacial temperature from step 4 and obtain an updated modelled velocity.

  6. We run the forward model again using the ratio of basal sliding to column-averaged velocity in Elmer/Ice from step 5 to constrain the slip ratio and obtain a further updated basal temperature.

  7. We run the inverse model again in Elmer/Ice with the improved englacial temperature from step 6 and obtain an updated modelled velocity and stress.

  8. We analyse the modelled results in step 7 and calculate basal friction heat and basal melt rate.

We perform the above procedure for all six sets of GHF to produce six different results for the basal thermal conditions.

4 Simulation results

4.1 Ice velocity

In the inverse model, the misfit between the modelled and the observed surface velocity is minimized. Therefore, we obtain very similar distributions of modelled surface velocity field using different GHF maps. Figure 6 shows the modelled velocity in the experiment, using the Martos et al. (2017) GHF as an example. The modelled surface velocity shows spatial similarities to the observed surface velocity (Fig. 6a, b). In total, three fast-flowing outlet glaciers (Lambert, Lepekhin, and Kronshtadtskiy glaciers) deliver ice to the ice shelf. The velocity of the Lambert Glacier exceeds 800 m yr−1 at the grounding line. The Lepekhin Glacier and the Kronshtadtskiy Glacier have maximum flow velocities of about 200 and 400 m yr−1 at their grounding lines, respectively. Regions with large differences between modelled and observed surface velocity occupy a small fraction of the whole area (Fig. 6c) and are associated with high-velocity gradients. Ice velocity decreases with depth. Figure 6c shows the modelled basal ice velocity. The maximum basal velocity on the Lambert Glacier exceeds 500 m yr−1 near the grounding line, and maximum basal velocities on Lepekhin Glacier and the Kronshtadtskiy Glacier reach about 150 and 200 m yr−1 at the grounding line.

Figure 6(a) Observed surface velocity, (b) modelled surface velocity in the experiment, using the Martos et al. (2017) GHF, and (c) modelled basal velocity. The white solid lines in panels (a), (b), and (c) represent speed contours of 30, 50, 100, and 200 m yr−1, respectively. The three fast-flowing outlet glaciers in panel (a) from left to right are the Lambert, Lepekhin, and Kronshtadtskiy glaciers.

4.2 Basal ice temperature and heat conduction

In Fig. 7, we show the modelled basal temperature from the six experiments. The modelled ice basal temperatures in the fast-flowing regions are all at the pressure melting point (warm). However, there are significant differences in the modelled distribution of warm-based conditions in the slow-flowing region using different GHF maps. The basal temperature is highly dependent on the GHF. In the experiment using Li et al. (2021) GHF (Fig. 7f), which has the highest GHF within the domain, the basal temperature is at the melting point over most of the domain, with extensive cold-based regions confined to the southern part. The experiment using the Martos et al. (2017) GHF (Fig. 7a), which has the second highest GHF, yields the second-largest area of a warm base, and the experiment using Purucker (2012) GHF (Fig. 7e), with the lowest GHF, gives the smallest warm-based area which is concentrated around the fast-flowing ice. All experiments display cold basal temperatures to the southwest of the Lambert Glacier basin, which is associated with thin ice over subglacial mountains (Fig. 1c).

Figure 7Modelled basal temperature relative to pressure melting point, in panels (a) to (f), corresponding to the GHF in panels (a) to (f) in Fig. 2. The ice bottom at the pressure melting point is delineated by a white contour.

Figure 8Modelled heat change in basal ice by upward englacial heat conduction (mW m−2). The negative sign means that the upward englacial heat conduction causes heat loss from the basal ice, as defined by the colour bar with cooler colours representing more intense heat loss by conduction. Panels (a) to (f) correspond to the GHF in panels (a) to (f) in Fig. 2. The white solid curves represent modelled speed contours of 30, 50, 100, and 200 m yr−1, which is the same as in Fig. 6b.

Figure 8 shows the modelled heat change of basal ice by upward englacial heat conduction in the six experiments. In most regions of the fast-flowing tributaries with velocity higher than 30 m yr−1, the heat loss caused by upward basal heat conduction is lower than 30 mW m−2 in all experiments, reflecting the development of a temperate basal layer that limits the basal thermal gradient. For the vast inland areas, experiments yield heat loss by upward heat conduction in the range of 45–60 mW m−2, except for the experiment driven by the Purucker (2012) GHF, which has lower values around 30–45 mW m−2. This is because the upward heat conduction equals GHF where the basal temperature is below the pressure melting point, and the Purucker (2012) GHF is lower than the others.

4.3 Basal friction heat

There is no significant difference in modelled basal friction heat across these six experiments, reflecting the fact that all of them have been tuned to match the surface velocity observations. So, we show only the modelled basal friction driven by the Martos et al. (2017) GHF (Fig. 5b). As expected, basal friction heat is high in fast-flowing regions. The three fast-flowing tributaries have friction heat amounting to more than 50 mW m−2, with the Lambert and Kronshtadtskiy glaciers having 2000 mW m−2 at the grounding line.

4.4 Basal melt rate

We obtain the basal melt rate using the thermal balance equation (Eq. 10). Figure 9 shows the modelled basal melt rate in the six experiments using different GHF maps. Regions with the basal melt rate coincide with a warm base where basal temperatures reach the pressure melting point. There are significant differences in the area of basal melting among the six experiments due to the large variability in GHF. The experiments using the Li et al. (2021) and Martos et al. (2017) GHF yield the largest area with basal melting. In contrast, the experiment using the Purucker (2012) GHF gives the least area with basal melting (Fig. 9).

The modelled basal melt rate is below 5 mm yr−1 in the parts of the vast inland region that are warm based. Higher basal melt rates occur in fast-flowing regions (Fig. 9) where frictional heat is high (Fig. 5b), despite the differences in GHF (Fig. 2). Basal melt rate is above 10 mm yr−1 near the grounding line, reaching 500 mm yr−1 at the grounding line of the central flowline running onto the Amery Ice Shelf. Thus, in fast-flowing regions, frictional heat is the dominant factor rather than GHF, consistent with Larour et al. (2012), who noted that slower-flowing ice in the interior of the ice sheet will be more sensitive to the GHF, but frictional heat dominates GHF in regions of fast ice flow. We use the positions of observed subglacial lakes to validate simulated regions with basal melting (Fig. 9). The modelled warm base in the experiment using Li et al. (2021) GHF covers all the observed subglacial lakes in the domain (Fig. 9f), including the recently discovered second-largest subglacial lake in Antarctica (Cui et al., 2020b). The warm base in the experiment using the Martos et al. (2017) GHF covers the second greatest number of observed subglacial lakes (Fig. 9a), and the experiment using the An et al. (2015) GHF covers the third (Fig. 9c). The experiment using the Shen et al. (2020) GHF captures two subglacial lakes in the southwest of the domain (Fig. 9b), while the experiment using the Shapiro and Ritzwoller (2004) GHF missed many known subglacial lakes in the southwest of the domain but successfully captures the recently discovered second-largest subglacial lake (Fig. 9d). The experiment using the Purucker (2012) GHF performs worst in recovering subglacial lake locations (Fig. 9e).

There are localized negative values of basal melt rate, indicating basal refreezing at three locations (Fig. 9). The modelled refreezing locations are generally characterized by large gradients in ice thickness, typically thinning by 700 m across a distance of 2 km. Radar surveys have not yet been carried out to confirm these freeze-on locations.

Figure 9Modelled basal melt rate (mm yr−1) in panels (a) to (f) correspond to the GHF in panels (a) to (f) in Fig. 2. The ice bottom at pressure melting point is surrounded by a red contour. The stars denote the locations of observed subglacial lakes, and the area surrounded by the black line is the likely second largest subglacial lake in Antarctica. There is modelled basal refreezing at three local places painted in black.

5 Discussion

Uncertainties and bias in our simulations can come from several sources. We expect that the present-day accumulation rate field in our modelling will be higher than the long-term average because of a lower accumulation rate during glacial periods (Watanabe et al., 2003; Van Ommen et al., 2004). This will tend to increase the downward advection of cold ice in our model, lowering the basal temperature in comparison to reality. On the other hand, we also expect that the modern-day surface temperature in our modelling will be higher than the long-term average temperature, again because of lower temperatures during glacial periods. This will tend to increase our modelled basal temperature in comparison with reality. It is unclear which of these competing biases is stronger.

Subglacial topography has an influence on geothermal heat at kilometre scales. Typically, it has been assumed that subglacial ridges receive less heat flow and subglacial valleys receive more heat flow, in comparison to the regional average (e.g. van der Veen et al., 2007; Colgan et al., 2021). However, the effect depends on subglacial rock type. Heat tends to follow the path of least resistance to the surface. The thermal conductivity of rock varies with lithology and can be either greater or smaller than the thermal conductivity of ice (Willcocks and Hasterok, 2019); thus, the sign of the topographic effect on GHF can be either negative or positive. Without knowing a priori whether the topographic effect will be positive or negative, it is hard to apply a topographic correction field to the GHF input field.

GHF distribution largely governs basal thermal conditions. Many previous studies (Larour et al., 2012; Pattyn, 2010; Pittard et al., 2016; Van Liefferinge and Pattyn, 2013; Van Liefferinge et al., 2018) on basal temperature and basal melt have used the Shapiro and Ritzwoller (2004), Fox Maule et al. (2005), Purucker (2012), and An et al. (2015) GHF datasets, with a few making use of the more recent Martos et al. (2017) and Li et al. (2021) GHF datasets. In this study, we find that the Li et al. (2021) and Martos et al. (2017) GHF datasets have higher GHF than the earlier datasets in the Lambert–Amery domain and consequently have the largest area with a warm base. The warmer basal conditions best match the observed distribution of subglacial lakes. However, it should be noted that observations of subglacial lakes are a one-sided constraint. A model result that does not predict basal melt at the location of the observed lakes is clearly too cold at that location. But if the model result shows basal melt at a place with no observed lakes, it is not clear whether this is because the model is too warm, the subglacial water exists in a form other than ponded lakes, or that lakes are present, but we do not have the data to detect them.

A lake complex beneath the Devon Island ice cap in Canada exists at temperatures well below pressure melting point due to large concentrations of dissolved salts (Rutishauser et al., 2018), and no similar ones are known to exist beneath the Antarctic ice sheet. Furthermore, relatively high electrical conductivity beds such as clay-rich sediments surrounded by bedrock can give rise to false positives in radar detections of subglacial water bodies (Tulaczyk and Foley, 2020​​​​​​​).

Our simulations make improvements on previous approaches. We use the full Stokes flow model in the inversion of basal friction field rather than a simplified physics model as in Wolovick et al. (2021). We also improve on the treatment of the basal friction field by imposing a larger basal friction where the ice bottom is colder than the pressure melting point and which increases with the temperature difference from freezing point. These modifications produce more physically meaningful results, since we expect frozen beds to have high basal friction. Hence, the basal friction field is constrained by simulated temperatures in addition to producing the best fitting match of simulated and observed surface velocities.

Van Liefferinge and Pattyn (2013) estimated basal temperature for the Antarctic ice sheet using three GHF datasets (Fox Maule et al., 2005; Shapiro and Ritzwoller, 2004; Purucker, 2012), and each of the datasets were improved by the method in Pattyn (2010). Their modelled temperatures show spatial similarities to our experiment field using the Purucker (2012) GHF. Pittard et al. (2016) did sensitivity experiments of the Lambert–Amery glacial system based on three GHF fields (Fox Maule et al., 2005; An et al., 2015; Shapiro and Ritzwoller, 2004), using the Parallel Ice Sheet Model (PISM), and found that modelled basal temperature reached the pressure melting point only under the fast-flowing ice, with maximum melting rates of 500 mm yr−1 at places very close to the grounding line of the central flowline onto the Amery Ice Shelf. We also model maximum basal melt at similar locations in the six GHF experiments. However, the Pittard et al. (2016) region of basal melt is mainly confined to the Lambert Glacier tributary and matches only that of our experiment using Purucker (2012) GHF.

We analyse the contribution of GHF and frictional heat to basal melt. The basal friction is a significant heat sources only under fast-flowing ice. Most GHF distributions (except Martos et al., 2017, and Li et al., 2021) in the grounded ice sheet near the ice shelf are homogeneous, but frictional heating in the fast-flowing ice is more than 10 times higher than that in the slow-flowing ice. Thus slower-flowing ice in the interior of the ice sheet is more sensitive to the GHF than fast-flowing ice (Larour et al., 2012).

GHF has its largest impact on the basal melt of the inland ice sheet. There are two principle ways to constrain GHF, i.e. (1) direct measurement and (2) inversion by multiple geophysical methods. The GHFs used in this study are based on the inversion of satellite or aeromagnetic data and seismic tomography. Direct observations of heat flux are difficult to obtain in Antarctica, and satellite data are low resolution. The most efficient method is to invert the heat flux through aerial geomagnetic observation such as for the Martos and Li GHF fields (Martos et al., 2017; Li et al., 2021). However, there are still large data gaps in remote regions, especially in PEL, leaving just inversion using satellite magnetic data with a lower resolution. The Li et al. (2021) field uses the latest aeromagnetic data to estimate the GHF in the PEL region, and this gives higher values than derived previously.

To validate the modelled basal melt, we use the locations of detected subglacial lakes. There may be many other undiscovered subglacial lakes beneath the study area, and further discoveries would help us validate the model results and possibly refine GHF maps. In addition, further observational constraints with a two-sided sensitivity to ice temperature, such as observations of subglacial freeze-on or measurements of englacial attenuation, would help us to identify areas in which the GHF maps are too warm, in addition to those areas in which they are too cold.

6 Conclusions

In this paper, we estimate the basal thermal conditions of the Lambert–Amery system by coupling a forward model and an inverse model, based on six different GHF datasets. We analyse the contribution of GHF, heat conduction, and basal friction to the modelled basal melt rate. We verify the result using the locations of all known subglacial lakes and evaluate the reliability of six GHF datasets in our study domain.

Our approach is distinct from that used to find GHF fields employed by Wolovick et al. (2021); in particular, the use of a full Stokes model allows the method to be extended to fast-flowing ice streams and ice shelf domains where neither the shallow ice nor shallow shelf approximations are valid. We also improve the basal friction calculation to include information on the basal ice temperature relative to its pressure melting point. This procedure results in removal of unrealistic noise manifested as local spikes in modelled basal friction heat.

We find significant differences in the spatial extent of temperate ice in the slow-flowing areas among the six experiments due to large variability in GHF. The experiments using the Li et al. (2021) and the Martos et al. (2017) GHF yield the largest area with basal melting and match the subglacial lake locations best. In contrast, the experiments using the Purucker (2012) GHF give the least area with basal melting and the worst match with subglacial lake locations. We suggest that the GHF datasets from Li et al. (2021) and Martos et al. (2017) are the most suitable choice for this study region. We cannot make our own GHF map from our analysis. While we can pick the GHF in places where the Li and Martos geothermal heat flow maps (Li et al., 2021; Martos et al., 2017) are consistent and both agree with the observations, we do not know which (if either) are correct where the Li and Martos GHF datasets disagree, and there are no observations. In order to make this determination, we would need additional observational constraints on the basal thermal state, such as measured basal temperatures from deep ice cores or observed refreeze-on, but neither are available in the region.

The fast-flowing region has fast basal velocities and high frictional heat, but there are large differences in basal melting rates between the six GHF datasets. The fast-flowing tributaries have frictional heating in the range of 50–2000 mW m−2. In the vast inland areas, our experiments generally yield high upward heat conduction in the range of 45–60 mW m−2, which means that GHF dominates the heat content of the basal ice in the slow flow regions. The modelled basal melt rate reaches 50–500 mm yr−1 locally in three very fast flow tributaries (Lambert, Lepekhin, and Kronshtadtskiy glaciers) feeding the Amery Ice Shelf and is in the range of 0–5 mm yr−1 in the inland region.

Data availability

MEaSUREs BedMachine Antarctica, version 2, is available at (Morlighem, 2020). MEaSUREs InSAR-based Antarctic ice velocity Map, version 2, is available at (Rignot et al., 2017). The fourth inventory of Antarctic subglacial lake is available at (Wright and Siegert, 2012). MEaSUREs Antarctic Boundaries for IPY 2007–2009 from Satellite Radar, version 2, is available at (Mouginot et al., 2017). ALBMAP v1 and the GHF dataset of Shapiro and Ritzwoller (2004) are available at (Le Brocq et al., 2010b). The GHF dataset of An et al. (2015) is available at​​​​​​​ (last access: 5 August 2022). The GHF dataset of Shen et al. (2020) is available at (last access: 5 August 2022). The GHF dataset of Martos (2017) is available at The GHF dataset of Purucker (2012) is available at (last access: 5 August 2022). The GHF dataset of Li et al. (2021) is available upon request from Lin Li. The modelled basal melt rate in this paper is available at (Zhao et al., 2022).​​​​​​​

Author contributions

LZ and JCM conceived the study. LZ, MW, and JCM designed the methodology. HK and LZ carried out the inverse model and produced the estimates and most figures. MW carried out the forward model and produced one figure. LZ wrote the original draft, and all the authors revised the paper.

Competing interests

The contact author has declared that none of the authors has any competing interests.


Publisher's note: Copernicus Publications remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.


This work has been supported by the National Natural Science Foundation of China (grant no. 41941006), National Key Research and Development Program of China (grant no. 2021YFB3900105), and State Key Laboratory of Earth Surface Processes and Resource Ecology (grant no. 2022-ZD-05).

Financial support

This research has been supported by the National Natural Science Foundation of China (grant no. 41941006), National Key Research and Development Program of China (grant no. 2021YFB3900105), and State Key Laboratory of Earth Surface Processes and Resource Ecology (grant no. 2022-ZD-05).

Review statement

This paper was edited by Pippa Whitehouse and reviewed by William Colgan, Chen Zhao, and one anonymous referee.


An, M., Wiens, D. A., Zhao, Y., Feng, M., Nyblade, A. A., Kanao, M., Li, Y., Maggi, A., and Lévêque, J.: Temperature, lithosphere-asthenosphere boundary, and heat flux beneath the Antarctic Plate inferred from seismic velocities, J. Geophys. Res.-Sol. Ea., 120, 359–383,, 2015 (data available at:​​​​​​​, last access: 5 August 2022). 

Arthern, R. J., Winebrenner, D. P., and Vaughan, D. G. Antarctic snow accumulation mapped using polarization of 4.3-cm wavelength microwave emission, J. Geophys. Res.-Atmos., 111, D06107,, 2006 

Budd, W. F., Warner, R. C., Jacka, T., Li, J., and Treverrow, A.: Ice flow relations for stress and strain-rate components from combined shear and compression laboratory experiments, J. Glaciol., 59, 374–392,, 2013. 

Colgan, W., MacGregor, J. A., Mankoff, K. D., Haagenson, R., Rajaram, H., Martos, Y. M., Morlighem, M., Fahnestock, M. A., and Kjeldsen, K. K.: Topographic correction of geothermal heat flux in Greenland and Antarctica. J. Geophys. Res.-Earth, 126, e2020JF005598,, 2021. 

Cuffey, K. M. and Paterson, W. S. B.: The physics of glaciers, fourth edition, Elsevier, Burlington, ISBN 978-0-12-369461-4, 2010. 

Cui, X., Jeofry, H., Greenbaum, J. S., Guo, J., Li, L., Lindzey, L. E., Habbal, F. A., Wei, W., Young, D. A., Ross, N., Morlighem, M., Jong, L. M., Roberts, J. L., Blankenship, D. D., Bo, S., and Siegert, M. J.: Bed topography of Princess Elizabeth Land in East Antarctica, Earth Syst. Sci. Data, 12, 2765–2774,, 2020a. 

Cui, X., Lang, S., Guo, J., and Sun, B.: Detecting and Searching for subglacial lakes through airborne radio-echo sounding in Princess Elizabeth Land (PEL), Antarctica, E3S Web Conf., 163, 04002,, 2020b. 

Diez, A., Matsuoka, K., Ferraccioli, F., Jordan, T. A., Corr, H. F., Kohler, J., Olesen, A. V., and Forsberg, R.: Basal Setings Control Fast Ice Flow in the Recovery/Slessor/Bailey Region, East Antarctica, Geophys. Res. Lett., 45, 2706–2715,, 2018. 

Dziadek, R., Ferraccioli, F., and Gohl, K.: High geothermal heat flow beneath Thwaites Glacier in West Antarctica inferred from aeromagnetic data, Commun. Earth Environ., 2, 162,, 2021. 

Fox Maule, C., Purucker, M. E., Olsen, N., and Mosegaard, K.: Heat flux anomalies in Antarctica revealed by satellite magnetic data, Science, 309, 464–467,, 2005. 

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

Fricker, H. A., Siegfried, M. R., Carter, S. P., and Scambos, T. A.: A decade of progress in observing and modelling Antarctic subglacial water systems, Philos. T. R. Soc. A, 374, 20140294,, 2016. 

Gagliardini, O., Zwinger, T., Gillet-Chaulet, F., Durand, G., Favier, L., de Fleurian, B., Greve, R., Malinen, M., Martín, C., Råback, P., Ruokolainen, J., Sacchettini, M., Schäfer, M., Seddik, H., and Thies, J.: Capabilities and performance of Elmer/Ice, a new-generation ice sheet model, Geosci. Model Dev., 6, 1299–1318,, 2013. 

Geuzaine, C. and Remacle, J. F.: Gmsh: A 3-D finite element mesh generator with built-in pre- and post-processing facilities, Int. J. Numer. Meth. Eng., 79, 1309–1331,, 2009. 

Greve, R. and Blatter, H.: Dynamics of Ice Sheets and Glaciers, Advances in Geophysical and Environmental Mechanics and Mathematics, Series Editor: Hutter, K., Springer, ISBN 978-3-642-03414-5, 2009. 

Gillet-Chaulet, F., Gagliardini, O., Seddik, H., Nodet, M., Durand, G., Ritz, C., Zwinger, T., Greve, R., and Vaughan, D. G.: Greenland ice sheet contribution to sea-level rise from a new-generation ice-sheet model, The Cryosphere, 6, 1561–1576,, 2012. 

Gudlaugsson, E., Humbert, A., Andreassen, K., Clason, C. C., Kleiner, T., and Beyer, S.: Eurasian ice-sheet dynamics and sensitivity to subglacial hydrology, J. Glaciol., 63, 556–564,, 2017. 

Hansen, P. C.: The L-curve and its use in the numerical treatment of inverse problems, in: Computational inverse problems in electrocardiology, edited by: Johnstone, P., WIT Press, Southampton, pp. 119–142, ISBN 1-85312-614-4, 2001. 

Jamieson, S. S., Ross, N., Greenbaum, J. S., Young, D. A., Aitken, A. R., Roberts, J. L., Blankenship, D. D., Bo, S., and Siegert, M. J.: An extensive subglacial lake and canyon system in Princess Elizabeth Land, East Antarctica, Geology, 44, 87–90,, 2016. 

King, M. A., Coleman, R., Morgan, P. J., and Hurd, R. S.: Velocity change of the Amery Ice Shelf, East Antarctica, during the period 1968–1999, J. Geophys. Res.-Earth, 112, F01013,, 2007.​​​​​​​ 

Larour, E., Morlighem, M., Seroussi, H., Schiermeier, J., and Rignot, E.: Ice flow sensitivity to geothermal heat flux of Pine Island Glacier, Antarctica, J. Geophys. Res.-Earth, 117, F04023,, 2012.​​​​​​​ 

Le Brocq, A. M., Payne, A. J., and Vieli, A.: An improved Antarctic dataset for high resolution numerical ice sheet models (ALBMAP v1), Earth Syst. Sci. Data, 2, 247–260,, 2010a. 

Le Brocq, A. M., Payne, A. J., and Vieli, A.: Antarctic dataset in NetCDF format, PANGAEA [data set],, 2010b. 

Le Brocq, A. M., Ross, N., Griggs, J. A., Bingham, R. G., Corr, H. F. J., Ferraccioli, F., Jenkins, A., Jordan, T. A., Payne, A. J., Rippin, D. M., and Siegert, M. J.: Evidence from ice shelves for channelized meltwater flow beneath the Antarctic Ice Sheet, Nat. Geosci., 6, 945–948,, 2013. 

Li, L., Tang, X., Guo, J., Cui, X., Xiao, E., Latif, K., Sun, B., Zhang, Q., and Shi, X.: Inversion of Geothermal Heat Flux under the Ice Sheet of Princess Elizabeth Land, East Antarctica, Remote Sensing, 13, 2760,, 2021. 

Malczyk, G., Gourmelen, N., Goldberg, D., Wuite, J., and Nagler, T.: Repeat subglacial lake drainage and filling beneath Thwaites Glacier. Geophys. Res. Lett., 47, e2020GL089658,, 2020. 

Martos, Y. M.: Antarctic geothermal heat flux distribution and estimated Curie Depths, links to gridded files, PANGAEA [data set],, 2017. 

Martos, Y. M., Catalán, M., Jordan, T. A., Golynsky, A., Golynsky, D., Eagles, G., and Vaughan, D. G.: Heat flux distribution of Antarctica unveiled, Geophys. Res. Lett., 44, 11417–11426,, 2017. 

Maule, C. F., Purucker, M. E., Olsen, N., and Mosegaard, K.: Heat flux anomalies in Antarctica revealed by satellite magnetic data, Science​​​​​​​, 309, 464–467,, 2005. 

Morlighem, M.: MEaSUREs BedMachine Antarctica, Version 2, Boulder, Colorado USA, NASA National Snow and Ice Data Center Distributed Active Archive Center [data set],, 2020. 

Morlighem, M., Rignot, E., Seroussi, H., Larour, E., Ben Dhia, H., and Aubry, D.: Spatial patterns of basal drag inferred using control methods from a full-Stokes and simpler models for Pine Island Glacier, West Antarctica, Geophys. Res. Lett., 37, L14502,, 2010. 

Morlighem, M., Rignot, E., Binder, T., Blankenship, D., Drews, R., Eagles, G., Eisen, O., Ferraccioli, F., Forsberg, R., Fretwell,P., Goel,V., Greenbaum, J. S., Gudmundsson, H., Guo, J., Helm,V., Hofstede, C., Howat, I., Humbert, A., Jokat, W., Karlsson, N. B., Lee, W., Matsuoka, K., Millan, R., Mouginot, J., Paden,J., Pattyn, F., Roberts, J., Rosier, S., Ruppel, A., Seroussi, H., Smith, E. C., Steinhage, D., Sun, B., Van den Broeke, M. R., Van Ommen, T. D., Van Wessem, M., and Young D. A.: Deep glacial troughs and stabilizing ridges unveiled beneath the margins of the Antarctic ice sheet, Nat. Geosci., 13, 132–137,, 2020. 

Mouginot, J., Scheuchl, B., and Rignot, E.: MEaSUREs Antarctic Boundaries for IPY 2007–2009 from Satellite Radar, Version 2, National Snow and Ice Data Center [data set],, 2017. 

Näslund, J.-O., Jansson, P., Fastook, J. L., Johnson, J., and Andersson, L.: Detailed spatially distributed geothermal heat-flow data for modeling of basal temperatures and meltwater production beneath the Fennoscandian ice sheet, Ann. Glaciol., 40, 95–101,, 2005. 

Pattyn, F.: Antarctic subglacial conditions inferred from a hybrid ice sheet/ice stream model, Earth Planet. Sc. Lett., 295, 451–461,, 2010. 

Pittard, M., Roberts, J., Galton-Fenzi, B., and Watson, C.: Sensitivity of the Lambert-Amery glacial system to geothermal heat flux, Ann. Glaciol., 57, 56–68,, 2016. 

Purucker, M.: Geothermal heat flux data set based on low resolution observations collected by the CHAMP satellite between 2000 and 2010, and produced from the MF-6 model following the technique described in Fox Maule et al. (2005), Interactive System for Ice sheet Simulation [data set], (last access: 5 August 2022), 2012. 

Rezvanbehbahani, S., Stearns, L. A., Van der Veen, C. J., Oswald, G. K. A., and Greve, R.: Constraining the geothermal heat flux in Greenland at regions of radar-detected basal water, J. Glaciol., 65, 1023–1034,, 2019. 

Rignot, E., Mouginot, J., and Scheuchl, B.: MEaSUREs InSAR-based Antarctica ice velocity map, version 2, Boulder, Colorado USA, NASA National Snow and Ice Data Center Distributed Active Archive Center [data set],, 2017. 

Rignot, E., Mouginot, J., Scheuchl, B., Van Den Broeke, M., Van Wessem, M. J., and Morlighem, M.: Four decades of Antarctic Ice Sheet mass balance from 1979–2017, P. Natl. Acad. Sci. USA, 116, 1095–1103,, 2019. 

Rutishauser A., Blankenship, D. D., Sharp, M., Skidmore, M. L., Greenbaum, J. S., Grima, C., Schroeder, D. M., Dowdeswell, J. A., Young, D. A., Discovery of a hypersaline subglacial lake complex beneath Devon Ice Cap, Canadian Arctic, Sci. Adv., 4, eaar4353,, 2018. 

Shapiro, N. M. and Ritzwoller, M. H.: Inferring surface heat flux distributions guided by a global seismic model: particular application to Antarctica, Earth Planet. Sc. Lett., 223, 213–224,, 2004. 

Shen, W., Wiens, D. A., Lloyd, A. J., and Nyblade, A. A.: A geothermal heat flux map of Antarctica empirically constrained by seismic structure, Geophys. Res. Lett., 47, e2020GL086955,, 2020 (data available at:, last access: 5 August 2022). 

Stearns, L. A., Smith, B. E., and Hamilton, G. S.: Increased flow speed on a large East Antarctic outlet glacier caused by subglacial floods, Nat. Geosci., 1, 827–831​​​​​​​,, 2008. 

Tulaczyk, S. M. and Foley, N. T.: The role of electrical conductivity in radar wave reflection from glacier beds, The Cryosphere, 14, 4495–4506,, 2020. 

Van der Veen, C. J., Leftwich, T., von Frese, R., Csatho, B. M., and Li, J. Subglacial topography and geothermal heat flux: Potential interactions with drainage of the Greenland ice sheet, Geophys. Res. Lett., 34, L12501,, 2007. 

Van de Berg, W. J., Van den Broeke, M. R., Reijmer, C. H., and Van Meijgaard, E.: Characteristics of the Antarctic surface mass balance, 1958–2002, using a regional atmospheric climate model, Ann. Glaciol., 41, 97–104,, 205. 

Van Liefferinge, B. and Pattyn, F.: Using ice-flow models to evaluate potential sites of million year-old ice in Antarctica, Clim. Past, 9, 2335–2345,, 2013. 

Van Liefferinge, B., Pattyn, F., Cavitte, M. G. P., Karlsson, N. B., Young, D. A., Sutter, J., and Eisen, O.: Promising Oldest Ice sites in East Antarctica based on thermodynamical modelling, The Cryosphere, 12, 2773–2787,, 2018. 

Van Ommen, T. D., Morgan, V., and Curran, M. A. J.: Deglacial and Holocene changes in accumulation at Law Dome, East Antarctic, Ann. Glaciol., 39, 359–365,, 2004. 

Watanabe, O., Shoji, H., Satow, K., Motoyama, H., Fujii, Y., Narita, H., and Aoki, S.: Dating of the Dome Fuji Antarctica deep ice core, Mem. Natl. Inst. Polar Res. Spec. Iss., 57, 25–37, 2003. 

Willcocks, S. and Hasterok, D.: Thermal refraction: Impactions for subglacial heat flux, ASEG Extended Abstracts, 2019​​​​​​​​​​​​​​, 1–4​​​​​​​,, 2019. 

Wolovick, M. J., Moore, J. C., and Zhao, L.: Joint inversion for surface accumulation rate and geothermal heat flow from ice-penetrating radar observations at Dome A, East Antarctica. Part I: model description, data constraints, and inversion results, J. Geophys. Res.-Earth, 126, e2020JF005937,, 2021. 

Wright, A. and Siegert, M.: A fourth inventory of Antarctic subglacial lakes, Antarct. Sci., 24, 659–664,, 2012. 

Zhao, C., Gladstone, R. M., Warner, R. C., King, M. A., Zwinger, T., and Morlighem, M.: Basal friction of Fleming Glacier, Antarctica – Part 1: Sensitivity of inversion to temperature and bedrock uncertainty, The Cryosphere, 12, 2637–2652,, 2018.  

Zhao, L., Kang, H. Wolovick, M., and Moore, J. C.: Modelled basal melt rate underneath the Antarctic Lambert–Amery glacial system, Zenodo [data set],, 2022. 

Short summary
Basal thermal conditions are important to ice dynamics and sensitive to geothermal heat flux (GHF). We estimate basal thermal conditions of the Lambert–Amery Glacier system with six GHF maps. Recent GHFs inverted from aerial geomagnetic observations produce a larger warm-based area and match the observed subglacial lakes better than the other GHFs. The modelled basal melt rate is 10 to hundreds of millimetres per year in fast-flowing glaciers feeding the Amery Ice Shelf and smaller inland.