Crustal Heterogeneity of Antarctica Signals Spatially Variable Radiogenic Heat Production

Geothermal heat flow (GHF) is a key basal boundary condition for Antarctic ice‐sheet flow. Large‐scale variations are resolved by several recent models but knowledge of the smaller‐scale variations, crucial for ice sheet dynamics, is limited by unresolved variations in crustal radiogenic heat production. To define this at continent‐scale we use 3D gravity inversion constrained by seismic Moho estimates to identify variations in crustal composition and geometry beneath thick ice. Geochemically‐defined empirical relationships between density and heat production capture the global average trend and its variability, and allow to estimate from upper‐crust density spatial variations in radiogenic heat production. Significant variations are observed typically 1.2–1.6 μW/m3, and as high as 2 μW/m3 in West Antarctica. The contribution to GHF from these heat‐production variations is similarly variable, typically 16–24 mW/m2 and up to 60 mW/m2. The mapped variations are significant for correctly representing GHF in Antarctica.

, and may vary significantly on small length-scales.Interpolations derived from large-scale models from univariate data will fail to represent local variations, and yet these are critical for controlling ice-sheet dynamics (McCormack et al., 2022;Reading et al., 2022).Therefore, a way forward is to generate high-resolution crustal models that can capture the variations in crustal properties and structure.Such models could then be used to understand GHF variations using forward or statistical approaches (Lösing et al., 2020).
In Antarctica, both seismic and gravity applications resolve substantial regional variations in crust and mantle structure (Abrehdary & Sjöberg, 2021;Haeger & Kaban, 2019;Haeger et al., 2019;Lloyd et al., 2020;Pappa et al., 2019aPappa et al., , 2019b;;Stål et al., 2019).Studies of Antarctic crustal structure have mainly focused on the interface between crust and mantle (marked by Moho discontinuity).This interface is primarily constrained by receiver function analysis (Baranov & Morelli, 2013;Chaput et al., 2014;Dunham et al., 2020;Hansen et al., 2009).Recent applications have also begun to resolve structure within the crust (Zhou et al., 2022).However, due to the challenging fieldwork environment, seismic observations are sparse and have uneven spatial coverage, which limits their ability to resolve crustal structures in detail.
Benefiting from the more consistent spatial coverage of airborne and satellite gravity measurement, gravity inversion is an alternative approach to estimate the Moho undulations, variably constrained by seismic information (Chisenga et al., 2019;Llubes et al., 2018;Pappa et al., 2019a;Zhang et al., 2020).These studies assume that long-wavelength Bouguer gravity reflects Moho undulations and use geometry inversion techniques with a single density contrast to recover the Moho.However, this assumption is problematic in continental-scale applications due to lateral variations in crust and mantle density structure that can generate gravity anomalies of comparable magnitude to the Moho as well as affecting density contrast across the interface (Aitken, 2010;Aitken et al., 2013).An additional challenge is raised due to the large uncertainties in the topography and ice thickness data in Antarctica (Morlighem et al., 2020), which would directly impact the resolved Bouguer gravity.Free-air gravity is unaffected by this uncertainty, although there is still a requirement for topographic gravity effects to be modeled.
An alternative approach involves incorporating all critical surfaces and density structures as updatable features within a single model domain.A 3D gravity inversion is then utilized to optimize the density and surface geometries to fit the gravity data (Aitken, 2010;Aitken et al., 2013;Alghamdi et al., 2018).This approach accounts for a laterally heterogeneous crust and mantle with the capability to generate seismic-constrained Moho and density structure in a continental-scale application (Aitken, 2010;Aitken et al., 2013).
Here, we build a 3D model domain including sedimentary basin, upper and lower crust layers and mantle.During the inversion, the density and thickness of each layer are modified within seismic constraints and density bounds, and with differing constraints applied to the solution.An ensemble of results is assembled that fits both the gravity observations and the prior seismic constraints.The ensemble result presents the most representative model of Antarctic crustal structure possible within the capacity of the method and current data constraints.From the ensemble result and empirical relationships between density and radiogenic heat production from a global geochemical database, we explore the implications of heterogeneity in crustal composition and structure for radiogenic heat production and GHF in Antarctica.

