the Creative Commons Attribution 3.0 License. Solid Earth Earth’s surface heat flux

Abstract. We present a revised estimate of Earth's surface heat flux that is based upon a heat flow data-set with 38 347 measurements, which is 55% more than used in previous estimates. Our methodology, like others, accounts for hydrothermal circulation in young oceanic crust by utilising a half-space cooling approximation. For the rest of Earth's surface, we estimate the average heat flow for different geologic domains as defined by global digital geology maps; and then produce the global estimate by multiplying it by the total global area of that geologic domain. The averaging is done on a polygon set which results from an intersection of a 1 degree equal area grid with the original geology polygons; this minimises the adverse influence of clustering. These operations and estimates are derived accurately using methodologies from Geographical Information Science. We consider the virtually un-sampled Antarctica separately and also make a small correction for hot-spots in young oceanic lithosphere. A range of analyses is presented. These, combined with statistical estimates of the error, provide a measure of robustness. Our final preferred estimate is 47±2 TW, which is greater than previous estimates.


Introduction
Heat flow measurements at Earth's surface contain integrated information regarding the thermal conductivity, heat productivity and mantle heat flux below the measurement point. By studying Earth's surface heat flux on a global scale, we are presenting ourselves with a window to the processes at work within Earth's interior, gaining direct information about the internal processes that characterize Earth's 'heat engine'. The magnitude of the heat loss is significant compared to Correspondence to: J. H. Davies (daviesjh2@cardiff.ac.uk) other solid Earth geophysical processes. Consequently, the study and interpretation of surface heat flow patterns has become a fundamental enterprise in global geophysics (Lee and Uyeda, 1965;Williams and Von Herzen, 1974;Pollack et al., 1993).
The global surface heat flux provides constraints on Earth's present day heat budget and thermal evolution models. Such constraints have been used to propose exciting new hypotheses on mantle dynamics, such as layered convection with a deep mantle interface (Kellogg et al., 1999), and that D" is the final remnant of a primordial magma ocean . One class of thermal evolution models are the so-called parameterised convection models (Sharpe and Peltier, 1978). While such models have been examined for several years, they have recently been reenergised, with work in fields including: (i) alternative thermal models (Nimmo et al., 2004;Korenaga, 2003); (ii) refined estimates of mantle radioactive heating (Lyubetskaya and Korenaga, 2007a); (iii) estimates of core-mantle heat flow, following from observations of double crossing of the perovskite/post-perovskite phase transition (Hernlund et al., 2005); and (iv) advances in models of the geodynamo (Christensen and Tilgner, 2004;Buffett, 2002). Earth's global surface heat flux plays a fundamental role in all of the above.
A comprehensive estimate of the global surface heat flux was undertaken by Pollack et al. (1993) (hereafter abbreviated to PHJ93), producing a value of 44.2 TW ± 1 TW, from a data-set of 24 774 observations at 20 201 sites. Until the recent work of  (hereafter abbreviated to JLM07) this was widely regarded as the best estimate. JLM07 revisited this topic, making alternative interpretations of the same heat flow data-set. One of their major contributions was in reassessing the corrections required for, and errors in, a reasonable estimate of Earth's total surface heat flux. Their final value is 46 TW ± 3 TW.
In this study, we use Geographical Information Science (GIS) techniques, coupled with recently developed highresolution, digital geology maps, to provide a revised estimate of Earth's surface heat flux. We employ a significantly (55%) larger data-set (38 347 data points) than previous work. While we focus on our preferred value, we also present a range of possible values based upon a variety of assumptions. This, combined with careful estimation of errors, provides a rigorous error estimate.
The structure of the paper is as follows. A brief introduction to the work of PHJ93 is first presented, along with some background to our related methods. The heat flow and geology data-sets used are then described. This is followed by a description of our methodology, including the corrections, which are guided by JLM07. The actual results and, importantly, the errors, are then presented and discussed. We conclude by discussing our preferred estimate of 47 TW, (rounded from 46.7 TW given that our error estimate is ±2 TW) which is greater than recent estimates. We note, however, that this value overlaps with JLM07, when considering the associated errors (±2 TW). Such error estimates are of great importance; our total error of around 2 TW is less than JLM07 (3 TW) but double the 1 TW of PHJ93.

