Topographic Correction of Geothermal Heat Flux in Greenland and Antarctica

We present a new approach to account for the influence of subglacial topography on geothermal heat flux beneath the Greenland and Antarctic ice sheets. We first establish a simple empirical proportionality between local geothermal flux and topographic relief within a given radius, based on a synthesis of existing observations of these properties elsewhere on Earth. This analysis essentially yields a high‐pass filter that can be readily applied to existing large‐scale geothermal heat flux fields to render them consistent with known subglacial topography. This empirical approach avoids both the geometric limitations of existing analytic models and the complex boundary conditions required by numerical heat flow models, yet it also produces results that are consistent with both of those methods, for example, increased heat flux within valleys and decreased heat flux along ridges. Comparison with borehole‐derived geothermal heat flux suggests that our topographic correction is also valid for non‐ice‐covered areas of Earth and that a borehole location uncertainty of >100 m can limit the value of its inferred heat flux. Ice‐sheet‐wide application of this approach indicates that the effect of local topography upon geothermal heat flux can be as important as choice of regional geothermal heat flux field across a small portion of Antarctica (2%) and a larger portion of Greenland (13%), where subglacial topography is best resolved. We suggest that spatial variability in geothermal heat flux due to topography is most consequential in slower‐flowing portions of the ice sheets, where there is no frictional heating due to basal sliding. We conclude that studies of interactions between ice sheets and geothermal heat flux must consider the effect of subglacial topography at sub‐kilometer horizontal scales.

While several studies have estimated geothermal heat flux across Greenland and Antarctica as functions of spatially varying crust and mantle properties Fox Maule et al., 2005;Greve, 2005;Greve & Herzfeld, 2013;Martos et al., 2017Martos et al., , 2018Petrunin et al., 2013;Rezvanbehbahani et al., 2017;Schroeder et al., 2014), no study has yet directly considered a topographic correction for geothermal heat flux across the entirety of either ice sheet. While any topographic dependence of geothermal heat flux is not expected to influence its regional value at larger spatial scales (≳ 10 km), because of the constrained self-affinity of subglacial topography , it may have significant implications for subglacial hydrology and basal conditions at smaller scales (≲ 10 km). Local reorganizations of water or ice flow can instigate larger-scale downstream changes (Fahnestock et al., 2001;Mantelli et al., 2019;Pittard et al., 2016). Spatial variability in geothermal heat flux can also be a consideration when reconstructing the subglacial hydrology of paleo-ice sheets (Näslund et al., 2005).
In this study, we develop and evaluate a new topographic correction for geothermal heat flux across Greenland and Antarctica. The topographic correction is a spatially variable but dimensionless value that accounts for the effect of local topographic relief on local geothermal heat flux. We do not present a new geothermal heat flux model per se, but rather we consider the implications of this topographic correction on existing large-scale geothermal heat flux fields. We describe our data product as a topographic "correction," as it modifies a large-scale geothermal heat flux field to account for small-scale topographic features, in a manner similar to downscaling of large-scale climate models via a topographic correction (Noël et al., 2017).

