Studies of geothermal heat provide essential insight into the internal structure and history of the Earth (Kelvin, 1863; Pollack and Chapman, 1977; Beardsmore and Cull, 2001; Artemieva, 2011; Davies, 2013; Hasterok, 2013; Jaupart et al., 2016; Podugu et al., 2017; Lucazeau, 2019). A range of mechanisms control the amount of heat observed: thermal properties of the upper mantle, ongoing or recent tectonism (e.g. Pasquale et al., 2014; Goes et al., 2020), crustal heat production (e.g. Jaupart et al., 2016; Hasterok et al., 2018), topographic focusing by refraction (e.g. Lees et al., 1910), erosion and sedimentation (e.g. Von Herzen and Uyeda, 1963; Fukahata and Matsu’ura, 2001), advection by groundwater (e.g. Mansure and Reiter, 1979), and preserved variations and anomalies of paleoclimatic conditions (e.g. Huang et al., 1997; Šafanda et al., 2004).

Insights from geothermal measurements are applied in mineral prospecting (e.g. Cull et al., 1988), and hydrocarbon exploration, for example, to constrain the oil and gas windows (e.g. Royden and Sclater, 1980; Shalaby et al., 2011), and geothermal energy (Dickson and Fanelli, 2013). One particular aspect of geothermal heat flow that has recently gained attention is the impact of heat transfer at the base of ice sheets in Greenland and Antarctica (Burton-Johnson et al., 2020; Karlsson et al., 2021; Colgan et al., 2022). Even moderate geothermal heat can generate basal melting at the pressure melting point, reducing the friction of ice over bedrock or sediment and changing the rheology of the ice (Greve and Hutter, 1995; Pattyn, 2010).

Heat flow is calculated from the thermal gradient in a borehole combined with measurements or assumptions regarding thermal conductivity. Factors that impact the uncertainty and reproducibility of the recorded heat flow include the depth of the borehole, integration time for bore fluid to equilibrate to the surrounding temperature, assumptions regarding thermal conductivity, and groundwater flow (Beardsmore and Cull, 2001). In marine settings, short bottom-penetrating thermistor-line probes, up to 10 m, are commonly used for measurements in soft sediments. Such sensors require relatively limited equipment and infrastructure, and are quick to deploy (Hyndman et al., 1979; Dziadek et al., 2017). The technique involves the same set of assumptions regarding thermal conductivity and groundwater circulation as for borehole measurements, but further includes uncertainties due to local shallow anomalies, frictional heat from the probe and stabilisation time, and periodic variations in the water temperature. Some deep-water probes include a core sample that is used for estimating thermal conductivity (Gerard et al., 1962; Hyndman et al., 1979; Beardsmore and Cull, 2001).

Most measurements have been conducted in the northern hemisphere, particularly in Western North America and Southern Europe (Figure 1B). Measurements are costly and often sparse in remote areas (Figures 1C,D). A few techniques have been established to generate continuous maps where in-situ measurements are unavailable. Forward models compute heat flow values from thermal gradients modelled from geophysical data (An et al., 2015; Martos et al., 2017; Gard and Hasterok, 2021), energy balance (Stål et al., 2020), geological association (Davies and Davies, 2010; Burton-Johnson et al., 2017), or isostasy (Hasterok and Gard, 2016). Each of these approaches is associated with assumptions, particularly regarding to the crustal heat production and the strength of association between the dataset used and observed heat flow (Ebbing et al., 2009; Haeger et al., 2019; Lösing et al., 2020). A different approach has been to interpolate and model heat flow empirically by linking heat flow measurements elsewhere to a target through the similarity of one or many observables (Goutorbe et al., 2011; Lucazeau, 2019; Shen et al., 2020; Li et al., 2021; Lösing and Ebbing, 2021; Stål et al., 2021). This empirical approach shows promising and converging results in the case of Antarctica; however, the choice of observables and how well they capture thermal properties have been discussed and challenged (Davies and Davies, 2010; Stål et al., 2021; Artemieva, 2022), and further analysis will likely refine the choices made.

One aspect of empirical heat flow studies that has attracted less debate, but has the potential to impact the results, is the representation of the reference heat flow catalogue. Stål et al. (2021) suggest that a focus on economic exploration, particularly, in the Gondwana continents (e.g. Africa and Australia), could lead to sampling bias. Another potential factor that might skew the representation of the catalogue is a tendency to drill in flat and accessible locations within mountainous regions. Heat flow values are sometimes corrected for topographic factors, but not always, and it can be difficult to determine if such corrections have been applied in older studies. Such clarifications are within the scope of the ongoing Global Heat Flow Data Assessment Project of the IHFC catalogue (Fuchs et al., 2021a).

Heat flow measurements have been collected in cumulatively growing databases (e.g. Chapman and Pollack, 1975; Pollack et al., 1993; Hasterok, 2019; Lucazeau, 2019). From analyses of those catalogues, the total heat balance and average heat flow, of Earth have been calculated. The spatial bias has been recognised and causal relationships are used to integrate heat flow from lithospheric age in the oceans (e.g. Lucazeau, 2019) or geological setting (e.g. Davies and Davies, 2010).

