1. Introduction
The Antarctic ice sheet represents the largest source of uncertainty in projections of future sea level rise (e.g., Cornford and others, Reference Cornford, Martin, Lee, Payne and Ng2016; DeConto and Pollard, Reference DeConto and Pollard2016; DeConto and others, Reference DeConto2021; Edwards and others, Reference Edwards2021). These uncertainties essentially pertain to uncertainties in climate forcing and ice dynamics, parametric uncertainty (uncertainty due to the values assigned to model parameters) and the initialization method used (Seroussi and others, Reference Seroussi2019; Reference Seroussi2023; Reference Seroussi2024). The initial state of ice-sheet models is largely conditioned by the thermal conditions of the Antarctic ice sheet, particularly at the base. Because of the temperature-dependent ice rheology, temperate ice (ice at the pressure melting point that contains a small fraction of liquid water) deforms and flows faster, while temperate conditions at the base also control to a large extent subglacial processes, as the presence of subglacial meltwater lubricates the bed, facilitates basal sliding, and promotes the formation of fast-flowing glaciers. Fast-flowing ice streams drain about 90% of the grounded ice (Bennett, Reference Bennett2003), and their acceleration may lead to an increased grounded ice discharge (Bell, Reference Bell2008). While cold basal conditions help maintain the ice sheet configuration by inhibiting fast ice flow, the thawing of frozen bed patches could provoke an increased ice mass loss (Dawson and others, Reference Dawson, Schroeder, Chu, Mantelli and Seroussi2022).
Despite its critical role in ice sheet dynamics, direct temperature measurements are mainly limited to temperature profiles retrieved from deep boreholes (depth of ∼1000–3700 m) (e.g., Engelhardt, Reference Engelhardt2004). Radar surveys (e.g., Carter and others, Reference Carter, Blankenship, Young and Holt2009; Dawson and others, Reference Dawson, Schroeder, Chu, Mantelli and Seroussi2024) or data derived from satellite observations (e.g., Macelloni and others, Reference Macelloni, Leduc-Leballeur, Montomoli, Brogioni, Ritz and Picard2019) can help constrain the thermal state of the Antarctic ice sheet. Nevertheless, direct observations and measurements present a limited coverage of the continent, and numerical modeling remains necessary to estimate the distributions of basal temperatures and basal melt rates for the entire ice sheet. While recent work on the basal thermal state of the Greenland ice sheet is available (e.g., MacGregor and others, Reference MacGregor2016; Reference MacGregor2022; Karlsson and others, Reference Karlsson2021), studies focusing on the thermal state and basal melting of the Antarctic ice sheet remain sparse, outdated (e.g., Wilch and Hughes, Reference Wilch and Hughes2000; Llubes and others, Reference Llubes, Lanseau and Rémy2006; Pattyn, Reference Pattyn2010; Van Liefferinge and Pattyn, Reference Van Liefferinge and Pattyn2013), or focus on a regional scale (e.g., Zhao and others, Reference Zhao, Gladstone, Warner, King, Zwinger and Morlighem2018; Kang and others, Reference Kang, Zhao, Wolovick and Moore2022; Huang and others, Reference Huang, Zhao, Wolovick, Ma and Moore2024). Others do not necessarily take into account the uncertainty in initial and boundary conditions or the inherent uncertainty of the ice sheet model used (e.g., Dawson and others, Reference Dawson, Schroeder, Chu, Mantelli and Seroussi2022; Park and others, Reference Park, Jin, Morlighem and Lee2024).
Englacial temperatures result from a complex interplay between thermodynamic processes (heat conduction, advection) and heat sources (strain heating, geothermal heat flow, frictional heat from basal sliding, latent heat variations from melting and refreezing), where geothermal heat flow (GHF) probably represents the largest source of uncertainty for both Greenland (Rezvanbehbahani and others, Reference Rezvanbehbahani, Stearns, Van der Veen, Oswald and Greve2019; Zhang and others, Reference Zhang2024) and Antarctic (Pattyn, Reference Pattyn2010) ice sheets. Geothermal heat flow underneath ice sheets has been estimated using different methods, such as seismically-derived models (e.g., Shapiro and Ritzwoller, Reference Shapiro and Ritzwoller2004; An and others, Reference An2015; Shen and others, Reference Shen, Wiens, Lloyd and Nyblade2020), magnetically-derived models using airborne magnetic surveys (Martos and others, Reference Martos2017) or satellite geomagnetic data (Fox-Maule and others, Reference Fox-Maule, Purucker, Olsen and Mosegaard2005; Purucker, Reference Purucker2013), multivariate analysis (Stål and others, Reference Stål, Reading, Halpin and Whittaker2021), Machine learning with multiple geophysical and geological datasets (Lösing and Ebbing, Reference Lösing and Ebbing2021), and borehole measurements of temperature profiles (Mony and others, Reference Mony, Roberts and Halpin2020; Talalay and others, Reference Talalay2020). However, building large-scale and accurate geothermal heat flow datasets from field measurements is challenging, and results from geophysically or statistically derived models often diverge in magnitude and spatial distribution (Burton-Johnson and others, Reference Burton-Johnson, Dziadek and Martin2020; Reading and others, Reference Reading2022), though they can be understood as suited to distinct interdisciplinary use cases. Moreover, continental-scale geophysical models fail to capture small-scale anomalies in geothermal heat flow that may arise from various processes (Reading and others, Reference Reading2022), such as heterogeneous crustal heat production (e.g., Carson and others, Reference Carson, McLaren, Roberts, Boger and Blankenship2013) and thermal refraction due to subglacial topography and differences in conductivities between the bedrock and the ice (e.g., Willcocks and others, Reference Willcocks, Hasterok and Jennings2021).
While the geothermal heat flow defines basal temperature gradients and influences heat diffusion from the base, surface temperatures and accumulation rates control vertical advection from the surface (Pattyn, Reference Pattyn2010). The interactions between vertical diffusion and advection modulate the shape of vertical temperature profiles, particularly in slow-flowing regions (Pattyn, Reference Pattyn2010; Park and others, Reference Park, Jin, Morlighem and Lee2024), where vertical advection tends to cool down the ice, an effect that decreases with depth. Therefore, the geothermal heat flow required to reach the pressure melting point at the base depends on the ice thickness, surface temperature, and accumulation rate (Pattyn, Reference Pattyn2010).
Finally, frictional heat contributes significantly to basal melt, especially in regions that host fast-flowing glaciers (e.g., Larour and others, Reference Larour, Morlighem, Seroussi, Schiermeier and Rignot2012; Zhao and others, Reference Zhao, Gladstone, Warner, King, Zwinger and Morlighem2018; Karlsson and others, Reference Karlsson2021; Kang and others, Reference Kang, Zhao, Wolovick and Moore2022; Huang and others, Reference Huang, Zhao, Wolovick, Ma and Moore2024). Strong basal melt rates ranging from 10 to 500 mm a−1 are found beneath fast-flowing glaciers and tributaries feeding ice shelves, while slow-flowing regions experience lower basal melt rates ranging from 1 to 10 mm a−1 (e.g., Beem and others, Reference Beem, Jezek and van der Veen2010; Pittard and others, Reference Pittard, Roberts, Galton-Fenzi and Watson2016; Kang and others, Reference Kang, Zhao, Wolovick and Moore2022; Huang and others, Reference Huang, Zhao, Wolovick, Ma and Moore2024). Frictional heat reaches up to 2000 mW m−2 in the Lambert-Amery glacial system (Kang and others, Reference Kang, Zhao, Wolovick and Moore2022) and Totten Glacier (Huang and others, Reference Huang, Zhao, Wolovick, Ma and Moore2024) and is the primary contributor to the high basal melt of fast-flowing regions.
In this paper, we provide new estimates of basal thermal conditions and basal melting of the Antarctic ice sheet. We use different modeling approaches, various geothermal heat flow datasets, and the output of two regional climate models to conduct a sensitivity analysis of the basal thermal state of the Antarctic ice sheet. We further quantify the subglacial meltwater production of the grounded part of the ice sheet and evaluate the contribution of geothermal and frictional heat to basal melt. We compare our results with observational constraints, including temperature profiles measured in deep boreholes, englacial temperatures retrieved from the Soil Moisture and Ocean Salinity (SMOS) satellite observations (Macelloni and others, Reference Macelloni, Leduc-Leballeur, Montomoli, Brogioni, Ritz and Picard2019), and the distribution of subglacial lakes (Livingstone and others, Reference Livingstone2022). Finally, most likely boundary conditions are presented that best fit the observational constraints and model approaches.
2. Methods
2.1. Ice thermodynamics
Various methods exist to model the thermodynamics of ice sheets. As a first approach, we employ a conventional diffusion–advection equation to calculate the three-dimensional temperature field of the Antarctic ice sheet, following Huybrechts (Reference Huybrechts1992); Ritz (Reference Ritz1992), and Pattyn (Reference Pattyn2010), i.e.,
\begin{align}
\frac{\partial T}{\partial t} = \frac{K}{\rho_i c} \ \frac{\partial ^2 T} {\partial z^2} - w \ \frac{\partial T} {\partial z} - u \ \frac{\partial T}{\partial x} - v \ \frac{\partial T}{\partial y} + \frac{\Phi_d}{\rho_i c} \, ,
\end{align} where K is the ice thermal conductivity, c is the ice heat capacity, ρi = 917 kg m−3 is the ice density, 
$\mathrm{v}=(u,v,w)$ is the three-dimensional velocity field, and 
$\Phi_d$ corresponds to internal heat dissipation (strain heating). Depending on the approach (see below), (1) is solved either time-dependently or directly for the steady state. A numerical solution is obtained using a normalized sigma coordinate 
$\zeta = (h_s-z)/h$ (Pattyn, Reference Pattyn2010), where 
$h_s=b+h$ is the surface elevation (m), b is the bed elevation (m), h is the ice thickness (m), and z is the elevation in a Cartesian reference system (m). A common approximation in ice-sheet models consists of assuming constant thermodynamic parameters. However, thermal conductivity and heat capacity are both temperature-dependent. Some of our simulations use a commonly constant value for both parameters, i.e., K = 2.1 W m−1 K−1 and c = 2009 J kg−1 K−1, while others account for the temperature dependence according to 
$K = 3.101 \times 10^{8} \exp(-0.0057 \ T)$ J m−1K−1 a−1 (Ritz, Reference Ritz1987) and 
$c = 2115.3 + 7.79293 (T-273.15)$ J kg−1K−1 (Pounder, Reference Pounder1965). Finally, strain heating is given by:
\begin{align}
\Phi_d = \rho_i g h \nabla h_s \frac{\partial \mathrm{v_d}} {\partial z} \, ,
\end{align} where g = 9.81 m s−2 is the gravitational acceleration and 
$\partial \mathrm{v_d} / \partial z$ is the ice deformation rate. While the surface temperature defines the surface boundary condition, the basal boundary condition is given by a basal temperature gradient, depending on the basal heat sources, i.e.,
\begin{align}
\frac{\partial T_b}{\partial z} = \frac{G + F_b}{K} \, ,
\end{align} where G is the geothermal heat flow (W m−2), 
$F_b = \mathrm{v_b} \tau_b$ is the frictional heat with 
$\mathrm{v_b}$ the basal velocity and 
$\tau_b = \rho_i g h \nabla h_s$ the basal shear stress. Ice temperature is not allowed to exceed the pressure melting point, so that
\begin{align}
T_{pmp}= T_0 - \beta h \, ,
\end{align} where 
$T_0 = 273.15$ K is the absolute melting point of ice and 
$\beta = 8.66 \times 10^{-4}$ K m−1 is the Clausius–Clapeyron gradient. When the temperature reaches the pressure melting point, a temperate layer is allowed to form, and the excess heat is used for basal melting (Cuffey and Paterson, Reference Cuffey and Paterson2010), i.e.,
\begin{align}
\dot{b} = \frac{1}{\rho_i L} \left(G + F_b - K \frac{\partial T_{bc}}{\partial z}\right) \, ,
\end{align} where 
$L = 3.34 \times 10^{-5}$ J kg −1 is the latent heat of fusion, and 
$\partial T_{bc}/\partial z$ is the basal temperature gradient corrected for the pressure melting point.
However, the so-called cold-ice method fails to conserve energy in temperate ice, as variations in latent heat and water content are not accounted for. Alternatively, polythermal approaches (e.g., Greve, Reference Greve1997a; Reference Greve1997b) or enthalpy formulations (e.g., Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012; Kleiner and others, Reference Kleiner, Rückamp, Bondzio and Humbert2015; Hewitt and Schoof, Reference Hewitt and Schoof2017) have been developed to ensure energy conservation in the presence of polythermal conditions. Here, we compare the cold-ice method (solving equation (1)) with the simplified enthalpy gradient method (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012; Blatter and Greve, Reference Blatter and Greve2015; Kleiner and others, Reference Kleiner, Rückamp, Bondzio and Humbert2015; Greve and Blatter, Reference Greve and Blatter2016; Hewitt and Schoof, Reference Hewitt and Schoof2017). The three-dimensional enthalpy field E (J kg −1) can be computed from:
\begin{align}
\frac{\partial E}{\partial t} = \frac{\partial}{\partial z} \left(k_{c,t}\frac{\partial E}{\partial z} \right) - w \ \frac{\partial E} {\partial z} - u \ \frac{\partial E}{\partial x} - v \ \frac{\partial E}{\partial y} + \frac{\Phi_d}{\rho_i} - \frac{\rho_w}{\rho_i} L D(\omega) \, ,
\end{align} were 
$k_{c,t}=K_c=K/(c \rho_i)$ and 
$k_{c,t}=K_0=K_c \times 10^{-5}$ is the diffusivity in cold 
$(E \lt E_{pmp})$ and temperate ice 
$(E \geqslant E_{pmp})$, respectively, 
$E_{pmp}=c(T_{pmp}-T_{ref})$ is the enthalpy of ice at the pressure melting point without water content, and 
$T_{ref}=223.15$ K is a reference temperature. The term 
$- (\rho_w / \rho_i)L D(\omega)$, where 
$D(\omega)$ is a drainage function (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012), is added to prevent unrealistic high water content (Greve, Reference Greve1997b; Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012; Greve and Blatter, Reference Greve and Blatter2016; Hewitt and Schoof, Reference Hewitt and Schoof2017). We explicitly enforce the continuity condition on the cold-temperate transition surface (CTS), considering the melting conditions on the CTS only, following the ENTM scheme of Blatter and Greve (Reference Blatter and Greve2015) and Greve and Blatter (Reference Greve and Blatter2016). This is achieved through a two-step procedure. First, the enthalpy field is calculated for the entire vertical ice column, and the CTS is defined as the highest grid point of the temperate layers. Then, a corrector step is applied for the cold part of each ice column that possesses a CTS by applying the continuity condition as a boundary condition at the CTS. The final enthalpy profiles combine the corrected profiles for the cold part and the initial estimates for the temperate part. The temperature T and the water content ω are then diagnostically determined from the modeled enthalpy field:
\begin{align}
T =
\begin{cases}
E/c + T_{ref} \, , & E \lt E_{pmp}\, , \\
T_{pmp} \, , & E \geqslant E_{pmp}\, .
\end{cases}
\end{align}
\begin{align}
\omega =
\begin{cases}
0 \, , & E \lt E_{pmp} \, , \\
(E-E_{pmp})/L \, , & E \geqslant E_{pmp} \, .
\end{cases}
\end{align} The surface boundary condition is given by 
$E_s = c (T_s - T_{ref})$. The basal boundary conditions are determined from the decision chart proposed by Aschwanden and others Reference Aschwanden, Bueler, Khroulev and Blatter2012:
\begin{align}
\begin{cases}
K_c \partial E_b / \partial z = G + F_b \, , & E_b \lt E_{pmp} \text{\ and \ } H_w = 0 \, , \\
E_b = E_{pmp} \, , & E_b \lt E_{pmp} \text{\ and \ } H_w \gt 0 \, , \\
E_b = E_{pmp} \, , & E_b \geqslant E_{pmp} \text{\ and \ } H_w \gt 0 \text{\ and \ } H_t = 0 \, , \\
K_0 \partial E_b / \partial z = 0 \, , & E_b \geqslant E_{pmp} \text{\ and \ } H_w \gt 0 \text{\ and \ } H_t \gt 0 \, , \\
\end{cases}
\end{align} where 
$0 \leqslant H_w \leqslant 2$ is the thickness of the water layer and Ht is the thickness of the temperate layer. The water layer thickness is determined by:
\begin{align}
\frac {\partial H_w}{\partial t} = \dot{b} - C_d \, ,
\end{align} where 
$C_d = 0.001$ m a−1 is a constant drainage rate. The basal melt rate is computed from:
\begin{align}
\dot{b} = \frac{G + F_b - q_i}{(1-\omega) \rho_i L} \, ,
\end{align} where 
$q_i = K_c \partial E_{b} / \partial z$ is the heat flux in the ice (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012; Kleiner and others, Reference Kleiner, Rückamp, Bondzio and Humbert2015). In contrast to (5), (11) allows refreezing as long as 
$H_w \gt 0$. Furthermore, the basal melt increases with the drainage of excess englacial meltwater production (Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012).
2.2. Subglacial water routing
Subglacial water flow is calculated following the method of Le Brocq and others (Reference Le Brocq, Payne, Siegert and Alley2009) and the continuity equation for the subglacial water flow can be expressed as (Pattyn, Reference Pattyn2010; Kazmierczak and others, Reference Kazmierczak, Sun, Coulon and Pattyn2022):
\begin{align}
\frac{\partial d_w}{\partial t} = \dot{b} - \nabla (\mathrm{v}_w d_w) \, ,
\end{align} where dw is the thickness of the water film (m), 
$\dot{b}$ is the basal melt rate (positive for melting), 
$\mathrm{v}_w = d^{2}_w \nabla \phi / (12 \mu)$ is the vertically integrated velocity of water in the water film, 
$\mu = 1.8 \times 10^{-3}$ Pa s is the viscosity of water, 
$\phi = \rho_w g b + \rho_i g h - N$ is the hydraulic potential, 
$N = p_i - p_w$ is the effective pressure, 
$p_i = \rho_i g h$ is the pressure exerted by the overlying ice and pw is the subglacial water pressure, with 
$\rho_w = 1000$ kg m−3 the density of water. In steady state, the basal melt rate 
$\dot{b}$ must balance the water flux divergence 
$\nabla (\mathrm{v}_w d_w)$ (Pattyn and others, Reference Pattyn, De Brabander and Huyghe2005; Pattyn, Reference Pattyn2010). A balance flux approach is used to calculate the subglacial water flux 
$\Psi_w$ (e.g., Budd and Warner, Reference Budd and Warner1996; Le Brocq and others, Reference Le Brocq, Payne and Siegert2006; Pattyn, Reference Pattyn2010) by integrating the basal melt rate in the direction of the hydraulic potential gradient (Pattyn, Reference Pattyn2010). Assuming that the basal water pressure pw is equal to the pressure exerted by the overlying ice pi, the effective pressure vanishes (N = 0) and the gradient of the hydraulic potential can be written as:
\begin{align}
\nabla \phi = {\rho_i} g \nabla h_s + (\rho_w - \rho_i) g \nabla b\, .
\end{align}2.3. Ice flow models
The Kori-ULB ice flow model (follow-up of f.ETISh; Pattyn, Reference Pattyn2017) is a vertically integrated and thermomechanically coupled ice flow model of intermediate complexity and has been used in recent projections of the contribution of the Antarctic ice sheet to sea level rise (Edwards and others, Reference Edwards2021; Coulon and others, Reference Coulon2024). For the purpose of this paper, the model is initialized using a two-step procedure. First, the model is run over the grounded part of the ice sheet using the Shallow Ice Approximation (SIA; Hutter, Reference Hutter1983; Pattyn, Reference Pattyn2017). Thermomechanical coupling is introduced using an Arrhenius relationship, following Ritz (Reference Ritz1987, Reference Ritz1992) that accounts for the shape of the vertical profile of the horizontal velocities, i.e.,
\begin{align}
A(T) = \frac{1}{2} A^{*} \ \exp{\left[ \frac{Q}{R} \left( \frac{1}{T_{m}}-\frac{1}{T_b} \right) \right]}
\end{align} where A is the flow rate factor, 
$A^{*}= 1.66 \times 10^{-16}$ if 
$T^{*} \ge 266.65$ K and 
$A^{*}= 2.00 \times 10^{-16}$ if 
$T^{*} \lt 266.65$ K, 
$Q= 78.20 \times 10^{3}$ J mol−1 if 
$T^{*} \ge 266.65$ K and 
$Q= 95.45 \times 10^{3}$ J mol−1 if 
$T^{*} \lt 266.65$ K, R = 8.314 J kg−1 mol−1, 
$T^{*}=T-T_m$ is the homologous basal temperature, 
$T_m = \beta \zeta h$ is the pressure melting point correction, and Tb is the basal temperature. When using the enthalpy gradient method, we also consider the dependence of A on the water content ω (Duval, Reference Duval1977; Greve and Blatter, Reference Greve and Blatter2009; Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012):
\begin{align}
A(\omega) = A(T=T_{pmp})(1 + 181.25 \omega) \, .
\end{align}Since horizontal velocities are vertically integrated, shape functions are used to determine the vertical profiles of horizontal velocities. We take into account the temperature effects by adopting the method of Lliboutry (Reference Lliboutry1979) and Ritz (Reference Ritz1987, Reference Ritz1992), which consider increased strain rates due to warmer and more deformable ice at the base, i.e.,
\begin{align}
\mathrm{v}(\zeta)= \frac{p_r+2}{p_r+1} \times (1-\zeta^{p_r+1}) \mathrm{v} + \mathrm{v_b}
\end{align} where 
$p_r = n-1+ (Q G_0 h)/(R T_{b}^{2})$, 
$\mathrm{v} = (\dot{\gamma}_b h)/(p_{r}+2)$, and 
$\dot{\gamma}_b = B_0 \tau_d^{n} \exp \left[ \frac{Q}{R} \left(\frac{1}{T_{bc}}-\frac{1}{T_b} \right) \right]$.
Vertical velocities w are calculated with an analytical expression from horizontal velocities (Hindmarsh, 1999; Hindmarsh and others, Reference Hindmarsh, Gwendolyn and Parrenin2009), using the conservation of mass and the incompressibility of ice:
\begin{eqnarray}
w(\zeta) &=& -\dot{a} \left[1 - \frac{p_r+2}{p_r+1} (1 - \zeta) +\frac{1}{p_r+1} (1 - \zeta)^{p_{r}+2} \right] \nonumber\\
&& -\dot{b}+ \mathrm{v} \nabla b + (1-\zeta) \mathrm{v} \nabla h
\end{eqnarray} where 
$\dot{a}$ is the accumulation rate (negative for melting), and 
$\dot{b}$ is the basal melt rate (positive for melting). Basal velocities are estimated using a Weertman sliding law:
\begin{align}
\mathrm{v_b}= A^{'}_{b} \ |\tau_d|^{m-1} \ \tau_d
\end{align} where m = 3 is the sliding law exponent, and 
$A'_b$ are spatially varying basal sliding coefficients, which are iteratively determined by minimizing the difference between the observed and modeled ice thickness (Pollard and DeConto, Reference Pollard and DeConto2012).
The second step in the initialization consists of running the model in hybrid mode with evolving ice shelves (HySSA; Bueler and Brown, Reference Bueler and Brown2009; Winkelmann and others, Reference Winkelmann2011). In this second phase, melt and accretion rates under the floating ice shelves are determined following the method of Bernales and others (Reference Bernales, Rogozhina, Greve and Thomas2017). Alternatively, basal velocities are calculated from the difference between observed surface velocities and modeled deformational velocities, similar to Karlsson and others (Reference Karlsson2021), i.e.,
\begin{align}
\mathrm{v_b} = \mathrm{max}(0,\mathrm{v_{obs}} - \mathrm{v_{d,SIA}})
\end{align} where 
$\mathrm{v_b}$ is the basal velocity, 
$\mathrm{v_{obs}}$ is the observed surface velocity, and 
$\mathrm{v_{d,SIA}}$ is the SIA deformational velocity. This approach is used in combination with a fixed ice sheet geometry instead of the optimization with the two-step initialization approach. We refer to simulations using the optimization of basal slip coefficients as the ”Kori-ULB Opt approach” and simulations using the observed surface velocity to calculate basal sliding as the ”Kori-ULB Obs approach”.
We compare the results of the Kori-ULB ice flow model with the static hybrid ice sheet/ice stream model of Pattyn (Reference Pattyn2010) and Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) to infer the subglacial conditions of the Antarctic ice sheet for the observed ice sheet geometry. Contrary to the Kori-ULB ice flow model, the static model uses a balance flux approach (e.g., Budd and Warner, Reference Budd and Warner1996; Fricker and others, Reference Fricker, Warner and Allison2000; Le Brocq and others, Reference Le Brocq, Payne and Siegert2006) to determine the velocity field for a given surface mass balance and ice sheet geometry (Pattyn, Reference Pattyn2010). Another difference is that the shape functions for the vertical profiles of the velocity field are taken for an isothermal case, i.e., 
$p_r = n$. While the SIA is employed across the whole grounded ice domain, SSA is applied for ice flow speeds larger than 100 m a−1 and across large subglacial lakes (Pattyn, Reference Pattyn2010). Basal velocities are computed using a Weertman sliding law. Finally, geothermal heat flow, surface temperatures, and surface accumulation rates were locally adjusted using data from deep borehole measurements and geothermal heat flow was further calibrated with the presence of subglacial lakes that suggest the basal ice to be at pressure melting point. Here, we employ both the static corrected/calibrated and uncorrected/uncalibrated version of the model.
2.4. Ensemble of simulations
We carry out 180 simulations using different modeling approaches (Table 1) and various combinations of input datasets for the geothermal heat flow and surface climatic conditions. We use the Kori-ULB Opt approach to compare the enthalpy gradient method (solving Eq. (6); Kori-ULB Enth) with the cold-ice method (solving Eq. (1); Kori-ULB Opt) and employ the latter to assess the sensitivity to ice-dynamical approximations, ensuring comparable results between the Kori-ULB ice flow model and the hybrid ice sheet/ice stream model (Pattyn, Reference Pattyn2010). All simulations are carried out at a spatial resolution of 8 km on a regular spaced grid, using a nonlinear vertical discretization of 21 points with the thinner layers at the base to account for higher strain rates close to the bed (Huybrechts, Reference Huybrechts1992; Pattyn, Reference Pattyn2010). We use the data from MEaSUREs for the bed topography, surface elevation, and ice thickness (Bedmachine v3; Morlighem and others, Reference Morlighem2020; Morlighem, Reference Morlighem2022) and the observed surface velocities (Phase-Based Antarctica Ice Velocity Map v1; Mouginot and others, Reference Mouginot, Rignot and Scheuchl2019a; Reference Mouginot, Rignot and Scheuchl2019b). Datasets for geothermal heat flow are due to Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004), Fox-Maule and others (Reference Fox-Maule, Purucker, Olsen and Mosegaard2005), Purucker (Reference Purucker2013), An and others (Reference An2015), Martos and others (Reference Martos2017), Shen and others (Reference Shen, Wiens, Lloyd and Nyblade2020), Lösing and Ebbing (Reference Lösing and Ebbing2021), Stål and others (Reference Stål, Reading, Halpin and Whittaker2021), and Haeger and others (Reference Haeger, Petrunin and Kaban2022). Surface temperatures and accumulation rates are prescribed from the outputs of two regional climate models: MAR v3.11 (Kittel and others, Reference Kittel2021) and RACMO2.3p2 (van Wessem and others, Reference van Wessem2018). The variability among the GHF datasets and the differences between the estimates of MAR and RACMO are presented in Figures S1 and S2 of the supplementary material. In this study, we consider multiple geothermal heat flow datasets to cover a range of possible scenarios. However, the geophysical community is increasingly moving toward employing results from multivariate analysis with robust uncertainty bounds (e.g., Stål and others, Reference Stål, Reading, Halpin and Whittaker2021) for ensemble-based ice sheet studies (Reading and others, Reference Reading2022). Nevertheless, one limitation of our multi-model approach is the need to select representative geothermal heat flow maps for each ice sheet model, making it challenging to fully sample within the uncertainty bounds associated with each geothermal heat flow dataset.
Table 1. Summary of the modeling approaches, i.e., ice flow approximations, thermodynamics and thermomechanical coupling, basal sliding, ice sheet geometry, and calibration of input data.
2.5. Evaluation with observational constraints
To evaluate the ensemble of simulations, observational constraints are used and the misfit between the modeled and observed temperatures is quantified through the root mean square error (RMSE) for grid cells where observations are available. A first constraint is the borehole temperature profiles associated with deep ice cores in Antarctica. These are South Pole (Price and others, Reference Price2002); Vostok (Salamatin and others, Reference Salamatin, Lipenkov and Blinov1994; Parrenin and others, Reference Parrenin, Rémy, Ritz, Siegert and Jouzel2004); Kohnen (Wilhelms and others, Reference Wilhelms2014); Dome F (Fujii and others, Reference Fujii2002; Hondoh and others, Reference Hondoh, Shoji, Watanabe, Salamatin and Lipenkov2002); Dome C (Parrenin and others, Reference Parrenin2007); Law Dome (Dahl-Jensen and others, Reference Dahl-Jensen, Morgan and Elcheikh1999; Van Ommen and others, Reference Van Ommen, Morgan, Jacka, Woon and Elcheikh1999); Taylor Dome (G. Clow and E. Waddington, personal communication 2008); Siple Coast (MacGregor and others, Reference MacGregor, Winebrenner, Conway, Matsuoka, Mayewski and Clow2007); Byrd (Gow and others, Reference Gow, Ueda and Garfield1968); and WAIS divide (Cuffey and Clow, Reference Cuffey and Clow2014). Borehole temperature measurements are the most reliable direct observations of the temperature field, with an accuracy ranging from 
$0.0053^{\circ}\mathrm{C}$ (WAIS Divide) to 
$0.1^{\circ}\mathrm{C}$ (Byrd Station) (Talalay and others, Reference Talalay2020). However, these are local measurements whose representativeness of regional-scale thermal conditions is inherently limited.
A second constraint is given by the passive microwave L-band from the Soil Moisture and Ocean Salinity (SMOS) satellite, which has the advantage of exhibiting high penetration depths. Here, we use the englacial temperature profiles retrieved from SMOS observations with glaciological and microwave emission models by Macelloni and others (Reference Macelloni, Leduc-Leballeur, Montomoli, Brogioni, Ritz and Picard2019). The temperature profiles were derived using the Robin model (Robin, Reference GdQ1955), which provides reliable results as long as horizontal advection is negligible. The authors introduced a quality flag to evaluate the reliability of the retrieved temperatures based on two criteria: (i) the balance velocity (< 5 m a−1) and (ii) the minimization of the cost function. Here, we only use the temperature retrieved from SMOS observations where this quality flag is zero (Figure 6 in Macelloni and others, Reference Macelloni, Leduc-Leballeur, Montomoli, Brogioni, Ritz and Picard2019). However, the reliability of the temperature retrievals decreases in the lower part of the ice sheet, notably when the ice is thick, as the brightness temperature and L-band sensitivity decrease below 1500 m from the ice surface (Macelloni and others, Reference Macelloni, Leduc-Leballeur, Montomoli, Brogioni, Ritz and Picard2019).
The last measure is the spatial distribution of subglacial lakes, witnessing subglacial conditions at pressure melting point (e.g., Pattyn, Reference Pattyn2010; Van Liefferinge and Pattyn, Reference Van Liefferinge and Pattyn2013; Van Liefferinge and others, Reference Van Liefferinge2018; Karlsson and others, Reference Karlsson2021). We employ the inventory of Livingstone and others (Reference Livingstone2022) that counts 675 subglacial lakes. However, the presence of subglacial lakes as a quality measure is biased toward warmer subglacial conditions, as it is an end-member of basal temperatures.
3. Results
3.1. Subglacial conditions
Figure 1a presents the ensemble mean basal temperatures relative to the pressure melting point. The pie plots indicate the fraction of temperate basal conditions, i.e., the proportion of the ice sheet where basal temperatures reach the pressure melting point. Results show that temperate or near-temperate basal conditions prevail in Antarctica, with, on average, 60 (min: 36.5 – max: 71.1) % of basal temperatures reaching the pressure melting point over the grounded part of the ice sheet. Temperate basal conditions represent 61 (32–74) % of the subglacial conditions of the East Antarctic Ice Sheet, 64 (48–80) % of the West Antarctic Ice Sheet, and only 23 (12–40) % of the Antarctic Peninsula. Cold basal conditions prevail in regions where the ice is thin, such as mountain ranges (e.g., Transantarctic Mountains, Gamburtsev Subglacial Mountains, coastal mountain ranges of Dronning Maud Land, Antarctic Peninsula) and in the Wilkes Subglacial Basin. Figure 1c further shows the prevalence of temperate basal conditions for 27 subglacial basins delineated according to Zwally and others (Reference Zwally, Giovinetto, Beckley and Saba2012). Temperate basal conditions characterize 75 to 79 % of basal conditions for basins 18, 19, 21, 22 of the West Antarctic Ice Sheet, which includes Thwaites and Pine Island glaciers and the Siple Coast region. In contrast, the East Antarctic Ice Sheet exhibits more heterogeneous basal conditions. Here, a higher variability among the simulated basal temperatures is found (Figure 1b). Finally, Figure 1d presents the likely basal thermal state of the Antarctic ice sheet estimated from the level of agreement between the ensemble of simulations in predicting cold or temperate basal conditions. A likelihood of 100% indicates that all simulations predict either temperate (dark red) or cold (dark blue) basal conditions. The level of disagreement increases as the values deviate from 100%, revealing that some simulations predict temperate and others cold basal conditions for a given location. Simulations generally agree on whether the base is likely cold or temperate at a given location, even when the standard deviation is high.
Figure 1. The basal thermal state of the Antarctic ice sheet. (a) Ensemble mean basal temperatures relative to the pressure melting point (
${}^{\circ}\mathrm{C}$). Values below 
$-20^{\circ}\mathrm{C}$ are truncated. The pie plots show the fraction of temperate (red) and cold (blue) basal thermal conditions for the entire ice sheet, East Antarctic Ice Sheet (EAIS), West Antarctic Ice Sheet (WAIS), and Antarctic Peninsula (AP). (b) Standard deviation from the ensemble mean (
${}^{\circ}\mathrm{C}$). (c) Fraction of temperate basal conditions (%) at the scale of subglacial basins delineated and enumerated according to Zwally and others (Reference Zwally, Giovinetto, Beckley and Saba2012). (d) The likely basal thermal state of the Antarctic ice sheet showing the percentage of agreement between simulations predicting cold or temperate basal conditions at a given location.
Figure 2a shows the ensemble mean basal melt rates and the resultant subglacial meltwater volume produced over the grounded ice sheet. The pie plots display the relative contribution of geothermal and frictional heat, which are the two main contributors to basal melt in Antarctica. The mean basal melt rate is found to be 6.9 (4.1–9.1) mm a−1. However, basal melt rates > 20 mm a−1 occur in West Antarctica and over the ice sheet margins, coinciding with high geothermal heat flow or fast-flowing glaciers. The ensemble mean subglacial meltwater volume is 69.3 (45.1–95.6) Gt a−1 over the grounded part of the Antarctic ice sheet, with 51.4 (30.3–76.6) Gt a−1 coming from East Antarctica, 16.2 (11.3–22.8) Gt a−1 from West Antarctica, and 1.8 (1.2–2.6) Gt a−1 from the Antarctic Peninsula. The overall contribution of geothermal and frictional heat is of the same order, with 51 (16–67) % of basal melt being due to geothermal heat and 49 (33–84) % to frictional heat. Figure 2c shows that the relative contribution of geothermal heat is the highest in slow-flowing regions inland of East Antarctica. The relative contribution of frictional heat increases toward the margins of the ice sheet and coincides with the occurrence of high basal melt rates. The subglacial basins contributing the most to the subglacial meltwater production are the basins with the largest surface area located in East Antarctica, i.e., basins 13 (9 Gt a−1), 17 (7.7 Gt a−1), and 3 (6.2 Gt a−1; Figure 2c). The standard deviation (Figure 2b) is higher over the ice sheet margins, where more variability among the estimates results from using different calculations of basal sliding in the model approaches, which modulate the contribution of frictional heat to basal melt. Finally, Figure 2d shows the ensemble mean subglacial water flux, computed from the basal melt rate estimates, following the method of Le Brocq and others (Reference Le Brocq, Payne, Siegert and Alley2009). Subglacial water flux follows the general ice flow pattern, with significant water flow being concentrated over the margins of the major subglacial drainage systems, as shown in Pattyn (Reference Pattyn2010). All simulations reproduce the same pattern.
Figure 2. The subglacial meltwater production beneath the Antarctic ice sheet. (a) Ensemble mean basal melt rate (mm a−1). Values above 20 mm a−1 and below 0 mm a−1 are truncated. The pie plots show the relative contribution of geothermal heat (dark blue) and frictional heat (light blue) to basal melt for the entire ice sheet, East Antarctic Ice Sheet (EAIS), West Antarctic Ice Sheet (WAIS), and Antarctic Peninsula (AP). (b) Standard deviation from the ensemble mean (mm a−1). (c) Ensemble mean relative contribution of geothermal and frictional heat to basal melting (%). Numbers indicate the subglacial meltwater volume (Gt a−1) integrated over each subglacial basin delineated as in Figure 1c. (d) Ensemble mean subglacial water flux (103 m2 a−1) computed following the method of Le Brocq and others (Reference Le Brocq, Payne, Siegert and Alley2009). The contribution of the englacial meltwater drained to the bed is neglected in panels a–c to ensure consistency between simulations using enthalpy and cold-ice methods but is accounted for in the calculation of the subglacial water flux shown in panel d.
While Figure 2a–c shows the basal melt resulting from the contribution of geothermal and frictional heat, simulations conducted with the enthalpy gradient method provide additional information about refreezing and englacial meltwater production. Figure 3 displays the mean basal melt and refreezing rates of simulations using the enthalpy method (a–b) and the excess englacial meltwater drained to the bed (c). Although refreezing areas exist, the average refreezing rate is low (less than 1 Gt a−1). Nevertheless, the subglacial meltwater production is 59.5 (46.1–71.7) Gt a−1, which is lower than the ensemble mean value. The contribution of frictional heat to basal melt is 66 (53–84)% and prevails over the contribution of geothermal heat, which is only 34 (16–47)%. However, the resultant lower subglacial meltwater volume is compensated by the drainage of the excess englacial meltwater, reaching 30.2 (27.7–32.1) Gt a−1 and leading to a total volume of subglacial water of 89.7 (73.7–103.8) Gt a−1. The temperate ice layer, within which the englacial meltwater production occurs, presents an average thickness of 17 m but is thinner than 8.9 m in 80% of the temperate regions (Figure 3d).
Figure 3. Mean subglacial and englacial meltwater production over the Antarctic ice sheet for simulations conducted with the enthalpy gradient method. (a) Mean basal melt rate (mm a−1). Values above 20 mm a−1 are truncated. The pie plots show the relative contribution of geothermal heat (dark blue) and frictional heat (light blue) to basal melt for the grounded part of the ice sheet. (b) Mean refreezing rate (mm a−1). Values below -0.1 mm a−1 are truncated. (c) Englacial meltwater drained to the bed (mm a−1) and the resultant additional meltwater volume (Gt a−1). (d) Thickness of the temperate ice layer (m), within which englacial meltwater is produced. The numeric values in panels a and c indicate the volume of subglacial water produced by basal melting (a) and the volume of englacial meltwater drained to the bed (c) for each subglacial basin delineated as in Figure 1c.
3.2. Sensitivity analysis
Modeled basal temperatures and melt rates vary significantly depending on the model approach as well as boundary conditions. The detailed results from the ensemble of simulations are presented in Figures S3–S9 of the supplementary material, while Figure 4 provides an overview of the ensemble sensitivity.
Figure 4. Sensitivity of basal thermal conditions to the geothermal heat flow and modeling approaches. The figure shows (a,b) the fraction of temperate basal conditions (%), (c,d) the subglacial meltwater production (Gt a−1), and (e-h) the contribution of (e,f) geothermal and (g,h) frictional heat to basal melt (Gt a−1). The contribution of the englacial meltwater drained to the bed (30.2 ± 2 Gt a−1) is not included in the subglacial meltwater volume for the Kori-ULB Enth approach to ensure consistency with simulations conducted with the cold-ice method. Geothermal heat flow datasets are ordered from lowest to highest averaged GHF values, which are provided for convenience but should not be considered statistically robust averages (Figure S1).
Geothermal heat flow uncertainty has the largest impact on basal temperatures and basal melt rates, and explains 73% of the total variability in the fraction of temperate basal conditions and 42% in the subglacial meltwater fluxes (Figures 4; S3). The geothermal heat flow from Purucker (Reference Purucker2013) is associated with the coldest basal thermal conditions, with, on average, only 42 (36–51) % of temperate basal conditions and an integrated basal melt of 57 (45–73) Gt a−1 (Figures 4; S3–S5). The warmest basal conditions are obtained with the geothermal heat flow of Fox-Maule and others (Reference Fox-Maule, Purucker, Olsen and Mosegaard2005) and Martos and others (Reference Martos2017) (Figures 4, S3–S5). The sensitivity to the geothermal heat flow is more pronounced in slow-flowing regions of East Antarctica (Figures S4–S5; S8), but basal melt rates are influenced by the magnitude of the geothermal heat across the whole ice sheet (Figure S5). On average, the contribution of geothermal heat to basal melt ranges from 23 (9–35) to 43 (24–54) Gt a−1, corresponding to 41–56% of the total subglacial meltwater production (Figure 4).
The second source of uncertainty pertains to ice-sheet model approximations, which account for 23% of the total variability in the fraction of temperate basal conditions and 50% in the subglacial meltwater production (Figures 4; S3). In terms of thermodynamics, the distributions of basal temperatures and the estimates of temperate basal conditions are similar between the Kori-ULB Enth (solving Eq. 6) and the Kori-ULB Opt (solving Eq. 1) approaches (Figures 4; S3; S6). However, the enthalpy model generates lower basal melt rates inland of East Antarctica (Figure S7), leading to a decreased subglacial meltwater production of 20 Gt a−1 on average compared to the Kori-ULB Opt approach (Figures 4; S3). This difference is associated with a lower contribution of geothermal heat to basal melt (Figure S8) and can be ascribed to the distinct basal boundary conditions in both approaches (Eq. 3 and 9). Moreover, simulations using constant thermodynamic parameters systematically predict slightly colder basal thermal conditions (Figures 4; S3; S6–S7).
For simulations conducted with different ice flow models using the cold-ice method (solving Eq. 1), the variability in the subglacial meltwater production mainly results from different contributions of frictional heating (Figures 4; S3; S8). On average, the contribution of frictional heat to basal melt ranges from 24 to 45 Gt a−1, corresponding to 39–56% of the total subglacial meltwater production (Figure 4). The Kori-ULB Opt approach predicts the highest subglacial meltwater volume (80 Gt a−1 on average) while the uncalibrated model of Pattyn (Reference Pattyn2010) estimates the lowest one (64 Gt a−1 on average). Basal thermal conditions are also more homogeneous with the hybrid ice sheet/ice stream model than the Kori-ULB ice flow model, with notable differences in East Antarctica (Figures S6–S7). This can be assigned to the optimization method used for the latter approach. The calibrated model of Pattyn (Reference Pattyn2010) estimates a higher subglacial meltwater volume (70 Gt a−1 on average) than the uncalibrated model, resulting from a higher contribution of geothermal heat to basal melt (Figure S3). This is essentially due to the correction of GHF with subglacial lakes, allowing for some colder basal areas to become temperate. The subglacial meltwater volumes are similar for the Kori-ULB Obs and the calibrated model of Pattyn (Reference Pattyn2010) (Figure 4). However, the spatial distributions of basal melt rates and the relative contributions of frictional and geothermal heat to basal melt are different for the reasons mentioned above (Figures S3; S7–S8). While all approaches adequately reproduce the velocity field over the grounded part of the ice sheet (Figure S9), the Kori-ULB ice flow model presents a better fit with the observed surface velocities (RMSE between 35 and 63 m a−1) than the hybrid ice sheet/ice stream model (RMSE between 81 and 100 m a−1).
The remaining variability can be assigned to the differences between the surface climatic conditions from the regional climate models MAR and RACMO (Figure S2). Simulations using MAR predict slightly higher basal temperatures and basal melt rates, but overall, the sensitivity is very low between MAR and RACMO, except for the Kori-ULB Opt approach, where simulations using MAR predict higher contributions of frictional heat to basal melt (Figure S3).
3.3. Comparison with observational constraints
Evaluation of the modeled results can be done to a certain extent by using (i) temperature profiles measured in deep boreholes, (ii) temperature fields retrieved from SMOS satellite observations in slow-flowing areas (Macelloni and others, Reference Macelloni, Leduc-Leballeur, Montomoli, Brogioni, Ritz and Picard2019), and (iii) the distribution of 675 subglacial lakes (Livingstone and others, Reference Livingstone2022) considered as an indicator of temperate basal conditions. Figure 5 locates the observational constraints and compares the modeled temperature profiles with borehole measurements. Figure 6 summarizes the performance of the ensemble of simulations with respect to these observational constraints. The RMSE is taken here as a spatially averaged value across the different observed variables. The detailed RMSE values and spatially distributed misfits are given in Figures S10–S12 of the supplementary material.
Figure 5. Observational constraints used for the evaluation of the ensemble of simulations, which include the temperature profiles from borehole measurements, the presence of subglacial lakes, and the englacial temperature field derived from SMOS data. The figure further compares the modeled temperature profiles with borehole measurements. The ensemble mean RMSE for a given temperature profile is displayed at the top of the corresponding panel. Symbols represent simulations conducted with different modeling approaches. Colors differentiate simulations employing different geothermal heat flow datasets. Filled and empty symbols distinguish simulations using MAR and RACMO. The detailed legend is shown in Figure 6. The temperature profiles are also presented separately for each model approach in Figures S13–S17 of the supplementary material to highlight individual model behaviors.
Figure 6. Evaluation of the ensemble of simulations with the observational constraints. The figure shows the RMSE (
${}^{\circ}\mathrm{C}$) of the modeled temperatures for each ensemble member with respect to borehole temperature profiles, SMOS-derived englacial temperatures, and the presence of subglacial lakes. The lines inside the panels indicate the ensemble mean RMSE for the corresponding observational constraint. Symbols represent simulations conducted with different modeling approaches. Colors differentiate simulations employing different geothermal heat flow datasets. Filled and empty symbols distinguish simulations using MAR and RACMO.
Based on the borehole temperature measurements, overall good results (average RMSE of 
$2.8^{\circ}\mathrm{C}$) are obtained with the calibrated model of Pattyn (Reference Pattyn2010). This is quite obvious, as the observational constraint has been used to tune some of the parameters in the model. The lowest performance (RMSE of 
$7.8^{\circ}\mathrm{C}$) is obtained with the uncalibrated model of Pattyn (Reference Pattyn2010) using temperature-dependent thermodynamic parameters, the geothermal heat flow of Martos and others (Reference Martos2017), and the surface climatology of RACMO. As for non-calibrated approaches, both Kori-ULB Opt and Kori-ULB Obs score well, with a mean RMSE of 2.9 and 
$2.7^{\circ}\mathrm{C}$ respectively. The performance of the Kori-ULB Enth approach is similar to Kori-ULB Opt when using constant thermodynamic parameters but decreases when accounting for their temperature dependence. Overall, simulations using constant thermodynamic parameters perform better, leading to a lower RMSE of 
$0.7^{\circ}\mathrm{C}$ on average. Regarding boundary conditions, simulations using geothermal heat flow datasets associated with colder thermal conditions (e.g., Purucker, Reference Purucker2013; An and others, Reference An2015; Shen and others, Reference Shen, Wiens, Lloyd and Nyblade2020) present a better fit with measured temperature profiles (mean RMSE of 3.3–
$3.5^{\circ}\mathrm{C}$). Simulations using MAR present a lower average RMSE (
$3.5^{\circ}\mathrm{C}$) than those using RACMO (
$3.8^{\circ}\mathrm{C}$).
The ensemble mean RMSE with borehole measurements is 3.7 (2.3–7.8)
$^{\circ}\mathrm{C}$ (Figure 6). However, we find considerable variability depending on the location of the temperature profiles (Figure 5). For instance, the mean RMSE for the ensemble of simulations ranges from 1.6 (0.6–3.2)
$^{\circ}\mathrm{C}$ (Law Dome) to 6.6 (1.9–14.4)
$^{\circ}\mathrm{C}$ (WAIS divide), and misfits are generally more pronounced in West Antarctica (Figures S11–S12). While most of the ice flow models satisfactorily simulate the temperature profiles in slow-flowing regions of East Antarctica, simulations conducted with the hybrid ice sheet/ice stream model (Pattyn, Reference Pattyn2010) fail to reproduce the shape of the temperature profiles and estimate unrealistic high thickness of temperate ice layers, in particular at WAIS Divide and Byrd Station (Figures 5; S16–S17). Such overestimation of the temperate ice layers is not observed in simulations conducted with the Kori-ULB model (Figure S13–S15), which explains their better performance (Figure 6). Furthermore, surface temperatures from the regional climate models can be higher than those from borehole measurements, potentially leading to biased temperature profiles at some locations (e.g., WAIS Divide, Byrd Station, Taylor Dome).
The comparison with temperature profiles derived from SMOS observations leads to a mean RMSE of 4.8 (3.6–6.8)
$^{\circ}\mathrm{C}$ (Figure 6). The best fit (RMSE of 
$3.6^{\circ}\mathrm{C}$) is found with the Kori-ULB Enth approach using constant thermodynamic parameters, the geothermal heat flow of Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004), and the surface climatology of RACMO. Overall, the approaches using the Kori-ULB ice flow model fit better with the SMOS observations (mean RMSE of 4.3–
$5^{\circ}\mathrm{C}$). The uncalibrated model of Pattyn (Reference Pattyn2010) results in the poorest fit, especially in combination with the highest heat flow datasets (RMSE of 
$6.8^{\circ}\mathrm{C}$). However, the mean RMSE for the uncalibrated and calibrated model of Pattyn (Reference Pattyn2010) is very similar, so, the calibration is not necessarily leading to an overall improved fit of the temperature profiles. The largest differences between modeled temperatures and SMOS observations are near the Gamburtsev subglacial mountains and Lake Vostok (Figures S11–S12). In terms of geothermal heat flow, simulations using the datasets of Shapiro and Ritzwoller (Reference Shapiro and Ritzwoller2004), Purucker (Reference Purucker2013), An and others (Reference An2015), and Shen and others (Reference Shen, Wiens, Lloyd and Nyblade2020) lead to a better fit with SMOS data (mean RMSE of 4.5–
$4.8^{\circ}\mathrm{C}$). However, most of these datasets (e.g., Shapiro and Ritzwoller, Reference Shapiro and Ritzwoller2004) have low spatial resolution and represent only a broad envelope of GHF values, in contrast to more recent datasets that display much higher variability. Additionally, simulations using the surface climatic conditions from RACMO systematically present an RMSE lower of 
$\sim 0.5^{\circ}\mathrm{C}$. Although SMOS data provides insights into the englacial temperature field of the Antarctic ice sheet, one should note that its limited sensitivity and reduced accuracy in the deepest parts of the ice sheet may hinder our evaluation of ice sheet model performances.
The use of subglacial lakes as a goodness of fit is probably the least adequate evaluation method, as it is biased toward the pressure melting point limit of the ice. The best-performing simulation results from using the uncalibrated model of Pattyn (Reference Pattyn2010) with temperature-dependent thermodynamic parameters, the geothermal heat flow of Haeger and others (Reference Haeger, Petrunin and Kaban2022), and the regional climate model RACMO (RMSE of 
$2.2^{\circ}\mathrm{C}$). As expected, the choice of the geothermal heat flow has the most influence. Simulations predicting warmer basal conditions, i.e., using the geothermal heat flow from Fox-Maule and others (Reference Fox-Maule, Purucker, Olsen and Mosegaard2005), Martos and others (Reference Martos2017), and Haeger and others (Reference Haeger, Petrunin and Kaban2022) better fit with the observed subglacial lake distribution (mean RMSE 
$ \lt 3.3^{\circ}\mathrm{C}$), since a larger area of the bed is at pressure melting point. However, simulations conducted with the geothermal heat flow of Shen and others (Reference Shen, Wiens, Lloyd and Nyblade2020) and Haeger and others (Reference Haeger, Petrunin and Kaban2022) perform reasonably well regarding both subglacial lakes (average RMSE of 3.3 and 
$3^{\circ}\mathrm{C}$, respectively) and borehole measurements (average RMSE of 3.5 and 
$3.6^{\circ}\mathrm{C}$).
4. Discussion
Only a few previous studies focused on basal thermal conditions of the Antarctic ice sheet, such as Llubes and others (Reference Llubes, Lanseau and Rémy2006) and Pattyn (Reference Pattyn2010). Van Liefferinge and Pattyn (Reference Van Liefferinge and Pattyn2013) further updated the latter study, and Park and others (Reference Park, Jin, Morlighem and Lee2024) provided new estimates using the Ice-sheet and Sea-level System Model (ISSM) that includes an enthalpy formulation and assessed the sensitivity to geothermal heat flow and vertical velocity, compared against borehole measurements. In addition, other continent-scale modeling also calculated basal thermal conditions. For instance, Sutter and others (Reference Sutter2019) explored the influence of geothermal heat flow on the evolution of the Antarctic ice sheet across the mid-Pleistocene transition and the implications for Oldest Ice with the Parallel Ice Sheet Model (PISM). Our study combines the approaches of Pattyn (Reference Pattyn2010) and Park and others (Reference Park, Jin, Morlighem and Lee2024), using the Kori-ULB model with both enthalpy and cold-ice temperature calculations for present-day conditions, compared with borehole temperature measurements, SMOS satellite observations, and the distribution of subglacial lakes.
Overall, our basal temperature distributions are consistent with those of previous studies (e.g., Llubes and others, Reference Llubes, Lanseau and Rémy2006; Pattyn, Reference Pattyn2010; Van Liefferinge and Pattyn, Reference Van Liefferinge and Pattyn2013; Dawson and others, Reference Dawson, Schroeder, Chu, Mantelli and Seroussi2022; Park and others, Reference Park, Jin, Morlighem and Lee2024). While temperate basal conditions prevail in West Antarctica and beneath fast-flowing glaciers, cold basal conditions are mainly found in slow-flowing regions of East Antarctica with low geothermal heat flow (e.g., Wilkes Subglacial Basin) or when the ice is thin enough (e.g., subglacial mountain ranges). Our likely basal thermal state also corroborates the estimates of Pattyn (Reference Pattyn2010), except for the Wilkes and Victoria Lands, where our results indicate likely temperate basal conditions. Greater variability in basal temperature estimates occurs in cold-based regions, as the temperature is constrained by the pressure melting point of ice, preventing significant variations in temperate areas (Pattyn, Reference Pattyn2010). Because thicker ice favors temperate basal conditions through its insulating effect and influence on the pressure melting point of ice, regions with large standard deviations also correspond to those with a rough bed and lower ice thickness. In this study, we used the topography from Bedmachine v3 (Morlighem and others, Reference Morlighem2020; Morlighem, Reference Morlighem2022), although Bedmap3 (Pritchard and others, Reference Pritchard2025) was recently released, showing more defined troughs and resolved mountain features. Local topographic features could significantly impact high-resolution ice flow and temperature modeling but would likely not be captured at the 8 km resolution of our simulations.
Our basal melt rate distributions also align with previous estimates (Pattyn, Reference Pattyn2010; Park and others, Reference Park, Jin, Morlighem and Lee2024), with low basal melt rates inland of East Antarctica and high basal melt rates in West Antarctica and over the ice sheet margins, where fast-flowing glaciers generate significant frictional heat. Whereas frictional heat significantly contributes to basal melt in fast-flowing regions, geothermal heat has more impact on basal melt in slow-flowing regions, in agreement with the results of, e.g., Kang and others (Reference Kang, Zhao, Wolovick and Moore2022) and Huang and others (Reference Huang, Zhao, Wolovick, Ma and Moore2024). However, our ensemble mean basal melt rate of 6.9 mm a−1 is higher than estimated by Llubes and others (Reference Llubes, Lanseau and Rémy2006) and Pattyn (Reference Pattyn2010). While close to the result of Pattyn (Reference Pattyn2010), our ensemble mean subglacial meltwater volume of 69.3 Gt a−1 for the grounded ice sheet is more than twice the estimate of Park and others (Reference Park, Jin, Morlighem and Lee2024). These differences can be attributed to using different input datasets, notably for the geothermal heat flow, and the thermodynamical model approaches. Specifically, Park and others (Reference Park, Jin, Morlighem and Lee2024) use an enthalpy method to solve the thermodynamics of the ice sheet, whereas we employ both enthalpy and cold-ice thermodynamical approaches. Contrary to cold-ice methods, enthalpy (e.g., Aschwanden and others, Reference Aschwanden, Bueler, Khroulev and Blatter2012) and polythermal (e.g., Greve, Reference Greve1997b) methods are energy-conserving and, therefore, better represent the phase transition between cold and temperate ice and the processes occurring therein. Moreover, enthalpy methods account for refreezing rates, which may explain the lower integrated subglacial meltwater volume estimates in Park and others (Reference Park, Jin, Morlighem and Lee2024). Nevertheless, our simulations conducted with the enthalpy gradient method exhibit low refreezing rates and a mean subglacial meltwater production of 59.5 Gt a−1, which is lower than the values obtained with the cold-ice method but is still higher than the estimates of Park and others (Reference Park, Jin, Morlighem and Lee2024).
For any given ice flow model, the factor that most influences the distributions of basal temperatures and basal melt rates is the geothermal heat flow. Overall, simulations using geothermal heat flow datasets leading to colder thermal conditions (e.g., Purucker, Reference Purucker2013; An and others, Reference An2015; Shen and others, Reference Shen, Wiens, Lloyd and Nyblade2020) better fit the borehole measurements and SMOS satellite observations. Conversely, simulations using geothermal heat flow datasets leading to warmer thermal conditions (e.g., Fox-Maule and others, Reference Fox-Maule, Purucker, Olsen and Mosegaard2005; Martos and others, Reference Martos2017; Haeger and others, Reference Haeger, Petrunin and Kaban2022) better fit the subglacial lake distribution as they predict higher fractions of temperate basal conditions. Based on our ensemble modeling, simulations conducted using the geothermal heat flow of Shen and others (Reference Shen, Wiens, Lloyd and Nyblade2020) and Haeger and others (Reference Haeger, Petrunin and Kaban2022) present a reasonably good fit with both borehole temperature measurements and the distribution of subglacial lakes and may, therefore, represent a good compromise. However, other geothermal heat flow datasets could also perform well within their respective uncertainty (e.g., Supplementary Figure 3 in Reading and others, Reference Reading2022). In particular, multivariate empirical models (e.g., Lösing and Ebbing, Reference Lösing and Ebbing2021; Stål and others, Reference Stål, Reading, Halpin and Whittaker2021) provide robust uncertainty bounds or consider most of the available information related to geothermal heat flow (Reading and others, Reference Reading2022). To assess whether the geothermal heat flow models compare well within their uncertainty bounds, a substantially larger ensemble, complemented with emulators to explore the full uncertainty range, would be required and could be the focus of future work. Furthermore, we do not account for small-scale anomalies in the geothermal heat flow (e.g., Willcocks and others, Reference Willcocks, Hasterok and Jennings2021), which could be associated with hotspots and locally enhanced basal melt (e.g., McCormack and others, Reference McCormack2022), potentially corresponding to subglacial lakes presence. Accounting for these anomalies in higher-resolution simulations might thus improve the performance of simulations predicting colder basal thermal conditions. Future work could also benefit from exploiting radar observations and radiostratigraphically informed internal architecture of the Antarctic ice sheet as an additional constraint on the spatial variations in geothermal heat flow and basal melt estimates (Matsuoka and others, Reference Matsuoka, MacGregor and Pattyn2012; Bingham and others, Reference Bingham2024).
Englacial and basal thermal conditions are also affected by ice-thermodynamical approximations. While enthalpy and cold-ice methods simulate similar distributions of basal temperatures, differences in basal melt estimates and modeled temperature profiles are non-negligible in areas where temperate ice is found. Simulations conducted with the enthalpy model estimate thinner temperate ice layers compared to the cold-ice models (Figures S13–S17), in agreement with the results of Greve and Blatter (Reference Greve and Blatter2016). Cold-ice models may overestimate the amount of temperate ice owing to heat dissipation and inadequate treatment of the boundary between cold and temperate ice. Furthermore, the modeled temperature profiles are sensitive to the choice of the ice flow model, even when the same thermodynamics method is used. Contrary to the Pattyn (Reference Pattyn2010) model, the Kori-ULB ice flow model includes thermomechanical coupling, which accounts for the increased deformation in layers closest to the bed and allows for accurately modeling the shape of the temperature profiles despite the absence of calibration of input datasets. Thermomechanical coupling influences to a large extent the vertical velocity profiles, which have been shown to significantly affect the modeled temperature profiles (Park and others, Reference Park, Jin, Morlighem and Lee2024). In addition, the hybrid ice sheet/ice stream model, which uses uniform basal slip coefficients in the sliding law, simulates more homogeneous basal thermal conditions and a lower contribution of frictional heat to basal melt. This explains why the Kori-ULB Opt approach leads to higher meltwater production due to locally higher estimates of frictional heating.
Finally, our results suggest that simulations leading to colder basal thermal conditions show better agreement with borehole measurements and SMOS observations. However, this may be biased by overestimated surface temperatures in the regional climate models and the lack of thermal memory in our model simulations, as we assume a thermal steady state for the present-day conditions. The englacial temperature field reacts slowly to climatic changes and is likely still adjusting to transient effects from previous glacial and interglacial periods (Pattyn, Reference Pattyn2010). Ritz (Reference Ritz1987) suggested that accounting for transient effects could lead to basal temperatures that are 2K lower than those estimated under steady-state conditions. Accounting for the climatic history with a spin-up over glacial-interglacial cycles would likely lead to more realistic englacial temperature distributions (Huybrechts, Reference Huybrechts2002).
5. Conclusion
We provide new estimates of thermal conditions and subglacial meltwater production of the Antarctic ice sheet. These estimates are based on an ensemble of simulations accounting for uncertainties arising from the geothermal heat flow, present-day surface climatic conditions, ice thermodynamics, and ice-dynamical approximations. The modeled temperature fields are compared with borehole temperature measurements, SMOS satellite observations, and the subglacial lake distribution.
Our results reveal warmer basal thermal conditions and higher basal melt than previously reported, with a fraction of temperate basal conditions of 60 (36.5–71.1)%, a mean basal melt rate of 6.9 (4.1–9.1) mm a−1, and a subglacial meltwater volume of 69 (45.1–95.6) Gt a−1 over the grounded part of the ice sheet. We find, in agreement with previous studies, that likely temperate basal conditions and high basal melt rates prevail in West Antarctica and beneath fast-flowing glaciers, while likely cold basal conditions are mainly found in slow-flowing regions or where the ice is thin (e.g., subglacial mountain ranges). The relative contribution of geothermal and frictional heat to basal melt is 51 (16–67)% and 49 (33–84)%, respectively. In addition, the englacial meltwater drained to the bed reaches 30.2 (27.7–32.1) Gt a−1.
The geothermal heat flow primarily controls the distributions of basal temperatures and basal melt rates in slow-flowing regions. While it remains the largest source of uncertainty, simulations conducted with geothermal heat flow datasets leading to overall colder thermal conditions (e.g., Purucker, Reference Purucker2013; An and others, Reference An2015) compare better with englacial temperatures derived from borehole measurements and SMOS satellite observations but present larger misfits with respect to the subglacial lake distribution. The opposite is true for simulations using geothermal heat flow datasets leading to warmer thermal conditions (e.g., Fox-Maule and others, Reference Fox-Maule, Purucker, Olsen and Mosegaard2005; Martos and others, Reference Martos2017). The geothermal heat flow of Shen and others (Reference Shen, Wiens, Lloyd and Nyblade2020) and Haeger and others (Reference Haeger, Petrunin and Kaban2022) are associated with simulations presenting a reasonably good fit with respect to both borehole temperature measurements and the distribution of subglacial lakes, but other geothermal heat flow datasets could also perform well within their uncertainty ranges (e.g., Lösing and Ebbing, Reference Lösing and Ebbing2021; Stål and others, Reference Stål, Reading, Halpin and Whittaker2021). Assessing the geothermal heat flow models within their uncertainty bounds would require larger ensembles and emulators and could be addressed in future work.
Nevertheless, the performance and sensitivity of the ensemble of simulations also vary depending on the ice-sheet model, in which the thermodynamic approach and the representation of vertical velocities strongly affect the shape of the modeled temperature profiles. In particular, cold-ice models could overestimate the thickness of the temperate ice layers due to inadequate treatment of the transition between cold and temperate ice. Furthermore, the approximations considered in ice flow models, especially the representation of basal slip, considerably influence the estimates of basal melt by modulating the contribution of frictional heat in fast-flowing regions.
Finally, our results indicate that simulations estimating colder thermal conditions better fit with englacial temperature measurements. However, this might be biased by a potential overestimation of thermal conditions based on present-day boundary conditions and the lack of thermal memory in our model simulations.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/jog.2025.10087.
Acknowledgements
This publication was generated in the frame of Beyond EPICA. The project has received funding from the European Union’s Horizon 2020 research and innovation programme under grant agreement no. 815384 (Oldest Ice Core). It is supported by national partners and funding agencies in Belgium, Denmark, France, Germany, Italy, Norway, Sweden, Switzerland, The Netherlands and the United Kingdom. Logistic support is mainly provided by ENEA and IPEV through the Concordia Station system. The SMOS data used in this study (Macelloni and others, Reference Macelloni, Leduc-Leballeur, Montomoli, Brogioni, Ritz and Picard2019) have been developed in the ESA Project - “Cryosmos. STSE SMOS+Cryosphere.” ESA contract No.4000112262/14/I-NB. Computational resources have been provided by the Consortium des Équipements de Calcul Intensif (CÉCI), funded by the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS) under Grant No. 2.5020.11 and by the Walloon Region. Olivia Raspoet is a FRIA grantee of the Fonds de la Recherche Scientifique de Belgique (F.R.S.-FNRS).
Code and data availability
The codes used to run the Kori-ULB model and produce the figures, as well as data supporting the results presented in this study, are accessible on Zenodo: https://doi.org/10.5281/zenodo.15664172 (Raspoet and Pattyn, Reference Raspoet and Pattyn2025). The code and reference manual of the Kori-ULB ice-sheet model are publicly available on https://github.com/FrankPat/Kori-ULB. The borehole measurements of temperature profiles are available through their original references or the Antarctic ice sheet paleo-constraint database (Lecavalier and others, Reference Lecavalier2022; Reference Lecavalier2023). All datasets used in this study are freely accessible through their original references. Additional information can be provided upon request from the authors.
Competing interests
The authors declare that they have no conflicts of interest.
Open access
