Groundwater production from geothermal heating on early Mars and implication for early martian habitability

See allHide authors and affiliations

Science Advances  02 Dec 2020:
Vol. 6, no. 49, eabb1669
DOI: 10.1126/sciadv.abb1669


In explaining extensive evidence for past liquid water, the debate on whether Mars was primarily warm and wet or cold and arid 4 billion years (Ga) ago has continued for decades. The Sun’s luminosity was ~30% lower 4 Ga ago; thus, most martian climate models struggle to elevate the mean surface temperature past the melting point of water. Basal melting of ice sheets may help resolve that paradox. We modeled the thermophysical evolution of ice and estimate the geothermal heat flux required to produce meltwater on a cold, arid Mars. We then analyzed geophysical and geochemical data, showing that basal melting would have been feasible on Mars 4 Ga ago. If Mars were warm and wet 4 Ga ago, then the geothermal flux would have even sustained hydrothermal activity. Regardless of the actual nature of the ancient martian climate, the subsurface would have been the most habitable region on Mars.


The habitability of Mars over geologic time, inextricably linked to the availability of liquid water and energy, remains a key theme in planetary exploration. Much of Earth’s microbial biomass resides within its crust, where water is readily available. Substantial biological diversity exists throughout the huge volume of subsurface habitable environments, which may reach >5-km depth (1). The thermophilic, chemoautotrophic nature of the last universal common ancestor of modern life and the deep rooting within the phylogenetic tree of organisms from hydrothermal systems point to the importance of subsurface water-rock reactions in the potential development of early terrestrial life and maintaining a habitable subsurface environment. Therefore, the subsurface could have been the most viable habitat for ancient simple life forms on early Earth and possibly Mars.

Much morphological and compositional evidence on Mars suggests abundant surface water sustained by subsurface hydrology during the Noachian eon [4.1 to 3.7 billion years (Ga) ago]. Those include spectroscopic detections of clays and other hydrous minerals (2), abundant valley networks (VNs) in the mid to low latitudes (3), lakes (4), putative shorelines (5), and in situ observation of conglomerates (6). However, whether that points to surface precipitation or groundwater is hotly debated. The uncertainty arises because stellar evolution models indicate that the Sun’s luminosity was approximately 30% lower 4 Ga ago (7). Mars receives only 43% of the solar flux incident on Earth; thus, provided that the orbit of Mars has not changed notably in the past 4 Ga, the climate of early Mars should have been extremely cold. The resulting faint young Sun paradox, between climate models that struggle to elevate surface temperature above 273 K and geological evidence suggesting the presence of abundant liquid water during the Noachian, is an outstanding question in Mars science, with major implications for early martian climate, hydrology, and habitability.

Isotopic data (8) and the contemporary rate of gas loss to space (9) indicate that Mars had a substantially thicker atmosphere during the Noachian eon. However, even with a thick atmosphere where surface pressure exceeds the triple point of H2O, a long-term, warm climate is unattainable with a CO2- and H2O- rich atmosphere alone (10, 11). In the past few decades, various greenhouse gases and clouds have been proposed as solutions to the faint young Sun paradox. Those include transient warm conditions driven by atmospheric injection of SO2, H2S, CH4, and H2 into the martian atmosphere by impacts, volcanism, or seasonal effects [e.g., (12)]. However, most of the proposed climate models have recently been discounted upon closer inspection (12), exacerbating the problem of creating and sustaining liquid water on the surface. Consequently, the view of long-term warm and wet surface conditions on Mars is problematic, if not incompatible [e.g., (13)], with our present-day understanding of the early martian climate powered by a faint young Sun.

An alternative view posits that Mars was mainly frozen, with aquifer discharge or intermittent snow and ice melt as the main source of liquid water (Fig. 1) [e.g., (14)]. In that scenario, liquid water forms primarily through basal melting of thick ice deposits and is maintained at temperatures above the freezing point in the martian regolith by either high geothermal heat flow, heat provided by impact, or volcanism (14). Multiple lines of evidence provide support for groundwater on Mars. Hydrothermal mineral phases such as serpentine and phyllosilicates are detected within the materials excavated from the deep subsurface by large impact craters (15, 16). Clear evidence of groundwater diagenesis can be found at Meridiani Planum (17) and at Gale crater (18). The widespread detection of crustal Fe/Mg clays over the Noachian regions of Mars implies that aqueous alteration mostly occurred deep in the crust and implicates the existence of widespread low-grade metamorphism/diagenesis or a hydrothermal system [e.g., (16)].

Fig. 1 Schematics of two possible end-member aqueous environments on Mars during the Noachian.

(Left) Warm and wet view of early Mars, where liquid water can be stable on the surface. (Right) Cold, arid view of early Mars, where liquid water is mostly generated in the subsurface by intermittent melting of snow and ice. The curly red arrows show where hydrothermal reactions could take place.