Methods
Heat flow observations are sparse and non-uniformly distributed across the globe. Pollack et al. (1993) (PHJ93) show that even on a 5×5 degree grid, observations from their dataset cover only 62% of Earth's surface. Therefore, to obtain a global estimate of surface heat flow, PHJ93 derived empirical estimators from the observations by referencing the heat flow measurements to geological units. The underlying assumption is of a correlation between heat flow and surface geology. In this way, it is possible to produce estimates of surface heat flow for regions of the globe with no observations (this could also be thought of as an area weighted average based upon geological category). PHJ93 did this by first attributing every 1 × 1 degree grid cell to a specific geology. The heat flow data in each cell was then averaged and resulting cell values were used to estimate an average heat flow for each geology. An estimate was produced globally for each geological unit by evaluating its area in terms of 1 × 1 degree grid-cells and multiplying by the estimated mean heat flux derived for that geology. The total heat flux was evaluated by summing the contribution of each geological unit. Finally, estimates were made for hydrothermal circulation in young oceanic crust, using the model of Stein and Stein (1992) (hereafter abbreviated to SS).
There has been a major revolution in the handling of spatial data in the intervening years, with the growth of GIS. GIS allows geological units to be defined by high-resolution irregular polygons in digital maps. The highest resolution geology data-set (a combination of two data-sets) utilised in this study has over 93 000 polygons. GIS allows us to evaluate the areas of these geology units exactly (to mapping accuracy) and to also match heat flow measurements with specific local geology. PHJ93 had to estimate the area of geological units by dividing Earth's surface into 1 × 1 degree equal longitude, equal latitude cells. They then hand-selected the predominant geology of each cell and summed the number of cells. Such a methodology could potentially generate errors in the estimates of geological unit areas. In addition, in cells with greater than one geological unit, this could lead to heat flow measurements being associated with the incorrect geology. Further advantages of GIS methods include the ability to: 1. Easily match heat flow measurements with individual geology polygons; 2. Undertake accurate re-calculations with different grids and for different geology data-sets; 3. Use equal area grids (rather than equal longitude grids); 4. Make robust error estimates, weighted by area.

