High geothermal heat flow beneath Thwaites Glacier in West Antarctica inferred from aeromagnetic data

Geothermal heat flow in the polar regions plays a crucial role in understanding ice-sheet dynamics and predictions of sea level rise. Continental-scale indirect estimates often have a low spatial resolution and yield largest discrepancies in West Antarctica. Here we analyse geophysical data to estimate geothermal heat flow in the Amundsen Sea Sector of West Antarctica. With Curie depth analysis based on a new magnetic anomaly grid compilation, we reveal variations in lithospheric thermal gradients. We show that the rapidly retreating Thwaites and Pope glaciers in particular are underlain by areas of largely elevated geothermal heat flow, which relates to the tectonic and magmatic history of the West Antarctic Rift System in this region. Our results imply that the behavior of this vulnerable sector of the West Antarctic Ice Sheet is strongly coupled to the dynamics of the underlying lithosphere. The Amundsen Sea sector, Antarctica, is underlain by shallow Curie depths – where the magnetic properties of rocks change – according to airborne magnetic data. This suggests high geothermal heat flow in this region of the West Antarctic Rift System.

T he Amundsen Sea sector of West Antarctica (Fig. 1a, b) is the Earth's prime region to test hypotheses on the interacting relationship and coupling between the solid earth and the cryosphere. In this study, we examine the distribution of geothermal heat flow (GHF) as a key factor to reveal past and present geodynamic, tectonic, and lithospheric processes and their correlation to the presently observed enormous ice mass loss in this region.
Geothermal heat flow of the Antarctic continent is one of the essential geophysical parameters for both delineating and identifying tectonic and geodynamic features and processes, and improving the parameterization of basal conditions in ice-flow dynamic models. Apart from a few locations 1 , thick ice cover and logistical challenges have largely prevented the installation of a dense network of direct temperature gradient measurement sites in drill holes, and will likely do so in the future. Indirect determination of heat flow from other geophysical parameters such as magnetic anomalies 2,3 , seismological models [4][5][6] , and interpretation of ice-penetrating radar images 7 as well as the estimation of crustal heat production from exposed geological provinces 8 have improved with greater coverage of data and accuracy of the models, and increasingly help relate tectonic and magmatic features to ice sheet dynamics.
The Amundsen Sea Embayment is underlain by the eastern branches of the West Antarctic Rift System, which extends from the Ross Sea to the Amundsen Sea and Bellingshausen Sea sectors of West Antarctica 9-11 . Lithospheric and crustal properties, such as thin crust, thin elastic lithospheric thickness, and high uplift rates, as well as interpreted tectonic and magmatic features, such as individual rift basins, fault systems, and young volcanic activity, coincide with the presently observed largest ice mass loss in Antarctica (outside the northern Antarctic Peninsula) by rapid thinning and retreat of the Pine Island, Thwaites, Pope, and other glaciers in the embayment [12][13][14] . This tectonic setting of the ASE is likely to favor relative high geothermal heat flow. Continent-wide studies of geothermal heat flow derived from magnetic and seismological data consistently indicate a substantially higher geothermal heat flow of West Antarctica compared to that of East Antarctica 2-5 . However, these continent-scales models still suffer from large uncertainties and lack the regional detail for the Amundsen Sea and other sectors.
Results and discussion Curie depth distribution. The Curie depth is considered as a proxy for geothermal heat flow, because it represents the depthto-bottom of a magnetic source in a first approximation of the temperature dependence of magnetization in crustal rocks (details are in Supplementary Method 1.1 to 1.3). We find that Curie depths (Fig. 2b) range between 10 and 18 km with the 580°C isotherm dipping down towards Thurston Island. In the central part of the Thwaites and Pope glacier catchment systems, the Curie isotherm exhibits shallow depths of 12-16 km, propagating onto the Amundsen Sea Embayment shelf. Shallower 580°C isotherms are located in the Byrd Subglacial Basin region and the southwestern part of the Bentley Subglacial Trench. The Pine Island Rift structure does not show prominent Curie depth anomalies, but indicates shallow thermal anomalies in the eastern part. Our estimates for the standard deviation of the Curie depth ( Supplementary Fig. 3) is ±1 km for most of the study area and peaks with ±2.7 km in the magnetic basin between Thurston Island and Ellsworth-Whitmore Mountains, and towards Marie Byrd Land, where magnetic data coverage is poor. Further considerations are reported in Supplementary Method 1.3 and Supplementary Discussion 2.1. We note, that the shallow Curie depth distribution correlates well with volcanic centers (e.g., Mount Takahe), which are further indicative for presumed elevated crustal heat flow. Geothermal heat flow data from continental rifts, such as the Baikal, Basin and Range, East African, Rhine and Rio Grande rifts, indicate that young rifts are characterized by high heat flow 18 . Mean values from the rift graben systems within rifts are typically in the range of 70-125 mW/m 2 . In non-volcanic sections of the East African Rift System, heat flow is normal to low, but this system shows strong volcanic heat transport in its young rift grabens 19 .
The Amundsen Sea sector of the West Antarctic Rift System is characterized by recent or ongoing magmatic activity, which has been hypothesized for the Pine Island Glacier region 20 and the Executive Committee mountain range in Marie Byrd Land 21 . Further volcanic edifices are found to cluster in the southwestern flank of the BST, and beneath Thwaites glacier 22 . The absence of relevant seismicity 23 indicates that there is no active tectonism of the rift zone. Narrow regions of low seismic velocities in the upper mantle, for instance beneath the BST, may however infer a recent localized phase of extension 24 .
High geothermal heat flow beneath Thwaites and Pope glaciers. Cumulative evidence, particularly the results from our study, points towards strongly elevated heat flow in the Amundsen Sea sector of the West Antarctic Rift System (Fig. 3a). The Moho is shallow with depths mostly between 17 and 32 km 25,26 , and its distribution correlates well with our Curie depth estimates. The thinned crust was attributed to Late Cretaceous/Early Cenozoic extensional phases 25 , and is underlain by an upper mantle with presumably high temperatures in the depth range from 75 to 250 km 27 . The lithosphere is relatively thin (<100 km), and density variations here predominantly stem from high upper mantle temperatures 28 . Our shallow Curie depth estimates and inferred high heat flow correlate with areas of low lithospheric effective elastic thickness 29 T e (Fig. 3b), indicative of a weak rifted lithosphere. The relatively high T e estimates across the continental shelf are likely an effect of the up to 7 km thick sedimentary sequence of inferred early Cretaceous to Cenozoic age 30 . The region of high T e propagates east towards Thurston Island Block, which was part of a vast region along the formerly convergent paleo-Pacific margin stretching as a magmatic arc from the eastern Amundsen Sea sector to and along the western Antarctic Peninsula 31 . Regions of low T e are prone to react more sensitively to changes in crustal loads. The observed glacial isostatic adjustment from ice mass loss in the Amundsen Sea sector exceed the uplift rates predicted from modeled glacial isostatic adjustment by more than 20 mm/year 32 , suggesting a rapid upper mantle flow response. The highest measured uplift rates were found just west of Thwaites Glacier and near Pope Glacier in regions of low upper mantle viscosity 32 , which notably coincide with the region of our estimated high geothermal heat flow and low T e (Fig. 3a, b).
Results of previous studies of Antarctic geothermal heat flow (Fig. 3c) vary strongly in their heat flow distribution and resolution for the Amundsen Sea sector. Heat flow based on continent-wide Curie depths 3 was derived from the first edition of the Antarctic Digital Magnetic Anomaly Project (ADMAP) 33 . The regional magnetic anomaly dataset presented in this study was not available for that compilation. Hence, their Curie depths for our study area mostly relied on satellite magnetic data that have lower spatial resolution, and therefore cannot resolve more local geothermal heat flow highs that we infer from aeromagnetic data in the Amundsen Sea sector. Other studies used radar echograms to interpret areas of basal melting for the estimation of geothermal heat flow in the Thwaites Glacier catchment 7 . Their GHF estimates are comparable to those of the continental scale models 3,4 . GHF estimates from insitu temperature gradient measurements lie within the lower end of that range 15,16,34,35 . Our GHF results presented in this study yield a higher resolution estimate for this highly dynamic part of the West Antarctic Ice Sheet.
The elevated geothermal heat flow band (Fig. 3a) is interpreted as caused by an anomalously thin crust 25,36 , underlain by a hot mantle 26,27 and possible tectonic reactivation along older fault systems (e.g., Mount Murphy Rift). The magnetic anomaly grid reveals NNE-SSW trending lineaments located between the Thurston Island and Marie Byrd Land crustal blocks (Fig. 3d). These lineaments are regionally confined and attributed to the magnetic signature of Cenozoic rifting in the eastern branch of the West Antarctic Rift System. Crustal heterogeneity and magmatic intrusions are other likely contributors to the observed heat flow variability, for instance in the vicinity of Byrd Subglacial Basin. High seismic shear-wave velocities (>5 km/s) at 20 km depth may relate to lower crustal mafic intrusions 37 that likely relate to subglacial volcanism 22 . Previous studies suggested that large, deep crustal plutons seated in extended continental crust may have formed rapidly but induced a prolonged perturbation of the deep-crustal geotherm 38 .
The thermal anomalies, attributed to a thin and laterally heterogeneous rifted crust, magmatism and inferred fault reactivation, are likely to cause a heat-advective effect on the deep hydrological system and, therefore, exert a profound influence on the flow dynamics of the West Antarctic Ice Sheet in the Amundsen Sea sector. The direct transfer of heat can facilitate basal melting and control the ice rheology and basal sliding, and thus erosion 39 . High geothermal heat flow beneath Thwaites and Pope glaciers could further contribute to rapid past and future changes in the glacier system. Our results in the Amundsen Sea region provide a new base for discussing the location and extent of crustal-scale thermal anomalies. This is a key finding to better characterize basal sliding properties and subglacial hydrology, as well as refine thermal boundary conditions for studies of ice sheet dynamics in the most rapidly changing sector of the West Antarctic Ice Sheet.