In this contribution, we analyse the most recent and extensive heat flow catalogue (Fuchs et al., 2021b) for factors that might bias the distribution and hence impact the integrated heat flow maps. We also discuss qualitative reasons for the uneven distribution and suggest statistical weights of individual samples for use when the catalogue finds ongoing use.

Spatial characteristics of the data are analysed to quantify clustering and misrepresentation. The analysis is carried out using Python libraries: geopandas (Jordahl, 2014), rasterio (Gillies, 2019), and numpy (Harris et al., 2020). Methods are implemented from agrid (Stål and Reading, 2020), a python-based grid for representing multidimensional geophysical data. All code is made available to ensure reproducibility, and all datasets used are provided in open repositories.

The map in Figure 1A shows the records in the IHFC. The spatial distribution is also shown as a kernel density estimation in Figure 1B. This smoothed distribution highlights that the highest density of measurements is in the Western United States and Southern Europe. The inset to the bottom-left of Figure 1B shows centered Ripley’s K function. For context, the expected value of the K function for spatially uniform sampling is K ̂ ( t ) = 0 . The measured K function is positive on all scales, indicating clustering.

Figure 1C shows the distance to the nearest IHFC heat flow record, measured from grid cell centres. Central Africa, the Amazon Basin, and parts of the Middle East, have extensive areas with a distance of over 500 km to the nearest record. In parts of interior Africa there are areas with over 1000 km, and up to 1359 km in South Pacific. Figure 1D shows the mean distance to the closest ten measurements. The overall distribution is similar; however notably, the Southern Ocean is highlighted as having only a few measurements representing large areas. For both metrics, Antarctica is exceptionally sparsely surveyed.

In Figure 2, we show the reference and sampled distributions, and the heat flow associated with each class for the eight categorical auxiliary variables. The calculated weights are listed. The horizontal bar charts show the heat distribution within each class. Generally, the mean heat flow values tend to be much higher than the median due to extremely high measurements in active geothermal settings. The difference between the median and mean values of IHFC and the weighted average heat flow indicates the magnitude of the impact on heat flow models from sampling bias. We also calculate the robust mean, excluding 1 and 10% upper and lower percentiles (Supplementary material Table S1).

To better understand the origin of the sampling bias, we extract the economic setting for the measurements: 17.8% of the measurements in IHFC are within the polygons defined as oil and gas fields (Rose et al., 2018), in relation to only 8.7% of the landmass. Meanwhile, 7.3% of the measurements are within 0.5° from a mine, as mapped (Maus et al., 2020), in relation to 12.3% of the global landmass, excluding oceans. Prospective regions are only slightly over-represented on a global average; however strongly pronounced in sparsely populated Australia, where 26.2% of the measurements in IHFC are located in mining areas (c.f. 22.1% by landmass), as defined above, and 38.3% of the measurements are within oil and gas fields (c.f. 5.6% by landmass, Rose et al. 2018). Moreover, many of the remaining measurements are in regions targeted for geothermal heat extraction, e.g. North West Tasmania and South-Western Victoria (Holgate et al., 2010; Bahadori et al., 2013). Figure 3C shows the Australian records in IHFC and the polygons used to define mining and oil and gas fields.

Figure 4 shows the weights calculated. Figure 4A shows the weights derived from seismic tomography regionalisation (Schaeffer and Lebedev, 2015), highlighting the general under-representation of the oceanic crust. Figure 4B shows the smaller weights from geomorphometrics (Amatulli et al., 2020). For this analysis, local rather than global and reg ional distribution impacts the weight. Figure 4C shows the weights based on the province (Hasterok et al., 2022). Figure 4D shows the weights from lithologies (Hartmann and Moosdorf, 2012). Oceans are set to a weight of 1. All ten maps are provided in Supplementary material Figure S4.

We now consider selected weights derived for the classes in each auxiliary variable using IPF, as described in Section 2.5.

The IPF calculated weighting from TYPE and TOPO produce a weighted average of 84 mW m ⁻², this is closest to the global mean (Pollack et al., 1993; Lucazeau, 2019), and also the lowest weighted mean for any fitted weighting. The weighted median is 58 mW m ⁻², which is lower than the median of IHFC, 64 mW m ⁻² (Figure 4E). IPF calculated weighting from the three variables CRUST, TYPE, and PLATE yield a weighted mean of 93 mW m ⁻², which is also close to expected global average. The weighted median is 60 mW m ⁻² (Figure 4F). An IPF calculated weighting from a combination with the potential to capture both tectonic and geological misrepresentation, TYPE and GLiM, gives a weighted mean of 95 mW m ⁻², and a weighted median is 58 mW m ⁻² (Figure 4G). An IPF calculated weighting from TYPE, GLiM and TOPO represents tectonic, geological and topographic settings. The weighted mean for this combination is 145 mW m ⁻², and the weighted median is 59 mW m ⁻² (Figure 4H).