Heat flow data-set
The heat flow data-set utilised in this study (which we abbreviate to DD10) was provided by Gabi Laske and Guy Masters (Scripps Oceanographic Institution, La Jolla, California, USA) in the autumn of 2003 (personal communication). It subsumes the data-set of PHJ93 and has been supplemented by a large number of observations made in the intervening years, giving a total of 38 374 data-points (this is a 55% increase from the 24 774 points of PHJ93  In this study, we use geology for this purpose, following the work of Pollack et al. (1993). (b) Focussed on the African continent -note the sparsity of data points. (c) Focussed on Europe, where good coverage is apparent, particularly in areas such as the Central North Sea, Black Sea and Tyrrhenian Sea, where interest in surface heat flow has been extensive due to exploration and tectonism. Figure 1 illustrates the global distribution of the dataset, clearly showing the inhomogeneous spread of measurements. Figure 2 includes a breakdown of the new data-set into new points, those points included in PHJ93, and those modified from PHJ93, while histograms of the heat flow measurements are given in Fig. 3. We stress, like PHJ93, that while the histograms of the ocean and continental heat flow measurements look similar, this is misleading. The oceanic region is dominated by sites with sediment cover and these are known to be biased systematically downwards by hydrothermal circulation (Davis and Elderfield, 2004).
In addition, the uneven geographical distribution noted in Fig. 1 should make one cautious in making global estimates directly from the raw data. We follow PHJ93 in methods to address the issue of hydrothermal circulation and un-sampled regions, though our implementation is different.

Geology data-sets
We utilise two geology data-sets: 1. A global data-set, CCGM/CGMW (Commission for the Geological Map of the World, 2000), abbreviated to CCGM (the French initials for the data-set -Commission de la Carte Géologique du Monde). It ascribes every point on Earth's surface to a geological unit (see Fig. 4).

2.
A data-set of continental geology (Hearn et al., 2003), abbreviated to GG -for Global GIS. It includes virtually all land above sea-level, excluding Antarctica and Greenland (see Fig. 5).
The CCGM has 14 202 geology polygons, while the GG data-set has 91 964 polygons. Therefore, GG has a much higher level of detail, especially in the USA. When using the GG data-set, the CCGM data-set is used for areas with no coverage (i.e. when the GG data-set is used, we also use 1066 geology polygons from the CCGM data-set to represent the absent areas with heat flow observations -e.g. oceans).  assigned to one of three classes (either igneous, sedimentary or other (endogeneous -plutonic or metamorphic)). For GG, the continental rocks are divided by geologic period, with no further division according to rock type. The GG classification therefore has more periods than PHJ93 but does not subdivide them according to rock type. Figure 6 shows the geology, together with three different types of grid. One can see that even at the 1 × 1 degree scale many cells contain more than one geological category. PHJ93 used this scale to define the geology. While we have already commented on the fact that there is a strong variation in the density of heat flow observations we should also note that there is a strong spatial variation in the detail of geological classification provided by the digital data-sets. This reflects the varying geological mapping that has been undertaken in different regions of the world, in addition to the intrinsic variability of global geology.

Grids
PHJ93 used a 1 × 1 degree equal longitude grid (64 800 grid cells) (see Fig. 6a for an example) in their analysis. We have undertaken our preferred analysis using a 1 × 1 degree equal area grid, with 41 252 cells in total. These cells are 1 × 1 degree at the equator, but at pole-ward latitudes the cell longitudes increase to approximate equal area. In this way, each cell has the same and equal weight (see Fig. 6b for an illustration of this grid in the North Atlantic). The difference between an equal area and equal longitude grid is greatest at high latitudes. Since there are limited heat flow observations at high latitudes, we expect that the improvement of using an equal area grid might be limited in this study. Amongst our range of investigations we also used a 5 × 5 degree equal longitude grid (see Fig. 6c).

Analysis
To better illustrate the impact of various aspects of the methodology, we have undertaken a series of alternative analyses, giving us a handle on the level of uncertainty in our estimate. We shall next describe, in detail, the methodology used in our preferred analysis, since this is the most complex. This will allow us to more easily describe the other analyses, without having to repeat the complete description of each stage: 3. The resulting polygons were spatially joined to the heat flow data. In the spatial join, the attributes of the geology polygon containing the heat flow point is added to the table of attributes of the heat flow points (i.e. each heat flow point has its geology associated with it -see Fig. 8).
4. The mean heat flow for each geology polygon (unioned with the grid as described above) was calculated, producing a new output (summary) table (see Fig. 9). The geology label of each polygon and the average heat flow was stored in the summary table.
5. The mean heat flow value was calculated for each geology class using the summary table of the previous step. This, again, was a straightforward arithmetic mean of all polygons that had non-zero heat flow polygons for each specific geology.
6. The global area of each geological unit was calculated.
7. The final estimate of the global average heat flow was evaluated by assuming that the average heat flow found for each geology could be assigned to all similar geology (i.e. for each geological category, its average heat flow was multiplied by its global area to find its contribution to the global heat flux).
We follow JLM07 and PHJ93 and make an estimate for hydrothermal circulation, assuming that the heat flow in young oceanic crust is best described by a half-space model. Unlike PHJ93 who used the parameters of SS, we have used the value suggested by JLM07, but have also repeated the analysis using the parameters of Parsons and Sclater (1977) (hereafter abbreviated to PS) and SS, to examine the differences and uncertainties arising from these models. Thus, for all geology younger than 66.5 Ma (66.4 Ma for 1983 timescale) we have replaced the heat flow obtained from the raw data with a value obtained from the equation: where q is surface heat flux (mW m −2 ), C is a constant (mW m −2 Myr 0.5 ), and t is the age of the oceanic lithosphere in millions of years. The value of C preferred by JLM07 is 490 ± 20 mW m −2 Myr 0.5 ; while the values of SS and PS are 510 and 430 mW m −2 Myr 0.5 , respectively. The error of JLM07 is small, partly because they use the constraint that at infinite age the half-space model should predict zero heat flux. While that is certainly correct, it might be over optimistic to believe that a half-space model with a single constant is the correct model, at least as fit through all the data selected over such a wide range. We note that JLM07 ignore data at old age, since the half-space model is known not to fit that data well (that miss-fit led to the development of plate models), and at very young age, since the high variability reduce their usefulness in constraining the parameters. As a result we take a slightly more conservative estimate of the error in C and assume errors 50% greater (i.e. ±30 mW m −2 Myr 0.5 at 2 standard deviations). Table 3 shows the results for the 56 individual geology units (the 20 units on the continents from GG, and the 36 units from CCGM that represent the rest of the globe not included in the GG data-set) for this case. As described above, the final result is produced by multiplying the average heat flow calculated for each geology by the global area for that geology. Notice that some geological units have no (or very few) heat flow measurements. However, these make up only a very small proportion of the total area (1.5%, for less than 50 readings (excluding the Glacier category, which is discussed below and makes up ∼3% of area)). For cases of geology with insufficient readings for the area to be included, we have shared the area among similar regions.
Taking the standard deviation of the heat flow values contributing to the estimate of a mean for a geological unit as the basis for estimating the error in the mean, we find that the resulting errors, weighted by area, are small (i.e. ±2.2 and 3.0 mW m −2 , for the GG and CCGM data-sets respectively, 2 std. dev.). There is, however, additional uncertainty, arising from: (i) the inaccuracy of the area; (ii) the fact that the extrapolating method is not perfect (i.e. the geology class  Fig. 8. Schematic illustrating the spatial join between a layer with three points and a second layer of two polygons. The process adds the attributes of the appropriate polygon that contains the point of interest. This is how we assign a geology to each heat flow measurement.

Summarise on Geology
Original  Fig. 9. Schematic illustrating the process of producing a summary table. A target field is selected, in this case geology. The process "summarises" selected fields in defined manners e.g. mean, sum, standard deviation, that correspond to each unique entry in the selection field. In the example illustrated, the mean heat flow has been summarised for each geology category. The summary table will also contain the number of rows (heat flow measurements) contributing to each summary value.
is not a perfect predictor of the heat flow); and (iii) the fact that large areas of Africa, Antarctica and Greenland are unsampled. Errors in the area could arise due to poor definition of geological boundaries in the digital geology maps. We do not estimate such errors here. However, we find that our estimate of continent (ocean) area is slightly greater (less) than JLM07, but as JLM07 point out, since the total area is fixed, these differences have only a small effect. The issue of the inadequate extrapolation and the poor coverage is more significant as a source of error and is later discussed in detail.
We have decided to include a value for the Glacier category in our preferred analysis. The Glacier category of the CCGM geology covers 3% of Earth's surface area, primarily across Greenland and Antarctica. Depending upon the exact methodology used, we get between 2 and 5 final readings (grid, geology), and a mean heat flow value of between 105 and 120 mW m −2 , which is probably far too high a value. An estimate of 65 mW m −2 was recently made for the heat flux in Antarctica, from estimates of the depth to the Curie temperature, based upon magnetic field measurements (Maule et al., 2005). The work of Shapiro and Ritzwoller (2004) also refers to the problem of estimating the heat flow of Antarc-tica, in this case using seismic measurements as a proxy. They suggest that the heat flow in West Antarctica is almost a factor of three higher, and more variable (more like the small number of actual heat flow observations in our data-set), than in East Antarctica, where the heat flow has an estimated "local mean" of 57 mW m −2 . Since East Antarctica has a larger area than West Antarctica (by a factor of ∼3), the work of Shapiro and Ritzwoller (2004) would also argue for an overall value for Antarctica closer to 65 mW m −2 rather than the ∼105-120 mW m −2 given by the raw measurements. While this is similar to current predictions of the average heat flow through continents, it must be viewed as an uncertain estimate. However, our preference is to use this estimate and (105-65) mW m −2 as the two standard deviation estimate of the error (105 mW m −2 being the minimum direct estimate from the data).
In Table 4 we list the heat flow and geology data-sets, the grids, and the methods used for the various alternative analyses undertaken. Each case is next described in detail: 1. While we have not repeated the work of PHJ93, in our first analysis, we used their heat flow data-set (taken from Gosnold and Panda, 2002), a 5 × 5 degree equal longitude grid and their methodology of selecting geology for the whole of a cell, based on the majority geology of that cell. However, we used the CCGM geology data-set.
2. As in Case 1, but using the new heat flow data-set (DD10).
3. As in Case 2, but a spatial join was undertaken between the heat flow data and the underlying geological polygons.
4. As in Case 3, but, in addition, we undertook a union between the geology and a 1 degree equal area grid.
5. As in Case 3, but with the combination of GG and CCGM geology data-sets.
6. The preferred analysis described above.
Solid Earth, 1, 5-24, 2010 www.solid-earth.net/1/5/2010/  Table 4. Description of the various analyses undertaken. 5 deg eq lon; a grid with 5 degree spacing in longitude. 1 deg eq area; a grid with 1 degree spacing in longitude at the equator, but varying longitude at other latitudes to maintain an equal area. Majority Geology; the geology which takes up the greatest area inside the grid cell is ascribed to the whole grid cell. These and other aspects related to this table, including union, are described further in the text.

Results
In Fig. 10, we plot the heat flow map of the world from our preferred analysis. The standard error (2 std. dev.) for the heat flow ascribed to each geological unit is presented in Fig. 11. We note that the error is highest for the young ocean estimates and the Glacier domains. Of course, these estimates of error are not useful locally; for example, some parts of Africa have no measurements. Therefore, like the value on the Global heat flow map being only indicative for that geology unit, likewise the error. Results from the various cases examined are listed in Table 5. One can see that the raw global heat flow in each case varies from 35.8 to 36.7 TW, with our preferred analysis giving a value of 36.7 TW. After the Stein and Stein (1992) SSbased young ocean estimate, one gets values ranging from 44.1 to 47.2 TW, with 47.2 TW for our preferred analysis; or 40.7 to 43.5 TW using the Parsons and Sclater (1977) PSbased estimate (43.4 TW for our preferred analysis). We note that the difference between estimates using only the raw data, and those using the SS half space model can vary between 7.8 and 11.1 TW, with a difference of 10.5 TW for our preferred analysis. We also note that using the SS and PS models leads to a difference of 3.4 to 3.8 TW (3.8 TW preferred analysis). Our preferred correction is between those of PS and SS, giving a total heat flow value of 46.7 TW.
The preferred analysis can be divided into four components (see Table 6): (i) young oceans; (ii) the rest of the oceans and continents; (iii) the Glacier category; and (iv) the contribution from hot-spots. Each is next described in detail. and (iii) the fact that ours covers oceanic crust out to 66.5 Ma while JLM07 extends out to 80 Ma), while Wei and Sandwell (2006) give a value of 20.4 TW. The PS model, with our area, leads to a young ocean heat flow estimate that is 2.8 TW less than that adopted here, while the SS model gives a value 1 TW greater. Considering the PS-model estimate might be too low, it is clear that an uncertainty of ∼2 TW is suggested by the alternative analyses listed above; our assumption that the error in C is ±30 mW m −2 Myr 0.5 at 2 standard deviations, gives this component an error of ±1.3 TW.
2. Estimates obtained for the rest of the oceans and continental components depend upon: (i) the fundamental assumption of a correlation between heat flow and geology; and (ii) the use of a 1 degree equal area grid to reduce the influence of clustering. Our preferred analysis predicts values of 13.8 TW/73 mW m −2 (continents), 7.8 TW/66 mW m −2 (rest of oceans); for this component of the heat flow. This compares to 13.2 TW/65 mW m −2 (continents), 7.6 TW/56.4 mW m −2 (rest of oceans) from PHJ93 and 14 TW/65 mW m −2 (continents), 4.4 TW/48 mW m −2 (>80 Ma) from JLM07. A significant percentage of the increased global heat flux in this study therefore arises due to this component, suggesting that heat flow values recently added to the global heat flow data-set have been, on average, slightly higher than earlier values (as illustrated in the histogram of the raw data of Fig. 3b). This continues the slight upward trend in the estimate for global heat flux over recent years (see Table 4 PHJ93, and JLM07).
3. As described above, for the Glacier category we have used a value of 65 mW m −2 based on the depth to the Curie temperature found from undertaking a spectral analysis of aeromagnetic data (Maule et al., 2005), rather than the very small (2 to 5) measurements that fall within this category in the various analyses (the raw data gives heat flow estimates of 105 to 120 mW m −2 for the various analyses -113.5 mW m −2 in the preferred analysis). This gives a value of 0.9 ± 0.3 TW. We note that using this alternative value for Antarctica (and Greenland), rather than the raw heat flow measurements, can make a difference of between 0.7 and 1.0 TW (0.7 TW in our preferred analysis). PHJ93 do not really address this issue while JLM07 include it in their total continental area and effectively use the global average, which is 65.3 mW m −2 . JLM07 include the error from this component in their total error estimate (for continents), which is 1 TW.
4. Hot-spot regions typically do not exhibit higher heat flow values compared to that expected for their crustal age, excluding a few sites at or near volcanic features (Stein and Von Herzen, 2007). Nonetheless, , utilising the underlying estimates for each geology category. Note that in regions with no data (see Fig. 1) the estimate is based solely on the assumed correlation between geology and heat flux. As a result, locally, the values presented could provide a poor estimate. Estimates in better-covered areas of the globe, however, suggest that on average, the current estimates are robust; higher resolution plots, focussing upon Europe (b) and Japan (c) are also shown. JLM07 include a contribution for hot-spots in their analysis, which they argue are not accounted for in their method. They estimate that this contribution is between 2 and 4 TW globally, based on Davies (1988) and Sleep (1990); with an error of ±1 TW. In our estimate for heat flow across young oceanic crust, we assume that the calibration measurements for the models have been selected to avoid hot-spots and, therefore, the effect of hot-spots is not included (thus a correction is needed). In contrast, we feel that hot-spot anomalies are included in the rest of the measurements on the ocean floor and continent. As a result, we only include a hot-spot correction for the young oceanic domains, which is proportional to the surface area included in the young ocean estimate. We take a contribution of 1 ± 0.33 TW.
As mentioned above, some geology classes have very few heat flow measurements. We have looked at limiting our estimates to only geology classes with: (i) at least 5; and (ii) at least 50 readings. Such a change has only a small effect on the final global estimate, of order ±0.3 TW (see Table 5). However, the random errors are reduced substantially by restricting the analysis to geology classes with at least 50 readings, although the errors in the young ocean estimates are much higher. This is because geology classes with less than 50 readings cover only a small percentage of Earth's surface. Indeed, the remaining classes (i.e. with >50 readings) cover >96% of Earth's surface. Our preference is to ignore geology classes with <50 readings. Our final preferred value is 46.7 TW, with an error of 2 TW (2 standard deviations). To be consistent with our error estimate of ±2 TW we round 46.7 to 47 TW. This final preferred estimate is separated into the oceanic and continental realms in Table 7.