Our understanding of aqueous and hydrothermal alteration on Mars has mainly relied on spectroscopic detection of altered phases (e.g., Fe/Mg phyllosilicates, prehnite, serpentine, and chlorite) in crustal rocks excavated from depth by large impact craters [e.g., (16)]. An alternative approach lies in first identifying regions on Mars with elevated geothermal heat flow where hydrothermal alteration could have likely occurred. Here, we use the latest geochemical and geophysical data of Mars to provide a first-order estimate of the Noachian surface heat flow, where, as in prior works (19), the current regional bulk chemistry to decimeter depths represents its ancient counterpart. Under the assumption that the Noachian Mars was primarily cold and arid [e.g., (11)], we first estimate the volume of liquid water that would form by basal melting of globally deposited ice sheets by our estimated Noachian geothermal heat flow, providing a possible solution to the faint young Sun paradox. We also consider the alternative of a Noachian Mars characterized by warm and wet conditions and the associated effects of geothermal heating on subsurface hydrothermal activity.


Feasibility of basal melting on Noachian Mars depends on the surface temperature, the thickness of the ice sheets, and the surface heat flow (fig. S1). We first investigate the thermophysical evolution of the ice by using the approach outlined in Materials and Methods. Specifically, we seek to ascertain the thickness of the ice and surface heat flow required for basal melting. Once these parameters are estimated, we assess whether surface heat flow necessary for basal melting would have been available in Noachian Mars. Modeling a CO2-dominated atmosphere of 1-bar pressure and a solar flux of 75% of the current value, (20) estimated the mean annual surface temperature for the Noachian Mars for various obliquities using three-dimensional climate models. As expected, the mean annual surface temperature is too low (maximum of ~230 K) for liquid H2O to be stable anywhere on the surface. The model also shows that elevated regions of Mars such as the southern hemisphere would have experienced substantial deposition of ice and snow (10, 11) (colloquially referred to as the icy-highland hypothesis) (Fig. 1). The thickness of snow and ice in the southern highlands would have depended on the availability of atmospherically circulated water during the Noachian. We conservatively model that the mean thickness of ice deposits on the southern highlands (above the predicted ice stability line of +1 km) did not exceed 2 km (text S1). That corresponds to roughly ~14% of the 640-m globally equivalent layer (GEL) estimate and ~1% of the 5-km GEL estimate of water considered necessary to form VNs (text S1). The ice sheet thickness in Noachian highlands could have significantly deviated from our estimate; however, the goal here is to assess whether even a 2-km-thick ice cap could have undergone melting with our estimated geothermal heat constraints. Ice sheets thicker than 2 km during the Noachian could have melted with even lower geothermal heat flow than the estimates that we provide here. Similarly, higher surface heat flux would have been required to melt ice sheets thinner than 2 km.

Considering a mean annual surface temperature of 230 K (20), we find that over 0.5 to 1 million years (Ma), the melt column from a 2-km-thick ice deposit ranges between 150 and 800 m m−2 for surface heat flow between 60 and 85 mW m−2 (Fig. 2). The basal melting of ice deposits thinner than 1 km requires surface heat flow in excess of 100 mW m−2, irrespective of the snow accumulation rate or the time elapsed since their deposition (Fig. 2). In addition, the surface temperature could be considerably lower than 230 K if a substantial portion of Mars was covered by ice. This is because ice’s high albedo reduces the insolation available to warm the surface (text S2). For a correspondingly lower temperature of 200 K, no melt is generated over 1 Ma, regardless of the snow deposition rates considered here (text S3), even with surface heat flow exceeding 100 mW m−2. Nevertheless, any thick martian ice sheets may have been covered by thin layers of dust or volcanic ash as also observed for recently detected ice sheets, thus reducing total albedo and alleviating the problem associated with the icy-highlands melt scenario (text S2). Accordingly, we use 60 mW m−2 as the minimum surface heat flow for basal melting in the subsequent discussion.

Fig. 2 Heat flow required for basal melting and the crustal heat flow budget during the Noachian.

(A and B) Basal melt expected as a function of the geothermal heat flux for a variety of ice accumulation rates and total thickness of the ice sheet. ka, thousand years. Two-dimensional density plots (C and D) show crustal heat flow as a function of the crustal heat production rate and crustal thickness for crustal densities of 2500 and 3000 kg m−3. The black contour shows the isochron for a heat flux of 60 mW m−2. As expected, only thick crust with a high heat production rate could generate surface heat flow that exceeds 60 mW m−2. (E and F) Mantle heat flow contribution required for surface heat flow to exceed 60 mW m−2 during the Noachian for regions in (A) and (B) where the crustal heat flow was less than 60. Areas in gray are for situations where mantle heat flow above 40 mW m−2 is required for the surface heat flow to exceed 60 mW m−2.

