Because global coverage of heat flow measurements is still poor in many areas, empirical estimators based on the geology, the thermotectonic age or the velocity structure of the upper mantle have often been used to affect an estimate to regions where such measurements are absent. On the basis of the assumption that heat flow is strongly related to its geodynamic environment, one may integrate multiple proxies derived from a large body of global geo- logical and geophysical data sets assembled during the past decades; these should help to better capture the variety of present-day settings. This idea is illustrated through two simple empirical methods: both of them are based on a set of examples, where heat flow measure- ments are associated with relevant terrestrial observables such as surface heat production, upper-mantle velocity structure, tectono-thermal age, on a 1◦ × 1◦ grid. To a given target point owning a number of observables, the methods associate a heat flow distribution rather than a deterministic value to account for intrinsic variability and uncertainty within a defined geodynamic environment. The 'best combination method' seeks the particular combination of observables that minimizes the dispersion of the heat flow distribution generated from the set of examples. The 'similarity method' attributes a weight to each example depending on its degree of similarity with the target point. The methods are transparent and are able to handle sets of observables that are not available over the whole Earth (e.g. heat production). The resulting trends of the mean heat flow deduced from the two methods do not differ strongly, but the similarity method shows a better accuracy in cross-validation tests. These tests suggest that the selected proxies have the potential to recover at least partly medium- to large-scale features of surface heat flow. The methods depict the main global trends of low heat flow in stable and ancient regions, and thermal high in active orogens and rift zones. Broad thermal anomalies are outlined in the Sahara and in the tectonically active eastern part of Antarctica. The similarity method estimates a continental heat loss of 13.6 ± 0.8 TW (2σ uncertainty), which is consistent with previous estimates.