Discussion
In this section we primarily focus on various components that lead to the uncertainty (±2.0 TW) in our final estimate. Aspects considered include: (i) the fundamental correlation assumed between geology and heat flow; (ii) potential issues with the heat flow data, including the raw measurements and the areal coverage of the data-set; (iii) the implications of the various methods undertaken; (iv) issues related to hydrothermal corrections in oceanic crust; (v) hydrothermal circulation on continents; (vi) the glacier category; and (vii) combining statistical errors. We end by discussing the implications of our preferred value of 47 TW. Table 5. The total global surface heat flux from the various analyses undertaken. Note that the first row of results is from Pollack et al. (1993). Pollack is used as an abbreviation for the heat flow data-set, the geology and the method of using geology, used in Pollack et al. (1993). SS - Stein and Stein (1992); PS - Parsons and Sclater (1977); S. Join -spatial join; 1 × 1 long -1 degree by 1 degree equal longitude grid; 1 × 1 area -1 degree by 1 degree equal area grid; CCGM -Commission Geology Map of World (2000); GG -Global GIS (2003); DD10 -the data-set presented in this paper. The original version was provided by Laske and Masters (pers. comm.) and DD10 has only minor changes. See text for further explanation of terms and abbreviations used in table.
Case number from Table 4 Source of Heat Flow Dataset Geology  Table 6. The preferred analysis of Earth's surface heat flux, broken down into five categories: (i) continents from data; (ii) oceanic crust contribution by calculation (age <66.5 Myr); (iii) oceanic crust from data (age ≥66.5 Myr); (iv) the glacier category; and (v) the hot-spot contribution. Columns of alternatives are also presented. "+1" signifies that this answer uses a 1-degree equal-area grid through a "union", to reduce the influence of clustering. "50" means that only geology categories with at least 50 heat flow measurements were used in this answer. GG -Global GIS geology (Hearn et al., 2003); CCGM -CCGM geology (Commission for the Geological Map of the World, 2000); PHJ93 -result from Pollack et al. (1993) (breakdown calculated by combining information in their Tables 2 and 3); JLM07 -result from ). C510, C490, C430 -estimation using C values in Eq. (1) of 510, 490 and 430 mW m −2 Myr 0.5 , respectively (SS - Stein and Stein (1992), where C was 510 mW m −2 Myr 0.5 ; PS - Parsons and Sclater (1977), where C was 430 mW m −2 Myr 0.5 ); WS -result for oceanic crust heat flow younger than 66 Myr from Wei and Sandwell (2006); 1983 T.S. -the 1983 Geological Time Scale (Palmer, 1983). Maule et al. (2005) is the estimate of Antarctica's heat flow from Maule et al. (2005). Preferred total = 13.77 + 23.14 + 7.84 + 0.92 + 1.00 = 46.67 TW.   Fig. 12. A plot of the average heat flow for each geology category as a function of the square root of its mean age. There is a trend in the data, of decreasing mean heat flow with increasing age. If one excludes the very oldest datum (the Archean -which stretches over a very long geologic age) and the very youngest data (these categories have very few measurements), there is a strong correlation. While the correlation of heat flow with geology could apply without a correlation between mean heat flow and age; the fact that such a correlation exists does lend some support that this fundamental assumption might have some merit. The correlation is strong enough (R 2 = 0.75) that it certainly seems better to take it into account than ignoring it.