Free-Air Gravity Data and Initial Model
In our inversion, we invert free-air gravity data.For the data set, airborne gravity data sets (the AntGG compilation (Scheinert et al., 2016) and regional-scale free-air gravity data, Text S1 in Supporting Information S1) are combined with the GOCO06s satellite gravity model (Kvas et al., 2019) calculated at 5 km ellipsoid height (Figure S1 in Supporting Information S1).After merging, the long-wavelength gravity information predominantly caused by the mantle convection or sub-lithospheric mantle density variations was removed using the spherical harmonics 0-10° (>2,000 km half-wavelength) from GOCO06s model.
The initial model incorporates the current understanding of Antarctic crustal structure, with a horizontal grid resolution of 20 km and extending to a depth of 75 km.It is divided into nine layers, including ice, lake water, seawater, sedimentary basins (with three layers accounting for sediment compaction with depth), crystalline crust (upper and lower crust), and mantle (Text S2 in Supporting Information S1).The topographic components of the 10.1029/2023GL106201 3 of 12 model are represented by BedMachine Antarctica (Morlighem et al., 2020).Initial sedimentary basin thickness is scaled from the relationship between sedimentary basin distribution likelihood model (Li et al., 2022) with thickness from individual estimations (Text S3 in Supporting Information S1).Initial Moho surface and its uncertainty are based on the kriging prediction and error using individual Moho estimates (see An et al. (2015a) and Table S2 in Supporting Information S1).The prediction error from kriging is used to represent the bounds of permissible Moho surface (±2 standard deviation).Within the crust, a two-layer equal thickness upper and lower crust is used in the continental crust to provide essential density contrasts at surface and Moho and a single-layer crust is used in the oceanic crust.The initial and permitted density ranges for each layer are shown in Table S1 in Supporting Information S1.

Gravity Inversion
3D gravity inversions were undertaken using VPMG software (Fullagar et al., 2008), which divides the model domain with prism-elements of pre-defined area but with arbitrary thickness.VPMG has two main gravity inversion styles: density-style, where the densities of prism-elements are varied, and geometry-style, where the elevations of the prism-element boundaries are modified.Each style works independently, and the other property cannot change.
This study followed the alternating inversion method used by Aitken et al. (2013), Alghamdi et al. (2018), andMoro et al. (2023).Each of our model runs ran three cycles of four single-iteration inversions, respectively one iteration for mantle density, then one iteration each for the density within the crust layers, the geometry of the Moho, and the sedimentary basin thickness.The overall inversion process is regulated by the maximum permitted change in each style: for a density inversion, the permitted change is the maximum density change in the current iteration; for a geometry style inversion, the permitted change is set as the maximum percentage of the interface elevation change.Increasing permitted changes allows more change to be accommodated in that phase, which increases the potential to resolve misfit through that style of inversion.
In total 35 inversions were run to generate a model-ensemble each containing results constrained by different permitted density or geometry changes.The permitted density change ranged from 5 kg/m 3 to 40 kg/m 3 by 5 kg/m 3 intervals while geometry change ranged from 5% to 13% by 2% intervals (Figure S5 in Supporting Information S1).Six models that do not sufficiently converge with gravity observation (misfit >10 mGal) are rejected and excluded from the model-ensemble.We use the mean of the model-ensemble and its standard deviation to show a representative crustal structure (Text S4 in Supporting Information S1).

Estimating Radiogenic Heat Production
For mapping of density to heat production, we estimate the statistical relationship between rock density and heat production computed from the global whole-rock geochemical database (Gard et al., 2019a) (Figures 1a and 1b, Text S6 in Supporting Information S1).For the global database, there is a large-scale trend that heat production decreases with increasing density value, albeit with non-linear behavior at the tails.It is clear from the scatter plot (Figure 1a) that the correlation between heat production and density is generally weak due to large natural variability between samples, which has been well documented in previous studies (e.g., heat production vs. silica content (Gard et al., 2019b), heat production versus Vp or density (Hasterok & Webb, 2017)).Despite the absence of a physically-based model to link heat production with density, the consistency between average density and heat production values across continents suggests a robust predictive power in using this empirical relationship (Hasterok & Webb, 2017).
It is of note that our geophysical model elements comprise prisms with areas of 400 km 2 and several kilometers to tens of kilometers thick.This volume comprises a mixture of lithologies and compositions, including variable proportions of mafic, felsic and intermediate composition rocks, as well as trace amounts of dense minerals such as garnet, olivine and iron oxides.The natural variations in heat production for individual rocks and lithologies will be mitigated by the very large volume sampled; we will detect in our model only systematic trends in the composition of rocks and factors such as metamorphic grade, that are also well represented in the geochemical database.
Although we cannot precisely estimate a heat production value from density for individual samples, the global median heat production and density are strongly correlated in density range between 2,650 and 2,950 kg/m −3 (R 2 = 0.993), which captures most igneous and metamorphic rocks.A similar relationship is evident for the 25th and 75th percentiles of the data set.Therefore, for a given density, we can estimate the most central heat production value in the global database and the associated interquartile range.Although the potential exists for systemic biases in the global geochemical database or anomalous crust in parts of Antarctica, we suggest that the density-heat production trends for large rock volumes should sit within the boundaries of the globally defined relationship.