Whether basal melting by geothermal heat flux occurred depends on the Noachian surface heat flow of Mars. The surface heat flow of a planet has contributions from both the crust and the mantle, and, to a lower degree, from the deeper interior. First, we consider a nonphysical scenario of the crust as the only thermal source, excluding the mantle and core. That reveals the upper limit on the crustal heat flow from a range of plausible values for thickness, density, and heat production rate of the crust. Using the Mars Odyssey Gamma Ray Spectrometer (GRS)–derived mass abundance of key heat-producing elements (HPEs; Th and K) (21), cosmochemically estimated U abundance from Th (fig. S2), and crustal thickness models of Mars (fig. S3), we find the present-day crustal heat production rate to vary between 2.5 × 10−11 and 10 × 10−11 W kg−1, similar to previous estimates (fig. S4) (22). As in prior works (18) (justified in the methodology), we simplify both the current spatial variation of HPE in the martian crust as representative of the Noachian and the decimeter-deep bulk regolith chemistry as representative of crustal depths. The implied Noachian heat production rate (using Eq. 6 in Materials and Methods) is ~1 × 10−10 to 3 × 10−10 W kg−1 (fig. S4). For a crustal bulk density of 2500 kg m−3, only thick crust (>80 km), with a high heat production rate of at least 2.4 × 10−10 W kg−1 could have generated crustal heat flow in excess of 60 mW m−2 (Fig. 2). The average crustal thickness of the southern highlands has been constrained to 57 ± 24 km using geoid-topography techniques (23), and almost 80% of the thickness of the crust is estimated to have been emplaced by the Noachian (24). We inverted the gravity data of Mars (25) to estimate the crustal thickness and find the average southern highlands crust to be 55 km for a crustal density of 2900 kg m−3. Combining the heat production rate from GRS data with the crustal thickness estimates, we find that the crustal heat flow would not have exceeded 60 mW m−2 anywhere on Mars (Fig. 3). Exceeding crustal heat flow of 60 mW m−2 becomes more unlikely with decreasing density as the crust becomes thinner (fig. S3). Combining the gravity-derived crustal thickness estimates (fig. S3) with HPE-based crustal heat production rate estimates (fig. S4), we find that the crustal heat flow in Noachian Mars could have ranged between 5 and 45 mW m−2, approximately four times the present-day estimate of crustal heat flow (22), insufficient to independently induce basal melting of even thick ice deposits (Fig. 3).

Fig. 3 Surface heat flow estimates for Noachian Mars.

(A) Crustal heat flow based on GRS chemical maps (fig. S2) and crustal thickness map of Mars (fig. S3) during the Noachian. (B and C) Surface heat flow arising from the combined crustal (A) and laterally constant mantle (denoted as qm) contributions. An H mask has been applied to exclude areas with high H content.

As expected, for basal melting to occur, our crust-only model suggests that a notable contribution from mantle heat flow is required. The lack of lithospheric flexure in the north polar region of Mars suggests that the current mantle heat flow likely does not exceed ~10 mW m−2 (26). Using that as the upper limit on the present martian mantle heat flow, the Noachian mantle heat flow would have been at least four times the current value, i.e., ~40 mW m−2. The fourfold higher heat flow during the Noachian is a natural consequence of radioactivity of HPE with billion-year half-lives. Consequently, 40 mW m−2 can be considered a conservative upper limit for Noachian mantle heat flow. Given that upper limit in the Noachian, we ask how much thermal contribution from the mantle would allow the surface heat flow to reach 60 mW m−2 necessary for basalt melting. From our crust-only model results, we find that a mantle heat flow between 20 and 40 mW m−2 would suffice across the southern highlands (Fig. 3). The range of modeled Noachian surface heat flow in Fig. 3 is consistent with previous estimates of surface heat flow derived from elastic loading models (27), crater relaxation models (28), and mantle convection models (29). Furthermore, mean surface heat flow estimates for current Mars range between 20 and 25 mW m−2 (30, 31). Heat flow during the Noachian would have been at least four times the present-day value; thus, even the lowest estimate of the mean surface heat flow value of modern Mars suggests Noachian heat flow in excess of 80 mW m−2. On Earth, the radioactive decay of HPE contributes only about half of the total heat flux, and the primordial accretionary heat supply has not yet been exhausted (32). Thus, it is likely that primordial heat would potentially increase the effective mantle heat flux. Consequently, our estimates are unlikely to overestimate the Noachian surface heat flow.

The surface heat flow maps from Fig. 3 provide the boundary conditions needed to estimate how much meltwater could be produced from ice sheets of various thicknesses. No melt is produced from sheets thinner than 1 km. For thicker units, over 0.5 to 1 Ma, a considerable amount of meltwater forms by basal melting from the surface heat flow range that we have constrained here (Fig. 4). For example, per square meter of a column of thick ice, melt exceeding 800 m is predicted in some parts of Mars over 1 Ma. Within the limits of the models described here, melt is only observed in locations with high surface heat flow, which correspond to areas on Mars with a high concentration of HPE, thick crust, or both. Additional heat from impacts, volcanism, and mantle plumes could have substantially magnified the geothermal heat.

