Present-day heat flow model of Mars

Until the acquisition of in-situ measurements, the study of the present-day heat flow of Mars must rely on indirect methods, mainly based on the relation between the thermal state of the lithosphere and its mechanical strength, or on theoretical models of internal evolution. Here, we present a first-order global model for the present-day surface heat flow for Mars, based on the radiogenic heat production of the crust and mantle, on scaling of heat flow variations arising from crustal thickness and topography variations, and on the heat flow derived from the effective elastic thickness of the lithosphere beneath the North Polar Region. Our preferred model finds heat flows varying between 14 and 25 mW m−2, with an average value of 19 mW m−2. Similar results (although about ten percent higher) are obtained if we use heat flow based on the lithospheric strength of the South Polar Region. Moreover, expressing our results in terms of the Urey ratio (the ratio between total internal heat production and total heat loss through the surface), we estimate values close to 0.7–0.75, which indicates a moderate contribution of secular cooling to the heat flow of Mars (consistent with the low heat flow values deduced from lithosphere strength), unless heat-producing elements abundances for Mars are subchondritic.

In order to have an understanding of the internal structure and the composition of Mars, as well as the thermal state of its subsurface and lithosphere, accurate knowledge of the present-day pattern of heat flow of the planet is required. Indeed, the heat flow should vary across the surface of Mars (as well as in time), depending on the abundances and distribution of the crust and mantle heat-producing elements, and on the thermal state of the planet's interior and ability to dissipate heat. On Mars, such variations of heat flow can probably affect, among many other things, the internal geodynamics of the planet, the distribution of groundwater, ice and clathrates, and processes such as subsurface melting and fluid circulation. It is essential, therefore, in deriving global heat flow models to predict those variations throughout the planet geographically. Furthermore, the determination of average planetary heat flow provides restrictions on radioactive element content, which in turn, is important to constrain cosmochemical models and the total radioactive heat production in a planet.
Moreover, the Urey ratio (Ur), the ratio between the total radioactive heat production and the total surface heat loss, provides knowledge about the thermal state of a planet and whether the planet is geologically active nowadays. Thus, Ur constitutes an excellent means to visualize the present-day internal heat budget of the red planet, because it describes possible interior cooling or warming. Therefore, if we compared the available information on the Urey ratio throughout the history of Mars, we should be able to constraint more accurately its thermal history 1 .
Until the arrival of InSight (Interior Exploration using Seismic Investigations, Geodesy and Heat Transport), a NASA Discovery Program mission (http://solarsystem.nasa.gov/insight/overview.cfm), which includes the HP 3 (Heat Flow and Physical Properties Probe) instrument, direct measurements of heat flow on Mars are non-existent. Foreseeably in 2018, InSight will land on the surface of Mars with the aim of determining the heat flow at the landing site in Elysium Planitia (the favored site is centered at 136°E 4°N), from the regolith thermal conductivity and the subsurface thermal gradient. In order to do this, the HP 3 instrument, a probe ("the mole") will penetrate up to 5 m into the martian regolith and work for one martian year 2 .
Until these in-situ measurements are acquired, the information currently available on past and present Martian heat flows must be deduced from indirect methods. A commonly used indirect method is based on the relation between the thermal state of lithospheric rocks and their mechanical strength, usually analyzed through the effective elastic thickness of the lithosphere (e.g., refs 3-8) or from the depth to the brittle-ductile transition beneath large thrust faults (e.g., refs 9-12). The so-obtained heat flows are valid for the time when the lithosphere was loaded or faulted, and therefore, the information deduced from regions deformed at different ages provides knowledge on the thermal evolution of Mars 8 .
Although numerous paleo-heat flow estimates have been derived using this approach, information on the present-day heat flow of this planet is limited by the difficulty in finding young structures on the surface of Mars with the polar regions currently providing the only estimates. The very limited (or non-existent) flexure caused by the north polar cap load indicates a very thick (> 300 km) effective elastic thickness of the lithosphere (T e ) on the North Polar Region (hereafter NPR) 13 . Similarly, admittance modelling of the South Polar Region (hereafter SPR) found a best fit for T e of 161 km 14 . These effective elastic thicknesses can be used to estimate the present-day heat flow in those regions (refs 11,13,15; see also below), because polar caps (and the associated loading) are a recent phenomenon 16 .
On the other hand, it is recognized that large crustal thickness variations should have an expression on the surface heat flow pattern, due to the contribution of radiogenic heat production originated by the heat-producing elements (HPEs) placed in the crust [17][18][19][20] . The significant homogeneity of elemental abundances measured by Mars Odyssey GRS suggests that the martian crust is much less geochemically variable than the Earth's crust 21 , which is consistent with a strong mixing by impact cratering (e.g. ref. 21), and the absence of large-scale crustal recycling, at least since the early part of the history of Mars (e.g., ref. 22), which could induce compositional differences between provinces. The influence of crustal thickness and heat production on the pattern of surface heat flow can therefore be reliably addressed, at least in a general sense.
In this work we use the currently available information on heat flow deduced from the present-day effective elastic thickness of the lithosphere, along with crustal differences in vertically integrated heat production arising from crustal thickness variations, in order to produce a preliminary global heat flow model for the present-day Mars. Previous global models were based on thermal history models 17,20 , whereas the present model is more empirically constrained by heat flow estimates based on lithospheric strength. Our aim is to gain understanding on the present heat flow patterns, which may be compared with our previous knowledge on the heat flow history of Mars, with thermal evolution and geochemical models, and, in the future, with the data obtained by InSight. Thus, this model should help us to advance another step in the study of the thermal state and evolution of the terrestrial planets.

Results
Crustal heat flows. If lateral heat transfer is neglected, the surface heat flow of Mars may be considered the sum of the heat generated in the crust and the heat flow from the mantle. In turn, the heat flow from the mantle is a consequence of the heat radioactively produced in the mantle lithosphere (or, more generally, in the stagnant lid), and of the heat coming-up from the deep interior (through mantle convection). The component of the heat flow arising from crustal radioactive heat sources is the sum of the respective contributions from all the HPEs. Global maps of elemental concentrations indicate that there is variation across the surface of Mars 21,23,24 , which is related to rock type or to extent of aqueous alteration 21 .
Previously, Hahn et al. 18 , obtained a crustal heat flow map based on local surface abundances of K and Th (and estimated abundances of U), considered as representative of the entire local crust, for latitudes roughly between 60°N and 60°S, because influence of volatile in higher latitudes obscures the GRS signal. However, the variations in heat production across the surface are relatively little; also, as above indicated, there are no reliable GRS data for latitudes above 60° and below − 60°, where HPEs measurements are heavily affected by the presence of water molecules near the surface 18,24 . Thus, we here prefer to consider an average value for the crustal heat production in our model, which additionally permit modeling crustal heat flows at high latitudes (see Methods Section).
We first derived a crustal heat flow map from crustal heat production and a crustal thickness model (for the description of the crustal thickness model see Methods Section). The crustal thickness model was derived from topography and gravity following the procedure of potential theory 25 by assuming an average thickness of 50 km (Fig. 1). This mean crustal thickness is slightly higher than previous models [26][27][28] , but in line with geophysical and geochemical evidence 5,12,29 , and consistent with the average thickness usually used in the modeling of lithospheric strength (e.g., refs 5 and 8). In agreement with previous geophysical modeling, our crustal thickness model assumes densities of 2900 and 3500 kg m −3 , respectively, for the crust and lithospheric mantle (see Methods Section). A crustal density of 2900 kg m −3 can be consistent with a basaltic crust somewhat fractured, or consistent with a basaltic crust including a felsic component. Indeed, some recent works suggest that the Martian crust could contain a substantial amount of felsic rocks [30][31][32][33] , although other works consider as non-solids the orbital-based evidences for a felsic component (see ref. 34). On the other hand, Baratoux et al. 32 proposed a high density (> 3200 kg m −3 ) for the martian basalts, although a crust so dense is incompatible with geophysical studies 32 , maybe favoring the presence of a felsic component. In any case, note that the HPEs abundances used in our model were derived from Mars Odyssey GRS observations, not assumed as a function of composition.
The crustal component of heat flow varies between 1.3 and 13.5 mW m −2 , with an average value of 7.0 mW m −2 (Fig. 2). The regional pattern shown by our model is similar to that observed in Hahn et al. 18 , albeit the absolute values are higher due to the higher crustal thicknesses of our crustal model. There is a strong correlation between the crustal thickness and crustal heat flow obtained, which implies that regional variations of crustal heat flow are strongly dependent on crustal thickness, which produces higher (lower) heat flow in zones with thickened (thinned) crust.
Present-day heat flow model. We can scale heat flow differences across the martian surface arising from crustal thickness and topographic differences across the planet. We complete the model by assuming a constant value for the heat production in the mantle lithosphere (see Methods Section) and a constant (sub-lithospheric) heat flow from the deep interior, which would include the contribution from the convective mantle; although mantle convection should produce some regional variations of convective heat flow, we assume a constant sub-lithospheric heat flow by necessity, because we cannot specify the magnitude of these variations.
For anchoring our nominal model, we use a heat flow of 17 mW m −2 derived (see Methods Section) for an effective elastic thickness of the lithosphere of 300 km at the NPR, which was derived from the limited response to loading in this region 13 ; lithosphere strength at the NPR currently constitutes the more reliable heat flow indicator for the present-day Mars 15 .
Although we do not adopt any specific compositional assumption for our crustal model or for crustal HPEs abundances, the composition of the crust might influence their thermal conductivity, which in turn influences the calculation of heat flows at the NPR (see ref. 8 and Methods section). The thermal conductivity used for our nominal model is 2 W m −1 K −1 , a value appropriate for basalts (whose thermal conductivity is usually between 1.5 and 2.0 W m −1 K −1 35 ). For the felsic terrestrial continental upper crust a thermal conductivity of 2.5 W m −1 K −1 is frequently used, although a recent revision 36 has found that the thermal conductivity of felsic rocks in the terrestrial crust usually is not significantly higher than 2 W m −1 K −1 . Because a thermal conductivity of 2 W m −1 K −1 is an upper limit for basaltic rocks, the use of this value in our nominal model is also roughly consistent with a basaltic crust including a felsic component.
The obtained nominal heat flow model is shown in Fig. 3. The average present-day surface heat flow is 19 mW m −2 , and it varies between 14 and 25 mW m −2 . The lowest values are found in the northern lowlands (e.g., Acidalia or Utopia Planitiae) and the giant impact basins (e.g., Hellas, Argyre and Isidis Planitiae) where the  Present-day Urey ratio of Mars. The Urey ratio characterizes the degree of internal cooling (or heating) of a planetary body, and it can be defined as the ratio between the total radioactive heat production and the total heat loss through the surface for a given time. If Ur > 1 the planet interior, as an average, is heating-up; conversely, if Ur < 1 the planet interior is cooling-down. For Mars, if we take into account a present-day radioactive heat production equivalent to a surface heat flow of 14.3 mW m −2 , according to the compositional model of Wänke and Dreibus 37 , and the average surface heat loss deduced from our model (19 mW m −2 ), then an Urey ratio of around 0.75 is obtained for the present-day Mars. This value is higher than the Ur (0.594 ± 0.024) proposed from some recent thermal evolution models 38 , but consistent with a more limited cooling deduced from lithospheric strength analysis 1 .

Discussion
The terrestrial planets have evolved from the same primitive elements and have a similar inner structure (with differentiated core, mantle and crust), although, in addition to distance to the Sun and size, they show chemical and geological differences that make them unique within our planetary system. Major events occurred in the history of a planet (large impacts, global dynamics and environmental changes, etc.), have left signatures which define the current appearance of, and differences between, those planets. Thus, the present-day surface heat flow is very important for understanding the thermal and geologic evolution of a planetary body. Therefore, in order to improve our knowledge of the evolution of Mars, including the prospect for crustal water circulation and availability, we need to understand more fully the present thermal state of Mars, as well as the regional pattern of surface heat loss.
In absence of direct heat flow measurements for Mars, the loading of the lithosphere of the polar regions by the northern and southern polar caps (which are recent features; refs 13,16) has allowed the estimation of their respective effective elastic thickness 13,14 , and hence to derive estimates of the present-day heat flow 8,15 . In this work, we have obtained a global heat flow model (the F NPR -based model) for present-day Mars from (i) scaling heat flow variations related to variations in crustal thickness and topography, (ii) assuming constant mantle heat flow, and (iii) using the heat flow obtained for the NPR from the effective elastic thickness as an anchoring value. Also, for comparative and control purposes, we derive two alternative models (see Methods Section.): a model (the F SPR -based model) similar to the nominal model, but in this case constructed using as anchoring value the heat flow obtained for the SPR from the effective elastic thickness; and a model (the F NPRfelsic -based model) using as anchoring value the heat flow calculated for the NPR from a crustal thermal conductivity of 2.5 W m −1 K −1 (a very generous upper limit thermal conductivity if the crust contains a felsic component).
Our nominal (F NPR -based) model (Fig. 3)  We also obtain regional heat flow variations on Mars of, at most, ± 6 mW m −2 with respect to the mean value. In order to verify that the results of both models are Geographical heat flow variations in our models are related to variations in the HPE content of the crust and mantle and, especially, to crustal thickness variations. There are other processes that could modify the surface heat flow pattern, as hydrothermal circulation 39 or the presence of mantle plumes below the lithosphere (e.g., refs 17 and 40). There is some evidence for recent water-related geologic activity, as for example water channels flowing from Cerberus Fossae, which could be related to magma intrusion 41 , although it is not clear how this activity would affect to the background heat flow. On the other hand, the existence of mantle plumes below volcanic regions should result in a locally increased heat flow. The analysis of gravity and topography does not find evidence for substantial dynamical support for Tharsis, the largest igneous complex on Mars, and therefore for large mantle plumes active beneath this region 42,43 . Thus, whereas mantle heat flow variations due to mantle plumes could exist, we have ignored their possible contributions in our present model, although future evidence could require local adjustments.
Our heat flow results can also be interpreted in terms of the Urey ratio. The total radioactive heat production obtained from the bulk compositional model of Wänke and Dreibus 37 (with chondritic HPEs abundances) can be converted to an equivalent heat flow, giving 14.3 mW m −2 for the present-day. Thus, for the average heat flow deduced from both the F SPR -and F NPR -based models, we obtain Ur of, respectively, 0.68 and 0.75; and Ur = 0.715 for an average heat flow of 20 mW m −2 , the value intermediate between both models. Similarly, the Urey ratio calculated for F NPRfelsic -based model, Ur = 0.73, is consistent with our other results. These values indicate a relatively small difference between the total radioactive heat production and heat loss through the surface, and therefore a moderate contribution from secular cooling to the heat flow of Mars.
Our Ur results can be compared with predictions calculated by convective history models (Fig. 4). Our Ur values are clearly higher than those obtained from the 2D and 3D convective history models of Plesa and Breuer 44 and Plesa et al. 38 , which usually vary between 0.52 and 0.62 as a consequence of higher predicted heat flows; these models would suggest therefore a substantial amount of secular interior cooling. Similarly, earlier 1D parameterized convective models predicted present-day heat flows around 23-25 mW m −2 (e.g., refs 45 and 46), equivalent to Ur ≈ 0.6 for the composition model of Wänke and Dreibus 37 ; more recent and sophisticated (including melting effects on mantle viscosity and hydration) 1D convective models obtained lower heat flows, in the range between 19 and 21 mW m −2 (refs 17, 47 and 48), which correspond to a Ur values between 0.68 and 0.77 (see Fig. 4), similar to those obtained by the present work.
On the other hand, Urey ratios of ~0.6-0.7 have been proposed from the melting conditions, deduced from elemental abundances, of magmas covering some Amazonian-aged (< 3 Ga) terrains 49 . However, the so-calculated mantle potential temperatures informs us on volcanic (non-average) regions, and Ruiz et al. 8 show that the Ur values obtained by these authors are therefore lower limits. However, volcanic regions could be associated with higher than average heat flows, making any Ur value based on those regions a lower limit; taking this into account, these results would be consistent with our own results.  48 were calculated from the present-day heat flow obtained by these models and the total heat production from the composition model of Wänke and Dreibus 37 . Ur = 1 separates the fields where the planet's interior is cooling-down (Ur < 1) or heating-up (Ur > 1). A limited influence of secular cooling through post-Noachian times is consistent with evidence for a liquid core existing at present 50 , recent abundant volcanism 51 and limited global contraction 52 ; limited secular cooling might also be consistent with high-temperature shergottite 53,54 , although this class of meteorites might have originated in mantle-plume settings unrepresentative of background mantle 55 . (For a review of the evidence for limited secular cooling see ref. 1). Alternatively, a low average heat flow could be explained by a sub-chondritic heat-producing elements abundances 13 ; in this case, our heat flow model would give proportionally lower Urey ratios for the current Mars. Figure 5 shows the Urey ratio as a function of HPE Mars /HPE Chondritic (where HPE Chondritic represents the standard chondritic model for Mars; Wänke and Dreibus 37 ). The blue and red lines represent, respectively, Ur values obtained from the average heat flows of the F NPR -and F SPR -based models. For example, if HPE Mars / HPE Chondritic = 0.9, the total radioactive heat production is equivalent to 12.8 mW m −2 , and the corresponding value of Ur is 0.68 and 0.61 for, respectively, the F NPR -and F SPR -based models. Thus, for obtaining Ur values similar to values found by Plesa et al. 38 (between 0.5 and 0.6), requires an HPE Mars /HPE Chondritic ratio between 0.75 and 0.85, and a total radioactive heat production equivalent to ~10-12 mW m −2 . Moreover, an Earth-like Ur < 0.5 36 would require an excessively low HPE Mars /HPE Chondritic ratio around ~0.7 or less, and we can conclude therefore that heat loss is less efficient for Mars than Earth.
Our results on the present-day thermal state of Mars, and the corresponding construction of heat flow models and maps, are preliminary and a first step to characterizing the heat flow structure of the planet, and provide an understanding of the thermal history and dynamics of Mars. These results will also be helpful in evaluating specific landing zones for future missions to Mars; for example, contributing to the development of further studies on the distribution of groundwater or clathrates 56,57 , subsurface geological/biological production and transport of methane 58 , ice melting processes, localizing areas with higher thermal gradient, or other potential geological, geophysical and geochemical issues on the red planet. Also, our model can be improved with the constraints provided by additional observations, measurements and theoretical advances. The most immediate example is the InSight mission, which will provide heat flow measurements at a single location that according to our models are predicted to be 18 ± 1 mW m −2 at the proposed landing site at Elysium Planitia. The in-situ heat flow data can be implemented easily in this methodology, and will provide further anchoring for our model.

Methods
Crustal and mantle heat production. K and Th abundances on the surface have been estimated from measurements by the Gamma Ray Spectrometer (GRS) instrument on board the 2001 Mars Odyssey spacecraft 21,24 , and then corrected for volatile presence 18 , whereas U abundances are determined by assuming a Th/U ratio of 3.8 (e.g., ref. 18). As we use surface average HPEs abundances as representative for the whole crust (3652, 0.69 and 0.18 ppm for, respectively, K, Th and U), the corresponding average heat production of Mars is currently 4.86 × 10 −11 W kg −1 (see also ref. 18). Concurrently, with a crustal density of 2900 kg m −3 (see below) this is equivalent to 0.1408 μ W m −3 , a contribution of ≈ 0.14 mW m −2 to the surface heat flow per each kilometer of crustal column. We assume HPEs homogeneously distributed in the crust, because the significant homogeneity of elemental abundances measure by Mars Odyssey GRS suggests that the martian crust is much less geochemically varied than the Earth's crust 21 , which is consistent with a strong mixing by impact cratering 21 , and with the absence of large-scale crustal recycling, at least since the early part of the history of Mars 22 . Although there are some evidences suggesting that the crust of Mars could be stratified, with heat-producing elements concentrated in a layer thinner than the whole crust 12,59 , currently these evidences lack enough detail, and do not justify specifying whichever given crustal stratification structure in our modeling. Also, HPEs abundances in the mantle lithosphere are poorly constrained. Here we use HPEs mantle lithosphere abundances 0.1 times the average value for the martian crust (see ref.  the thickness of the Martian crust, we assume (1) that the observed gravitational anomalies arises only from relief along the surface and crust-mantle interface, and (2) constant crustal and mantle densities to overcome the non-uniqueness associated with potential modeling. We use the spherical harmonic model MarsTopo2600 of the shape of Mars from Wieczorek 60 , and the gravitational potential model JGMRO_110 C from Konopliv et al. 61 . Topography data are useful only up to the resolution of gravity data, so both gravity and topography coefficients are truncated beyond degree and order 90 in our analysis given the dramatic decrease in spectral correlation that is observed between gravity and topography beyond this degree (e.g., ref. 32). It is further required either to assume a mean crustal thickness or to anchor the inverted crustal thickness to a given value at a specific location 60 . Here we assume a mean crustal thickness of 50 km (that satisfy the condition that the inverted crustal thickness is not negative anywhere on the planet). For crustal density we use 2900 kg m −3 for several reasons: this value is the usually used by works calculating crustal thickness from gravity and topography 60,62 and the effective elastic thickness of the lithosphere (e.g., refs 4, 5, 13 and 26); a crustal density of 2900 ± 50 kg m −3 produces the best fit for the admittance/coherence signal of most of regions of Mars 4,5 ; an upper limit of 3020 ± 70 kg m −3 has been calculated for the crustal density at the highlands of Mars 63 from Geoid-to-topography ratios and Bouguer inversion. For the mantle lithosphere density, we assume 3500 kg m −3 , a value widely used for Mars (e.g., refs 5, 8, 26 and 27). Under these assumptions, we first calculate the Bouguer gravity anomaly from surface topography and free air anomaly, and then calculate by downward continuation the shape of the crust-mantle interface necessary to minimize the difference between the observed and predicted Bouguer anomalies (for reviews see refs 25 and 60). In order to mitigate errors in downward continuing the Bouguer anomaly, we applied a minimum amplitude filter (see ref. 25) for the Moho relief at degree l = 50. Finally, we obtain the crustal thickness by subtracting the relief on the Moho from surface topography.
Scaling of present-day surface heat flow. We scale heat flow differences across the martian surface due to crustal thickness and topographic differences in the planet using the expression where F local is the local heat flow for a given location on the surface of Mars, F ref is a reference value used for anchoring the model, H c and H m are, respectively, the crustal and mantle heat production by volume unit, Δ t is the difference between the local elevation and the elevation at the reference location, and Δ b is the difference between the local crustal thickness and the crustal thickness in the reference location. The term H c Δt corresponds to the crustal contribution to the surface heat flow due to topographic differences across the surface; the term (H c − H m )(Δb − Δt) corresponds to the effect from Moho topography and differences in heat production between crust and mantle. Surface topography and crustal thicknesses are derived, respectively, from Mars Orbiter Laser Altimeter (MOLA) topography 64 , and from our crustal thickness model described above. For reference location we take the North Polar Region (NPR), and we average crustal thickness and topography values in the region including 180°W to 180°E and 75°N to 90°N; in this case, for elevation determination it is necessary to subtract the average thickness of the north polar ice cap (2 km; refs 13, 14 and 65). The effective elastic thickness of the lithosphere (T e ), which is a measure of the total strength of the lithosphere, can be converted to heat flow following the equivalent strength envelope procedure described by McNutt 66 . This methodology is based on the condition that the bending moment of the mechanical lithosphere must be equal to the bending moment of the equivalent elastic layer of thickness T e . The link between the mechanical structure and heat flow comes from the dependence of the ductile strength on temperature. Here we use the implementation of this methodology by Ruiz et al. 8,15 , which takes into account the strength contributions from the crust and the lithospheric mantle, and the effect of the polar cap above (but not included in) the lithosphere, in order to calculate an upper limit heat flow for the NPR. We calculate the temperature profile following the procedure described by Ruiz et al. 8,15 , and using crustal heat-producing elements (HPE) abundances based on Mars Odyssey GRS measurements (see above). We perform this calculation using a Cartesian geometry, because it is sufficient for deriving upper limit heat flows. Indeed, for a spherical geometry (which could be considered as appropriate due to the small radius of Mars): i) the temperature as a function of depth is somewhat increased, for a given heat flow, with respect to the Cartesian geometry, which in turn makes the temperature profile accordingly somewhat hotter; ii) increased temperatures in the lithosphere reduce the lithosphere strength; iii) because a fixed T e implies necessarily a fixed lithosphere strength, in order to fit a given T e derived from the observations, we must have a lower surface heat flow than in the Cartesian geometry. For F ref , we use a value of 17 mW m −2 , which correspond to a heat flow calculated following this procedure, and taking for the NPR an effective elastic thickness of the lithosphere (T e ) of 300 km (ref. 13), a crustal thickness of 33 km (after our crustal thickness model), an average elevation of − 6 km, a thermal conductivity of 2 W m −1 K −1 for the crust, and a temperature-dependent thermal conductivity for the lithospheric mantle 8 . For comparative purposes, we have carried out in the same way an alternative heat flow model, based on scaling a F ref of 23.5 mW m −2 deduced for the South Polar Region (180°W to 180°E and 80°S to 90°S) from T e = 161 km 14 , and average crustal thickness and elevation of, respectively, 71 km and 1 km. We base our nominal model on the lithosphere strength in the NPR because T e is better constrained than in the SPR (see ref. 15). For a complete explanation and discussion of the values used in the calculations see Ruiz et al. 8,15 . We also have performed a surface heat flow model from a F ref of 17.5 mW m −2 at the NPR, based on an upper limit crustal thermal conductivity of 2.5 W m −1 K −1 .

Urey ratio calculations.
We derive present-day Urey ratios by calculating the ratio between the total radioactive heat production obtained from the compositional model of Mars of Wänke and Dreibus 37 and standard decay constants, and the average of present-day surface heat flow derived from our model. Thus, Urey ratio is given by Ur = H M M mars /(4π R mars 2 F), where H M is the mean heat production rate by mass unit, M mars and R mars are, respectively, the mass and the mean radius of Mars, and F is the average surface heat flows.