Data
For both Greenland and Antarctica, we employ the BedMachine data sets (v3 for Greenland, v1 for Antarctica), for which bed topography is derived from multiple interpolation methods of existing ice-thickness measurements, mostly from radar sounding. These methods include ordinary kriging in the ice-sheet interior and mass conservation in fast-flowing regions (Morlighem et al., 2017(Morlighem et al., , 2019. BedMachine has spatial resolutions of 150 m for Greenland and 500 m for Antarctica. The spatial resolutions of these input bed-topography data sets determine the 150-and 500-m grid resolutions of the topographic correction fields that we generate for each ice-sheet domain. These output fields represent first-order estimates of the topographic correction for geothermal heat flux, and they are self-consistent with their respective BedMachine topographies.
For both Greenland and Antarctica, we apply these topographic corrections to existing background geothermal heat flux data sets that are consistent with available observations of geology, aeromagnetic and gravimetric anomalies and associated Curie depth estimates (Martos et al., 2017. Both of these geothermal heat flux data sets have been evaluated against sparse available measurements. We bilinearly interpolate these 15-km resolution geothermal heat flux data sets to the 150-m resolution of our geothermal heat flux anomaly field for Greenland and 500 m for Antarctica. This interpolation is necessary to reconcile the geothermal heat flux and BedMachine spatial grids, and it results in nearly identical background geothermal heat flux values at adjacent grid nodes.

Topographic Correction for Geothermal Heat Flux
Several approaches have been developed to account for the influence of topographic relief on local geothermal heat flux in ice-free terrain. The commonly applied Jeffreys-Bullard method is a two-dimensional (2-D) cross-sectional solution of the heat equation (Bullard, 1938;Jeffreys, 1938). In this approach, the far-field lower boundary is the geothermal heat flux, while the upper boundary is the mean annual surface temperature. In the absence of an ice sheet, this upper boundary temperature can be estimated locally using an atmospheric lapse rate (Gruber et al., 2004). However, applying the Jeffreys-Bullard approach beneath an ice sheet requires well-constrained temperatures at the ice-bed interface to serve as the upper boundary condition. In practice, it is exceptionally challenging to reliably constrain ice-bed interface temperatures across large scales (Liefferinge & Pattyn, 2013;MacGregor et al., 2016;Pattyn, 2010).
Presently, the most common approach for estimating topographic correction of geothermal heat flux is the one-dimensional (1-D) horizontal-profile analytical method developed by Lachenbruch (1968) (e.g., Ayunov & Duchkov, 2008;Ganguly et al., 2000;Shi et al., 1988;van der Veen et al., 2007). This method decomposes a given topographic profile into a set of idealized plane slopes that bracket the observed subaerial or submarine topography. The geothermal heat flux correction is then estimated from solution tables as a function of toe and brink angles and slope height. Because this method solves the topographic correction of slope segments independently, it cannot be easily applied to model the interaction between two or more topographic features, such as opposing valley walls. Thus, the Lachenbruch (1968) method is not well-suited for implementation to the complex 2-D topography beneath ice sheets (Lee, 1991).
Here, we develop a simple geostatistical approach for estimating a 2-D topographic correction for geothermal heat flux. Our objective in developing this method is to produce the simplest possible ice-sheet-wide estimate of geothermal heat flux that is consistent with the known effects of subglacial topography. This approach consists of three steps. First, we empirically estimate the linear proportionality between local geothermal heat flux and topographic relief within a given radius by compiling a new synthesis of terrestrial observations of these properties. Second, we evaluate the topographic relief in Greenland and Antarctica using a 2-D filter within an empirically determined radius. Finally, we combine the empirically derived geothermal heat flux dependence on topographic relief with the calculated topographic relief. This approach allows us to express the topographic influence on geothermal heat flux as relative anomalies-or a dimensionless correction-rather than absolute values. This approach effectively yields a high-pass filter to apply to an existing geothermal heat flux field, to render it more physically consistent with local 2-D topography.
We conceptualize the influence of topography upon geothermal heat flux as the perturbation of a given geothermal flux (G) by an anomaly (ΔG G / ): where G denotes the local geothermal flux after topographic correction. This approach is broadly analogous to the formulation of topographic correction in local gravimetry studies (Hammer, 1939). We estimate the local topographic correction for geothermal heat flux (ΔG/G)-a non-dimensional quantity-as a function of local topographic relief as: characteristic height  can be interpreted as the local topographic relief (within r) necessary to induce a 100% change in local geothermal heat flux. While this linear model provides a first-order assessment of topographic correction, we acknowledge that there may be a physical basis for employing higher-order models to describe the relation between topographic relief and geothermal heat flux.
Following Equation 1, we calculate a first-order estimate of topographic correction for geothermal heat flux across Greenland and Antarctica by perturbing existing geothermal heat flux estimates for both regions with our self-consistent topographic correction for geothermal heat flux to yield topography-corrected geothermal flux. We emphasize that this topographic correction is independent from the geothermal heat flux estimated in previous studies. For simplicity, in this study, we only apply our topographic correction to one geothermal heat flux data set per region, but we emphasize that it can be readily applied to other existing geothermal heat flux data sets.
We estimate the local uncertainty associated with topographically corrected geothermal heat flux ( G  ) as the quadratic sum of the independent uncertainties associated with both the large-scale regional geothermal heat flux estimate ( G  ) and the topographic correction for geothermal heat flux (  

Empirical Model Parameters
We empirically determine optimal values of the characteristic height () and averaging radius (r) based on application of Equation 2 to previous studies that explicitly attributed measured or modeled local geothermal heat flux anomalies to topographic relief. These studies were selected for their deliberate investigations of the influence of topography on geothermal heat flux at the kilometer scale, which we consider most suitable for consideration of the effect of subglacial topography, in particular for deeply incised troughs where we expect the topographic influence to be greatest. Previous studies variously present observed topographically corrected geothermal heat flux (Westaway & Younger, 2013), modeled topographically corrected geothermal heat flux (Ayunov & Duchkov, 2008;Blackwell et al., 1980;Safanda, 1987;van der Veen et al., 2007) or both (Ganguly et al., 2000;Shi et al., 1988) in a total of 12 continental, three lacustrine and three oceanic settings (Table 1).
At each of these 18 sites, we applied averaging radii between 1,000 and 10,000 m and characteristic heights between 400 and 1,400 m to identify the combination of these two parameters that minimizes the rootmean-squared error (RMSE) between our model and previous studies (Figures S1.1-S1.18). Based on this synthesis, we adopt an averaging radius (r) of 5,000 m and a characteristic height () of 950 m. This parameter combination yields an average RMSE of 9.0% between our model and previously observed or modeled topographically corrected geothermal heat fluxes (Table 1). Given multiple additional factors that are poorly constrained, but which likely further influence this uncertainty, such as subglacial geology or groundwater flow, we conservatively assign a relative uncertainty of twice this RMSE (18%) to the geothermal heat flux anomaly inferred by our approach. This uncertainty estimate is independent of BedMachine uncertainty, as the spatially variable bed elevation uncertainty in BedMachine is independent between adjacent cells (Morlighem et al., 2017(Morlighem et al., , 2019. While BedMachine uncertainty provides a good estimate of the absolute elevation uncertainty of a given cell, it does not provide a good estimate of the relative elevation uncertainty between adjacent cells. Ideally, this uncertainty could be characterized by calculating subglacial topographic relief across an ensemble of simulations, but that is, beyond the scope of our present study. With the empirical parameters we prescribe, our geostatistical model reasonably reproduces previously studied 1-D horizontal profiles that highlight topographic influence on geothermal flux. The topographic correction for geothermal heat flux associated with an idealized 1-D cross-sectional topography of Jakobshavn Isbrae estimated by our statistical approach yields an RMSE of 22.3% compared to that modeled independently by van der Veen et al. (2007) using the Lachenbruch (1968) analytical method ( Figure 1b). Notably, our method generally infers geothermal heat flux anomalies of smaller amplitude. A sub-kilometer horizontal offset in minimum values likely reflects the inability of the analytical method to fully capture the interacting effects of opposing valley walls. At Cascadia Basin (L21), our model reproduces the observations of Ganguly et al. (2000) with an RMSE of 4.5% (Figure 1d), which is substantially lower than using the Lachenbruch (1968) analytical approach (8.1%). The RMSE between the two model methods- Lachenbruch (1968) and our model-is 7.0%. Similar to Jakobshavn Isbrae, the topographic correction for geothermal heat flux estimated by our model at Cascadia Basin is lower in amplitude than that inferred by the Lachenbruch (1968) analytical model. However, for this field site, the observed values suggest that the topographic effect is indeed less extreme than that predicted by the Lachenbruch (1968) analytical model.

Comparison with 3-D Numerical Solution
To provide an evaluation target for our geostatistical approach, we also estimate geothermal heat flux at the ground surface with a three-dimensional (3-D) numerical model of heat conduction. This first-principles approach is effectively a numerical solution of a 3-D variant of the Jeffreys-Bullard solution for estimating the topographic influence on geothermal heat flux (Bullard, 1938;Jaeger, 1965;Jeffreys, 1938), and it follows similar recent studies (Petrunin et al., 2013;Rogozhina et al., 2016). The governing equation for steadystate heat conduction in isotropic solids in the absence of internal heat sources is given by: Note. The RMSE associated with the best-fit mean characteristic height (950 m) and averaging radius (5 km) adopted for this study is also shown. "Data type" indicates the target data used in this study.
The bold values denote the mean (or average) of each column. where T is the temperature (Carslaw & Jeager, 1959). Heat flux through the solid is given by Fourier's Law of thermal conduction: where q is the heat flux and k is the thermal conductivity of the solid (Carslaw & Jeager, 1959). For homogenous materials, Equation 4 reduces to Laplace's Equation.
We solve these equations numerically over a portion of Disko Island, West Greenland, using FEniCS, an open-source and automated computing platform based on the finite element method (Alneas et al., 2015;Logg et al., 2012). This Disko Island model domain encompasses a 900 km 2 area in the northwest of the island at a spatial resolution of 150 m.
The surface elevation of the Disko Island model domain is derived from the same BedMachine v3 data set that we use for the geostatistical approach (Morlighem et al., 2017, Figure 2a). The mean annual surface temperature boundary condition for the Disko Island model domain is derived from topographic-downscaling of simulations from the MAR v3.5.2 regional climate model forced by ERA-20C climate forcing over the period of 1961-1990(Fettweis et al., 2017. This Disko Island model domain was deliberately selected as predominantly ice-free terrain, which allows the surface boundary condition for geothermal heat diffusion to be constrained as the climatological mean annual air temperature, and thus avoids assuming ice-bed temperatures (Bullard, 1938;Jeffreys, 1938). The bottom boundary condition is a regional geothermal flux of 60 mW m −2 applied uniformly at a depth of 10 km below sea level ). The lateral model edges are specified as insulated boundaries, which is a reasonable approximation for large model areas. This approach of Type 1 (Dirichlet) surface boundary condition and Type 2 (Neumann) bottom and lateral boundary conditions is common for modeling near-surface geotherms (Jaeger, 1965;Noetzli et al., 2007;Petrunin et al., 2013). We assume the domain has a homogenous thermal conductivity of 2.25 W m −1 K −1 , which is characteristic of the basalts that comprise c. 90% of Disko Island's surficial geology (Fuchs et al., 2013;Jennings et al., 2019;Larsen & Pedersen, 2009 Lachenbruch (1968) approach and this study (with uncertainty envelope). (c) Idealized submarine topography at Cascadia Basin (L21) considered by Ganguly et al. (2000). (d) Corresponding topographically corrected relative geothermal heat flux estimated by Ganguly et al. (2000) using the Lachenbruch (1968) approach and this study, with uncertainty envelope as light red fill. The normalized geothermal heat fluxes observed by Ganguly et al. (2000) are shown for reference.
The model domain is discretized into tetrahedral elements, with horizontal coordinates of the ground surface vertices aligning with the locations of geothermal heat flux estimates from our geostatistical approach. The temperature field is represented in a continuous, piecewise-linear function space. Estimates of geothermal heat flux are found in a post-processing step using Equation 5. Since the geothermal gradient derived from the discrete temperature solution is discontinuous and piecewise-constant, geothermal heat flux is calculated at the centers of the mesh cells. This cell-centered flux field is projected linearly onto the mesh vertices to facilitate a direct comparison between our first-principles and geostatistical approaches. The 30 × 30 × 10 km Disko Island model domain contains ∼1.4 million elements. The computational time required to estimate topographic correction for geothermal heat flux within this 900 km 2 domain from numerical solution of the 3D steady-state heat equation is comparable to the computational time required to apply the geostatistical approach over the entire Greenland BedMachine domain (4 × 10 6 km 2 ). This highlights that the geostatistical approach is ∼10 4 times faster than the numerical approach per unit area.
Topographic correction for geothermal heat flux is estimated from numerical results for the Disko Island model domain following Equation 1. These 3-D numerical results indicate that the geostatistical approach generally characterizes the magnitude and spatial distribution of the topographic correction for geothermal heat flux ( Figure 3). While no systematic biases are apparent between the geostatistical and numerical approaches, the numerical approach clearly captures substantially more kilometer-scale spatial variability than the geostatistical approach. Agreement between both methods appears to be best along valley bottoms.
The geostatistical approach appears to have the most positive disagreement along steep valley walls, and the most negative disagreement across elevated highlands. The RMSE between the numerical and geostatistical approaches across the entire Disko Island model domain is 9.4%. This suggests that the ±18% relative uncertainty that we assign to geothermal anomalies inferred by our approach is reasonable.
Comparison against a numerical model is not ideal because the numerical solution does not fully capture reality. It does not capture known spatial differences in thermal conductivity associated with rock type, systematic biases in mean annual temperature associated with solar heating or persistent microclimates, and unknown spatial variability in subsurface processes that influence geothermal flux like groundwater flow or heat production. These known differences occur both horizontally and vertically. A 3-D cross-sectional view of the numerical model suggests that geothermal heat flux anomalies >10% can propagate to several kilometers depth (Figure 4). We also acknowledge that the Disko Island numerical domain features geothermal heat flux anomalies across what is mostly an air-rock boundary, while we apply the geostatistical model to COLGAN ET AL.
10.1029/2020JF005598 7 of 26 the subglacial ice-rock boundary where spatial variability in the thermal conductivities of both rock and ice likely play a significant role.
Following conservation of energy, the total heat flux through the topographic surface must be equal to the total heat flux prescribed along the bottom model boundary. Therefore, while the prescribed geothermal heat flux at 10 km depth will influence the absolute magnitude of geothermal heat flux at the surface, it does not influence the normalized geothermal heat flux that is our main focus. Likewise, a uniform change in mean annual surface temperature field will not alter the model's estimate of either the total geothermal heat flux or the normalized geothermal heat flux variation across the topographic surface. However, spatial variability in the mean annual surface temperatures (i.e., the surface boundary condition) will impact the model's estimate of the normalized geothermal heat flux variations. The results presented in Figure 3 thus represent a general topography-dependent normalized heat flux variation, independent of the heat flux stipulated at 10 km depth, the spatially uniform thermal conductivity value, and uniform changes to the mean annual surface temperature.
COLGAN ET AL.
10.1029/2020JF005598 8 of 26  In reality, of course, we must acknowledge that thermal conductivity is spatially variable. Globally, the thermal conductivity of common rock types is generally recognized to range between 0.5 and 7.0 W m −1 K −1 ), although the variation among individual rock types is much smaller (Čermák & Rybach, 1982). At Disko Island, however, the vast majority (c. 90%) of surficial geology is flood basalts (Maligât formation) and picrite basalts (Vaigat formation), with only a fringe of Cretaceous-Tertiary sediments covering the eastern coast (Larsen & Pedersen, 2009). Spatial variations in thermal conductivity due to spatial variations in geology are therefore very limited around the value of 2.25 W m −1 K −1 that we employ for Disko Island. At other locations, however, it is conceivable to find geologic settings in which neighboring rock types have significantly different thermal conductivities. In these settings, geotherms close to the geological boundary will be influenced by the contrast between relatively high and low conductivity rock types. Layered geology with thermal conductivity contrasts can lead to refraction of heat flux at interfaces. It is theoretically conceivable that sharp spatial contrasts in thermal conductivity associated with rock type can result in geothermal corrections of a similar magnitude to those associated with topography. Jaeger (1965) presents a discussion of the influence of thermal conductivity variations on geothermal heat flow. In principle, our numerical model can incorporate the influence of thermal conductivity variations if the geological formations underlying a site are well constrained. Improving knowledge of subglacial ice-sheet geology presents the tantalizing possibility to constrain geothermal corrections associated with local subglacial geology (Dawes, 2009). This type of local geology correction would represent an additional layer of refinement for local geothermal heat fluxes-over and above the local topographic correction presented here-and is an interesting topic for future research.

Topographically Corrected Geothermal Heat Flux
Across Greenland and Antarctica, our estimates of the topographic correction for geothermal heat flux highlight spatial patterns that are consistent with increased geothermal heat flux within valleys and decreased geothermal heat flux along ridges or mountains ( Figures 5 and 6). We predict that geothermal heat flux in several regions-most notably central East Greenland and the Antarctic Peninsula-is routinely enhanced by ∼50% within deeply incised glacier valleys and correspondingly reduced by a similar magnitude along adjacent ridges. In Greenland, BedMachine resolves ∼100 outlet glaciers that are both sufficiently narrow and deeply incised to result in a topographic correction of >150%. This pattern is consistent with the hypothesis first presented by van der Veen et al. (2007)-that subglacial topography likely has a significant influence on geothermal heat flux beneath ice sheets. While we only show areas on land and beneath ice sheets here, this pattern also holds between fjords and their adjacent walls. Ocean areas are included in the data product that accompanies this study.
In regions with less topographic relief, geothermal heat flux has correspondingly less topographic dependence. While geothermal heat flux within the interiors of both ice sheets appears to have limited topographic dependence in our model, we note that this pattern is primarily a function of poorly resolved subglacial topography there. Where interior ice-sheet bed topography is relatively well-resolved, for example, inland of Jakobshavn Isbrae in Greenland (Figure 5g) or the subglacial Gamburtsev Mountains in East Antarctica (Figure 6g), we infer appreciable spatial variability in geothermal heat flux (±50%). Within ice-sheet interiors, even a single individual airborne radar transect can sufficiently resolve bed topography to infer non-trivial topographic corrections for geothermal heat flux. This is especially apparent along the coast of central East Greenland, where the topographic dependence of geothermal heat flux is clearly evident along widely spaced airborne transects (Figure 7a).

Spatial Variability in Geothermal Heat Flux
We next place our modeled topographic correction for geothermal heat flux field in the context of the range of existing subglacial geothermal heat flux distributions. We characterize uncertainty in geothermal heat flux as the fractional ensemble spread across previously published and publicly available geothermal heat flux models for both Greenland and Antarctica. We calculate fractional ensemble spread as local ensemble spread divided by local ensemble mean across both domains. This field highlights regions of maximum and minimum agreement between existing subglacial geothermal heat flux data sets. In Greenland, we consider an ensemble of five geothermal heat flux models (Greve & Herzfeld, 2013;Martos et al., 2018;Rezvanbehbahani et al., 2017)-treating the three realizations of Rezvanbehbahani et al. (2017) as independent samples ( Figure S2). In Antarctica, we use an ensemble of four geothermal heat flux models ( Figure S3; An et al., 2015;Fox Maule et al., 2005;Martos et al., 2017;Shapiro & Ritzwoller, 2004). Most of the diversity within both ensembles reflects differences in modeling approaches between studies. This analysis highlights that the relative magnitude of the topographic correction that we infer exceeds that of inter-model range in geothermal heat flux around most of the periphery of Greenland, where pronounced bed topography is well resolved (Figure 8). In Antarctica-despite the coarser spatial resolution COLGAN ET AL.  of BedMachine there-there are also extensive subglacial regions where the topographic corrections for geothermal heat flux is more significant than the inter-model range in geothermal heat flux (Figure 9). Simply put, acknowledging the effect of topography upon geothermal heat flux appears to be just as important as choice of geothermal heat flux model across substantial sectors of both Greenland and Antarctica. We emphasize that the effect of this topographic correction is localized, and the ultimate underlying regional geothermal heat flux value is prescribed from an independent geothermal heat flux model.
The peripheries of the Greenland and Antarctic ice sheets not only inherently sample a wide range of regional geothermal heat fluxes, but also have relatively well-resolved bed topography in comparison to the ice-sheet interiors. We therefore further highlight the significance of the spatial variability in geothermal heat flux due to topography around both ice-sheet margins. Accounting for topography effectively increases the spatial variability of geothermal heat flux characteristic of both ice-sheet margins. In Greenland, where there is less apparent spatial variability in regional geothermal heat flux, this effect is especially apparent. There, including our topographic correction increases the 99th percentile-or unusually high-geothermal flux value from 75 to 100 mW m −2 , and decreases the 1st percentile value from 51 to 26 mW m −2 (Figure 10). Simply put, in regions of extreme regional geothermal heat fluxes, acknowledging topographic correction can result in even more extreme local values.

Comparison with Borehole-Derived Measurements
There are insufficient borehole-derived geothermal heat flux measurements in Greenland or Antarctica to robustly demonstrate that our topographic correction for geothermal heat flux significantly decreases discrepancy between modeled and measured geothermal heat fluxes. We therefore turn to the continental  (Rysgaard et al., 2018) highlights how deep incision of a fjord can contribute to anomalously high local geothermal heat fluxes. Here, the topographic correction is shown, rather than masked, in ocean areas. (c) The South Greenland Ilimaussaq sites (Sass et al., 1972) highlight, with black boxes, the ± 1.8 km of positional uncertainty associated with borehole latitudes and longitudes each reported to the nearest arcminute. United States, where thousands of borehole-derived measurements of geothermal heat flux are publicly available, to explore the utility of the topographic correction we present. We calculate our topographic correction for geothermal heat flux at 150-m spatial resolution over a domain extending between 30°-45 °N and 90°-120 °W (Figure 11). We derive this 150 m spatial resolution topography from bilinearly downscaling a 100-m spatial resolution digital elevation model (USGS, 2013). Aside from the Yellowstone Hotspot, there are no extreme geothermal heat flux anomalies associated with mantle or tectonic processes within this domain.
COLGAN ET AL.

10.1029/2020JF005598
14 of 26 Within this domain, there are 5,160 borehole-derived geothermal heat flux measurements available with a stated positional accuracy of ±0.0001 decimal degrees, or approximately ±10 m (Hasterok, 2010). However, among borehole positions reported to four decimal places, there likely remains a non-trivial quantity of decimal sequences that suggest some borehole positions were converted from the nearest arcminute to degree decimal, leading to false precision. For example, positions ending in 0.3333 likely reflect the conversion of positions initially reported as 20 arcminutes into decimal degrees. Because we cannot presently be certain where arcminutes have been converted into degree decimals of false precision, we cannot systematically filter out boreholes where the true positional uncertainty is one arcminute, or ∼±1.8 km, from this comparison.
To demonstrate the value of our topographic correction, we calculate the RMSE between available borehole geothermal heat flux measurements and the geothermal heat flux field modeled by Davis (2013)-both with and without the topographic correction applied. We bi-linearly interpolate the ∼1° (or ∼110 km) Davis (2013) global model to the 150-m resolution of the continental United States domain. We refer to comparing borehole-derived fluxes to the Davis (2013) model as the "standard" comparison, while comparing borehole-derived fluxes to the same model "including topographic correction" as the improvement yielded from a systematic topographic correction. We make these comparisons several times, each time for smaller subsets of the borehole measurements restricted to higher minimum topographic corrections associated with greater topographic relief ( Figure 12). The RMSE between the borehole and modeled geothermal heat fluxes "with topographic correction" clearly decreases relative to the "standard" comparison with increasing topographic relief. In settings of greater relief, applying topographic corrections of >25% decreases RMSE to only ∼18% of the RMSE when topographic correction is omitted (9 vs. 49 mW m −2 ). Because the borehole sample size decreases as topographic correction increases, smaller borehole subsets correspond to more limited regions over which the Davis (2013) model and Hasterok (2010) database are being compared.
If we consider only the 419 observations with topographic corrections >0.05-the largest subset in which topography may be expected to influence geothermal heat flux-we find that the Davis (2013)  distributions span several orders of magnitude, we calculated these linear regressions using only values within one standard deviation of the mean (89 ± 57 mW m −2 ). This case also reveals a slight negative bias in the topographic correction (−4 mW m −2 ), which highlights that accounting for topography more frequently down-corrects geothermal fluxes measured in valleys, than up-corrects geothermal fluxes measured on ridges. This pattern is consistent with the notion that geothermal flux measurements have historically been collected disproportionately in valley bottoms, due to their relative ease of access in comparison to ridge tops (Westaway & Younger, 2013   While this improvement in correlation (r) and RMSE are insignificant, we interpret this continental United States case study as indicating that topographic correction does statistically improve comparisons between globally modeled and locally measured geothermal heat fluxes. While this comparison also highlights the previously recognized challenges of comparing global geothermal heat flux models with local boreholes measurements, it also highlights the previously unrecognized challenge presented by the non-trivial positional uncertainties presently embedded in borehole databases. Additionally, as the continental United States is nearly entirely ice-free, this case study inherently excludes the complex interactions between an ice sheet and bedrock that also influence geothermal heat flux (e.g., Alley et al., 2019;Stevens et al., 2016).

Contextualizing Existing Ice-Sheet Measurements
As sparsely distributed geothermal heat flux observations are critical for constraining regional geothermal heat flux models in both Greenland and Antarctica, it is desirable to apply a systematic topographic correction to these observations. Simply put, a borehole topographic correction of even a few percent can be important to the regional heat budget when it is the only borehole measurement constraining geothermal heat flux within hundreds of kilometers (Greve, 2019). Geothermal heat flux measurements also tend to be preferentially located within topographic lows, such as accessible valley bottoms, which can result in a systematic warm bias within observational data sets if a topographic correction has not been applied (Westaway & Younger, 2013).
For example, at Young Fjord's Dybet site (74.46307°N, 21.19075°W), Rysgaard et al. (2018) measured a geothermal heat flux of 93 ± 21 mW m −2 . At this site, we calculate a topographic correction of 0.54 ± 0.10, meaning that the fjord's relatively narrow and steep geometry enhances geothermal heat flux by 54% ± 10% relative to the regional mean geothermal flux (Figure 7). Applying our inferred topographic-geothermal anomaly to the Dybet observation would yield a topographically corrected regional geothermal heat flux of 60 ± 23 mW m −2 . Our topographic correction therefore suggests that the relatively high local geothermal flux measured at the Dybet site is primarily due to local topography, rather than site proximity to the present-day onset location of the Northeast Greenland Ice Stream or the Iceland hotspot track (Rysgaard et al., 2018).
An additional practical challenge in applying topographic corrections to historical geothermal heat flux measurements in Greenland is knowing the precise location of those measurements. For example, Sass et al. (1972) report geothermal borehole locations at Ilimaussaq in South Greenland to the nearest arcminute, equivalent to a positional uncertainty of ∼1.8 km. Given the complex terrain in the vicinity of Ilimaussaq, the topographic corrections we calculate within ± 1.8 km of the reported Ilimaussaq boreholes range between −50% and +30% (Figure 7). In other cases, such as evaluating even the most recent geothermal heat flux models (Martos et al., 2017, the reported positional uncertainty of available geothermal heat flux measurements may be the nearest tenth of a decimal degree, equivalent to ∼11 km of positional uncertainty. We suggest that a positional uncertainty of <100 m is required for confident calculation of a topographic correction for geothermal heat flux measurements. At present, borehole-temperature profiles from deep ice-core sites provide our only in situ information about subglacial geothermal heat flux. We therefore explore the topographic correction we calculate in the vicinity of six sites in Greenland and 12 sites in Antarctica where subglacial geothermal heat flux has been determined in situ (Table 2). While we provide the coordinates of these sites in polar stereographic projection, we note that the previously published positional accuracies we employ have only been reported to a tenth of a decimal degree for many sites (Martos et al., 2017. In Greenland, topographic corrections are <3% at all deep-core sites ( Figure S4). These minor topographic corrections may partially reflect the deliberate placement of deep-cores at flow divide sites with minimal topographic relief. Spatial variability in topographic correction is greatest at Dye-3, where it reaches ±6% within 3 km of the core site, which is the only site not deliberately selected as a deep-core location. The subglacial topography beneath Greenland's three southernmost deep-core sites (Dye-3, Greenland Ice Core Project [GRIP] and Greenland Ice Sheet Project 2 [GISP2]) may act to very slightly reduce geothermal flux by 2%-3%.
In Antarctica, where the spatial resolution of subglacial topography is coarser, topographic corrections are <4% at all borehole sites ( Figure S5). Our estimated topographic correction, however, does reach +46 ± 27% at the Dyer Plateau site on the Antarctic Peninsula, where Nicholls and Paren (1993) report a geothermal heat flux of 100 ± 5 mW m −2 . This suggests that, similar to the Dybet site in northeastern Greenland, subglacial topography may play a role in the relatively high local geothermal heat flux measured at Dyer Plateau. Applying our inferred topographic-geothermal anomaly to the Dyer Plateau observation would yield a topographically corrected regional geothermal heat flux of 68 ± 27 mW m −2 . At Whillans Ice Stream, however, vastly different geothermal heat fluxes have been measured at Lake Whillans (285 ± 80 mW m −2 ; Fisher et al., 2015) and Whillans Grounding Zone (88 ± 7 mW m −2 ; Begeman et al., 2017), which are only separated by ∼110 km. With presently available BedMachine topography (Morlighem et al., 2019), we assess negligible topographic corrections for geothermal heat fluxes at both sites (<2%). This supports existing interpretations that the extreme geothermal heat flux observed at Lake Whillans may be characteristic of a larger region or reflect dynamic subglacial hydrology, rather than local bed topography as it is presently resolved (Gooch et al., 2016;Mikucki et al., 2016).

Subglacial Refreezing
Acknowledging increased spatial variability in geothermal heat flux may have counterintuitive implications for the interpretation of contemporary ice-sheet velocity. Accounting for topographically increased geothermal heat flux could influence the apparent thickness of deeply incised glaciers inferred from mass conversation (Morlighem et al., 2017(Morlighem et al., , 2019. Higher local geothermal heat flux, and warmer basal ice temperatures, can allow more internal ice deformation in regions with both frozen and thawed ice-sheet beds, as well as a greater fraction of surface ice velocity to be explained by basal motion where the bed is thawed. However, within major outlet glaciers, basal frictional heating (∼1,000 mW m −2 ) is generally an order of magnitude greater than geothermal heat flux (∼100 mW m −2 ). The basal ice of most major outlet glaciers in Greenland is expected to be at the pressure melting point, or thawed (MacGregor et al., 2016). Accounting for the topographic influence of geothermal heat flux is therefore likely to have trivial implications on COLGAN ET AL. Note. Site location is given in native decimal degree precision, as well as projected coordinates. mass-conversation derived ice thickness in settings with active subglacial hydrology where ice flows rapidly via basal sliding (Smith-Johnsen et al., 2020).
Instead, we suggest that the spatial variability in geothermal heat flux due to topography can potentially be important in slow-flowing inland ice-sheet areas. There, where frictional heat flux is closer in magnitude to the local geothermal heat flux (∼100 mW m −2 ; Marshall, 2005), the decreased geothermal heat flux associated with prominent subglacial ridges will make them colder relative to their surroundings. Subglacial ridges are also cooled by thinner overlying ice, which increases the pressure-melting-point temperature relative to warm valleys overlaid by thick ice. The combination of decreased local geothermal heat flux and increased local pressure-melting-point temperature is expected to result in sharp spatial contrasts in basal thermal conditions that promote preferential refreezing of subglacial water at cold subglacial ridges (Rezvanbehbahani et al., 2019).
As a case study, we assess spatial variations in topographically corrected geothermal heat flux and basal pressure-melting-point temperature along the Eqip Sermia ice-sheet profile presented in Bell et al. (2014). This profile depicts a massive refrozen basal ice unit extending downstream of a prominent bedrock ridge at km 22 ( Figure 13). We calculate basal pressure-melting-point temperature (T pmp ) as: where T 0 is the 273.15 K melting point of ice, β is a melting-point depression factor of 0.87 K km −1 , and H is the ice thickness (Marshall, 2005). Between km 4 and 22, the upstream onset zone of the basal ice unit, pressure-melting-point temperature increases 0.7 °C, from −1.3 °C to −0.6 °C. Over this same onset distance, topographically corrected geothermal flux decreases by >25%, from 84 to 62 mW m −2 . We suggest that this coincident increase in pressure-melting-point temperature and decrease in geothermal heat flux precondition the water and ice flowing over the subglacial ridge at km 22 to form a refrozen basal ice unit. We therefore speculate that some massive refrozen basal ice units not only require a local meltwater supply, but also cold ridges upon which this meltwater can refreeze (Bell et al., 2014).
The apparent spatial discontinuity of basal water identifications beneath the Greenland ice sheet represents an unresolved challenge to larger-scale but smoother inferences of that ice sheet's basal thermal state (Chu et al., 2018;Jordan et al., 2018;MacGregor et al., 2016;Oswald et al., 2018;Rezvanbehbahani et al., 2019). Similarly, the distribution of subglacial lakes beneath the Antarctic and Greenland ice sheets may partly reflect such basal thermal state transitions (e.g., Bowling et al., 2019;Smith et al., 2009). In this context, it is plausible that topographic variations in geothermal heat flux can influence the spatial distribution of frozen and thawed basal thermal states. Topographic variations in geothermal heat flux likely influence the basal thermal state of the ice-sheet interior at small scales (≲ 10 km) where both the bed is near the pressure-melting point and the relative topographic correction is significant (≳ 10%).
For Antarctic subglacial lakes in particular, we note that radar sounding only rarely images the flanks of subglacial lakes, so there are few constraints on their depths, except where seismic observations are available (Horgan et al., 2012). It is unlikely that these subglacial lakes mask exceptionally deep and narrow troughs needed to augment geothermal heat flux to exceptional levels, as subglacial troughs typically shallow substantially toward the ice-sheet interior (Morlighem et al., 2017(Morlighem et al., , 2019. Thus, identifying regions where topographic variations in geothermal heat flux appear to exert an influence on ice-sheet basal thermal state more likely serves to identify regions where basal temperature is near the pressure-melting-point, rather than to identify regions with unrecognized high-relief basal topography (Chu et al., 2018).

Model Limitations
The primary limitation of our model is the use of single values of the characteristic height and radius to capture the sensitivity of geothermal heat flux to topographic relief. The local topographic relief necessary to induce a 100% change in geothermal heat flux clearly depends on additional processes. These processes include both horizontal and vertical variabilities in the thermal conductivity associated with different rock or sediment types (Clauser & Huenges, 1995), variable rates of radiogenic heat production (Roy et al., 1968), spatial variability in groundwater saturation and hydraulic permeability (Saar, 2011), and spatial variability in paleoclimate (Westaway & Younger, 2013). For example, where geothermal heat flux is locally elevated due to one of these properties or processes (e.g., presence of radiogenic igneous bedrock, bedrock with anomalously high thermal conductivity, or no groundwater flow), the characteristic height we determined empirically would underestimate the topographic influence on geothermal heat flux. Similarly, as the contrast in thermal conductivities between underlying rock and overlying air, ice or water influences the magnitudes of thermal deflection, best-fit parameters can be expected to vary between subaerial, subglacial and submarine domains.
There are well-known geological differences between oceanic and continental crust types. Restricting the empirical data to just the 12 continental sites of previously observed or simulated topographically corrected COLGAN ET AL. geothermal heat fluxes, and thus excluding the three lacustrine and three oceanic sites, yields a similar characteristic height of 910 m with an averaging radius of 4 km. This best-fit, continental-only parameter combination also satisfactorily reproduces the 12 previously observed or simulated topographically corrected geothermal heat fluxes with an identical RMSE (Table 1). This makes the continental-only parameter combination statistically indistinguishable from the global (continental, lacustrine and oceanic) parameter combination we employ. We therefore conclude that our first-order geostatistical approach is relatively insensitive to geological subsets, although the presently available number of studies assessing kilometer-scale variations in geothermal flux is admittedly small. We expect future applications can derive more optimal length-scales, either by weighted-analysis of the historical observations that we have assimilated, or by evaluating against new observations. As different averaging radii yield different topographic reliefs, the characteristic height and averaging radius are clearly co-dependent. The 950-m characteristic height we discuss here is therefore valid only for a 5-km averaging radius.
The model parameters we adopted satisfactorily reproduce : (1) 18 previously observed or simulated topographically-corrected geothermal fluxes and, (2) a first-principles numerical solution of the heat equation at Disko Island, and (3) variation in borehole-derived geothermal fluxes in high-relief portions of the United States. However, their universal application across all of Greenland and Antarctica-regardless of known geology-clearly represents an imperfect first attempt to constrain the influence of topography on geothermal flux. Using a single characteristic height is more defensible in Greenland, where the subglacial geology appears to be more uniform than in Antarctica, where the geology is generally believed to be more complex (Dawes, 2009). Implementing a spatially variable characteristic height, whereby characteristic height becomes a function of rock properties associated with known underlying geology and overlying material, could overcome this limitation. In practice, however, there are substantial challenges associated with implementing a geology-dependent characteristic height.
Finally, the Greenland and Antarctic topographic data sets that we employ have nominal grid resolutions of 150 and 500 m, respectively. Their true spatial resolution, however, is less than this nominal grid resolution in most interior ice-sheet regions. We acknowledge that this coarser true resolution may result in "lower highs" for topographic correction within valleys, but given the self-consistent nature of our approach, it would also lead to "higher lows" for topographic correction along ridges. Thus, while stated topographic resolution may be optimistic in certain areas, the state-of-the-art topographic data sets we employ provide the best present opportunity to highlight geothermal heat flux dependencies on local topography.

Improving Constraints on the Topographic Effect
Our approach for inferring kilometer-scale spatial variations of the topographic correction for geothermal heat flux could potentially be validated through more intensive sampling of the uppermost sediment temperatures in steep Greenlandic fjords. Geothermal heat flux can be inferred by measuring the steady-state temperature profile and thermal conductivity of the uppermost few decimeters of submarine sediment, using gravity-driven probes equipped with pulse heaters and rapidly equilibrating temperature sensors (Hyndman et al., 1979;Pfender & Villinger, 2002). By contrast, sampling geothermal heat flux on land requires deeper drilling, often into permafrost, to measure the temperature gradient well below the penetration depth of the annual temperature cycle (Gruber et al., 2004).
Variations in geothermal heat flux due to crust or mantle properties are commonly assumed to occur over much larger spatial scales than the distance between two adjacent fjords. We hypothesize that it should be feasible to discern the topographic influence on geothermal heat flux by sampling geothermal heat flux within adjacent narrow (<2r or <10 km) and wide (>2r or >10 km) fjords of similar depth. For example, the relatively narrow Tyroler Fjord (∼2 km wide) on the north side of Clavering Island (the "Dybet" site of Rysgaard et al., 2018) with the wider (∼15 km) Godthåb Gulf just south of Clavering Island (74.09°N, 22.03°W). Within both fjords, a significant along-fjord gradient in geothermal heat flux would be unlikely, but their across-fjord gradients should differ significantly. A fjord as wide as Godthåb Gulf should approach a uniform across-fjord geothermal heat flux near its center.
Our work highlights that better subglacial mapping of the ice-sheet interior is important not only for better constraints of ice-sheet form and flow, but also for understanding the kilometer-scale topographic relief that can influence their basal thermal state and erosive potential (e.g., Lai & Anders, 2020). The geostatistical approach we adopted means that our understanding of spatial variability in geothermal heat flux will improve in concert with knowledge of subglacial topography. However, the recent conclusion of NASA's Operation IceBridge, which collected numerous ice-penetrating radar profiles across both ice sheets, means that the rapid growth in new subglacial observations witnessed over the past decade is likely to decelerate.

Conclusions
This study was motivated by the need to better understand the spatial distribution of the geothermal heat flux beneath Earth's ice sheets. It is widely understood that this property is poorly known, but it is of particular value to understanding past, present, and future ice-sheet evolution. We developed a topographic correction for geothermal heat flux based on a simple 2-D model for estimating the effect of kilometer-scale topography upon larger-scale geothermal heat flux fields inferred by other means. This model overcomes previous limitations of analytic methods by drawing on an empirical synthesis of terrestrial measurements of geothermal heat flux that directly considered topographic relief. This model reproduces existing observations, generates physically self-consistent patterns of topography-induced anomalies in geothermal heat flux beneath the Greenland and Antarctic ice sheets and is consistent with a numerical model. We identified key regions in both Greenland and Antarctica where the effect of topography upon local geothermal heat flux is equivalent to the ensemble uncertainty in its regional value.
This topographic correction is an additional, but unavoidable, layer of complexity in understanding geothermal heat flux. This underappreciated effect carries multiple consequences for modeling ice-sheet flow and challenges interpretation of both older, poorly georeferenced measurements of geothermal heat flux and newer ones in areas of high relief. We conclude that future studies of subglacial geothermal heat flux cannot credibly represent spatial variability in this physical property without also considering the effect of subglacial topography first highlighted by van der Veen et al. (2007).
There is presently a mismatch between the spatial scales needed to model ice-sheet flow accurately (<1 km; Aschwanden et al., 2016) and those of present ice-sheet-wide geothermal heat flux estimates (>10 km; Martos et al., 2017). Our study highlights the importance of bridging this gap and presents an imperfect first approach for doing so. Because our model sidesteps the 3-D time-varying heat-flow problem indicated by both the geometry and history of our polar study regions, there remains a significant opportunity for methodological improvements in topographic corrections of subglacial geothermal heat flux.

Data Availability Statement
The topography-dependent geothermal flux anomaly fields we calculated for Greenland and Antarctica, as well as their associated uncertainties, are now available on the Program for Monitoring of the Greenland ice sheet (PROMICE) data portal at http://doi.org/10.22008/promice/data/topographiccorrectiongeothermalflux/v1. The Greenland BedMachine v3 topographic data set we use is available for download from the National Snow and Ice Data Center (NSIDC) at https://nsidc.org/data/IDBMG4. The Antarctica BedMachine v1 topographic data set we use is also available for download from NSIDC at http://nsidc.org/data/nsidc-0756.