Fig. 4 Meltwater from the estimated Noachian surface heat flow on Mars.

Meltwater produced (as column equivalent from a 2-km ice sheet) for surface heat flow for a 40 mW m−2 contribution from the mantle for 1 Ma (A) and 0.5 Ma (B). A Noachian Mask has been applied to exclude areas that are younger than Noachian in age. The red and black symbols show the location of hydrous minerals on Mars. (C and D) Same as (A) and (B), but with a mantle heat of 30 mW m−2.

The uncertainties, associated with the Noachian surface heat flow estimations that are notably agnostic of lateral variations in mantle heat, mostly preclude quantifiable insight from spatial correlations between heat flow and the distribution of hydrous minerals. The lack of concentration maps with global coverage for key hydrous minerals further limits the informativeness of these analyses. Nevertheless, we consider limited qualitative comparisons between our modeled heat flow and hydrous mineral locales. Our heat flow models, where HPE variance drives the thermal flux variation, indicate relatively low surface heat flow in Arabia Terra compared to the rest of the southern highlands (Fig. 4). Coincidentally, Arabia Terra also has relatively fewer sites with hydrous minerals compared to the rest of the southern highlands (Fig. 4). However, the lack of abundant hydrous minerals in Arabia Terra may reflect extensive burial and erosion that it is thought to have undergone (33). An extensive dust mantling could also obscure HPE content of the underlying crust.

The Terra Sirenum–Terra Cimmeria region has the highest concentration of Th and K and resulting crustal heat flow in the southern highlands of Mars (text S4, fig. S5, and Fig. 3). Terra Sirenum not only contains the largest cooccurrence of chlorides and phyllosilicates on Mars (Fig. 4) (34) but has also been interpreted to bear paleosubmarine hydrothermal units from the synthesis of chemical, magnetic, geomorphic, and mineralogic observations (35). The chlorides in Terra Sirenum have been hypothesized to form as a result of groundwater discharge and evaporation (34). Chlorides and other salts could have also notably contributed to the eutectic depression of the melting point of ice. The discovery of alunite in Cross crater (36) provides further evidence for the presence of groundwater in Terra Sirenum. Water in Cross crater was likely supplied by regionally upwelling groundwaters. Nevertheless, the Sirenum-Cimmeria region preserves the strong K-Th spatial correlation seen elsewhere on Mars, reducing the likelihood of temporally sustained aqueous activity (37).

An empirical cumulative distribution function (ECDF) plot of heat flow values from areas bearing various hydrous minerals shows that locations bearing chlorides have consistently higher heat flow than those of other mineral phases as well as Noachian terrane generally (text S5 and fig. S6). For example, while only ~40% of chloride sites correspond to surface heat flow not exceeding 70 mW m−2, that corresponds to nearly 60% of Noachian terrane and phyllosilicate sites as well as 80% of sulfate sites. We used a Kolmogorov-Smirnov (K-S) one-sided hypothesis test for further validation. The null hypothesis that the heat flow ECDF of chloride bearing regions and that of Noachian terrain represent the same distribution fails at better than 95% confidence (P ~ 0.001). The underlying signed difference in the distributions then supports the alternative hypothesis that the cumulative heat flow of chloride terrain exceeds that of the Noachian counterpart. Similarly, the null hypothesis that the heat flow ECDFs of sulfate-bearing regions and Noachian terrain are from the same distribution tests false (P = 0.0082). The sign indicates lower heat flow of sulfate terrain compared to the Noachian terrain. To first order, our case examples across Arabia and Sirenum-Cimmeria, along with ECDF and K-S test results over hydrated mineral locales, bolster the idea that locally higher geothermal heat flow in the Noachian enabled crustal conditions conducive to hydrous minerals. Furthermore, the detection of Fe/Mg smectites, chlorite, prehnite, illite/muscovite, serpentine, opaline silica, and analcime in the walls and far-flung ejecta of craters suggests that these minerals were excavated from substantial depths, where the Noachian crust experienced notable subsurface hydrothermal alterations at temperatures ranging between 50° and 400°C (16). Even with the most conservative geothermal estimates from our models, the temperature at kilometer depth would readily exceed that threshold (Fig. 5).

Fig. 5 Crustal conductive geotherm for a variety of heat flow estimates (qs) and thermal conductivity of the crust (k).

The vertical blue line denotes temperature value of 50°C (323 K), past which hydrothermal alterations can occur.