The robustness of the fundamental correlation assumed
The fundamental correlation assumed by Pollack et al. (1993) (PHJ93) is that regions of similar geology have, on average, similar values for heat flow. From Table 3, one can see that the standard deviation on the mean values for various geology categories is very high, suggestive that this assumption is not useful. However, this is misleading. In Fig. 12, we plot the mean heat flow of different geology categories as a function of the square root of age, for all categories with greater than 50 measurements (the age for the Precambrian is hard to set; most of the measurements are likely to be in the Late Proterozoic, but we have no way of knowing the exact age distribution. As a result, the Precambrian category is excluded from our plot). Polyak and Smirnov (1968) were among the first to suggest such a relationship, while Morgan (1985) pointed out that the relationship is weak. This weak relationship can be easily understood; although there is evidence of some decrease in heat production with age, the range of heat production within each age group is greater than the differences between age groups Mareschal, 2003, 2007). Nevertheless we find that there is a clear correlation in the data, with the average heat flow decreasing with increasing age (R 2 = 0.75). This implies that there is some power in the correlation and that using such a methodology as an empirical estimator for un-sampled regions gives a better estimate (with slightly lower error) than a straightforward global average (as was done by  (JLM07)). We remind readers that our assumption is that there is a correlation on average between geological category (not necessarily age) and heat flow. It should also be noted that our final estimate and that of JLM07 are similar. The relationship between average age and heat flow for the continents with our data (excluding the Precambrian) is: where q represents heat flow (mW m −2 ) and t is the age (millions of years). We note that the thermo-tectonic age might lead to a further improved correlation. However, at the time of writing, such information was not available for the digital geologies used here. Using this fundamental correlation, we find that the standard error on each geology category is low, since our data-set contains a large number of individual measurements. As a consequence, the contribution of this to the final error is very small (note that the errors are area weighted, as is the global estimate). However, since the correlation (between geology and heat flow) is not perfect, the error in areas with very few measurements could potentially be higher. As a result, it is likely that our estimate for this component of the error is unrealistically low. Nonetheless, once one realises that: (i) the spread of continental and old oceanic heat flow values have a limited range; and (ii) the uncovered area is not that great (at least at a 5 degree sampling), it cannot be very large. The only way to confirm our prediction is to continue measuring heat flow, especially in currently un-sampled areas.

