Geothermal Heat Flow and Thermal Structure of the Antarctic Lithosphere

High‐quality maps of Geothermal heat flow (GHF) are crucial when modeling ice dynamics, shape, and mass loss of the Antarctic Ice Sheet, which is one of the largest potential contributors to sea level rise. The determination of GHF remains challenging, as in situ data are sparse and geophysical models exhibit large discrepancies in amplitude and resolution, especially on regional scales. Using a novel approach implementing a joint inversion of gravity and seismic tomography data with various geophysical and mineral physics information, we estimate the 3D thermal lithospheric structure and present a new GHF map. The resulting surface heat flow correlates with the location of subglacial volcanism and can represent a boundary condition for accurate ice dynamics models that can explain observed acceleration in the ongoing ice mass loss. Absolute values are within the range of other seismology‐based methods and are much lower than those obtained using for example, magnetic data. High uncertainties remain in the parametrization of the upper crustal structure and thermal parameters.

slope drive ice retreat in these glaciers . The results of our study suggest that the elevated GHF along with the corresponding bedrock topography is also a destabilizing factor in the ice sheet dynamics.
Of all the factors influencing the ice mass balance, GHF is one of the least constrained parameters. Improving our understanding of the conditions at the base of the ice sheet is crucial to explain the observed ice mass loss in Antarctica and the feedback mechanisms this induces. It might even provide a key for predicting future changes in the ice mass balance and consequently sea level. The determination of GHF remains challenging, as in situ data from boreholes are sparse and ill-distributed, leaving large parts of the continental interior unexplored (Burton-Johnson, . These borehole estimates can also coincide with local anomalies and would thus not be representative of regional to continental scales (Fisher et al., 2015). Additionally, most of the Antarctic boreholes do not reach the solid bedrock and instead represent heat flow estimates derived from a temperature gradient within the ice sheet or unconsolidated sediments, which can be influenced by many other factors such as climate forcing, hydrothermal circulation, or ice dynamics (Burton-Johnson, . Consequently, geophysical models are necessary to derive GHF on larger scales. Within the last few years, significant progress has been made in estimating GHF. Most studies focused on only one type of data, generally using results from either seismological (e.g., An et al., 2015;Shapiro & Ritzwoller, 2004;Shen et al., 2020) or magnetic (e.g., Dziadek et al., 2021;Fox Maule et al., 2005;Martos et al., 2017) studies. Recently, thermal isostasy has also been used to estimate GHF (Artemieva, 2022). Additionally, some novel, probabilistic methodologies were proposed, including a multivariate approach (Stål et al., 2021), Bayesian inversion (Lösing et al., 2020) of several data sets, and a machine learning approach (Lösing & Ebbing, 2021). Demonstrating general accordance in the main features related to EANT and WANT, these models still exhibit large discrepancies in both amplitude and lateral variation and show rather partial consistency with the observed volcanic activity. For a comprehensive review of the state of the art of Antarctic GHF see Burton-Johnson, . To overcome uncertainty and incompleteness of the available data, it is therefore necessary to use existing data sets complementarily.

Method
In this study, we derive GHF for Antarctica using several independent data sets in a three-step approach. These data sets contain information on seismic velocities, topography, mineral properties, and the gravity field. First, we calculate the initial temperature distribution in the upper mantle by iteratively combining seismic tomography and gravity data considering composition and density variations self-consistently (Haeger et al., 2019). In this inversion scheme, s-wave velocities from a global tomography model (Schaeffer & Lebedev, 2013) were converted to temperatures following Stixrude and Lithgow-Bertelloni (2005) taking into account anharmonicity as well as anelasticity (Cammarano et al., 2003) based on an initial composition. Thermal density variations are determined based on the temperature field and are subsequently used to correct the gravity field and the residual topography. These fields are jointly inverted to obtain compositionally induced densities and a new composition model that is in turn used to update the temperature field. This process is continued until convergence is reached. The absolute temperatures obtained in the final iteration can still be biased due to the input tomography reference model, requiring further calibration using independent data on geotherms from well-studied cratons that shared a tectonic history with EANT (Kaapvaal Craton and Western Australian Craton, (Artemieva, 2006)). The velocities for the coldest locations of each craton are picked, averaged for each depth layer and the difference to the coldest locations in Antarctica is taken. The entire Antarctic velocity model is then shifted according to this difference (compare Section S1.1 in Supporting Information S1).
Second, we define the lithosphere-asthenosphere boundary (LAB) in a thermal sense based on the resulting geotherm by assuming it corresponds to the 1300°C isotherm. We still observe an overestimation of the LAB depth compared to previous determinations (e.g., O'Donnell et al., 2019;Pappa et al., 2019;Steinberger & Becker, 2016), indicating that the previous calibration step insufficiently corrects the bias caused by the input model. This is especially noteworthy in WANT, where seismic LAB determinations are more accurate than in cratonic EANT. In order to correct for this, we compare the seismic LAB depth at distinct locations in WANT (O'Donnell et al., 2019) to our model, average the differences and shift the entire LAB model accordingly.
Third, we solve the one-dimensional steady-state heat equation independently for each grid point to obtain the temperature distribution in the lithosphere. T is the temperature, z is the depth, and A is the radiogenic heat production. Introduction of the temperature dependence of the thermal conductivity ( , ) leads to non-linearity in the inversion that is dealt with using an iterative approach. The model consists of a 3-layer lithosphere and includes the upper and lower crust and lithospheric mantle with the 1300°C isotherm representing the LAB as a lower boundary condition ( Figure S1 in Supporting Information S1).
One of the most uncertain parameters is the radiogenic heat production in the upper crust (Lösing et al., 2020), which is estimated to cause 26%-40% of the surface heat flow (e.g., Artemieva & Mooney, 2001;Hasterok & Chapman, 2011). Only limited in situ data are available on thermal properties from samples collected mostly in coastal areas or outcrops (Gard et al., 2019). On a regional scale, the most comprehensive study of radiogenic heat production has been conducted by Burton-Johnson et al. (2017) for the Antarctic Peninsula. Additionally, we compared the data to those obtained from direct measurements on continents that were adjacent to Antarctica in the supercontinent Gondwana (Pollett et al., 2019). As a result, we present several alternative models of GHF.
The main parameters used in the calculations are summarized in Table S1 in Supporting Information S1. While the different parameterizations do have a large impact on the resulting GHF, it should be noted that the main characteristics, that we will describe and analyze in detail, are persistent in all models, although their amplitude can vary. Therefore, we discuss only one preferred GHF model. For this preferred model, we defined a variable upper crustal thickness that combines the relative upper crustal thickness presented in Baranov et al. (2018) with the total crustal thickness obtained in Haeger et al. (2019). Radiogenic heat production is assumed to be vertically homogeneous within the upper crust. Laterally, the continent was divided into EANT and WANT and the available heat production determinations (Gard et al., 2019) were averaged independently for both regions, yielding values of 1.82 and 1.68 μW/m³ for WANT and EANT, respectively. The thermal conductivity within the mantle is kept constant at 4 W/m K (Petrunin et al., 2013) (compare Section S1.4 in Supporting Information S1). This was determined as a good compromise between including as much of the sparsely existing data as possible and limiting the introduction of additional model assumptions and uncertainties from these data sets. The model is presented on a 10 by 10 km grid. Actual resolution however is limited to the input tomography data. Quantification of this resolution is difficult to obtain, however it is unlikely to exceed 200 km.
Since we use a 1D approach for the thermal modeling, we estimated the possible errors of the method compared with the full 3D approach. We obtain variations between the 1D and 3D models for the surface heat flow of less than ±1 mW/m 2 over most of the continent. At the EANT WANT transition, higher differences of up to 2.8 mW/m 2 occur, which is less than the deviations expected due to uncertainties in weakly defined crustal parameters (See Section S3.2 in Supporting Information S1 for more details).
In contrast to previous studies, we also present a 3D temperature model of the Antarctic lithosphere. The latter is important for further studies, since the GHF at the ice-bedrock interface can be highly affected by ice dynamics down to depths of 3-5 km (e.g., Petrunin et al., 2013). By providing the temperature gradient through the entire lithosphere, our model allows for direct calculation of the GHF at any depth or horizon (see Figure S5 in Supporting Information S1). Validation of the final surface heat flow model is done through comparison with independent data such as bore hole estimates, the location of subglacial volcanoes, ice mass loss and ice flow velocities.

Results and Discussion
The obtained LAB map ( Figure 1b) exhibits a clear division of the continent in two parts, with depths of over 200 km in EANT and shallower than 90 km in most of WANT. Minimal thickness is found in Marie Byrd Land (⁓75 km). The resulting temperature profile from bedrock surface to the LAB is shown in Figure 1c. The general trend of our preferred surface heat flow model (Figure 2a) correlates closely with the thickness of the thermal LAB: high heat flow dominates in WANT while EANT is characterized by generally lower GHF. Although the trend corresponds to previous models, vast differences occur on smaller scales. In WANT, the maximal GHF exceeds 80 mW/m 2 with local maxima in Marie Byrd Land, along the shore of the Antarctic Peninsula, northern Victoria Land and in proximity to the South Pole. Most of EANT displays GHF between 50 and 60 mW/m 2 . Three regions diverge from this general trend and exhibit GHF above 60 mW/m 2 . Elevated GHF is found in Wilkes Land surrounding the Totten Glacier, in the Tonian Ocean Arc Super Terrane (TOAST) and surrounding the Lambert Glacier with an emphasis on the Gamburtsev Subglacial Mountains.
Both lithospheric thickness and GHF mirror the dichotomy of Antarctica that is caused by the fundamentally different tectonic origin and age of its eastern and western sections (e.g., Boger, 2011;Veevers, 2012). The deep lithospheric roots and the GHF of 50-60 mW/m 2 observed in EANT are within the normal range for cratonic regions of Precambrian origin (e.g., Artemieva & Mooney, 2001). In WANT, the elevated GHF and thin lithosphere agree well with the younger tectonic age of this region that only reached its current configuration in the Cenozoic (Dalziel & Elliot, 1982). A closer look on regional variations of the new GHF model shows an improved correspondence with recent tectonic and volcanic activity as well as with elevated ice discharge and sparsely existing bore hole estimates of high confidence (for more detail compare Section S3.1 and Figure S6 in Supporting Information S1).
Within WANT, elevated GHF is found along the West Antarctic Rift System (WARS), but the highest values for GHF are found in Marie Byrd Land and along the Antarctic Peninsula which is shown in Figure 3a. The WARS has exhibited periodic tectonic activity from the Early Cretaceous, though the end of this activity is still under debate (e.g., Behrendt, 1999;Müller et al., 2007;Siddoway et al., 2004;Winberry & Anandakrishnan, 2004). Volcanic activity also plays a major role for the thermal state of WANT. Subglacial volcanoes have been found to . This agrees very well with the aforementioned regions exhibiting the highest heat flow within Antarctica. Compared to GHF estimated for regions with volcanic activity on other, better studied continents, the determined GHF (up to 80 mW/m 2 ) is relatively low (e.g., Flóvenz & Saemundsson, 1993;Tanaka et al., 2004). It is noteworthy though, that we do not consider possible ongoing or recent volcanism within the crust in our model, so the displayed heat flow values are mainly caused by the deeper GHF component. Such elevated surface heat flow predominantly originating from mantle heat sources could support the existence of a mantle plume below Marie Byrd Land that has been hypothesized for the past 30 years (e.g., Bredow et al., 2021;LeMasurier & Rex, 1989;Seroussi et al., 2017).
It is common to attribute changes in the ice dynamics and subsequent ice loss to atmospheric and oceanic forcing. However, recent studies suggest a direct link between the location of the origin of ice streams and zones of increased heat flow (Petrunin et al., 2013;Rogozhina et al., 2016;Smith-Johnsen et al., 2020) that allows to consider GHF as an important factor in ice dynamics. Our results also show a good spatial correspondence with the map of ice mass loss between 2003 and 2019 ( Figure 2c, , which warrants a detailed analysis. Most ice mass loss is located along the coast with a particular focus on WANT and the Amundsen and Bellingshausen Sea Sectors. Here, streams originating from the Antarctic Circumpolar Current reach the shelf break bringing comparably warm and salty Circumpolar Deep Water to the shelf, leading to rapid basal melting (e.g., Nakayama et al., 2013;Silvano et al., 2016;Wåhlin et al., 2010). The elevated GHF prevalent in WANT strengthens the instability, with basal melt lubricating the ice-bedrock interface, which can result in increased ice-slip  Table S1 in Supporting Information S1) (c) Ice mass change per year from 2003 to 2019 . Light gray lines correspond to the 0 m/yr of ice mass change isolines. rates (Goeller et al., 2013). A similar situation can be found in Wilkes Land surrounding the Totten Glacier. This region has been shown to house not only the highest ice mass loss (Rignot et al., 2019; but also the fastest acceleration of mass loss within EANT (Velicogna et al., 2020). First explanations have focused on similar mechanisms as observed in the Amundsen and Bellingshausen Sea Sectors. Combined with a retrograde slope, this could explain the far-reaching ice mass loss surrounding the Glacier. This coincides with a region of elevated GHF, which can further destabilize the ice sheet and could drive unstable retreat, which underlines the importance of reliable GHF estimates in accurately modeling ice dynamics.
Similarly high GHF is predicted for Dronning Maud Land with an emphasis on the TOAST, without a reflection in the ice mass loss. This region is even experiencing mass gain . The high GHF here is caused by the shallower LAB without significant crustal thinning. The stable ice sheet in this region is likely connected to the underlying bedrock topography (Figure 1a). The lower surface temperatures related to the higher altitude compete with the elevated GHF in this area. The net effect of these factors leads to basal temperatures below the pressure melting point (Van Liefferinge & Pattyn, 2013) impeding fast ice flow, which agrees with the observed ice flow velocities (Rignot et al., 2011). Except for the shoreline, the topography here has a negative gradient in the direction normal to the coast line (Morlighem et al., 2020), which further impedes ice flowing toward the coast. This is supported by numerical simulations of the paleo ice sheet evolution in the Oligocene (Siegert, 2008) and Miocene (Gasson et al., 2016), which have shown that Dronning Maud Land is among the regions housing the most stable ice sheet of the entire AIS.
Elevated GHF also surrounds the Lambert Glacier, especially in the Gamburtsev Subglacial Mountains. As it was shown for similar glaciers in Greenland (Smith-Johnsen et al., 2020), high GHF in combination with steep topography in the mountains can be one of the main reasons for the initiation of ice streams. A similar streaming onset could be argued for the Lambert Glacier that originates in the Gamburtsev Subglacial Mountains. At the same time, we find reduced GHF directly below the Lambert Glacier. Elevated GHF has two major origins in our model: radiogenic heat production in the upper crust, which is enriched in radiogenic elements and high temperatures in the uppermost mantle leading to shallow LAB. Surrounding the Graben, elevated GHF is governed by the thick upper crust with elevated radiogenic heat production. Crustal thickness and hence radiogenic heat production is reduced in the Graben structure itself. Additionally, lithospheric thinning, which should cause high GHF in such rifts, is not apparent in the LAB map. If crustal thinning is caused by glacial erosion rather than crustal stretching in the rift zone, no compensatory asthenospheric uplift is expected due to the small width of the structure. However, relatively low resolution of the tomographic data, used for the thermal model (Schaeffer & Lebedev, 2013), does not allow us to unambiguously determine the cause of this anomaly.
When comparing different models of GHF, a general difference in amplitude is visible between those obtained based on seismic determinations (e.g., An et al., 2015;Shen et al., 2020) and those based on magnetics (e.g., Dziadek et al., 2021;Fox Maule et al., 2005;Martos et al., 2017) with the latter producing much larger amplitudes especially in WANT. With maximal amplitudes of ∼85 mW/m 2 , our model falls into the same range as previous models based on seismology. A possible cause of the difference in amplitudes is that processes such as lateral heat transfer from recent volcanism or hydrothermal circulation are not considered in models as ours. The lateral distribution of GHF varies strongly between all models, especially in the WANT. Again, our model is in better agreement with the seismological models attributing maximal GHF to Marie Byrd Land while continental models based on magnetics (Martos et al., 2017) and based on machine learning (Lösing & Ebbing, 2021) associate the highest GHF toward the WARS and central WANT and have a lower correlation to recent volcanic activity. Dziadek et al. (2021) present a regional model for the Amundsen Sea Sector with GHF exceeding 110 mW/m 2 below the Thwaites and Pope glaciers. We also find elevated GHF in this area with a continuation of the maximum toward Marie Byrd Land, though for the reasons mentioned above, we cannot reproduce their amplitude. While our model cannot reach the lateral resolution of Burton-Johnson et al. (2017) for the Antarctic peninsula, a comparison still yields good agreement between the studies. They determine RHP values between 1 and 2 μW/ m³ for large parts of the peninsula which is in agreement with our chosen value of 1.82 μW/m³ for this area. Their resulting GHF spans values between ∼70 and 90 mW/m 2 for the majority of the study area matching our determinations of 70-85 mW/m 2 . Values exceeding 120 mW/m 2 are only found in very confined, small-scale structures and cannot be resolved with our model. The thermal isostasy-based model by Artemieva (2022) finds homogeneously elevated GHF in most of WANT with a distinct minimum in the WARS. They argue, that steady-state heat conduction approaches such as the one employed in this study are not valid within WANT which is characterized by recent tectono-magmatic activity. We agree that the nonstationary thermal state can noticeably affect GHF. However, we also suggest that uncertainties in the input data needed to calculate the unsteady thermal conductivity problem will increase the uncertainty in the calculated heat flux. Although we cannot avoid this problem with the method used, we must point out the existence of such a problem.
Many input parameters, especially those describing the upper crust, remain highly uncertain. Since it is challenging to estimate their individual numerical uncertainty, we apply a statistical approach. We estimate the uncertainties of our model assumptions by calculating the standard deviation of models 2-31 ( Figure 2b, description of models in Table S1 in Supporting Information S1). The simplified reference model 1 is excluded from the calculation. In general, the standard deviation seems to be directly correlated to the crustal thickness, underlining the previous assumption that the highest uncertainties are induced by the largely unknown thermal parameters of the (upper) crust. The highest uncertainties of around +19 mW/m 2 can thus be found in the Gamburtsev Subglacial Mountains in EANT. In WANT, the highest uncertainties are found in Marie Byrd Land (∼+14.4 mW/m 2 ) and in the Antarctic Peninsula (∼+15.6 mW/m 2 ). The uncertainties induced by mantle parameters appear to be small in comparison. In order to estimate the influence of the input tomography model and the different calibration steps, we calculate the simplified test case with a shifted LAB depth of ±20 km. Largest uncertainties are associated with Marie Byrd Land and reach −6.8 mW/m 2 for a downward shifted and +11.7 mW/m 2 for an upward shifted LAB. Deviations over EANT remain very small and do not exceed +2 mW/m 2 for the majority of the craton.

Conclusions
The model presented in this study provides not only a novel surface heat flow map, but also a three-dimensional thermal model of the entire Antarctic lithosphere. These data can help to further refine surface GHF values as new data on the crustal structure and its thermal properties are obtained, and can also be used as a reference model providing boundary/initial conditions for further geophysical modeling. We estimate uncertainties of our model through statistical analysis. High values are associated with regions of thick crust, especially in the Gamburtsev Subglacial Mountains where they reach +19 mW/m 2 . Uncertainties related to the seismological tomography input model and the calibration thereof are lower and reach their maximum in Marie Byrd Land, where the lithosphere is thinnest. Absolute values are within the range of other seismology-based methods and are much lower than those obtained using for example, magnetic data. We assume high reliability of our model for the deeper lithosphere, which is confirmed by independent data, such as the correlation of predicted thermal anomalies with the main tectonic components of the Antarctic continent and the identified centers of volcanic activity. This as well as the spatial correspondence of positive GHF anomalies with zones of maximum change in ice sheet dynamics indicates a direct connection of deep lithospheric and surface processes. This is in agreement with elevated GHF found below Thwaites and Pope glaciers in a regional study of the Amundsen Sea Sector (Dziadek et al., 2021), though we cannot reproduce their amplitudes. Therefore, it is crucial to include a consistent GHF distribution to accurately model ice dynamics, ice mass balance, and sea level change to better predict possible scenarios for the Earth system dynamics.