The faint young Sun paradox also challenges the formation hypotheses of VNs on Mars by surface water. Solving it likely requires endogenous parameters controlling the stability of liquid water on Mars (and Earth), as alternative stellar evolution models cause inconsistencies with helioseismology (38, 39). Groundwater production from geothermal heating, as we show here, may at least partially explain the formation of VNs on Mars. VNs could have then formed via sapping, by mass wasting aided by groundwater seepage, or by floods (40). Likewise, basal melting from high geothermal heat flow would slowly accumulate meltwater at the base of ice sheets and either infiltrate into the warm ground, form subglacial streams, or sustain periglacial surface water. Ice and snow packs can continually recharge by the migration of ice toward high-altitude regions due to adiabatic surface cooling, which would recharge VN sources after melting events in a predominantly cold climate (11). The rapid decline in the rate of valley formation at the end of the Noachian is also readily explained by a declining heat flow.

Even if Mars was once a warm and wet planet (41), over time, with the loss of the magnetic field, atmospheric thinning, and subsequent drop in global temperatures, liquid water would have been stable primarily in the subsurface (Fig. 5). Liquid water from basal melting in cold-icy highlands or from insolation-driven processes under a warm, wet Mars would have infiltrated to great depths owing to Mars’s fractured crust and lower gravity. Previously, Clifford (42) computed the thickness of groundwater that can reside in the martian subsurface by computing compaction curves to calculate the depth at which the porosity on Mars is less than 1%. The self-compaction depth predicted by that model estimated the depth of zero porosity to approximately 8.5 km and the total pore volume of the martian crust to be sufficient to store a GEL approximately ~500 m deep.

Previously, it has been shown that complex organic molecules would likely not be preserved in the shallow subsurface (<10 cm) of Mars over a time scale of few hundred million years (43). Thus, if life ever originated on Mars, then it may have followed the groundwater table to progressively greater depths where stable liquid water could persist. In addition, the deep subsurface would have guarded early life from the Late Heavy Bombardment (44). Any relict biomarkers from early martian life are also likely much better preserved at depth given the shielding from harmful radiation. At these depths, life could have been sustained by hydrothermal activity and rock-water reactions provided by the elevated Noachian geothermal heat flow, reminiscent of kilometers-deep ecosystems of Earth (45). In summary, we demonstrate that the geothermal heat flux would have played a key role in martian habitability and hydrology during the Noachian.

The faint young Sun paradox remains an outstanding problem in Mars science, especially with respect to the ancient hydrology and climate. Here, we modeled the specific conditions required for the basal melting of ice sheets on Mars revealing that basal melting provides a simple solution to the paradox. We show that geothermal heat flux in excess of 60 mW m−2 could have melted the base of ice sheets thicker than 1.5 km 4 Ga ago. Using the available geophysical and geochemical data of Mars, we show that prerequisite conditions for basal melting are readily satisfied. In addition, we show that geothermal heat flow in the shallow subsurface of Mars could have sustained hydrothermal alteration regardless of the Noachian climate. Improved constraints on the surface heat flow on Mars, as may result from InSight and future derivations of additional HPE abundances from gamma spectra, will allow us to better assess the role of geothermal heat in the habitability of Mars during the Noachian.


Thermal model