Measurement techniques
Heat flow readings are most robust when taken in deep boreholes. However, such boreholes are rare and their construction is expensive. Since the surface heat flux is small compared to the solar heat flux and the resulting advected heat by fluids through the near surface is large, it is always a challenge to obtain accurate estimates. The issues related to the collection of heat flow data are discussed in detail by Beardsmore and Cull (2001) for example, while Slagstad et al. (2009) illustrate a recent example of a leading edge study that undertakes correction of heat flow measurements based upon climate and topography (note that such modern corrections are likely very rare in the data-set used here). While obtaining more accurate heat flow measurements is essential in constraining the global surface heat flux, we must accept that this will be a slow and expensive process. In the meantime, as much as possible must be extracted from the measurements made to date.

Areal coverage of heat flow measurements
If we assume that the young ocean is well covered, since it is estimated from a model independent of our heat flow data, then at a 5 degree equal area grid sampling we find nearly 84% coverage; at 1 degree equal area spacing we have 53% coverage. The coverage is not very much more than PHJ93 (who had 62% coverage on a 5×5 degree grid), even though our data-set contains around 55% more heat flow measurements. This demonstrates that recent measurements added to the data-set were made in similar geographical locations to those in the data-set of PHJ93. Further measurements must be made in un-sampled regions to improve the reliability of global heat flow estimates.