From all variables considered, some extreme weights are suggested. In order to compensate for the sparse measurements from GLiM class glaciers and ice sheets, measurements for those two classes are calculated to have a weight of 813.9 (Figure 2). Other underrepresented tectono-geographical regions include, as a most recent orogen, the Birimian Orogen, Kuunga Orogen and particularly the Scotian Orogen, where no heat flow measurements are catalogued. Generally, Gondwana and the oceanic crust are underrepresented.

We have shown that the distribution of heat flow measurements is spatially biased and does not fully represent the Earth’s geometric, tectonic, geological or geomorphometric disposition. This conclusion is readily seen from the clustering of measurements (Figure 1B), and in the context of eight auxiliary variables, whereby we found that some settings are overrepresented whilst others are correspondingly underrepresented. The selection of those variables is somewhat arbitrary; however, together, they represent a meaningful range of parameters that should be expected to impact the distribution of geothermal heat. The scope of this contribution is not to analyse the relation between observables and geothermal heat, as has been done in previous studies (e.g. Davies and Davies, 2010; Goutorbe et al., 2011), but to investigate how well the heat flow catalogue represents the Earth’s surface.

The calculated weights from each variable compensate for the bias, and the combined fitted weights provide individually optimised weights for each record in IHFC. The somewhat diverging results highlight the need for caution when weighting is applied in empirical models. The choice of weights depends on the scale of the model and relevance of the variables used. It would be tempting to evaluate the weights based on their performance to reduce prediction misfits in existing empiric models; however, such an appraisal will necessarily contain the same bias. A useful test case could be cross-validation of a sufficiently large subset of the heat flow catalogue that can be shown to truly represent a random sampling of the Earth, and of good quality. Unfortunately, given that some settings are strongly underrepresented, such a sample distribution is not yet possible. With empirically driven heat flow models, there is potential to further explore how the weighting and processing of reference measurements can impact the results. Some observables, such as the distance to volcanoes, have been shown to have large impact (Lösing and Ebbing, 2021); however, the sensitivity changes if weights are applied, and in many cases extreme geothermal settings are weighted down.

One potential shortcoming in this study is that we have assumed that the coordinates of the measurements are given correctly and accurately to match the auxiliary variables; this is often not the case. The surprisingly low bias in the geomorphometrics variable might be, at least partially, explained by imprecise positions that sample random geomorphological forms in the vicinity rather than the actual topographic shape right at the borehole.

From an assumption that weighted mean and median values from the catalogue should approach the qualitative estimates of the global heat flow distribution, the closest weighted mean is given precedence in our interpretation. The total Earth heat loss is estimated to be 40–42 TW, or 80 mW m ⁻² (Lucazeau, 2019). Earlier studies come to similar results (e.g. Pollack et al., 1993). Most weighted mean values are closer to this average, suggesting that weight applied to metrics of the catalogue can improve the prediction if applied carefully.

Weights derived from one auxiliary variable, that generate weighted averages close to the expected global value: Last orogeny (OROGEN) average is +12 mW m ⁻² compared to global average, shallow geology (GLiM), +13 mW m ⁻², geomorphometric shapes (TOPO), +26 mW m ⁻², and tectonic plate, +26 mWm⁻². The fitted weight from plate type (TYPE) and topography (TOPO) (Figure 4E) produces a close value of +4 mWm⁻², however we have concerns regarding the validity of the geomorphometrics (TOPO), as the coordinates given in the catalogue might not be precise enough. Moreover, marine data is note weighted to topography. Plate type (TYPE) and shallow geology (GLiM) also generate a weighting that is similar to the global average at +16 mWm⁻² (Figure 4G). A qualitative reasoning of what processes would impact heat flow on different scales, is represented by weighting for tectonic setting (TYPE), shallow geology (GLiM) and topography (TOPO). This combination gives a weighted mean difference of +65 mWm⁻² (Figure 4H).

The spatial sampling bias is consistent with the challenges and lack of incentives to conduct investigations in remote regions or developing countries. More measurements of the thermal gradient have likely been conducted for prospecting and exploration reasons in some of such sparse areas; however, the results are not included in heat flow compilations and the understanding of some areas suffers from very sparse records (Figures 1C,D, (Discussed by e.g. Brigaud et al., 1985; Lesquer and Vasseur 1992; Noorollahi et al., 2009; Yousefi et al., 2010)). We also suggest that the quality of the measurements varies spatially. Some assumptions might be outdated in areas where older measurements dominate (Figure 1A). If one were to use borehole data alone, the marine portions of the map would be even more sparsely sampled, except near oil and gas provinces.

Acknowledging the achievement of the cumulative heat flow catalogues, the scope of this contribution is to support and add value to the ongoing, substantial undertaking to coordinate a representative compilation of measurements.

We make six recommendations from the results of this study:

1) Empirical heat flow model developers should consider applying weights when using the reference database and investigate how this can reduce uncertainty and misfit.