Antarctic Crustal Structure
Our model-ensemble results indicate spatial variations in sedimentary layer thickness beneath the ice sheet and ice shelf (Figure 2a and Figure S6a in Supporting Information S1).We have resolved a sedimentary layer with average 1 km thickness in onshore West Antarctica where sedimentary basin is present.These results are consistent with passive seismic estimations (Chaput et al., 2014;Dunham et al., 2020;Zhou et al., 2022).In contrast, the East Antarctica section preserves thicker sedimentary packages, with an average thickness of 2 km.Major and thick sedimentary basins (3-4 km) are preserved in the Aurora Subglacial Basin and Pensacola Pole Basin.The Wilkes Subglacial Basin is thinner (1-2 km) in the north, with a thicker basin in the south (3-4 km).Thick sedimentary basins are identified beneath major ice shelves, with a thickness of 8-12 km resolved in the Filchner-Ronne Ice Shelf (FRIS).The basin is thickest closest to the edge of the ice shelf and becomes thinner toward the western side of the inner ice shelf.For the Ross Ice Shelf (RIS), we have resolved near 1-3 km of sediment that is thicker from the grounding line to the edge of the ice shelf into the Ross Sea.The western side of the RIS is generally thicker than the eastern side, separated by a mid-shelf high (Tankersley et al., 2022).The sedimentary basin beneath the Amery Ice Shelf is narrower, with a thickness of 3 km.
The thinner crust in West Antarctica (25 km) and thicker crust in East Antarctica (35 km) are captured by the mean Moho as shown in Figure 2b and Figure S6b in Supporting Information S1.The thinnest crust occurs in the Ross Sea, Byrd Subglacial Basin, and Byrd Subglacial Trench as part of the West Antarctic Rift System (WARS), while the thickest crust occurs beneath the Gamburtsev Subglacial Mountains (GSM).The standard deviation of the Moho surface shows variability within the model-ensemble, as seen in Figure S6 in Supporting Information S1.This variability reflects both the consistency of the initial model with the gravity field and the influence of the applied constraints.High standard deviation occurs where constraints are weak, and the initial model does not explain gravity observations.Conversely, low standard deviation occurs where either constraint is strong or the initial model is already explaining gravity observations.For instance, in the GSM region, we see that the deep Moho, constrained by seismic data (Hansen et al., 2010) shows a large variance after inversion.This large Moho 10.1029/2023GL106201 6 of 12 variation indicates that the initial Moho does not match gravity observation, which requires either reducing the crust thickness (Pappa et al., 2019a) or including a dense mafic underplate structure (Ferraccioli et al., 2011).In the Wilkes Land region, the Moho shows low variance despite low seismic data constraints.This indicates that the initial isostatic assumption is close to matching the gravity observation.

Antarctic Crustal Heat Production Estimate
Applying the empirical density to heat production relationship to our modeled upper-crust density, we generated an upper-crust heat production map (Figure 1c).The map shows that West Antarctica generally has expected heat production higher than in East Antarctica, with broad areas exhibiting expected heat production higher than 1.8 μW/m 3 , except for the ridge of the WARS and the Haag block with expected heat production less than 1.4 μW/m 3 .East Antarctica shows lower expected heat production, with an average of 1.2 μW/m 3 , while local areas such as Law Dome, Enderby Land, Vostok Highland, and the Dronning Maud Land (DML) region have expected heat production higher than 2 μW/m 3 .On the other hand, GSM and the southern Wilkes Subglacial Basin exhibit low expected heat production, less than 1 μW/m 3 .