Various analysis choices
We have undertaken various analyses using different methods and data-sets (see Tables 4, 5 and 6). Such a range of analyses allows us to evaluate that: (i) the effect of including alternative geology is between 0.1 and 0.2 TW (see Table 5, 5th versus 7th row); (ii) not including the 1 degree grid increases clustering and leads to a decreased heat flow of ∼0.15 TW (Table 5, last two rows); (iii) increasing the threshold of points required to include a geology, as discussed in Sect. 3, only changes the global heat flow slightly, but reduces that (small) component of the error; and (iv) the geological time scales used give heat flow estimates that differ by 0.2 TW (see Table 6, 2nd column, 23.1-22.9 TW), which is very small (we utilised the 2004 Geological Time Scale (Gradstein et al., 2005) for our preferred estimate. The 1983 Geologic Time Scale (Palmer, 1983) was used for all other estimates). These analyses demonstrate that the final result is not sensitive to the details of these choices. In contrast, the various analyses undertaken show that the young ocean estimate is a major contributor to the final error.

Hydrothermal correction in oceanic crust
It is well known that measurements of heat flow on young oceanic crust grossly underestimate the actual heat flow (Davis and Elderfield, 2004). This is due to relatively shallow hydrothermal circulation. Theoretically, heat flow is expected to decline as the inverse square root of oceanic crustal age (Eq. 1). Indeed, this trend is observed in heat flow data over old oceanic crust. Therefore, in obtaining a complete estimate of Earth's surface heat flux, one must correct for hydrothermal circulation in young oceanic crust, using Eq. (1). A value for the constant C must therefore be selected and one must know the distribution of ocean floor area as a function of age. We note also that the theoretical expression cannot apply at the very youngest age. We discuss these aspects in turn.