To investigate the thermophysical evolution and basal melting of hypothetical martian ice sheets, we constructed a one-dimensional finite-difference model capable of simulating their deposition, densification, and conduction of geothermal heat (fig. S1). Deposition of the ice sheet is assumed to occur uniformly over the martian surface at a steady rate, b·, until the ice reaches a critical thickness, Htot, here used as a proxy for the extent of Mars’s glaciation. During this steady-state accretionary period, the density profile of the ice sheet is given by equation 2 of (46)dρdz=f0ρ2ρi1b·exp(QRT)[[ρiρ1]ρi(gρ(z)dz)ρ]n(1)where ρ is the density of the ice sheet at depth z, f0 is a constant coefficient [see (47)], ρi is the density of ice (zero void space), b· is the accumulation rate, Q is the activation energy, R is the gas constant, T is the temperature of the ice sheet at depth z, g is the acceleration due to gravity, and n is the ice deformation constant (power law of applied stress). Assuming a surface density, ρ0 = 350 kg/m3, Eq. 1 can be numerically integrated to produce density profiles within the growing ice sheet. Once the ice sheet reaches its critical thickness, Htot, accumulation ceases and the temporal evolution of the density profile is governed by equation 1 of (46)dρdt=f0ρ exp(QRT)[[ρiρ1]ρi(gρ(z)dz)ρ]n(2)where t is time. Equation 2 is solved using an explicit finite-difference discretization in time. In either case (steady-state accretion or postcritical thickness densification), the thermal and phase evolution of the ice sheet will be governed by the conduction of geothermal heat through the ice, latent heat associated with phase change (basal melting), and densification (Eqs. 1 and 2). We implement a volume-averaged conservation of energy coupled with the enthalpy method to simulate the evolution of the ice sheet (48, 49)ρc¯Tt=z(k¯Tz)ρiLϕwt(3)H=ciceT+Lϕw(4)ϕw={0H<Hs=ciceTm(HHs)/Lif HsHHs+L1 H>Hs+L(5)where c is the specific heat capacity, k is the thermal conductivity, L is the latent heat of fusion for the water to ice phase transformation, ϕw is the liquid fraction, H is the enthalpy, Hs is the enthalpy of a discretization cell consisting of only solid ice less any void space (ϕv = 1 − ρ/ρi, ρ from Eq. 1 or 2) that has not been eliminated by compaction, and Tm is melting/freezing temperature. Subscripts ice, w, and v refer to characteristics of the ice, melt, and void space components of the three-phase system, respectively, and variables carrying an over bar are volumetrically averaged quantities [i.e., y¯=ϕwyw+ϕvyv+(1ϕwϕv)yice]. Equation 3 ensures conservation of heat, and Eqs. 4 and 5, combined, make up the enthalpy method. The advantage of the enthalpy method is its ability to accurately simulate the active formation and evolution of a basal melt phase, in contrast to predictions of basal temperatures that exceed the melting point but do not include a melt phase (46, 50), allowing for volumetric estimates of subglacial melt production.

We simulate ice sheets as being in contact with an atmosphere of constant temperature, TS (Dirichlet boundary condition), and which are subject to the constant basal heat flux (Neumann boundary condition) in fig. S1. Given a density profile (output of Eq. 1 or 2), Eq. 3 is solved using an implicit finite-difference discretization in time and a centered spatial discretization. Equations 3 to 5 are iterated until both the temperature profile and the melt fraction profile stabilize, upon which these profiles are passed back to the densification calculation of Eq. 1 or 2, and the process repeats itself for the next time step. Given that the triple point of H2O is around 650 Pa, assuming a low-density surface firn of 350 kg m−3, less than 1 m of firn is needed to exceed the triple-point pressure. Thus, the basal pressure of all appreciably thick ice sheets will be greater than that of the triple point.

In its current state, our one-dimensional model simulates the accretion, densification, and thermal evolution of martian ice sheets, including the ability to predict the presence and quantity of basal melt. Two simplifications are introduced regarding any basal melt that is produced. First, as this is a one-dimensional model, there is no lateral transport of melt into or out of the domain, meaning that any basal melt produced accumulates as a layer of water beneath the remaining ice sheet (Fig. 2). In reality, hydraulic gradients would likely lead to the flow of meltwater through basal hydrological systems (51), akin to their terrestrial counterparts (52). Here, we have focused on the regional production of basal melt to assess its relationship with hydrous mineral distribution. Second, while there is no convection term in our conservation of energy equation (Eq. 3), it is likely to be the dominant form of heat transport in any basal melt. Given the low viscosity and thermal diffusivity of water, a basal melt layer heated from below will eventually become gravitationally unstable and undergo convection, mixing the liquid layer. Given that the convective time scale for the basal melt is always much less than the temporal discretization implemented in the finite-difference model (100 years), we assume that any geothermal heat conducted into the basal melt is evenly distributed throughout the liquid phase.

The primary input parameter being investigated in this application of the forward model is the critical thickness, Htot, of the martian ice sheet. This quantity determines the relationship between the regional production of basal melt and the extent of glaciation, providing a method to assess the spatial distribution of hydrous minerals and meltwater production under disparate ice sheet conditions. The model’s sensitivity to variations in a number of parameters was explored (accumulation rate, surface temperature, and void space composition), but all negligibly affected the results when compared to ice sheet critical thickness. A description of all variables, their units, and the values used during simulations can be found in table S1.

Noachian geothermal heat flux

In addition to the ice thickness, whether basal melting occurred during the Noachian depends heavily on the surface heat flow. We do not have in situ estimates of laterally varying surface heat flow of contemporary Mars; thus, any estimate of the Noachian surface heat flow will be riddled with considerable caveats and uncertainties. Our goal is thus not to attempt to estimate the Noachian surface heat flow but to assess whether the currently constrained geophysical and geochemical properties of Mars, when extrapolated to the Noachian, could have allowed basal melting.

The three major sources of heat in planetary interiors are due to (i) accretion of the planet by impacts, (ii) core formation/differentiation and the associated release of gravitational potential energy, and (iii) the radioactive decay of unstable isotopes. These processes were most intense during the first few million years of the planet’s history; thus, the internal heat energy of the terrestrial planets was exponentially greater in the early stages of their histories. Both accretion and core formation are rapid processes; thus, the heat associated with these processes is most substantial at the early time when accretion and core formation occur. Nevertheless, on the basis of geoneutrino measurements, it has been proposed that the radioactive decay of the HPE contributes only about half of Earth’s total heat flux and that Earth’s primordial heat supply has not yet been exhausted (32). If similar conditions apply for Mars, then our work only underestimates heat flow. After accretion and differentiation, the major source of heat in the interior of planets is the decay of long-lived HPEs with half-lives of billions of years (e.g., 238U, 235U, 232Th, and 40K). The surface heat flow (qs) includes contribution from the heat generated in the crust (qc) and the mantle (qm) by the radioactive decay of long-lived HPE. The relative contribution from the crust and mantle to the surface heat flow depends on the distribution of the HPE within the crust and the mantle.

Crustal heat flow. The crustal heat flow depends on the crustal heat production rate and the thickness of the crust. The crustal heat production rate (Qc) is given byQc=[0.9928 CUH238U exp(tln2τ12238U)+0.0071 CUH235Uexp(tln2τ12235U)+CThH232Thexp(tln2τ12232Th)+1.191×104CKH40Kexp(tln2τ1240K)](6)where C and H are the concentration and heat release constants of the radiogenic elements, t is time, and τ12 are the half-lives of the radioactive elements (53). Crustal heat flow (qc) is given by the product of the crustal heat production (Qc), crustal density (ρcr), and crustal thickness (Tcr).

We estimate the crustal heat production rate using the chemical mass fraction maps of K and Th derived from GRS data (21), with latest data shown in fig. S2. Chemical maps show lateral heterogeneity in the distribution of K and Th at a 5° × 5° pixel resolution; U mass fractions were calculated for a cosmochemically constant Th/U mass ratio of 3.8 (fig. S2) due to the lack of GRS-derived abundances. We calculated the radioactive 232Th, 235U, 238U, and 40K, assuming that 232Th is 100% of total Th abundance, 235U and 238U are 0.7204 and 99.2742% of total U abundance, respectively, and 40K is 0.012% of total K abundance, similar to previous work by Hahn et al. (22). Chemical maps represent the top tens of centimeters of the martian surface and, therefore, represent near-surface regolith, ice, and dust deposits. For such measurement to represent the bulk chemistry of the martian upper crust, it must be normalized to a volatile-free basis (22). That equates to a 7 to 14% increase in the K, Th, and U abundances (22), which we applied to the chemical maps by renormalizing to Cl, stoichiometric H2O, and S-free basis. At high latitudes, water ice saturating regolith pore space not only overwhelms the signatures of other elements in gamma spectra from a high–spectral resolution but low-efficiency instrument like the GRS but also requires compositional modeling different from lower latitudes (20, 54); therefore, the heat flow values at or near the poles are likely severely underestimated. Consequently, we use a cutoff based on stoichiometric H2O mass fractions, the H mask described in (21), to limit our heat flow estimates to martian northern mid-to-low latitudes mostly devoid of substantial shallow ice.

We explicitly simplify the crustal columns as homogeneous with the bulk regolith chemistry (55). That contrasts with Earth’s continental crust that typically exhibits an approximately exponential decrease in the abundance of the HPEs with depth. Such fractionation on Earth is due to intracrustal melting and tectonic processes that further concentrate the incompatible elements into the more evolved, terrestrial upper crust. However, these processes of intracrustal differentiation that so affect the terrestrial continental crust are unlikely to be active on Mars to any notable degree. Previous studies (21, 55) provide multiple other reasons in defense of vertically homogeneous distribution of HPE. For example, (i) the similar geochemical behavior of K and Th implies that the K/Th ratio at the surface should reflect the ratio in the mantle source region, especially in the absence of plate tectonics. (ii) The high rate of volcanic intrusions with respect to extrusions will act to average out a possible depth dependence of K and Th concentrations. (iii) Impact gardening will have mixed the crust to considerable depth, particularly for large basin-forming impacts. Although it has been suggested, there is no definitive evidence for current or past plate tectonic processes, nor have orbital remote sensing instruments detected any major provinces of the highly evolved crustal products that we see on Earth or the Moon. Chemical provinces do not suggest that ejecta from the largest (and deepest) impact basins have compositions that differ notably from surrounding areas. Thus, measured surface abundance can be reasonably inferred to be representative of concentrations from greater depths. Another simplification is that relative lateral variations in HPE concentrations in the martian crust have not changed substantially since the Noachian. That is reasonable for Mars as previous work has shown that both K and Th vary only subtly with apparent surface age (56).

Crustal heat production further depends on the density of the crust. The spectroscopic and gravity investigations of Mars have shown that the composition and the density of the martian crust can vary greatly. Gravity investigation of martian volcanoes constrain the crustal density to around 3200 ± 100 kg m−3 (57), whereas the density of the southern highlands is found to range between 2500 and 3000 kg m−3 (58). Given the large spread on the estimates of the martian crustal density, we adopt a range of 2700 to 2900 kg m−3 for our work.

Last, crustal heat flow is given by the product of the crustal heat production rate, crustal density, and the crustal thickness. Previously, the study by Wieczorek and Zuber (23) constrained the crustal thickness of the southern highlands via geoid-to-topography ratio analysis and found a zero-elevation thickness of 57 ± 24 km under the assumption of Airy isostasy. Others have solved for the crustal thickness of Mars with the aid of gravity data (59). An explicit assumption in gravity-derived crustal models is that the density of the martian crust either is constant or does not notably vary laterally. With this assumption, crustal thickness models are derived by assuming that the lateral heterogeneities in Bouguer anomaly are primarily due to differences in crustal thickness. Here, we use similar methods to compute the crustal thickness of Mars from gravity data. The gravitational potential U exterior to a planet’s surface can be represented as a sum of spherical harmonic functionsU(r)=GMrl=0m=ll(R0r)lClmYlm(Ω)(7)where G is the gravitational constant, M is the total mass of the object, Ylm is the spherical harmonic function of degree l and m, and Clm represents the spherical harmonic coefficient of the gravitational potential at a reference radius R0. The analyses of tracking data from the Mariner 9, Viking 1 and 2, MGS, Mars Odyssey, and Mars Reconnaissance Orbiter have allowed the construction of the gravitation field of Mars. The most recent gravity model of Mars (JGMRO120D) nominally contains the spherical harmonic coefficients up to degree and order 120 (spatial resolution ~177 km).

From the spherical harmonic potential coefficients, the radial component of the gravitational field (gr) is computed viagr=GMr2l=0m=ll(R0r)l(l+1)ClmYlm(8)where G is the gravitational constant and M is the total mass of the object.

Assuming a mantle density of 3500 kg/m3 and crustal density in the range of 2700 to 2900 kg/m3, we calculated the Bouguer gravity anomaly from the surface topography as followsΔgBA=grgBAandgBa=4πρcrustR2M(2l+1)TlmYlm(θ,Φ)(9)where Tlm are the coefficients of the surface topography and ρcrust is the density of the martian crust. A standard definition of topography H in planetary geophysics isH(θ,Φ)=S(θ,Φ)A(θ,Φ)(10)where S is the radius from the center of mass and A is the height of the reference equipotential surface. The reference equipotential surface is a hypothetical surface that has a specific value of the potential, equivalent to the terrestrial geoid. The Mars Orbiting Laser Altimeter onboard the MGS spacecraft has made more than 640 million ranges to the surface, providing a near-global topographic coverage of Mars. We use the spherical harmonic model of the topography of Mars called the MarsTopo719.shape, which was calculated from the gridded datasets available through the Planetary Data System (PDS).

The resulting Bouguer anomaly for Mars is shown in fig. S3. We use the method in (60) to solve for the relief on the crust-mantle interface (i.e., Moho), required to explain the Bouguer gravity anomaly. The crustal thickness is then calculated by subtracting the relief on the Moho from surface topography. We assume that the observed Bouguer anomaly on Mars arises because of relief along the surface and crust-mantle interface (i.e., Moho) and solve for the Moho topography that would produce the observed Bouguer anomaly using SHTOOLS (61). The resulting crustal thickness maps are shown in fig. S3 for various crustal densities.

Mantle heat flow. The surface heat flow of a planet also has a substantial contribution from the mantle. The lack of lithospheric flexure in the north polar region of Mars suggests that the present-day mantle heat flow likely does not exceed ~10 mW m−2 (26). Assuming this to be an upper limit on the present-day martian mantle heat flow, the Noachian mantle heat flow would have been at least four times the present-day value, which equates to ~40 mW m−2. Thus, 40 mW m−2 is a reasonable upper limit for the Noachian mantle heat flow. Here, we first estimate the crustal flow during the Noachian using the approach described above. If crustal heat flow alone is insufficient for basal melting (as will mostly be the case), then we add contribution from the mantle treating 40 mW m−2 as the absolute upper limit on the Noachian mantle heat flow. We also run a few thermal models in which we ascribe mantle heat flow of 20 to 40 mW m−2 globally to constrain the amount of meltwater produced (Fig. 4).


Supplementary material for this article is available at

This is an open-access article distributed under the terms of the Creative Commons Attribution-NonCommercial license, which permits use, distribution, and reproduction in any medium, so long as the resultant use is not for commercial advantage and provided the original work is properly cited.


Acknowledgments: We would like to thank S. McLennan whose early reviews were substantial. Reviews from two anonymous reviewers helped improve and clarify this manuscript. Funding: This work was funded by a start-up grant to L.O. by Rutgers University. S.K.’s participation was supported by MDAP grant 80NSSC18K1375. Author contributions: L.O. and J.B. contributed to all the thermal and heat flow models. S.K. and M.S. provided expertise with deriving chemistry from GRS data, statistical compositional analyses, and planetary heat flow modeling. L.O. wrote the paper with notable feedback from all other authors. Competing interests: The authors declare that they have no competing interests. Data and materials availability: All data needed to evaluate the conclusions in the paper are present in the paper and/or the Supplementary Materials or through NASA’s PDS (
View Abstract

Stay Connected to Science Advances

Navigate This Article