Methods
Magnetic anomaly dataset. The new magnetic anomaly grid for the Amundsen Sea sector of West Antarctica we compiled represents a further improvement with respect to the latest generation of the publicly available ADMAP 2.0 Antarctic-wide magnetic anomaly grid 40 . We integrated and processed existing offshore and onshore airborne magnetic anomaly data, unpublished onshore airborne data and data from new survey flights with detailed levelling to better constrain the lithospheric and crustal properties in the region. We expanded the AWI helicoptermagnetic anomaly grid 41 with 14 flights (~2880 km total survey line length) on the continental shelf of the eastern Amundsen Sea Embayment during RV Polarstern expedition PS104 42 (Supplementary Fig. 1). Furthermore, we added unpublished line data from the Pine Island Rift region and focused on cross-point calculations which exhibited a small offset (>50 nT). This step improved the data across the grounding lines of the Pine Island and Thwaites glaciers. Weather conditions did not allow for additional flights to link the onshore and offshore regions better. Additionally, a base station correction could not be applied to the helicopter-borne data offshore but was applied to all onshore data. Therefore, the data may still be affected by diurnal variations in the geomagnetic field. All magnetic data were DC shifted to the regional long-wavelength domain of the Magnetic Field Model MF7 (https://www.geomag.us/models/MF7.html) to help homogenize the longwavelength magnetic field across different surveys 43 . In areas, where no flight lines are available we also filled data gaps with the MF7 grid (Fig. 2a). A detailed description of the data and processing can be found in Supplementary Method 1.2.
Curie depth estimates. The depth-to-bottom of the magnetic source (Z b ) is in first approximation the Curie depth, which is considered as a proxy for geothermal heat flow, because of the temperature dependence of magnetization in crustal rocks. Regions found to have shallow Curie point depths are expected to have higher heat flow, and, therefore, higher average temperature gradients. Curie depth estimates mark a transition zone, rather than an exact boundary, and assume a main magnetic source from magnetite with a Curie temperature of 550-580°C 44,45 . This assumption neglects the compositional variability in plutonic rocks that lead to Curie temperature ranges between 300 and 680°C, and in cases of magnetic assemblages of Fe-Ni-Co-Cu metal alloys up to 620-1084°C 46 . Hence, we assume a temperature of 580°C at the lowest magnetic boundary. We applied the centroid method after Tanaka et al. 47 to window sizes of 200 × 200 km and 300 × 300 km ( Supplementary Figs. 2 and 3). The windows were extracted equidistantly from the newly compiled magnetic anomaly grid with a spacing of 50 km. We considered the fractal approach (power-law scaling) to calculate the Z b as not applicable. High fractal scaling could underestimate Curie depths due to overcorrection 48 and has been discussed previously 49,50 . From the slope of the power spectrum, the top boundary (Z t ) and the centroid of a magnetic layer (Z 0 ) composed of a horizontal (equivalent) layer are estimated ( Supplementary Fig. 2). The basal depth of the magnetic source is Z b = 2Z 0 − Z t . The obtained bottom depth Z b of the magnetic source is assumed to be the Curie point depth. The Curie depth error, which we calculated using the standard deviation from the spectral analysis, is less than ±2 km for most of the area (Supplementary Fig. 3).
Thermal modeling. A degree of uncertainty surrounds the mineralogical properties in the region. We therefore approximate the steady-state geothermal heat flow in a model of a homogenous material via where Q is the heat flow at the bed surface (in mW/m 2 ), T is the Curie temperature (580°C), z the Curie depth (in km), and k the thermal conductivity of the lithosphere (in W/mK). The resulting geothermal heat flow thus represents a linear transformation of the Curie depth map.

Data availability
The datasets collected for and generated in this study are available from PANGAEA Data Publisher (www.pangaea.de