The C constant
As noted previously, estimates of the constant C vary; 430, 490 and 510 mW m −2 Myr 0.5 , for PS, JLM07 and SS, respectively (note that these values are not strictly equivalent; JLM07 avoided hot-spots when calibrating data, but SS did not -in Table 5 we therefore do not include the hot-spot correction for the SS and PS cases. We have investigated the influence of this constant on the total young ocean estimate by using all three values; the difference, bounded by the PS and SS values, is 3.8 TW. This is significant compared to other sources of uncertainty. Our preference for C is the value of JLM07 (490 mW m −2 Myr 0.5 ), although, as discussed in Sect. 2.4, we specify a larger uncertainty (30 vs. 20 mW m −2 Myr 0.5 ). This increased uncertainty leads to a range of 2.6 TW, which is large, but less than the PS-SS range quoted above.

Area/age distribution
The estimate for heat flow in young oceanic crust is calculated by multiplying the average heat flow for a certain age range by the total area for that age. As a consequence, variations in the estimate of area for different age ranges could be significant. Table 1 lists the area of oceanic crust for CCGM geology. While using the listed area is the best approach, in this section we consider an alternative to get a sense of how significant this effect might be. JLM07, following Sclater et al. (1980) suggest that where A is the area of oceanic crust (km 2 ), τ is the age of the oceanic crust (years), and C A and τ m are constants. In Fig. 13, we plot cumulative area as a function of age, using the 2004 Geology Timescale, finding that it is well fit by the above expression (R 2 = 0.98). We find C A = 3.7 km 2 yr −1 while τ m = 150 Myr. This compares to estimates of C A = 3.45 km 2 yr −1 and τ m = 180 Myr (Sclater et al., 1980) and C A = 3.34 km 2 yr −1 (JLM07). If we use the area given by Eq.
(3), we get an estimate of 25.0 TW for the young oceans (for all ocean floor <66.5 Ma), compared to 22.3 TW from a calculation using the actual areas. Using the constants from JLM07 in Eq. (3), in an estimate of heat flux out to age of 66.5 Ma, yields a value of 23.4 TW. Therefore this source of error is ∼1-2 TW at most, with a smaller error likely with using actual area as done in this work.  Fig. 13. Cumulative area of oceanic crust (m 2 ) as a function of age (Ma). A least squares fit through the data for Eq.
(3) is shown as a dark curve. We note that the data are reasonably well fit by such an expression, as pointed out originally by Sclater et al. (1980).

Correction in very young oceanic crust
The half space formulation has a singularity (infinite value) at the ridge axis (see Eq. 1). This is an integrable singularity so it causes no numerical problems. Davies et al. (2008) modelled the surface heat flow across a spreading ridge in a twodimensional adaptive finite-element model. The model used a prescribed kinematic spreading upper surface to mimic the half space assumptions but avoided the unrealistic boundary condition at the ridge and had the resolution to numerically resolve the sharp variation. For a medium spreading rate of 5 cm yr −1 , Davies et al. (2008) found that the heat flow curve deviates from the half space curve at around 0.15 Myr and with an asymptotic value of just over 1 Wm −2 . If these conditions were the weighted average conditions for all spreading ridges, the half-space model would overestimate the heat flow by around 0.1 TW. The overestimate would be locally greater for slower spreading ridges and less for faster spreading ridges. While this effect is not negligible, it is less than the estimated error and, hence, we have not considered it in our final error estimate.

Hydrothermal circulation on continents
Since hydrothermal circulation is critical to estimates of oceanic heat flux, it is sensible to enquire what role it might play on continents. JLM07 suggest that the contribution is likely to be small, given that the estimate for the entire Yellowstone system is around 5 GW, with similar values predicted for the East African Rift. As stated by JLM07, it would require 200 "Yellowstones" to increase the continental heat loss by 1 TW. Our current understanding suggests that the contribution from this component is less than the estimated error.

Glacier category
The error arising from the Glacier category is potentially very large. As discussed earlier, our preference is a value of 65 mW m −2 (Maule et al., 2005). The error bounds of the authors (∼24 mW m −2 ) leads to an uncertainty of ∼0.3 TW (0.9 × 24/65), which is similar to our assumed error. We note that JLM07 estimate the error due to poor sampling to be around 1 TW; effectively this combines the error from: (i) the poor sampling discussed earlier; and (ii) the error arising from the Glacier category. Our weighted combined estimate for these sources of error is 0.8 TW.

Combining statistical errors
Our final estimate of the total area weighted error is 2 TW. This is slightly higher than the estimate of PHJ93 since we have assumed more uncertainty arising from the young ocean estimate (and not ignored the uncertainty arising from poorly sampled categories -in this case arising from the "Glacier" category). It is 1 TW less than the error estimate of JLM07, even though our estimates of the error for individual components are similar. We estimate a contribution of: 1.3 TW (oceanic correction); 0.3 TW (Glacier category, which is related to poor continental sampling in JLM07); 0.5 TW (rest of oceans); 0.3 TW (rest of continents -CCGM); 0.4 TW (rest of continents -GG); and 0.3 TW (hot-spots). Assuming that these error terms are independent and uncorrelated, the combined uncertainty is only 2 TW (not 4 TW if one incorrectly added the error contributions naively (e.g. for A = B + C, the errors should be combined A 2 = B 2 + C 2 , not A = B + C)). The biggest difference between our error and that of JLM07 is that we account for only 0.33 TW, compared to 1 TW, for the hot-spot correction.

Significance of result to thermal budget and thermal evolution models
Of the heat emerging at Earth's surface (estimated at 47 ± 2 TW herein), it is believed that 6-8 TW is generated within the crust (e.g. Rudnick and Fountain, 1995;Jaupart and Mareschal, 2003). It is argued that 5-13 TW originates within the core (e.g. Buffett, 2003;Lay et al., 2006), although this range remains highly uncertain. The remainder (24-38 TW) must be provided by heat generation within the mantle and by the secular cooling of the planet (plus minor contributions from other sources). Various models for the bulk silicate Earth (which includes the continental crust) lead to a total present-day heat production of ∼20 TW, with an uncertainty of ∼15% (e.g. McDonough and Sun, 1995;Palme and O'Neill, 2003;Lyubetskaya and Korenaga, 2007b). Removing the contributions from continental crust leaves a mantle heat production of 13±4 TW (JLM07). This would leave the balance (7-29 TW) attributable to secular cooling. Recent studies suggest that such rapid cooling is highly unlikely.
Indeed, JLM07 argue that the initial temperature of the solid Earth was only ∼200 K higher than the present-day. These figures therefore reveal an energy imbalance relating heat emerging at Earth's surface and heat generated within Earth's interior. There are, however, many hypotheses for how these diverging constraints can be satisfied, such as: (i) an increased CMB heat flux; and (ii) a delay between the generation of heat in Earth's interior and its arrival at the surface (see JLM07 for a discussion). Nonetheless, our revised estimate of Earth's surface heat flux only exacerbates this issue. Our result follows a recent trend of increasing global heat flow estimates, which makes the global "energy paradox" described above (Kellogg et al., 1999) more difficult to understand.

Conclusions
Our revised estimate of Earth's total surface heat flow is 47 ± 2 TW, which is larger than previous investigations. This estimate was derived from an improved heat flow data-set, with 38 347 heat flow measurements and the methodologies of Geographical Information Science. Given the sparse and inhomogeneous nature of heat flow measurements globally (poor sampling in Antarctica, Greenland, Africa, Canada, Australia, South America and parts of Asia), there remain uncertainties in our estimate. In addition, models for hydrothermal circulation in young oceanic domains exert a significant control on our final value. Our result follows a recent trend of increased estimates for Earth's surface heat flow, thus posing difficulties for simple interpretations of heat sources in the mantle. Nonetheless, our estimate will provide a concrete boundary condition for future investigations of Earth's thermal evolution.