Model Limitations and Value
At continental-scale, gravity inversion has low sensitivity to the depth of anomaly sources.Residual gravity anomalies longer than 400 km in wavelength (five times model vertical extension) are effectively indistinguishable from a 1D response, meaning that the depth of the sources of these anomalies is undefined.The low vertical sensitivity might account for the similar structure observed for both the mantle and crustal density.Our model sensitivity test indicates that, given the current constraints, changing the mantle density alone is insufficient to account for the residual gravity anomaly (Figure S4a in Supporting Information S1).On the other hand, crustal density has the power necessary to explain all onshore residual gravity anomalies (Figure S4c in Supporting Information S1).In this situation, the crustal anomaly could be up to 50% greater than the model-ensemble mean.However, a homogenous mantle is unlikely where seismic tomography models reveal large lateral variations in the upper mantle (An et al., 2015a;Lloyd et al., 2020;Shen et al., 2018).
The most robust constraint in our model is the Moho measurement.The model-ensemble mean Moho depth is close to the seismic Moho estimations with an RMS misfit of 3.8 km.This misfit falls within the uncertainty range of 3-5 km from receiver function studies (Figure S7 in Supporting Information S1), which is lower than the RMS misfit of other continental-scale models (e.g., 6.9 km in Pappa et al. (2019b), 6.2 km in Zhang et al. (2020)).Although constraints do not exist for the density distribution in the crust and mantle, the mantle density derived from model-ensemble mean exhibits a similar structure to the average mantle Versus (50 km beneath Moho) obtained from Bayesian inversion of Rayleigh wave and receiver functions (Shen et al., 2018).
Several models exist that resolve the large-scale variations of crustal structure (Afonso et al., 2019;Haeger & Kaban, 2019;Laske et al., 2013), and also the situation for focused regions (Capponi et al., 2022;Lösing et al., 2022), but these models cannot connect the finer-scale variations with the larger scale variations across the continent as showed in our model.

Heat Production Model Validation
A continental-scale heat production map had previously been estimated using inferred geological age and sparse outcrop (Stål et al., 2020).Due to the challenge of direct sampling from the subglacial environment, the geological age is largely unknown beneath the ice-sheet.Thus, in absence of ground-truth information in continental-scale, we compared our heat production map with the PetroChron Antarctica data set (Sanchez et al. (2021), Figure 1d), which is a geochemical estimation using outcrops.Our model shows a good match with the spatial distribution of heat production from the database, with elevated heat production in areas such as Antarctica Peninsula, TAM, Enderby Land, and Lambert Rift captured by our model.Our model also indicates low heat production values in the coastal region of East Antarctica.For a quantitative analysis, the misfit between the modeled heat production 10.1029/2023GL106201 7 of 12 is presented in the Figure S11 in Supporting Information S1 (mean = −0.14μW/m 3 , median = 0.26 μW/m 3 , standard deviation = 1.67 μW/m 3 ).Despite a general concordance, significant misfits are seen between our expected heat production values and those calculated from PetroChron Antarctica, due to the difference in both sampled volumes and sampling frequency.
The opportunity remains for parts of Antarctica to be anomalous with respect to their relationship between crustal density and heat production.East Antarctica includes several regions known to have extensive Archean crust including in DML and Enderby Land (Cox et al., 2023), as well as regions with high-grade metamorphic rocks, including eclogites (Brown et al., 2020) and high-heat-producing granite suites (Carson et al., 2014).West Antarctica meanwhile has regions with active subduction zones as well as volcanism and underlying plutonism (Jordan et al., 2020).These regions and their subglacial counterparts may exhibit atypical density-heat production relationships that may exacerbate or moderate the variability defined by our model, however, it is not possible reliably to express these relationships from currently available data.

Implication for Radiogenic Crust Contribution for GHF
We calculated the contribution of radiogenic heat production from the crust to the total GHF by using the estimated heat production value (Figure 3c) and the upper-crust thickness (Figure 3b).Our findings suggest that radiogenic heat production could contribute up to 50 mW/m 2 to GHF.Areas with thick crust and high heat production show a greater crustal contribution.In specific regions such as Lambert, DML, and Vostok Highlands, the crustal heat production could contribute to GHF exceeding 40 mW/m 2 (Figure 3c).Our map also highlights higher heat production values in West Antarctica.However, due to the thin crust in West Antarctica, the crustal contribution from this region averages 20 mW/m 2 .In some localized areas like the Antarctic Peninsula, Marie Byrd Land, and Ellsworth-Whitworth Mountains, the contribution from crustal heat production can be as high as 50 mW/m 2 .
In Antarctica, continental-scale GHF estimation usually assumes a laterally uniform crustal heat production value (e.g., 1 μW/m 3 at the upper-crust (An et al., 2015b), or 2.5 μW/m 3 at the surface and exponentially decreasing to 8 km (Martos et al., 2017)).To illustrate the contribution of the spatially variable heat production to the total GHF, we use a uniform 1 μW/m 3 as a reference heat flow model for comparison.We calculate the radiogenic crust heat flow contribution with a minimum (25th percentile), median (Figure 3c), and maximum heat production case (75th percentile).We then compare our model results with the reference model (1 μW/m 3 ).All models are calculated with the same upper-crust thickness (Figure 3b).
The heat flow variation compared to the reference model is within 10 mW/m 2 for the most conservative heat production estimation (25th percentile fitting), with East Antarctica showing negative variation and West Antarctica a positive heat flow variation (Figure 3f).For the median case, the heat flow change compared with the reference model is mostly positive, with an increase of 0.6 mW/m 2 at the 25th percentile and 6.8 mW/m 2 at the 75th percentile.For the maximum heat production case (75th percentile fitting, Figure 3e), areas including Antarctica Peninsula, Marie Byrd Land, TA, and the DML area show a positive heat flow variation up to 50 mW/m 2 above the reference.
Although our model exhibits significant variability depending on the chosen scenario, it is evident that, in all cases, highly variable heat flux is expected to be derived from radiogenic crust heat production variation.Further testing, which consider the uncertainties of resolved density structure, reveals a similar pattern for elevated heat flow variation (Figure S12 in Supporting Information S1).

Implication for Ice Sheet System
GHF is a key boundary condition of the ice sheet flow as elevated heat flow warms the ice sheet bed, affecting ice rheology and potentially increasing basal melt rate (McCormack et al., 2022).Regions with higher GHF also lead to the formation of temperate ice, which in turn results in faster surface flow (Pittard et al., 2016).GHF also impacts ice sheet dynamic related subglacial hydrology.In areas where the basal temperature exceeds the pressure melting point, subglacial meltwater is generated, forming an input to the hydrological system.Based on numerical ice sheet modeling, the basal melting rate increases linearly with GHF at a rate of approximately 0.1-0.4mm/yr/(mW/m 2 ) (Joughin et al., 2009;Seroussi et al., 2017).
When considering the contribution of heat production to GHF, the model suggests a strong spatial variation in heat flow (ranging from 6 to 60 mW/m 2 ) due to the incorporation of variations in crustal heat production.Considering that mantle heat flow contributes to a long-wavelength heat source beneath the ice sheet, the short-wavelength heat source resolved by our model (40 km) resulting from crustal heat production could result in a spatially heterogeneous basal melt (0.6-24 mm/yr).This spatial variability in basal conditions could be a crucial factor for the subglacial hydrology system and impact ice sheet dynamics.

Conclusions
We developed a model-ensemble of the crustal structure of Antarctica using gravity inversion, constrained by seismic information.Analyzing the resulting upper-crustal density we estimated expected volumetric heat production using an empirical relationship defined from a global geochemical database on the basis that the rock-volumes represented in our model elements are overall consistent with the global sample set.Incorporating also upper-crust thickness, we generated a continental-scale crustal radiogenic heat production estimation and its contribution to spatially variable GHF (ranging from 6 to 60 mW/m 2 ).
Our work offers crucial insights into the likely pattern and magnitude of heterogeneity in a key heat source beneath the ice sheet, and associated limits.The identified crustal heterogeneity and its uncertainties are key to

Figure 1 .
Figure 1.Estimation of radiogenic heat production in Antarctica.(a) 2D histogram shows empirical density and heat production relationship in log 10 space using a global data set (Gard et al., 2019a).The upper and lower triangles, and the cross, represent the 75th, 25th, and median heat production values in each density bin, respectively.The dash-dotted line indicates the linear regression; (b) linear regression of density and heat production, with the density extension is set with the dashed-line range in panel (a); (c) modeled heat production map based on the empirical relationship using median densityheat production relationship; (d) median heat production value from PetroChron Antarctica data set in 20 × 20 km cell bin.It shows spatially limited heat production estimates based on geological samples across Antarctica (Sanchez et al., 2021).EL, Enderby Land; LD, Law Dome; VH, Vostok Highland.Colormap is truncated at 3 μW/m 3 .

Figure 2 .
Figure 2. The estimated structures from model-ensemble mean includes: (a) sedimentary basin thickness; (b) Moho depth; (c) crystalline crust density; and (d) mantle density.The model suggests distinct structures when comparing East Antarctica (thick basin, thick and dense crust) with West Antarctica (thin basin, thin and less dense crust), and highlights small-scale variations within the crust.