Henri Samuel
email: samuel@ipgp.fr

Maxim D Ballmer

Sebastiano Padovan

Nicola Tosi

Attilio Rivoldini

Ana-Catalina Plesa

Mars' mantle early differentiation likely created an enriched layer at the coremantle boundary.

• A basal mantle layer delays planetary cooling and promotes deep melting through the present-day. • Massive deep mantle melting yields strong tradeoffs between core size and layer thickness when interpreting geodetic and seismic data.

Introduction

The present-day structure of Mars and other terrestrial planets results from billions of years of thermo-chemical evolution. It is known from geodetic data (gravity field, precession, and tides) that Mars is a differentiated planet with a liquid core [START_REF] Yoder | Martian precession and rotation from Viking lander range data[END_REF][START_REF] Yoder | Fluid Core Size of Mars from Detection of the Solar Tide[END_REF][START_REF] Van Hoolst | Chapter 18 -interior structure and evolution of mars[END_REF][START_REF] Smrekar | Pre-mission insights on the interior of mars[END_REF]. Such a large-scale differentiation is indirect evidence that the planet has experienced a global magma ocean stage during its early history. Indeed, as Mars was formed from the accretion of planetesimals containing both metallic iron and silicates, gravitational segregation combined with large-scale melting (and therefore low viscosities) of a mixture of both materials appears to be the only viable mechanism to efficiently separate metal from silicates at planetary scale. For example, both solidstate motion or diffusion would require time scales much longer than the age of the planet itself [START_REF] Stevenson | Models of the Earth's core[END_REF][START_REF] Karato | Core formation and chemical equilibrium in the Earth -I. Physical considerations[END_REF][START_REF] Rubie | Mechanisms of metal-silicate equilibration in the terrestrial magma ocean[END_REF]. In addition, the presence of an early magma ocean is also suggested by accretion scenarios. As planetary bodies reach sizes on the order of a few thousands of kilometers, the last stages of accretion become very energetic, and incoming impacts are likely to melt large fractions of the forming planet ( [START_REF] Senshu | Thermal evolution of a growing Mars[END_REF] and references therein). Moreover, the presence of short-lived radioactive heat-producing elements (HPE) such as 26 Al and 60 Fe contributed significantly to the occurrence of a deep and global magma ocean [START_REF] Nimmo | How rapidly did Mars accrete? Uncertainties in the Hf W timing of core formation[END_REF][START_REF] Dauphas | Hf-W-Th evidence for rapid growth of Mars and its status as a planetary embryo[END_REF][START_REF] Morishima | N-body simulations of oligarchic growth of Mars: Implications for Hf-W chronology[END_REF]. Finally, the very process of core formation generates large amounts of melting by converting gravitational potential energy into heat via viscous heating [START_REF] Senshu | Thermal evolution of a growing Mars[END_REF][START_REF] Samuel | Heat partitioning during core formation by negative diapirism in terrestrial planets[END_REF][START_REF] Rubie | Formation or Earth's core[END_REF]. Furthermore, the 182 W and 142 Nd isotopic anomalies that have been measured in Martian meteorites indicate that the Martian core formed within the first 13±2 Myr after the solar system formation [START_REF] Kleine | Rapid accretion and early core formation on asteroids and the terrestrial planets from hf-w chronometry[END_REF][START_REF] Foley | The early differentiation history of mars from 182w-142nd isotope systematics in the snc meteorites[END_REF]. All Martian meteorites show a strong depletion of highly siderophile elements [START_REF] Brandon | Evolution of the martian mantle inferred from the 187re-187os isotope and highly siderophile element abundance systematics of shergottite meteorites[END_REF], indicative of efficient metal-silicate separation. The latter can only be achieved if Mars has experienced a large-scale magma ocean during its early evolution [START_REF] Mezger | Core formation and mantle differentiation on mars[END_REF]. Thus, the presence of a silicate magma ocean at the end of core formation on Mars appears to be almost inevitable.

The solidification of a magma ocean during the progressive cooling of the planet is a complex process that involves significant chemical fractionation. The first solids that form are strongly depleted in incompatible elements, in particular HPE and iron oxides. During subsequent solidification, the newly formed cumulates become progressively more enriched (Elkins Tanton et al., 2003, 2005;Zeff & Williams, 2019). Since the crystallization of the Martian magma ocean is thought to occur from the bottom-up, the stacking of gradually iron-enriched and hence denser material results in a gravitationally unstable configuration. This gravitationally unstable stacking could lead to one or more episodes of Rayleigh-Taylor overturns of the cumulates (Maurice et al., 2017;Ballmer et al., 2017;Boukaré et al., 2018). Depending on a number of poorly constrained parameters (e.g., the solidification time of a magma ocean, the efficiency of melt-solid separation in the mushy freezing front [START_REF] Hier-Majumder | The origin of volatiles in the Earth's mantle[END_REF], or the crystal-melt density contrasts), the solidification of a Martian silicate magma ocean and the overturn of the resulting gravitationally unstable mantle stratification may ultimately lead to the presence of a significantly denser and enriched material at the bottom of the mantle compared to the overlying mantle. The enriched material could be either well-mixed and compositionally homogeneous (but distinct from the overlying mantle) or heterogeneous with a vertical compositional gradient (Maurice et al., 2017;Ballmer et al., 2017).

In both cases, if the iron enrichment is such that the induced compositional density contrast is significantly larger than thermal density contrasts, the compositionally distinct material will form a stable flat layer enveloping the Core-Mantle Boundary (CMB) [START_REF] Olson | An experimental approach to thermal convection in a two-layered mantle[END_REF][START_REF] Tackley | Strong heterogeneity caused by deep mantle layering[END_REF][START_REF] Lebars | Stability of thermal convection in two superimposed miscible viscous fluids[END_REF][START_REF] Samuel | Thermochemical convection and helium concentrations in mantle plumes[END_REF]Tosi, Plesa, & Breuer, 2013;[START_REF] Plesa | Can a fractionally crystallized magma ocean explain the thermo-chemical evolution of Mars? Earth Planet[END_REF][START_REF] Limare | Convection in an internally heated stratified heterogeneous reservoir[END_REF]. Such a stable basal layer can remain unmixed with the rest of the mantle for billions of years, with a negligible erosion [START_REF] Zhong | Entraiment of a dense layer by thermal plumes[END_REF]. The presence of such a mantle reservoir is also supported by the isotopic anomalies measured in Martian meteorites [START_REF] Harper | Rapid Accretion and Early Differentiation of Mars Indicated by 142Nd/144Nd in SNC Meteorites[END_REF][START_REF] Foley | The early differentiation history of mars from 182w-142nd isotope systematics in the snc meteorites[END_REF][START_REF] Debaille | Coupled 142Nd-143Nd evidence for a protracted magma ocean in Mars[END_REF]. Its long-term preservation can thus strongly influence the evolution and the present-day internal structure of Mars. However, the effects of such a strong mantle stratification on the long-term thermo-chemical evolution of Mars and on the interpretation of available and upcoming geophysical data have not yet been investigated in detail.

In this study, we quantify the consequences of a stable basal layer on the thermochemical evolution of Mars, and discuss the implications on the interpretation of available and upcoming geophysical data, with a focus on the ongoing InSight mission. The InSight lander touched down on the surface of Mars on November 2018 [START_REF] Banerdt | Initial results from the InSight mission on Mars[END_REF] and has since deployed short period and three-axis very broadband seismometers to record Martian seismic activity [START_REF] Lognonné | Seis: Insight's seismic experiment for internal structure of mars[END_REF]. The mission also aims at improving our knowledge of the core structure by precisely measuring the nutation of Mars with the radioscience experiment RISE [START_REF] Folkner | The Rotation and Interior Structure Experiment on the InSight Mission to Mars[END_REF]. Moreover, it features a heat-flow probe (HP 3 ) to measure the heat flux at the landing site [START_REF] Spohn | The Heat Flow and Physical Properties Package (HP 3 ) for the InSight Mission[END_REF]. Unfortunately, after almost one Martian year, the probe has failed to penetrate deep enough below the Martian surface, which prevents accurate measurments of tempearture time series and vertical thermal gradients.

The paper is organized as follows: in Section 2 we discuss the consequences of the presence of a deep enriched mantle layer resulting from the solidification of a Martian magma ocean on the composition of the mantle in HPE and iron. In Section 3, we model the thermal evolution of a Mars-like mantle in a stagnant lid convection regime with either no layering, or with a basal layer enriched in HPE and iron. In Section 4 we conduct a refined and systematic exploration of the parameter space, and extend the characterization of the influence of a basal layer on the thermo-chemical evolution of Mars, with Section 5 describing the obtained results. Section 6 discusses the implications of our results on Martian mantle layering on the interpretation of seismic, geodetic and heat-flow data, followed by a summary of the study in Section 7.

Mantle layering and enrichment

The crystallization of a Martian magma ocean may result in the formation of a stable basal layer enriched in iron and HPE. The enrichment of the basal layer depends on the style of crystallization and on its volume fraction relative to the entire silicate mantle.

Consider a differentiated planet of radius R p = 3389.5 km, with a metallic core of radius R c =1700 km, within the plausible range for Mars [START_REF] Van Hoolst | Chapter 18 -interior structure and evolution of mars[END_REF][START_REF] Smrekar | Pre-mission insights on the interior of mars[END_REF][START_REF] Rivoldini | Geodesy constraints on the interior structure and composition of Mars[END_REF], and a corresponding mantle volume V m = 4π(R 3 p -R 3 c )/3 and bulk volumetric heat production of radioactive elements, H m . In the presence of an enriched layer of thickness D d and volume V d above the core, the remaining overlying mantle volume is V m = V m -V d . The presence of the HPE-enriched layer, with associated heat production H d implies by mass balance that the overlying mantle heat production is reduced compared to the homogeneous case, namely:

H m = H m 1 - V d V m (Λ d -1) , (1)
where Λ d = H d /H m > 1 is the layer enrichment factor, the computation of which is detailed in the supporting information S1.

Figure 1 illustrates quantitatively the influence of deep mantle enrichment in HPE, assuming fractional crystallization and HPE abundances inferred by Wänke and Dreibus (1994) (i.e., U=16 ppb, Th=56 ppb, K=305 ppm).

The thickness of the basal layer cannot be constrained from magma ocean crystallization scenarios. It depends on the poorly-constrained physical conditions during overturn(s) of the gravitationally-unstable cumulate layers (Maurice et al., 2017;[START_REF] Ballmer | Compositional mantle layering revealed by slab stagnation at ˜1000-km depth[END_REF]Boukaré et al., 2018) (also see S1). A few scenarios could be ruled out. For example, a partial overturn promoted by the cold shallow temperatures, leaving highly enriched Fe cumulates at the surface is unlikely because the resulting mass distribution would not be compatible the moment of inertia of Mars [START_REF] Konopliv | An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data[END_REF][START_REF] Konopliv | Detection of the Chandler Wobble of Mars from Orbiting Spacecraft[END_REF] . Aside from such unlikely configurations, the a priori choice for the plausible range of layer thicknesses is arbitrary to some degree. We selected a range D d = 100 -500 km (or equivalently a range of layer volume fractions between 0.03 and 0.17) that explores values from thin layers (comparable to the thickness of boundary layers) to significant (yet smaller) thicknesses in comparison to that of the entire Martian silicate envelope. This range therefore allows one to assess the influence of the basal layer on Mars evolution. Within this range, the corresponding layer enrichment factor increases from about 5 to 19 with decreasing D d (Fig. 1a). This enrichment induces by definition a larger HPE content in the layer than in the overlying mantle, and therefore a larger radioactive heat production (Fig. 1b). The FeO content of the layer also increases with decreasing D d (Fig. 1c), leading to an increase in density ρ d relative to that of the overlying mantle, ρ m (see supporting information S1). This density contrast remains significant, as the layer buoyancy number

B d = ρ d -ρ m ρ m α∆T , (2)
which expresses the ratio of compositional to thermal density contrasts, is significantly larger than one (Fig. 1c, right axis) for the parameter range explored here (see Table 1). In the above equation, α is the thermal expansion and ∆T is the characteristic temperature scale, which is chosen as the superadiabatic temperature difference across the entire mantle. The large values of B d are sufficient to prevent convective mixing between the enriched layer and the overlying mantle despite its relatively large HPE content [START_REF] Olson | An experimental approach to thermal convection in a two-layered mantle[END_REF][START_REF] Samuel | Thermochemical convection and helium concentrations in mantle plumes[END_REF][START_REF] Lebars | Stability of thermal convection in two superimposed miscible viscous fluids[END_REF]Tosi, Plesa, & Breuer, 2013;[START_REF] Plesa | Can a fractionally crystallized magma ocean explain the thermo-chemical evolution of Mars? Earth Planet[END_REF][START_REF] Trim | The feedback between surface mobility and mantle compositional heterogeneity: Implications for the Earth and other terrestrial planets[END_REF][START_REF] Li | The stability and structure of primordial reservoirs in the lower mantle: Insights from models of thermochemical convection in three-dimensional spherical geometry[END_REF][START_REF] Nakagawa | Effects of thermo-chemical mantle convection on the thermal evolution of the Earth's core[END_REF][START_REF] Limare | Convection in an internally heated stratified heterogeneous reservoir[END_REF][START_REF] Langemeyer | The dynamics and impact of compositionally originating provinces in a mantle convection model featuring rheologically obtained plates[END_REF][START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF][START_REF] Li | The influence of deep mantle compositional heterogeneity on Earth's thermal evolution[END_REF].

3 Influence of the basal layer on the dynamic evolution of the Martian mantle: Finite-Volume modeling

To study the influence of iron and HPE enrichment in the basal layer quantified above on the dynamics of a Martian-like mantle, we consider the evolution of a solid-state slowly deforming mantle (i.e., with no inertia) in a stagnant-lid convection regime, under the Boussinesq approximation. Stagnant-lid convection is thought to currently occur inside Mars and most other terrestrial planets (the currently observed plate tectonics on Earth being an exception), and is essentially due to the strong dependence of mantle viscosity on temperature, leading to the presence of a very viscous lid in the coldest part of the shallow mantle. The viscosity η depends on temperature T and pressure P through the following Arrhenius relationship:

η = η 0 exp E * + P V * R T - E * + P ref V * R T ref , (3)
where E * = 200 kJ/mol is the effective activation energy, V * = 5 cm 3 /mol is the effective activation volume, R is the gas constant, and η 0 = 10 20 Pa s is the reference viscosity corresponding to the reference temperature T ref = 1600 K and reference pressure P ref = 3 GPa.

We carried out dynamic simulations in a 2D half-cylindrical domain of radial extent D = R p -R c = 1689.5 km in the (r,θ) space. To ensure temperature distributions resembling more closely those of a spherical geometry, we re-scaled the radius of the core to 1118.5 km (0.33 non-dimensional units) to keep the ratio of CMB-toplanet surface equal to that of a spherical body (Van Keken, 2001) (see also supporting information S2).

We considered two end-member cases. The first case corresponds to a compositionally homogeneous mantle. The other assumes compositional layering with a denser and enriched layer occupying the bottom 17% of the domain, whose iron and HPE enrichment decreases linearly with increasing height r -R c above the CMB. The domain is heated from below (T (r = R c ) = T c = 2000 K) and from within, and cooled from above (T (r = R p ) = T s = 220 K). All boundaries are free-slip and the side-walls are thermally insulating. The specific details of the modeling approach and equations are given in the supporting information S2. The enrichment in both iron and HPE is represented by a scalar, time-dependent compositional field, C(t, r, θ). In the homogeneous case, C = 0 everywhere. If the basal layer is present, the initial value of the compositional field decreases linearly from 1 at the bottom of the domain, to 0.2, at the top of the basal layer, and is set to zero elsewhere.

Four dimensionless numbers govern the dynamics of the system. The first is the thermal Rayleigh number that expresses the convective vigor:

Ra = ρ m gα∆T D 3 η 0 κ , (4)
where ρ m = 3500 kg/m 3 is the mantle density, g = 3.7 m/s -2 the gravitational acceleration at the surface of Mars, α = 2 10 -5 K -1 the thermal expansion coefficient, ∆T = T c -T s = 1780 K, and κ = 10 -6 m 2 /s is the thermal diffusivity. The second governing parameter is the buoyancy number defined in Equation 2. The third and fourth governing parameters are respectively the dimensionless internal heating parameter in the regular mantle and in the bulk enriched layer, when it applies. They correspond to the ratios of the whole-mantle Rayleigh number for internally heated convection to the whole-mantle thermal Rayleigh number defined above:

H m = ρ m D 2 H m k m ∆T , (5a)
H d = ρ m D 2 H d k m ∆T . ( 5b
)
The C-dependent non-dimensional heat production is:

H = H m + C d∆H dC , (6)
where d∆H/dC is set to 51.68 such that the volume-averaged internal heating in the basal layer amounts to H d . Similarly, the buoyancy number can more generally be expressed as follows:

B = C dB dC , (7)
where dB/dC is set to 6.8 such that the volume-averaged buoyancy number in the basal layer amounts to B d . Equations ( 6) and ( 7) are valid in the case of no stratification (C = 0), or in the case of a basal layer with a compositional gradient.

Given the value of physical parameters entering into the expression of the dimensionless quantities defined above, Ra = 2.2 10 7 . In the absence of compositional stratification, C = B = 0 and H d = H m . As done above, we assumed the HPE abundances of the [START_REF] Wänke | Chemistry and Accretion History of Mars[END_REF] Mars composition model for the mantle (i.e., U=16 ppb, Th=56 ppb, K=305 ppm), leading to a uniform heating with H = 5.8. When an enriched layer of constant composition is present, the corresponding mantle becomes heavily depleted in HPE and H m ∼ = 0.9 while there is a corresponding increase in H d ∼ = 29.8 (see Fig. 1b). The corresponding buoyancy number is B = B d = 3.8 (see Fig. 1c). In the case of a gradual enrichment of the basal layer, the non-dimensional heat production within the enriched layer increases from 11.2 at its top to 52.6 at its bottom. In these simplified experiments we do not account for radioactive decay, therefore the HPE input remains constant with time. This allows reaching a steady-state stage at which we can more easily compare the differences between the cases.

Each case starts with the same initial thermal condition: a uniform dimensionless temperature of 0.95 (i.e., 1691 K), with Thermal Boundary Layers (TBLs) of dimensionless thickness 0.05 (84.5 km) at the top and bottom boundaries, and a random perturbation of small amplitude to break the lateral symmetry. Each case is then evolved for 4.5 Gyr at which time evolution coincides closely with a statistical steady-state stage, where the averaged quantities (heat flux, temperature and velocities) do not evolve significantly or oscillate around a mean value. This corresponds to an elapsed time at which the mantle has 'forgotten' its initial thermal state.

Note that we do not impose a minimum viscosity for our finite-volume experiments. However, for computational efficiency purposes, we require that viscosity cannot vary beyond ten orders of magnitude. Such a requirement yields a maximum viscosity cutoff in the coldest part of the mantle, which however does not prevent the formation of a stagnant lid (see supporting information S2).

Figure 2 shows the model predictions at statistical steady-state. The temperature field in the homogeneous case (Fig. 2a) shows the presence of downwelling plumes originating at the base of a thermal boundary layer located underneath a thick cold, and therefore considerably more viscous layer (i.e., the stagnant lid, see also the horizontally averaged temperature profile, Fig. 2b). The magnitude of velocities in the lid is close to zero. In the mantle below, convective motions described above are vigorous, leading to an efficient homogenization of the temperature, as can be observed in the horizontally averaged temperature profile (Fig. 2b). The maximum viscosity contrast between the top of the lid and the CMB is 2.6 10 9 .

When considering a linear enrichment in iron and HPE with depth in the basal layer (Fig. 2f-g), the corresponding compositional density contrast stabilizes the layer against convection, and prevents significant mixing with the upper layer. Therefore, the compositional gradient remains preserved and efficiently prevents convection within the enriched layer, despite its higher temperature and corresponding lower viscosity. Consequently, no motion can develop and therefore heat is exclusively transported by conduction across it. This situation results in an increase of temperature in the layer and a vertically heterogeneous, purely diffusive temperature profile (Fig. 2d-e). In contrast, convective flow still occurs in the overlying mantle with upwelling and downwelling thermal plumes (Fig. 2d). The temperature in the convecting mantle region is lower than in the homogeneous case described above, which implies a smaller value of the corresponding effective Rayleigh number, and therefore a weaker convective vigor. In this case, the maximum viscosity contrast is 7.9 10 9 .

It should be noted that a possible alternative to the stably stratified basal layer exists, in which the enriched layer has a homogeneous iron and HPE content instead of the gradual increase of these quantities with depth considered above (Ballmer et al., 2017;Boukaré et al., 2018). A homogeneous basal layer opens up the possibility to double-layered convection as shown in the supporting information S3, which has a similar but considerably less pronounced influence on the thermal evolution than in the conductive basal layer case. However, the likelihood of occurrence of this scenario is small (supporting information S4), which is why we do not consider this possibility below.

Overall, the comparison of the two end-member cases described above has revealed the significant influence of an enriched denser basal layer at the bottom of a convecting mantle on the thermo-chemical evolution of the entire planet, and has shown that the stratification determines the efficiency of heat transfer, and shapes the thermal structure of the entire mantle, even after billions of years of evolution. In the following sections, we therefore systematically explore this influence in a more global and more exhaustive context of the thermo-chemical evolution of an entire Mars-like planet, i.e., including the buoyant crust and metallic core. The flat interface between the basal layer and the overlying mantle allows for a straightforward parameterization of the heat flux across this boundary because strong lateral variations in composition are absent.

Parameterized convection models: Approach

To model the thermo-chemical evolution of a Mars-sized planet, we used a parameterized approach [START_REF] Stevenson | Magnetism and thermal evolution of the terrestrial planets[END_REF][START_REF] Spohn | Mantle differentiation and thermal evolution of Mars, Mercury, and Venus[END_REF][START_REF] Hauck | Thermal and crustal evolution of Mars[END_REF][START_REF] Morschhauser | Crustal recycling, mantle dehydration, and the thermal evolution of Mars[END_REF]. This approach reproduces accurately the dynamic evolution of a stagnant-lid planet in spherical geometry with various complexities (e.g., temperature, melt-and pressure-dependent viscosity, heterogeneous heating), at a considerably smaller computational cost than modelling in 3D and even 2D geometries (Tosi et al., 2013;Plesa et al., 2015;Thiriet, Breuer, et al., 2018;[START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF]. The latter allows exploring a larger parameter space (tens of thousands of evolutions) compared with 3D or even 2D convection studies, as we did in this study. The parameterized approach computes the thermo-chemical evolution of a planet by considering the heat transfer between several concentric envelopes, as illustrated in Fig. 3a and Fig. 4a: an adiabatic convecting liquid iron-rich alloy, overlaid by a convecting silicate mantle that is covered by a conductive rigid lid. The latter includes a buoyant crust enriched in heat-producing elements.

The viscosity of the Martian mantle plays an important role, and is assumed to depend on temperature, T , hydrostatic pressure, P , and melt fraction φ following an Arrhenius relationship [START_REF] Karato | Rheology of the upper mantle: A synthesis[END_REF]:

η(T, P ) = max η 0 exp E * + P V * R T - E * + P ref V * R T ref -β * φ , 10 -2 , ( 8
)
where E * and V * are the effective activation energy and activation volume, R is the gas constant, and T ref and P ref are the reference temperature and pressure at which viscosity equals the reference viscosity, η 0 (in the absence of melt). The effective activation volume and energy can either directly account for viscous deformation in the diffusion creep regime, or mimic deformation in the dislocation creep regime (Plesa et al., 2015;[START_REF] Kiefer | Water undersaturated mantle plume volcanism on present-day Mars[END_REF]Thiriet et al., 2018;[START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF][START_REF] Schulz | Stagnant-lid convection with diffusion and dislocation creep rheology: Influence of a non-evolving grain size[END_REF]. In the first case, E * and V * correspond to the intrinsic values. In the latter case (dislocation creep), E * and V * correspond to the intrinsic values divided by the stress power-law index, whose value is close to 3.5 [START_REF] Christensen | Convection in a variable-viscosity fluid: Newtonian versus power-law rheology[END_REF]. The strong sensitivity of viscosity to temperature and the relatively small size of Mars imply that its mantle convects in the stagnant-lid regime. The above expression accounts for a dependence on melt fraction φ, and for a minimum threshold of 10 -2 Pa s based on experimental constraints on the viscosity of peridotitic liquids [START_REF] Liebske | Viscosity of peridotite liquid up to 13 GPa: Implications for magma ocean viscosities[END_REF]. The sensitivity of viscosity to melt fraction is expressed by β * , the value of which is set to 30 [START_REF] Hirth | Rheology of the upper mantle and the mantle wedge: A view from the experimentalists[END_REF]. The melt fraction is estimated as

φ(T, P ) = min max 0, T -T sol T liq -T sol , 1 , (9)
where we used the parameterization of solidus (T sol ) and liquidus (T liq ) curves appropriate for a Martian mantle [START_REF] Ruedas | On the relative importance of thermal and chemical buoyancy in regular and impact-induced melting in a Mars-like planet[END_REF][START_REF] Duncan | Extending the Solidus for a Model Iron-Rich Martian Mantle Composition to 25 GPa[END_REF]. Additional influences of shallow melt extraction and the influence of iron enrichment on the solidus and liquidus curves were taken into account and will be discussed later.

The crustal thickness evolves with time based on the occurrence of melting at shallow depths. The lithospheric thickness is determined by considering an energy balance between the convective heat flux at the top of the mantle, the conductive heat flux out of the lithosphere, and the energy consumed to transform a portion of convective mantle into additional viscous lithosphere material, and vice versa ( [START_REF] Schubert | Subsolidus convective cooling histories of terrestrial planets[END_REF][START_REF] Spohn | Mantle differentiation and thermal evolution of Mars, Mercury, and Venus[END_REF] and references therein). Similar to [START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF], the crust is enriched in HPE relative to the underlying mantle. Specifically, we assumed that the crust is ten times more enriched in HPE relative to the underlying mantle bulk mantle (Table 1). The latter corresponds to the entire silicate envelope in the homogeneous case, or the entire silicate envelope minus the basal layer if the basal layer is present.

The exact equations are given in [START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF] (note that there is a sign typo in the last term on the right hand side of Equation ( 20) in the previously quoted paper). As detailed below, heterogeneous internal heating is taken into account, as well as latent heat of melting-crystallization. In the case of a homogeneous mantle, model equations are those described in [START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF]. For completeness these are summarized below. Then, we detail below the approach we used to account for the presence of a denser and enriched layer.

Homogeneous mantle

In the absence of a basal layer, the parameterized approach consists of solving for energy balance for the convecting mantle and the core (see Fig. 3a for a schematic representation of the model), respectively:

ρ m C pm V m m (S t + 1) Ṫm = -q m + ρ cr Ḋcr [L m + C pcr (T m -T l )] A m + q c A c + H m V m , (10)
ρ c C pc V c c Ṫc = -q c A c , (11)
where dotted quantities indicate derivatives with respect to time t, T m is the uppermost convective mantle temperature, and T c is the temperature at the CMB; ρ m and ρ c are the mantle and the core densities; C pm , C pcr , and C pc are mantle, crust, and core specific heat at constant pressure; A m and A c are the surface areas of the convecting mantle and core; V m is the volume of the convective mantle, V c is the volume of the core, and H m is the volumetric internal heating rate due to the presence of heat-producing elements. L m is the latent heat of silicates melting; m expresses the ratio between the average temperature in the convecting mantle and the uppermost mantle temperature, T m [START_REF] Stevenson | Magnetism and thermal evolution of the terrestrial planets[END_REF], and is constantly updated, while c is the constant ratio between the average temperature of the core and T c , and is computed only once. These two quantities originate from the fact that Equation (10) and Equation ( 11) are written in terms of uppermost mantle and core temperatures, while the corresponding heat balance considers the average temperature for the convecting mantle and core. T l is the temperature at the base of the stagnant lid, at which viscosity has grown by one order of magnitude with respect to the convecting mantle, yielding [START_REF] Davaille | Transient High-Rayleigh-Number Thermal-Convection with Large Viscosity Variations[END_REF]:

T l = T m -a rh R T 2
m /E * , with a rh = 2.54, which provides a good correspondence between parameterized and 3D spherical modeling (Thiriet, Breuer, et al., 2018).

The (average) heat flux out of the surface of the convecting mantle, q m , is:

q m = k m T m -T l δ u , ( 12
)
where k m is the mantle thermal conductivity, δ u is the thickness of the upper thermal boundary layer of the convecting mantle. Similarly, the heat flux between the core and the mantle is:

q c = k m T c -T b δ c , ( 13
)
where T b is the temperature at the base of the convecting mantle. The latter is obtained by adding the adiabatic gradient to T m :

T b = T m + αgT m C pm ∆R, ( 14
)
where α is the mantle thermal expansion coefficient, g is the surface gravitational acceleration, and ∆R = R p -D l -R c -δ u -δ c is the thickness of the convecting mantle devoid of its thermal boundary layers.

The expression for the thickness of the upper thermal boundary layer stems from boundary layer theory:

δ u = (R l -R c ) Ra c Ra βu , (15)
where β u = 0.335 allows for the closest match between parameterized and 3D spherical modeling of stagnant lid convection with variable viscosity (Thiriet, Breuer, et al., 2018), and Ra c = 450 [START_REF] Choblet | 3D thermal convection with variable viscosity: Can transient cooling be described by a quasi-static scaling law?[END_REF] . The thermal Rayleigh number associated with the convecting mantle is defined as:

Ra = ρ m αg∆T (R p -D l -R c ) 3 η m κ , (16)
which expresses the mantle convective vigor. In the definition above, κ = k m /(ρ m C pm ) is the mantle thermal diffusivity, η m = η(T m , P m ) is the viscosity of the mantle below the stagnant lid, and ∆T = T m -T l + max(T c -T b , 0) is the sum of temperature differences across the upper and lower thermal boundary layers of the convecting mantle.

The thickness of the lower thermal boundary layer is:

δ b = κη c Ra δ b ρ m αg|T c -T b | 1/3 , ( 17
)
where η c = η((T b + T c )/2), P c ) is the mantle viscosity taken at temperature and pressure half-way across the lower thermal boundary layer. Ra δ b = 0.28 Ra 0.21 i is the bottom boundary layer Rayleigh number [START_REF] Deschamps | Inversion of two-dimensional numerical convection experiments for a fluid with a strongly temperature-dependent viscosity[END_REF], where

Ra i = ρ m αg∆T i (R p -R c ) 3 /(η m κ)
is the thermal Rayleigh number for the entire mantle thickness and ∆T i = T m -T s + max(T c -T b , 0) is the sum of the temperature contrasts across the entire rigid lid, and the basal thermal boundary layer.

The model accounts for crust formation, in which latent heat is consumed or released upon melting and crystallization at shallow depths, through the use of a time-dependent Stefan number that expresses the ratio of latent to specific heat: S t = L m (dφ m /dT m )/C pm , where φ m = (1/V m ) Vm φ(r)dV is the average melt fraction in the convecting mantle.

At depths where the pressure is below 7.4 GPa, the produced melt is buoyant [START_REF] Ohtani | Flotation of olivine and diamond in mantle melt at high pressure: Implications for fractionation in the deep mantle and ultradeep origin of diamond[END_REF][START_REF] Ohtani | Flotation of olivine and diamond in mantle melt at high pressure: Implications for fractionation in the deep mantle and ultradeep origin of diamond[END_REF], and is therefore assumed to be extracted upwards to contribute to the build-up of the crust. The model accounts for the fact that melt extraction alters the solidus and liquidus curves. For additional details related to shallow melt extraction and the crustal growth model we refer to [START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF] where the description for the growth rate Ḋcr is given. At depths where the pressure is above 7.4 GPa, no upward or downward melt extraction is assumed. The presence of melt in these deeper regions only influences viscosity, and the mantle energy balance through the consumption or the release of latent heat upon melting or crystallization, respectively (i.e., term containing S t Ṫm on the left-hand side of Equation 10).

The set of differential equations are integrated in time using a second-order Runge Kutta scheme with dynamic time-step, subject to the following initial temperature conditions: T c (0) = T c0 = 2100 K, T m (0) = T m0 = 1800 K, along with small values for D cr = 1 m and D l = 10 m. The values of the model parameters along with their meaning are listed in Table 1.

Layered mantle parameterization

When the enriched basal conductive layer is present (see Fig. 4a for a schematic representation of the model), an additional equation is required to describe the heat transfer outwards, inwards, and within the basal layer, along with additional modifications to the previously described energy balance. We neglect the possible erosion of the basal layer by plumes. Although such erosion exists (Figure 2f), it is relatively limited given the assumed density contrasts between the basal layer and the overlying mantle. Therefore, considering a constant thickness of the basal layer with time is a reasonable assumption. Under these conditions, the core energy balance (Equation ( 11)) remains unchanged. However, Equation (10) becomes:

ρ m C pm V m m (S t + 1) Ṫm = -q m + ρ cr Ḋcr [L m + C pcr (T m -T l )] A m + q d A d + H m V m , (18)
where V m = V m -V d is the volume of the convecting mantle without the basal layer (i.e., the volume contained within

r = R c + D d and r = R p -D l )
, and H m is the time-dependent output of HPE per unit volume within V m . Similarly, m is the timedependent ratio between the averaged temperature within V m and T m . A d is the surface area of the interface between the basal layer and the overlying mantle, q d is the heat flux across this boundary. The latter will be explicitly described further below.

As in Section 3, we set the initial content in HPE in the silicate envelope using the abundances listed in [START_REF] Wänke | Chemistry and Accretion History of Mars[END_REF]. The presence of the HPE-enriched layer implies by mass balance that the overlying mantle heat production H m is smaller than what is considered for the homogeneous case H m . Heat production in the depleted mantle H m is given by Equation ( 1). The thickness of the uppermost thermal boundary layer that enters for example in the expression of q m (Equation ( 12)) now becomes:

δ u = (R l -R c -D d ) Ra c Ra βu , (19)
where Ra , the Rayleigh number associated with the convecting volume V m is:

Ra = ρ m αg∆T (R p -D l -R c -D d ) 3 η m κ , (20)
and 2000) and T i is the temperature at top of the basal layer. The meaning of T b remains the same as in the homogeneous mantle case (i.e., the temperature at the bottom of the convecting mantle just above the lower thermal boundary layer). However, its expression now becomes:

∆T = T m -T l + ∆T b , where ∆T b = T i -T b = 1.43RT 2 m /E * (Deschamps & Sotin,
T b = T m + αgT m (R c + D d -D l -δ u -δ b )/C pm ,
where δ b is the thickness of the thermal boundary layer just above the interface between the basal layer and the overlying mantle. The latter is computed with the corresponding form of Equation ( 17).

For all cases shown in this study the Rayleigh numbers associated with the convecting mantle are always supercritical, which is consistent with our model assumptions and with the recent traces of volcanism observed at the surface of Mars [START_REF] Hartmann | Evidence for recent volcanism on Mars from crater counts[END_REF][START_REF] Neukum | Recent and episodic volcanic and glacial activity on Mars revealed by the High Resolution Stereo Camera[END_REF].

We considered the presence of a compositional gradient within the layer, such that the iron content linearly increases with depth [START_REF] Plesa | Can a fractionally crystallized magma ocean explain the thermo-chemical evolution of Mars? Earth Planet[END_REF]Maurice et al., 2017;Ballmer et al., 2017;Boukaré et al., 2018). As seen in Figure 2d-g, for a sufficiently large compositional gradient as we assumed here, thermal expansion cannot overcome the compositional gradients. This would result in a stratified layer that is stable against thermal convection. Additionally, even in the absence of compositional gradient the basal layer is likely to be too thin and too viscous to allow for convective motion (supporting information S4). Consequently, heat within the layer can only be transferred via conduction (Section 3). Therefore, similar to the thermal evolution within the stagnant lid, heat transfer across the stably stratified basal layer is described by the following time-dependent, spherically symmetric, diffusion equation:

ρ(r)C pm ∂T ∂t = k m r 2 ∂ ∂r r 2 ∂T ∂r + H(t, r) -ρ(r)L m ∂φ(t, r) ∂t ( 21
)
where T is the radially dependent temperature within the basal layer, r is the radius ranging here between R c and R c + D d , ρ is the density, H is the radially varying and time decaying heat production due to HPE. Note that the melt fraction φ is a function of r and t and that the basal layer remains diffusive regardless of the value of its viscosity. This remains true if the layer is entirely molten and has therefore very small viscosity (Equation ( 8)). Such stability against convection within the layer results from its increasing iron content with depth (supporting information S1). Since the layer enrichment increases linearly with depth, both the density and HPE content follow the same linear trend:

ρ(r) = ρ d f e (r), (22a)
H(t, r) = H d (t)f e (r), (22b)
where f e (r) is a linear function (see supporting information S5) that expresses the depth-dependence of the enrichment in incompatible elements within the basal layer.

Equation ( 21), with time varying boundary conditions (T (R c ) = T c and T (R c + D d ) = T i ) must be solved at each time step, and is discretized using finite-differences of second-order accuracy in space. The time integration can be performed using a firstorder implicit scheme. Higher-order explicit schemes (second-and third-order Runge-Kutta) were tested using sub-time-stepping. However, given the small step size, they did not result in a noticeable accuracy improvement. Therefore, we used the unconditionally stable and more efficient implicit scheme. With the knowledge of T (t, r), we can express the fluxes q c = -k m (∂T /∂r)| r=Rc and q d = -k m (∂T /∂r)| r=Rc+D d , where the temperature derivatives are obtained via second-order accurate finite differences.

Following Elkins-Tanton (2008) (Equation (2) therein), the influence of iron on both the solidus and the liquidus is accounted for by subtracting the term: 6 (Fe# d -Fe# m ), to T sol and T liq where Fe# m is the Fe-number (i.e., Fe# m =100 Fe/(Fe+Mg)) for the overlying mantle, and Fe# d is the Fe-number within the basal layer. The latter increases linearly with depth in the case of a stably stratified layer: Fe# d (r) = Fe# d f e (r) (as in Equation 22a), where Fe# d is the average iron number of the denser basal layer, whose computation is detailed in the supporting information S1 and S5.

Parameterized Convection Model Results

Using the approach described above, we extended in the following our comparison between the layered case and the homogeneous mantle case presented in Section 3, to the scale of a Mars-sized planet. Then, we explored systematically a wider parameter space defined by mantle rheological parameters to quantify the influence of the basal layer on various key quantities characterizing the evolution of Mars, and its resulting present-day structure. We checked that the parameterized convection model described above can reproduce the results displayed in Fig. 2 for the same conditions (i.e., no crust, Boussinesq approximation, fixed CMB temperature, no radioactive decay) both in the homogeneous and layered mantle cases.

Influence of the stably stratified basal layer

We first proceed with the comparison of the evolution of two selected cases. One without layering, and a second one with a 500-km thick basal layer in the mantle, analogous to the cases considered in Figures 2a-c andFigures 2d-g, with however additional complexities such as adiabatic heating, melting, crustal formation, radioactive decay, and core evolution. The values of the governing parameters are: R c = 1700 km, η 0 = 10 20 Pa s, E * = 200 kJ mol -1 , and V * = 5 cm 3 /mol. The values for all other quantities correspond to those given in Table 1.

Main evolutionary trends and present-day structure

The panels of Figures 3 and4 display the schematic view of the parameterized model along with the evolution of several key quantities (temperatures, crustal and lithospheric thicknesses) and the resulting present-day thermo-chemical structure, in the case of a homogeneous mantle and for a layered mantle with conductive basal layer, respectively.

The comparison of the present-day areotherms in the homogeneous (Fig. 3b) and layered case (Fig. 4b) shows that the layered case is globally hotter. The thermal evolution in the homogeneous mantle case is monotonous with a continuous cooling of the mantle and the core (Fig. 3c). On the contrary, the thermal evolution in the layered mantle case is more complex, with an increase in basal layer temperatures during the first 1 Gyr followed by a more steady decrease until the present-day (Fig. 4c). Unlike the homogeneous mantle case, the core temperature continuously and significantly increases, while the uppermost convecting mantle temperature decreases during the entire evolution. Importantly, the planet on average cools down more efficiently in the homogeneous case than in the layered case. The crustal and total lithospheric thicknesses evolve in a comparable way in the homogeneous (Fig. 3d) and layered cases (Fig. 4d). The timing for crustal formation occurring mostly during the first Gyr is similar in both cases, which is in line with photogeological estimates [START_REF] Greeley | Magma generation on Mars: amounts, rates, and comparisons with Earth, moon, and venus[END_REF][START_REF] Nimmo | Early Crustal Evolution of Mars[END_REF]. However, the predicted time evolution of crustal thicknesses is also sensitive to other model parameters, such as the initial thermal state or the value of the rheological parameters. Therefore the crustal evolution described above may change for different combinations of governing parameters (see e.g., Section 5.2).

The effect of the basal layer on the thermal evolution essentially originates from the reduced heat transfer between the layer and the overlying mantle (as noted in Section 3). This reduction limits deep mantle and core heat loss to space.

Comparative evolution of temperatures and heat fluxes across the planetary envelopes

Figure 5 displays the evolution of several additional quantities corresponding to the two cases described above and shown in Figures 3 and4. The presence of the basal layer induces a conductive flux, q d , at the interface (Fig. 5a). This heat essentially comes from the radioactive decay of HPE present in the layer. As time increases, the HPE content decreases due to the radioactive decay, which explains the observed 'bell-shaped' evolution of q d . The CMB heat flux becomes rapidly negative in the layered case (i.e., the enriched layer heats up the core, see Fig. 5b), while it remains always positive for the homogeneous mantle, leading to the continuous decrease of the CMB temperature displayed in Figure 3c. The presence of the conductive layer slightly delays mantle cooling (Fig. 5c), but only during the first few ∼ 100 Myr. Indeed, in the homogeneous case, the convecting mantle is more enriched in HPE than in the layered case. Accordingly, the mantle heat flux becomes larger in the layered case compared to the homogeneous case. This considerable difference results from the fact that the HPE in the layered case are concentrated in the deep mantle. Eventually, the HPE output is transferred to the overlying mantle via conduction at the interface, thereby enhancing q m , in the layered case. In contrast, in the homogeneous case, HPE are most abundant in the enriched buoyant crust above the mantle, and do not contribute to mantle thermal evolution (except in the shallowest part of the mantle where the enriched crust could slightly delay mantle heat loss by increasing the temperature locally). Consequently, the surface heat flux that accounts for mantle contribution and crustal heat production is larger in the homogeneous case than in the layered case (Fig. 5d). However, after 4.5 Gyr of evolution both cases show comparable heat flux values because the total HPE content is the same [START_REF] Wänke | Chemistry and Accretion History of Mars[END_REF], and only their distribution varies among the two cases. This prediction is highlighted in Figure 5e that displays the bulk (U r b ) and the convective (U r c ) Urey ratios, defined respectively as the ratio between the total heat HPE production, including and excluding that of the crust, to the total heat escape at the surface. Values larger than one for U r b and U r c are indicative of heating, and values smaller than one indicate cooling. However, by excluding the heat production contribution in the crust, the convective Urey ratio specifically expresses the efficiency of mantle cooling or heating, while the bulk Urey ratio is a proxy for the entire planet. In the homogeneous case, the crust grows rapidly, and the rest of the mantle is accordingly depleted of HPE. This leads to a very distinct evolution of bulk and convective Urey ratios (Fig. 5e). However, both U r b and U r c drop below unity early on, indicating that the mantle immediately looses heat, while the entire planet cools down during most of its history. On the contrary, for the layered case, U r b and U r c are very similar to each other, because most of the HPE are located in the deep mantle and contribute to the mantle heat balance. The Urey ratios are significantly above unity for the first half of the planet's history (Fig. 5e). At present-day, the convective Urey ratio for the layered case is considerably larger than that of the homogeneous case. On the contrary, the present-day surface heat flux is about 25% larger for the homogeneous case than in the layered case because of the different contributions of crustal heat production located just below the surface, and the HPE-enrichment in the deep layer.

The distribution of heat sources, and the evolution of the heat fluxes across the different planetary envelopes control the thermal evolution of each layer. As previously noted, the increase of q d followed by a decrease of q d , combined with the radioactive decay of heat sources in the basal layer yields the observed early increase in q m , followed by a decrease of the temperature in the basal layer (Fig. 4c). Similarly, the positive CMB heat flux (Fig. 5b) yields a rather slow and gradual cooling of the core in the homogeneous case (Fig. 3c), while the stronger and essentially negative q c for the layered case yields a significant increase in core temperatures during the entire planet history (Fig. 4c).

Influence of the basal layer thickness

The above description of the model results demonstrates that the presence of an enriched basal layer considerably influences the evolution of the planet, and impacts its resulting present-day thermo-chemical state. This is illustrated in Figure 6 that compares the present-day thermal profiles obtained for a homogeneous mantle (Figure 6a) with a case with a 500-km-thick basal layer (Figure 6b). The solidus and liquidus curves are also displayed and show that a large fraction of the basal layer is fully molten. This is due to both the hot temperature of the enriched layers, and the depression of the melting curves due to the iron enrichment in the deep mantle. Note that melting of the basal layer may lead to a reduction of the compositional density contrast between the basal layer and the overlying mantle. However, the iron enrichment is sufficiently large to maintain density of the layer considerably larger than that of the overlying mantle (supporting information S6). Therefore deep melting is unlikely to affect the stability of the basal layer for the cases considered in this study.

Figure 7 shows the present-day thermal state for different cases with either a homogeneous mantle or a basal layer of various thicknesses, D d . The presence of an enriched basal layer also prevents core heat loss because the layer acts as a heat buffer, and even heats up the core for a large part of the planet's history (Fig. 4c and Fig. 5b), leading to an increase in present-day core temperatures, T c . Figures 7b-c display the present-day crustal thicknesses and surface heat fluxes as a function of the layer thickness with other governing parameters being the same as in the layered case considered above. For relatively thin layers the crustal thicknesses are about 10% larger than that of the homogeneous case. We observe a decrease in crustal thickness with increasing D d , which becomes more pronounced for thicknesses above 400-km. This effect is due to the fact that thicker basal layers further deplete the shallow mantle in HPE and delay the heat transfer from the basal layer to the overlying mantle. This leads to a colder shallow mantle early on, and therefore to smaller associated crustal production rates with increasing values of D d . However, the corresponding range of D cr remains comparable to the value obtain for a homogeneous mantle for the range of D d explored. The present-day heat flux steadily decreases with increasing D d , (Figure 7c). The observed rather modest influence of D d on present-day surface heat flux may be due to the fact that for a given set of rheological parameters, the present-day surface heat flux is mostly governed by the HPE content of the bulk mantle, which does not change among the different cases considered.

Combined influences of the stably stratified basal layer and mantle rheology

The investigations described above did not explore the effects of mantle rheological parameters, which are key quantities to the thermo-chemical evolution of Mars. Therefore, in the following we consider their influences and focus on the case of a 200km-thick basal layer. As previously discussed, different values of layer thickness do not impact the qualitative behaviour described below. Figure 8 shows present-day values of several main quantities as a function of the reference mantle viscosity, η 0 , and the effective activation energy, E * . The average planet temperature (Fig. 8a), increases with increasing η 0 and E * . An increase in either of these two parameters implies an increase in mantle viscosity (Eq. ( 8)), which reduces the efficiency of convective heat transfer, and eventually diminishes planetary heat loss, leading to larger temperatures. This general trend is observed for the depleted mantle (T m , Fig. 8c), the enriched layer (T d , Fig. 8d), and the core (T c , Fig. 8e).

The temperature at the base of the crust, T cr (Figure 8b) shows a more complex dependence on η 0 and E * , but is essentially most sensitive to E * . In fact, T cr is more difficult to interpret within such a wide parameter space, because it is also sensitive to the thickness of the crust D cr (see Fig. 8f), which in turn depends on the early thermal history of the planet that may contain non-linear feedbacks with η 0 and E * . However, the ratio of T cr and D cr correlates well with the surface heat flux (Fig. 8g), which essentially follows a trend comparable to that predicted for temperature, and can be explained in the same way.

As previously discussed (Fig. 6), a fraction of the basal layer may be molten. Fig. 8h shows the thickness of the molten part of the basal layer, which has similar sensitivities to η 0 and E * as those of the predicted temperatures of the basal layer. The presence of a hot molten layer located at the top of the CMB affects the planet's reaction to tidal forcing and internal dissipation [START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF]. To quantify these effects, we computed the degree-two tidal Love number (k 2 ) and associated quality factor Q at Phobos semi-diurnal tidal frequency (5h33min), following the approach outlined in supporting information S7. These two quantities are displayed in Figure 8i and Figure 8j, respectively. The Love number is mostly sensitive to the core radius and to the rigidity of the mantle, which is why it follows closely the trend shown by the thickness of the molten layer displayed in Fig. 8h. In contrast, the tidal quality factor is essentially sensitive to the thermal state of the solid part of the planet, in particular in regions where its ability to deform is greater. In the present case, this corresponds to the hottest part of the solid-to-partially-molten mantle. The latter is defined as the mantle with a melt fraction smaller than φ c = 0.4, the critical value for the rheological transition from solid to liquid behavior of silicates [START_REF] Lejeune | Rheology of crystal-bearing silicate melts: An experimental study at high viscosities[END_REF][START_REF] Costa | A model for the rheology of particle-bearing suspensions and partially molten rocks[END_REF]. The trend observed in Figure 8j is consistent with Q scaling as the inverse of the viscosity in the hottest part of the solid denser layer, where most of the dissipation takes place. Therefore, the dependencies of the main present-day quantities displayed in Figure 8 are essentially governed by the rheology of the convecting mantle combined with the heat buffer effect of the enriched layer.

To better quantify the influence of the basal layer we define the relative difference function, ∆, such that:

∆(x) = x layered -x homogeneous x homogeneous , ( 23
)
where x is a given quantity associated with a given evolution (e.g., T m , k 2 , D cr ...). The subscript 'homogeneous' refers to a case without basal layer while the subscript 'layered' refers to the same case with a basal layer. Recall that in what follows, the latter always corresponds to a 200-km thick stably stratified layer with an initial layer temperature T d0 = T m0 (i.e., initially no temperature contrast between the basal layer and the overlying mantle). Large absolute values for ∆ indicate large differences with respect to the non-layered case, and vice versa.

Figure 9 displays ∆ field values for several key quantities after 4.5 Gyr of evolution, as a function of η 0 and E * . The first row (a-d) shows the quantities that are most significantly affected by the presence of the basal layer over a wider parameter space. The second row (e-h) shows the quantities that are less affected by the presence of the basal layer (or only affected over a relatively small region within the parameter space considered). Among the most affected quantities, the surface heat flux is reduced by 15% to 30%, even though present-day mantle temperatures are not necessarily very different compared to the case without layering (Fig. 7b). These differences mostly reflect the distinct thermo-chemical evolutions between the homogeneous and the layered cases. Systematic differences are also predicted for the field ∆(D cr ) (Fig. 9d), indicating a strong influence of the basal layer. However, as explained earlier, the present-day crustal thickness mostly results from the early history, which cannot be distinguished in present-day temperature fields. As seen above, the basal layer acts as an insulator and heat source for the core, resulting in a significant increase (∼50 to ∼70 percents) in T c , as illustrated in Figure 9b. Similar to the surface heat flux, this magnitude of ∆(T c ) is mostly sensitive to η 0 . In addition, the presence of a thick molten layer (Figure 8h) strongly affects the Love number k 2 (Figure 8c) resulting in a 30% to 50% increase. Other quantities, such as the planet average temperature (Fig. 9e), or the temperature of the convecting mantle (Fig. 9f), are less affected by the presence of the basal layer, because the core volume, hence its contribution to the average thermal state, is rather small, and because of the thermostat effect (i.e., the non-linearity induced by the temperature dependence of viscosity, which tends to attenuate thermal differences with time [START_REF] Schubert | Subsolidus convective cooling histories of terrestrial planets[END_REF]).

The presence of the basal layer also results in a reduction of the tidal quality factor Q. The influence is rather moderate compared to other quantities displayed in Fig. 9) (usually <30%, Fig. 9g), although it can be larger for specific combinations of small η 0 and large E * . This corresponds roughly to the region of the parameter space where the differences in mantle convective temperatures are the largest (Fig. 9f), which probably relates to the sensitivity of dissipation to mantle temperature mentioned above. In addition, this region also coincides with the region where the layer-induced decrease in lithospheric thickness are the largest (Fig. 9h), perhaps indicating a change in efficiency of stagnant-lid convection.

Consequently, the presence of a basal layer combined with the influence of mantle rheological parameters results essentially in a strong decrease in present-day surface heat flux and crustal thickness, associated with a significant increase in core temperature as well as k 2 . The impact of the presence of a basal layer on the tidal dissipation is moderate to large, and more modest for the overlying convecting mantle temperatures, due to the thermostat effect.

Implications

Our results have important implications for the long-term thermo-chemical evolution of Mars. For example, primordial heat can be stored efficiently in the deep interior of Mars if a stably stratified layer insulates the core, while most of the Martian mantle undergoes efficient cooling. This prediction may account for petrological evidence, which points to a protracted thermal history of Mars mantle [START_REF] Filiberto | Constraints on the depth and thermal vigor of melting in the martian mantle[END_REF]. Mantle source temperatures of igneous rocks estimated by [START_REF] Filiberto | Constraints on the depth and thermal vigor of melting in the martian mantle[END_REF] indicate the presence of a long-lived hot reservoir somewhere in the Martian mantle, but also of another reservoir that has steadily cooled with time, consistent with our model predictions. The effects of a denser basal layer on the long-term evolution moreover affect the interpretation of different types of geophysical data in terms of the present-day structure of Mars, including seismic wave arrival-times and geodetic measurements (k 2 , Q, moment of inertia factor). Indeed, as illustrated in Figures 9a,d andh, the presence of a basal layer leads to a dramatic change of the lithospheric and crustal thicknesses, as well as of the shallow thermal gradient. These quantities are expected to strongly affect the propagation of seismic waves that can be recorded by the InSight SEIS instruments. For instance, the shallow thermo-chemical profile of Mars can affect the presence of seismic shadow zones, and modulate their extent (e.g., [START_REF] Zheng | Seismological implications of a lithospheric low seismic velocity zone in Mars[END_REF], their Fig. 2).

One of the most dramatic influences of the basal layer is the presence of a molten silicate layer above the CMB. This prediction could influence the interpretation of seismic data in terms of core size. Indeed, the presence of a molten silicate layer could give the false impression of the CMB location to be shifted towards shallower depths if core phases are recorded by SEIS. This effect is illustrated in Figure 10 that displays seismic velocity profiles for P-and S-waves, along with the associated ray paths for deep reflected waves for a source located at 50 km depth, and an epicentral distance of 60 degrees. The case without (Fig. 10a-b) and with (Fig. 10c-d) a 300-km thick basal layer are considered. The thermal profiles for these two cases are displayed in Figure 7b and the ray paths were computed using the Tau p toolkit [START_REF] Crotwell | The TauP Toolkit: Flexible Seismic Travel-time and Ray-path Utilities[END_REF]. In absence of a basal layer the P-and S-wave seismic velocity profiles have comparable depth dependence in the entire mantle (Fig. 10a) leading to similar ray paths for compressional and shear waves in this region (Fig. 10b). This yields comparable reflections at the CMB for P-and S-waves. However, the presence of a molten basal layer changes P and S wavespeeds, V P = (K + 4µ/3)/ρ and V S = µ/ρ, in the deep mantle in two different ways. Since most of the layer is molten, the corresponding shear modulus µ (and hence V S ) is zero. The bulk modulus K also decreases due to melting and large temperatures, but its value remains significantly above zero (Fig. 10c). Consequently, the P and S ray paths become significantly different: S-waves are reflected at the interface (in this case located close to R c + D d ) where the mantle is molten, while P-waves can travel further towards deeper regions in the basal layer (however at considerably smaller wave speeds), and are reflected at the CMB (Fig. 10d). These different ray paths would lead to considerably smaller values of travel time difference between P-and S-waves, ∆t S-P , compared to the case where the basal layer is absent. Indeed, compared to the homogeneous mantle case the larger P-wave distances together with the smaller P-wavespeed in the molten mantle would yield larger P-waves travel times. On the other hand the smaller S-waves distance in the layered mantle case would contribute to the decrease of the travel time difference between P-and S-waves. Specifically, in the absence of a basal layer ∆t S-P = 474 s, while ∆t S-P = 287 s when the basal layer in present. This corresponds to a significant (i.e., ∼40%) relative difference in ∆t S-P between the homogeneous and the layered mantle cases. In the current example, the detection of a deep bouncing S-wave could lead to an erroneous interpretation of a core radius of ∼1940 km instead of 1700 km. Therefore, identifying such an 'anomalous'/unexpected 'S -P ' travel time difference would allow one to discriminate between the presence or the absence of a compositional stratification in the mantle of Mars.

An analogous influence of the presence of a basal layer can be expected on the extent of shadow zones due to deep mantle and core structures. In the absence of a basal layer, the angular extent of the shadow zone due to the liquid core is larger for S-waves than for P-waves. For a given core size, if a molten basal layer is present, the angular extent of the S-waves shadow zone would be even larger than in the absence of a basal layer, because shear waves are unable to propagate in the deep molten mantle. On the contrary, the P-wave shadow zone due to the deep regions of the planet would be less affected since compressional waves can still propagate in the molten mantle.

Additionally, the interpretation of geodetic data could be influenced by the possibility of a deep molten silicate layer in terms of core size and composition. This is shown in Figure 11, which illustrates the possible tradeoffs between Martian core size and the thickness of a basal layer. This figure was obtained by considering a set of prescribed temperature profiles sharing common characteristics, but also with distinct features. All the profiles have D cr =60 km and D l =300 km thick crust and lithosphere, respectively, with T m = 1600 K, T c = 2600 K, and T d = 2500 K. The denser layer is therefore completely molten, while the overlying mantle is entirely solid. The two main parameters that remain distinct among these models are the core radius R c and the molten layer thickness D d . These two parameters were varied systematically within 0-400 km for D d , and 1500-1850 km for R c . Following the approach detailed in [START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF] with the knowledge of the profiles, we computed the mantle elastic coefficients, density, and rigidity required to compute the corresponding degree-two Love number. The latter is obtained following the approach outlined in the supporting information S7. For the core, which is assumed to be composed of iron with a fraction x S of sulfur, we followed the approach outlined in [START_REF] Rivoldini | Geodesy constraints on the interior structure and composition of Mars[END_REF] to compute its physical properties. The obtained k 2 field shown in Figure 11a clearly illustrates the tradeoff between R c and D d : any given value of k 2 can be explained equally well by either a small core overlaid by a relatively thick molten silicate layer, or a larger core with a thinner or even no overlying molten silicate layer. For k 2 = 0.169 ± 0.012 (2σ) [START_REF] Konopliv | An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data[END_REF] (marked by the black lines in Fig. 11a), a significant range of compatible solutions in (D d , R d ) space exist. Such tradeoffs can be partially removed by considering additional constraints, such as the normalized moment of inertia factor, whose value in the 2σ range I/(M R 2 p ) = 0.36379 ± 0.0002, is known with a considerably better accuracy than that of the Love number [START_REF] Konopliv | An improved JPL Mars gravity field and orientation from Mars orbiter and lander tracking data[END_REF]. Figure 11b displays the normalized moment of inertia factor associated with R c and D d , assuming a core sulfur content x S = 0.11. The estimated 2σ range (black contours in Fig. 11b) yields a different and smaller set of compatible solutions. The k 2 -inferred and I-inferred sets of compatible solutions can be combined by considering their intersection in (D d , R d ) space to yield a considerably smaller set of solutions (Fig. 11d). However, the resulting set depends on to the assumed core composition, which was fixed to a sulfur content of 11%. If one instead adjusts the core sulfur content in order to match I/(M R 2 p ) = 0.36379, one can express the dependence of x S on R c and D d , as displayed in Figure 11c. The core sulfur content is more sensitive to R c than D d , which allows constraining the relationships between k 2 , I, x S , R c and D d . For instance, Figures 11d-f show the set of R c and D d values that satisfy both k 2 and I constraints for different values of core sulfur content. Therefore, for a given mantle composition, combining k 2 and I would allow one to constrain the value of R c , D d and the composition of the Martian mantle for a plausible content of light elements in the core. Additional constraints on the present-day Mars tidal quality factor (e.g., Q = 95 ± 10 ( [START_REF] Khan | A Geophysical Perspective on the Bulk Composition of Mars[END_REF]) may be considered. However, this constraint may be too loose to yield any improvement. In addition, Q likely depends on grain size, whose present-day value for Mars is not well known (supporting information S7). Nevertheless, the presence and the persistence of a denser, hotter and molten silicate layer overlying the Martian core will significantly affect the tidal dissipation of Mars. Over geological time scales, this could significantly influence the orbital evolution of Martian satellites, in particular Phobos, and could modulate the nature of the thermal-orbital constraints on the mantle rheology and the initial thermal state of Mars [START_REF] Samuel | The rheology and thermal history of mars revealed by the orbital evolution of Phobos[END_REF]. This aspect is worth investigating in the near future.

In this study we assumed that the basal layer results from fractional crystallization. We explored the influence of another style of crystallization : the intermediatebatch crystallization as in Ballmer et al. (2017). In this crystallization scenario the fractionation is due to compaction of the crystal mush at ∼ 50% melt fraction (instead of due to crystal settling at ∼ 100% melt fraction). We found that it does not significantly affect the model results.

Our modelling results rely on the assumption of a strong compositional stratification of the mantle, as discussed earlier and shown in Fig. S1. This implies that the interface between the basal layer and the overlying mantle is essentially flat and that no mixing occurs between the layers. In this case, the parameterized convection can reproduce accurately the evolution in curved geometry (supporting information S7). However, alternative scenarios are possible, such as the development of a strong topography at the interface and the partial erosion and mixing between the two mantle layers, which would occur in the case where the compositional density contrasts are smaller than what we considered [START_REF] Davaille | Two-layer thermal convection in miscible viscous fluids[END_REF][START_REF] Samuel | Thermochemical convection and helium concentrations in mantle plumes[END_REF]Tosi, Plesa, & Breuer, 2013;[START_REF] Mcnamara | Thermochemical structures beneath africa and the pacific ocean[END_REF][START_REF] Li | Chemical complexity of hotspots caused by cycling oceanic crust through mantle reservoirs[END_REF]Maurice et al., 2017). This would likely reduce the influence of the basal layer that we described in our study. Therefore, our modeling results for the (heterogeneous) layered and (homogeneous) non-layered cases can be considered as plausible end-members within a broad spectrum of intermediate scenarios.

We showed in the supporting information S6 that due to its iron enrichment, the basal layer remains denser than the overlying mantle, even if it is partially or entirely molten. However, our models do not account for possible melt transport within the basal layer. Indeed, if the layer is partially molten, melt-solid density differences may lead to episodes of melt segregation [START_REF] Boukaré | Modeling phase separation and phase change for magma ocean solidification dynamics[END_REF] and/or Rayleigh-Taylor overturns. Because both the enriched melt and solids are denser than the overlying mantle these processes would remain confined within the layer. If they occur, these episodes may result in smaller temperature gradients across the basal layer, but the layer's heat buffer effect we reported in this study would remain present. Hence, even in this case our main conclusions would not be altered by the occurrence of these complexities.

In this study, we considered the HPE abundances associated with the bulk composition model of [START_REF] Wänke | Chemistry and Accretion History of Mars[END_REF] while for other bulk composition models the HPE content can be significantly different [START_REF] Sanloup | A simple chondritic model of Mars[END_REF][START_REF] Treiman | Core formation in the Earth and Shergottite Parent Body (SPB): Chemical evidence from basalts[END_REF][START_REF] Lodders | An oxygen isotope model for the composition of Mars[END_REF]. Different HPE contents will affect the thermal evolution and the resulting present-day crustal thickness, surface heat flux, and deep thermal structure (Plesa et al., 2015). However, these changes are unlikely to qualitatively change the influence of the basal layer and our general conclusions.

Our main results did not account for the presence of water, which has been reported to be potentially important for the thermochemical history [START_REF] Kiefer | Water undersaturated mantle plume volcanism on present-day Mars[END_REF]. We have therefore considered cases where the influence of water on the thermochemical evolution of the planet is accounted for (supporting information S10). Water reduces mantle viscosity, which enhances heat transfer. Water also depresses the melting curves, which enhances melting. As a consequence, the presence of water enhances early cooling and crustal production, but the thermostat effect attenuates the long-term differences in thermal histories between the cases with and the cases without water. In addition, the effect of water is similar for cases with and without basal layer. Therefore, the influence of the basal layer in presence of water remains similar to what we observed when water is not present.

A possible limitation of our approach results from the assumed constant value of the thermal mantle conductivity, in particular in the basal layer. We evaluated the possible impact of this simplification in the supporting information S9 and found that larger values of thermal conductivity (possibly due to the larger temperatures in the basal layer (Hofmeister, 1999;[START_REF] Schumacher | Influence of a variable thermal conductivity on the thermochemical evolution of Mars[END_REF]) could notably reduce the temperature contrasts between the basal layer and the overlying mantle and generate thicker crusts. However, the iron enrichment may reduce the thermal conductivity of the basal layer (Zhang et al., 2019). Therefore, large uncertainties remain regarding the differences in thermal conductivity between the enriched basal layer and the overlying mantle. In spite of these uncertainties, the influence of the basal layer on the thermo-chemical evolution of Mars remains qualitatively comparable, and therefore does not affect the main conclusions of our study.

As pointed out in several previous studies [START_REF] Nakagawa | Effects of thermo-chemical mantle convection on the thermal evolution of the Earth's core[END_REF]Zeff & Williams, 2019), the insulating nature of the basal layer may prevent the sustainability of an early Martian dynamo suggested by magnetic data [START_REF] Acuña | Magnetic field and plasma observations at Mars: Initial results of the Mars global surveyor mission[END_REF][START_REF] Connerney | Mars crustal magnetism[END_REF][START_REF] Johnson | Crustal and time-varying magnetic fields at the InSight landing site on Mars[END_REF]. Indeed, unlike the homogeneous mantle case [START_REF] Breuer | Early plate tectonics versus single-plate tectonics on Mars: Evidence from magnetic field history and crust evolution[END_REF], the presence of a basal layer requires external sources in addition to favorable endogenous processes to power an early dynamo: (1) Soon after the emplacement of the basal layer, the initial super-heating of the core due to core formation processes could have been substantial [START_REF] Senshu | Thermal evolution of a growing Mars[END_REF][START_REF] Samuel | Heat partitioning during core formation by negative diapirism in terrestrial planets[END_REF][START_REF] Rubie | Formation or Earth's core[END_REF], and would have been sufficient to generate a strong CMB heat flux during a few tens of millions of years; (2) The overturn that led to the formation of the basal layer in the first place should have enhanced CMB heat flux by delivering cold material to the CMB (Elkins Tanton et al., 2005;[START_REF] Plesa | Can a fractionally crystallized magma ocean explain the thermo-chemical evolution of Mars? Earth Planet[END_REF], in particular if plate tectonics was operating on early Mars [START_REF] Breuer | Early plate tectonics versus single-plate tectonics on Mars: Evidence from magnetic field history and crust evolution[END_REF]; (3) Late giant impacts may have led to early dynamo episodes lasting for a few tens to a few hundreds of millions of years [START_REF] Reese | Early martian dynamo generation due to giant impacts[END_REF][START_REF] Monteux | Dynamics of core merging after a mega-impact with applications to Mars' early dynamo[END_REF] even though the heat associated with the impact would further reduce the heat flux at the CMB (Arkani-Hamed & Olson, 2010); (4) Elliptical instabilities in the Martian core lasting for several hundreds million years could be exited by early satellites orbiting the planet in retrograde fashion [START_REF] Arkani-Hamed | Did tidal deformation power the core dynamo of Mars?[END_REF]. This process has been shown to be sufficient to trigger and maintain an early dynamo for the first ∼ 500 Myrs of Mars evolution [START_REF] Arkani-Hamed | Did tidal deformation power the core dynamo of Mars?[END_REF][START_REF] Sauret | Tide-driven shear instability in planetary liquid cores[END_REF]. These mechanisms of internal or external origin, alone or in combination, would have led to the existence of an early Martian dynamo operating on the scale of a few hundred million years. Therefore, one can reasonably consider that the presence of deep mantle layering remains compatible with the existence of an ancient Martian dynamo.

Conclusions

We studied the influence of the presence of an iron-and HPE-enriched layer at the bottom of the Martian mantle, on the long-term thermo-chemical evolution of the entire planet. The presence of such a layer is a likely consequence of the solidification of an early silicate magma ocean that Mars most probably experienced. We conducted a restricted set of dynamic evolution calculations, and a broader systematic exploration using parameterized convection calculations to model the thermo-chemical evolution of Mars with or without the presence of a denser and enriched silicate layer overlying the Martian core. We focused on the most likely possibility of a non-convecting basal layer characterized by motionless, purely conductive heat transfer. The presence of the basal layer strongly influences long-term planetary thermal evolution, yielding a considerable temperature increase in the lower mantle and in the core. In turn, the rest of the mantle cools down more efficiently, affecting the crustal thickness and surface heat flux. The significant temperature increase due to the HPE enrichment in the denser basal layer systematically generates large amounts of melting in the enriched silicate regions. Consequently, the deep mantle enrichment in iron and HPE implies the presence of a molten silicate layer above the CMB. This hot and molten silicate material significantly increases the planet's Love number, and increases its tidal dissipation. These drastic changes induced by the deep Martian mantle layering are likely to alter our interpretation of seismic, geodetic and heat flux data that the Insight mission has started to collect since the deployment of its instruments. Therefore, the possibility of the presence of an enriched basal layer should be considered when interpreting available and upcoming geophysical data. In addition to constraining the present-day structure of Mars, this consideration will allow one to further connect the present-day structure of the planet with its early state, and to reconstruct its long-term thermo-chemical history. Mantle specific heat at constant pressure 1142.0

J kg -1 K -1 C pcr Crust specific heat at constant pressure 1000 J kg -1 K -1 C pc Core specific heat at constant pressure 840 J kg -1 K -1 k m Mantle thermal conductivity 4 W m -1 K -1 k cr Crust thermal conductivity 2.5 W m -1 K -1 α
Silicates thermal expansion coefficient 2 10 The resulting cumulate compositional profiles after MO freezing are gravitationally unstable (Figure S1a), promoting mantle-scale overturns. Overturns typically result in stable mantle stratification, and we consider that the bottom part of the post-overturn mantle (thick solid lines in Fig. S1a) turns into a denser basal layer, while the upper part is efficiently mixed. Depending on the value of the layer thickness, D d (see below for assumptions), and the initial Fe# of the MO (i.e., the Fe# of bulk-silicate Mars, Fe# BSM ), the average Fe# of the basal layer Fe# d is computed. We chose Fe# BSM such that the the average Fe# of the mantle above the basal layer is always Fe# U L . We consider Fe# U L = 25 (precisely 24.9884 (Bertka & Fei, 1997)), and as a lower bound, Fe# U L = 20. For both these end-members, Fe# BSM slightly increases with increasing D d . In any case, Fe# BSM ≤34.5 for any value of D d considered in this study (and even ≤29 for Fe# U L =20).

The average density of the basal layer is computed from the average Fe# of the layer Fe# d :

ρ d = ρ 0 + Fe# d (ρ Fe -ρ Mg ) (S3)
with the densities of FeO-rich and MgO-rich end-members, ρ Fe = 4192.2 kg/m3 and ρ Mg = 3266 kg/m 3 , respectively. These end-member rock densities are taken from Elkins Tanton et al. (2003), assuming an aggregate of 50% olivine and 40% clinopyroxene and 10% garnet as a rock analog (see their Table 4). ρ 0 is chosen such that the average density of the Martian mantle above the basal layer is 3500 kg/m 3 . For Fe# U L = 20, ρ 0 = 3314.76 kg/m 3 . For Fe# U L = 25, ρ 0 = 3268.56 kg/m 3 . The relevant density anomaly ∆ρ d is the difference between ρ d and the density of the rest of the mantle. ∆ρ d strongly depends on D d , as being controlled by Fe# d (Figure S2). Previous studies demonstrated that D d depends on the length-scale of overturns (i.e., small-scale incremental overturns vs. global-scale mega-overturn) (Maurice et al., 2017;Ballmer et al., 2017;Boukaré et al., 2018). Smaller-scale overturns tend to promote efficient mixing of a larger part of the mantle (i.e., above the deep layer), hence sustaining smaller D d . Since the scale of overturns during MO freezing, and thus D d , remain highly uncertain, being controlled by poorly-constrained rheological parameters in the mushy cumulate package and by the MO cooling history (Maurice et al., 2017;Ballmer et al., 2017;Boukaré et al., 2018), we varied D d as a free parameter in our study. Owing to convective mixing during and after gravitational overturn(s), we also assumed that the Fe# profile across the basal layer is gradual (dashed lines in Fig. S1a), consistent with numerical model predictions (Maurice et al., 2017;Ballmer et al., 2017), and that the overlying mantle is homogeneous in composition. In addition to FeO, heat-producing elements are enriched in the basal layer due to their highly incompatible behavior. We assume that all HPE remain in the MO during freezing, and hence eventually enter the deep layer, except those that are trapped in the solid as an interstitial porosity while the MO is freezing. We consider a conservative value of interstitial melt porosity of φ i = 5% (for example, the interstitial porosity beneath mid-ocean ridges on Earth is about 1% (Toomey et al., 1998)). During MO freezing, the concentration of HPE in the interstitial droplets, C i , evolves from that of the bulk Mars, C m , (during formation of first cumulates) to that of the denser layer, C d , (last cumulates). Assuming that this evolution can be approximated as a linear compositional trend due to the highly incompatible nature of HPE, the average concentration of HPE in the interstitial droplets is (C m + C d )/2. Thus, the relevant concentration of HPE in the (well-mixed) mantle above the deep layer is: ) with:

C m = φ i (C m + C d )/2, ( S4
C d = C m V m -φ i (V m -V d )/2 V d + φ i (V m -V d )/2 (S5)
by mass balance. V m and V d are the volumes of the entire mantle of Mars, and of the deep basal layer, respectively. We define the HPE layer enrichment as the ratio between C d and C m , which therefore writes:

Λ d = H d H m = C d C m = V m -φ i (V m -V d )/2 V d + φ i (V m -V d )/2 . ( S6
)
This quantity corresponds to the ratio of the basal layer to bulk mantle HPE heat productions.

S2 Finite-Volume modelling

To model the dynamic evolutions displayed in Figures 2 and Figures

DT Dt = ∇ 2 T + H, (S10)
In the above equations, p is the dimensionless dynamic pressure, η is the dimensionless dynamic viscosity, T is the dimensionless temperature, e r is a radial unit vector pointing upward, B is the buoyancy number, Ra is the thermal Rayleigh number, H is the C-dependent internal heating defined in Section 3. DX/Dt = ∂X/∂t + u • ∇X expresses the total/Lagrangian derivative of the scalar field X (i.e., T or C).

At low temperatures, the Arrhenius law (Eq. ( 8)) yields extremely large viscosity values that are difficult to numerically handle. We thus imposed a maximum viscosity contrast (∆η max ) of 10 10 , which is sufficiently large to ensure that a stagnant lid regime is established (Guerrero et al., 2018). While the minimum viscosity (η min ) obtained from Eq. ( 8) is free to vary according to the local conditions of temperature and pressure, the maximum viscosity (η max ) is limited at each time step to η max = η min ∆η max , which ensures that the viscosity obtained in the convecting part of the domain is not affected by artificial cutoffs.

We solved Eqs. (S7)-(S10) in a 2D half cylinder. All boundaries are free-slip. Lateral side-walls are thermally and compositionally insulating, while the upper and lower surfaces have imposed temperature and zero compositional flux. The domain was discretized using a uniform grid consisting of 122 equally-spaced radial shells with 318 grid cells for each shell. Resolution tests confirmed that this grid spacing was sufficient to yield accurate results. In order to obtain temperature distributions similar to those corresponding to a fully spherical geometry, we rescaled the core radius following the approach of Van Keken (2001) in such a way that the ratio of inner-to-outer surface in 2D is the same as in 3D. In formulas:

R c R p 2 = R c R p and R p -R c = 1, (S11)
where R c and R p are the true CMB and outer radius of the planet and R p and R c are the corresponding non-dimensional radii of the cylindrical shell (0.33 and 1.33, respectively with R p = 3389.5 km and R c = 1700 km).

The above equations are implemented in the computer code Gaia (Hüttig et al., 2013), which uses a co-located, primitive variable formulation. A three-level implicit scheme (Harder & Hansen, 2005) is used to solve the energy equation. The advection of the composition is treated with a particle-based approach (Plesa et al., 2012) using 40 particles per grid cell. Gaia has been extensively benchmarked for thermal and thermo-chemical convection against various numerical solutions in different domain geometries (Plesa et al., 2012;Hüttig et al., 2013;Tosi et al., 2015).

S3 Case of a convective basal layer

We considered the case of an enriched basal layer, whose composition is homogeneous (contrary to the case displayed in Figure 2d-g, where the iron and HPE enrichment increases with depth in the layer). Figure S3 displays the temperature field and several corresponding horizontally averaged quantities after 4.5 Gyr of evolution. Apart from the absence of a radial compositional gradient of the basal layer relative to the overlying mantle (Fig. S3d), all model parameters, boundary and initial conditions are identical to those corresponding to the stably stratified case displayed in Figure 2d-g and described in Section 3. As in the stably stratified case, the compositional density contrast between the basal layer and the overlying mantle prevents the complete mixing between the two layers. The non-zero values of the compositional field in the mantle above the layer (Fig. S3c) indicate that a limited amount of mixing between the two layers has occurred. The absence of a compositional gradient allows the development of convective motions within the basal layer. This results in a more homogeneous temperature profile than what can be observed for the stably stratified case (compare Fig. 2e and Fig. S3b). One can observe an intermediate thermal boundary layer located at the interface between the basal layer and the overlying mantle. Generally, the obtained thermal state is similar (but less hot) to that of the stably stratified case, with a hotter basal layer (whose average temperature is close to 2150 K while the temperature of the overlying convecting mantle is close to 1750 K) due to the HPE enrichment, and the fact that the intermediate thermal boundary layer at the interface reduces the heat transfer between the enriched basal layer and the more depleted overlying mantle. The value of the associated buoyancy number is B d = 3.8.

While the observed thermal evolution corresponds to solid-state convection one can expect a similar behavior if the basal layer is hot enough to be entirely molten. In this case, experiments and theory have shown that the dependence of heat transfer on the layer's Rayleigh number follows the same power-law dependence observed for the solid-state case (infinite Prandtl number limit) (Grossmann & Lohse, 2000;Ahlers et al., 2009;Stevens et al., 2013). Convective motions in the basal layer are also possible in this case because of the relatively low value of mantle viscosity (η 0 ) considered. For larger values of η 0 the basal layer is less prone to convection and will become purely all other constant quantities are listed in Table 1, and viscosity η d is computed using Equation ( 8), with V * = 5 cm 3 /mol and a temperature of 2150 K (Fig S3a -b). Due to the higher heat transfer efficiency of convective motions compared to conduction a convecting basal layer is expected to be colder than a purely conductive layer (Fig S3ab). Figures S4a-c show that a basal layer will most likely be sub-critical, except for a rather extreme combinations of low reference viscosities and large layer thicknesses (left side of Fig. S4c). Therefore, denser convective basal layers throughout the planet's history require rather large layer thicknesses and/or rather extreme temperatures that would have to be maintained for billions of years. While such combinations are not completely impossible, they appear to cover a very small fraction of the plausible parameter space.

In summary, an enriched basal layer of homogeneous composition is unlikely to be convective for the ranges of layer thicknesses considered. Consequently, conductive heat transfer is the most likely heat transfer mechanism within the basal layer, regardless of its compositional stratification.

S5 Depth-dependent enrichment factor

To account for a linear depth-dependent enrichment in the stratified basal layer case, we proceed as follows. We first express the radial dependence of the iron-number Fe# =100 Fe/(Fe+Mg) within the layer, assuming a linear increase with radius r (or alternatively with the depth beneath the interface:

h d = R d -r, where R d = R c + D d ): Fe# d (r) = Fe# di + (R d -r) dFe# d dh d , (S12)
where Fe# di = Fe# d (r = R d ) is the Fe-number at the top of the basal layer, whose volume-averaged Fe-number is

Fe# d = 1 V d V d Fe# d dV = 4π V d R d Rc Fe# d r 2 dr, (S13)
and ) is obtained by inserting Equation (S12) into Equation (S13).

dFe# d dh d = 4(R 3 d -R 3 c )(Fe# d -Fe# di ) R 4 d -R 3 c (4R d -3R c ) . ( S14
The value of Fe# di is chosen such that the associated compositional density contrast at the top of the basal layer is at least 0.11ρ m , which is sufficiently large to prevent significant erosion of the layer via entrainment into the overlying mantle by thermal plumes originating at the interface, and to keep the interface essentially flat (Fig. 2h). Then, with the knowledge of Fe# d (r) and Fe# d one can express the radius-dependent enrichment factor within the basal layer, being the ratio of these two quantities:

f e (r) = Fe# d (r) Fe# d if R d ≥ r > R c 0 if r > R d (S15)
The above radius-dependent enrichment factor accounts for iron enrichment and therefore also for the compositional density increase with depth in the basal layer (Eq. ( 22a)). We use the same dependence to account for depth-dependence of HPEs within the layer (Eq. ( 22b)).

S6 The stability of silicate liquids at the base of the Martian mantle

To assess the gravitational stability of a molten basal layer in the lowermost mantle of Mars, we estimated the density difference between silicate melts and solids at the relevant pressures and compositions. To calculate liquid densities, we considered the equation of state (EOS) by Suzuki et al. (1998), which is based on experimental constraints and valid for a similar (i.e., Mars-like) bulk composition (Morgan & Anders, 1979) as the one considered here (Bertka & Fei, 1997). Unfortunately, the experiments cannot constrain the K term (i.e., the pressure derivative of the bulk modulus) in the EOS, so we considered this uncertainty here within the thermodynamically allowed bounds of K =6 and K =10.1 (Suzuki et al., 1998). Finally, we compared these EOS to solid density profiles that are calculated for the same bulk composition (MA (Morgan & Anders, 1979)), each along a 1706.56 K isentrope, which corresponds to a CMB temperature of 2023.15 K. The solid density profiles and the isentropes for the MA composition have been calculated with Perple X (Connolly, 2005) using the thermodynamic formulation and mineralogy database of Stixrude and Lithgow-Bertelloni (2011). This approach provides us with the relevant bounds for the melt-solid density difference ∆ρ for Mars-like mantle compositions. In addition, we computed the effect of FeO-enrichment of the liquid on ∆ρ, considering the effects of the iron content on the EOS of liquids according to Jing and Karato (2009).

If the composition of the basal molten layer equals that of the overlying Martian mantle (Morgan & Anders, 1979)), ∆ρ is negative, regardless of the choice for K (black lines in Figure S5a). This implies that a molten mantle basal layer with the same composition as the solid mantle cannot remain stable at the base of the Martian mantle, as previously suggested (e.g., (Boukaré et al., 2015) and references therein). However, in the case where the basal molten layer is enriched in iron, it can be indeed stabilized in the deep Martian mantle. Within the uncertainty bounds for K (dashed and solid lines in Figure S5a), we find that ∆ρ for liquid compositions of Fe#=60 is positive over a wide pressure range (6<P<23 GPa), allowing the molten basal layer material to remain gravitationally stable (blue lines in Figure S5a). For Fe#=67, ∆ρ is clearly positive over the considered pressure range, by at least 60 kg/m 3 (red lines in Figure S5a). For iron contents of the mantle above the basal layer Fe# U L = 25 (Bertka & Fei, 1997), the average iron number of the basal layer is Fe#≥67 for D d ≤ 950 km. For our lower-bound Fe# U L = 20, Fe#≥67 for D d ≤ 690 km. Thus, even thick enriched molten basal layers (beyond the range we considered) can be well stabilized in the deep Martian mantle for these iron contents. The green lines in Figure S5b show the explicit ∆ρ profile for the case of D d = 500 km, i.e., the case corresponding to the thickest layer we considered for Fe# U L = 25 (Section 3), Fig. 1, and Fig S1). This case is compositionally stratified, with a minimum Fe# at the top of the layer of 71.8 (red dashed line in Figure S1a). If such a layer is liquid, it would always remain denser than the overlying solid with bulk Mars composition. For layers with D d < 500 km the density difference would be even larger, because the corresponding iron enrichment would be even more important than that of a 500-km thick basal layer (Figure S1). This conclusion on the stability of the basal molten layer holds for Fe# U L = 20 (orange lines in Figure S5b). These results are in good agreement with those inferred with another approach based on a solid-liquid thermodynamic modeling in the MgO-FeO-SiO 2 system ( (Boukaré et al., 2015) and references therein), which apply to the higher-end pressure range considered in our study (20)(21)(22).

∆ρ (kg/m³)

Fe#=67 Fe#=60 Fe#=22.9 S1a) and in Figure S1b, respectively.

S7 Computation of the degree-two Love number and tidal dissipation

The temperature and pressure profiles, along with the activation volume and energy, are used to compute the complex compliance, in the frame of the Andrade pseudoperiod model (Jackson & Faul, 2010;Padovan et al., 2014). First, the radiallydependent complex compliance is calculated with (subscripts R and I indicate the real and imaginary part, respectively):

J R = J U 1 + β * Γ (1 + n) ω -n a cos nπ 2 , ( S16
)

J I = J U β * Γ (1 + n) ω -n a sin nπ 2 + 1 ω v τ M , (S17) (S18)
where J U is the unrelaxed compliance (the inverse of the unrelaxed shear modulus G U ), Γ is the gamma function and n, β, and τ M are parameters related to the coefficients of the Andrade creep function J(t) = J U + βt n + t/η, by β * = β/J U and τ M = ηJ U . The subscripts 'a' and 'v' refer to anelastic and viscous components, respectively. The effects of temperature (T ), pressure (P ), and timescale (T 0 = 2π/ω) of the applied forcing enter the model through the frequency ω a,v , which is related to the pseudoperiod master variable X B = 2π/ω a,v , through

X B = 2π ω d d ref -ma,v exp - E * R 1 T - 1 T ref exp - V * R P T - P ref T ref . (S19)
Thus, ω a and ω v are ω obtained from the above equation with the corresponding grain size exponent, m a , and m v , for anelastic or viscous processes, respectively. The subscript 'ref' refers to the reference conditions, E * is the activation energy, V * the activation volume, and R the gas constant. Values for parameters that are not included in the convection model (e.g., β, n, m a,v , τ M , etc.) are based on the experimental results of (Jackson & Faul, 2010), which used olivine samples with an iron content lower than the Martian mantle. Thus, we assume that the Martian mantle is olivinedominated, and include the effect of the iron only by changing the unrelaxed shear modulus, with an approach similar to Plesa et al. (2018) and Nimmo and Faul (2013). However, we neglected the influence of iron enrichment between the basal layer and the overlying mantle. This simplification, however, does not affect the major conclusion of this work because the change in shear modulus in the basal layer is essentially due to the large temperature increase in this region. From the complex compliance, the shear modulus G and viscosity η at a given depth in the mantle are obtained as:

G = J 2 R + J 2 I -1 2 , (S20) η = 1 ωJ I . (S21)
The radial profiles of rigidity and viscosity are then used to calculate the degreetwo complex Love number k * 2 by solving the equations of motion inside the body with the matrix-propagator technique (Moore & Schubert, 2000;Padovan et al., 2014), to obtain in turn the elastic part k 2 = (k * 2 ) and the dissipation Q = k 2 / (k * 2 ). We assumed a grain size of 1 cm (Nimmo & Faul, 2013). In general, a smaller grain size makes the body more dissipative, with a corresponding k 2 increase and Q decrease (i.e., larger dissipation). However, the choice of the grain size has only a minor effect on the value of k 2 , as illustrated in Figure S6. In contrast, the value of Q is considerably more affected by the grain size. Given that 1 cm is likely an upper bound for the possible grain size of the mantle of Mars, and that for smaller grain sizes Q would decrease (Figure S6), the Q shown in the main text (Section 5 and Figure 8) can be seen as an upper bound value.

S8 Comparison between finite-volume and parameterized convection results for the stably stratified case

Several studies have proven the ability of the parameterized convection approach to reproduce spherical stagnant-lid convection, including when complexities such as pressure-dependent viscosity, adiabatic heating/cooling, crustal enrichment or the presence of time-decaying heterogeneous internal heating (Stamenković et al., 2012;Tosi et al., 2013;Plesa et al., 2015;Thiriet et al., 2018). To further extend this comparison to the stably stratified case displayed in Figure 2, we computed the evolution of a similar case with the parameterized model for a layered mantle described in Section 4.2). This implies the absence of adiabatic heating/cooling, a fixed CMB temperature, no crustal enrichment and radioactive heating without decay. The presence of the stably stratified basal layer does not rely on the use of additional scalings, since the heat transfer within the layer is computed by solving for the heat equation in spherical geometry with depth dependent heat sources (Eq. ( 21)). Figure S7 compares the obtained temperature profiles after 4.5 Gyr of evolution and shows that the parameterized convection model reasonably reproduces the thermal evolution obtained via finite-volume modeling in cylindrical geometry.

S9 Influence of the thermal conductivity of the basal layer

The iron enrichment of the layer together with the hot temperatures associated with this layer could affect the value of its thermal conductivity. This, in turn, could influence the thermo-chemical evolution of the planet. Experimental studies, suggest that the iron enrichment could reduce lattice thermal conductivity (Zhang et al., 2019), while hotter temperatures could enhance radiative conductivity (Hofmeister, 1999;Zhang et al., 2019). The uncertainties are large, because experimental works do not cover the same pressure and temperature ranges relevant to our study and are conducted on mono-crystals instead of mineralogical assemblages (Hsieh et al., 2018). Therefore, we studied the influence of thermal conductivity in the basal layer on our models results. Figure S8 displays the present-day areotherms (a) and crustal thicknesses (b) corresponding to different values of k d . Larger values of conductivity enhance the heat transfer across the basal layer, leading to smaller temperature contrasts between the basal layer and the overlying mantle, and lower core temperatures (Fig. S8a). The increase in heat transfer efficiency, implies hotter shallow temperatures early on, resulting in larger crustal production rates and therefore thicker crusts (Fig. S8b). In these tests k d does not explicitly depend on pressure, iron content or temperature. Therefore, possible feedback between temperature and thermal conductivity is not accounted for. However, the explored range of k d values should still allow to estimate the potential influence on this parameter on the thermo-chemical evolution of Mars with a layered mantle. While the effect displayed is significant, the influence of the layer remains qualitatively similar, and even in the coldest case corresponding to the most extreme value of k d , a significant fraction of the basal layer is above the solidus and liquidus curves, as seen in Fig. 9b and Fig. 8h.

S10 Influence of water

To account for the influence of water on the thermo-chemical evolution of the planet, we followed the approach detailed in Tosi et al. (2017), where water is assumed to behave as an incompatible element (i.e., it is removed from the mantle upon shallow melting). The presence of water reduces both viscosity and melting temperatures that will increase when water is removed during the evolution due to shallow partial melting. Therefore, high water contents in the mantle tend to enhance both mantle cooling and crustal production. This effect is illustrated in Figure S9a-b that compares the evolution for a case displayed in Figure 4 (where the effect of water is not accounted for), and another case, which only differs from the fact that the influence of water on melting curves and viscosity is taken into account. The case with water shows initially a faster cooling (during the first half of the planet's history). Then, cooling becomes less efficient than in the case without water, due to the thermostat effect. Overall, the differences between the two cases are relatively small (∼30-70 K, see Fig. S9a). The influence of water on the evolution of the crustal thickness is more pronounced (Fig. S9b) and leads to relative differences on the order of 30%. Note however that the crustal thickening due to the presence of water can be considerably reduced/compensated when considering different values for the rheological parameters, or slightly colder thermal state. For instance, a 50 K colder mantle would lead to a similar crustal evolution than in the case without water (compare thick plain and thin dotted curves in Fig. S9b). In addition, the influences described above are similar for cases without basal layer. Consequently, the effect of water is unlikely to significantly change the influence of the basal layer on the thermo-chemical evolution of Mars. This is shown in Figure S10 that displays the relative differences with respect to a case where the basal layer is absent, similar to the panels displayed in Figure 9, except that the influence of water is accounted for assuming an initial concentration of water in the mantle of 200 ppm. The trends observed in Figure S10 are similar to those seed in Figure 9. Figure S10. Relative differences with respect to a case where the basal layer is absent, expressed via the operator ∆ (defined in Eq. ( 23)). A 200-km thick stably stratified layer is considered, together with Rc =1700 km and V * = 5 cm 3 /mol, similar to the panels displayed in Figure 9, except that the influence of water is accounted for. See text for further details.

Supporting information references

Figure 1 .Figure 2 .

 12 Figure 1. Results of the magma ocean fractional crystallization process: present-day quantities as a function of the thickness (or volume fraction) of the enriched layer. (a) Deep layer bulk HPE enrichment factor. (b) Radioactive heat production in the bulk layer (orange) and in the depleted overlying mantle (blue). (c) Layer bulk compositional density contrast. The core radius is assumed to be 1700 km. See main text and supporting information S1 for further details.

Figure 3 .Figure 4 .

 34 Figure 3. Example of the thermo-chemical evolution of Mars without compositional layering in the mantle. The mantle rheological parameters are η0 = 10 20 Pa s, E * = 200 kJ mol -1 , and V * = 5 cm 3 /mol. (a) Present-day structure, and (b) areotherm resulting from 4.5 Gyr of evolution. (c) Evolution of CMB temperature (Tc), uppermost convecting mantle temperature (Tm) and average planet temperature. (d) Evolution of crustal thickness (Dcr), total lithospheric thickness (D l + δu), which includes the thickness of the uppermost convecting mantle thermal boundary layer. See text for further details and definition of the symbols.

Figure 5 .Figure 6 .Figure 7 .Figure 8 .

 5678 Figure 5. Time evolution of various (a-d) fluxes, and (e) bulk and convective Urey ratios corresponding to the layered case displayed in Figure 4 (blue curves), and the equivalent homogeneous case displayed in Figure 3 (black curves).

Figure 9 .

 9 Figure 9. Relative differences with respect to a case where the basal layer is absent, expressed via the operator ∆ (defined in Eq. (23)). A 200-km thick stably stratified layer is considered, together with Rc =1700 km and V * = 5 cm 3 /mol. (a-d) Quantities most affected by the presence of the basal layer. (e-h) Quantities, moderately to weakly affected by the presence of the basal layer.

Figure 10 .

 10 Figure10. Influence of mantle stratification on seismic velocities (a,c) and P-and S-wave ray paths (b,d). Top: Homogeneous mantle case (i.e., no basal layer). Bottom: case with a 300-km thick basal layer. The seismic velocity profiles in (a) and (b) correspond to the black and yellow thermal profiles displayed in Fig.7b, respectively. The depth of the quake is 50 km, and the epicentral distance is 60 degrees. See text for further details.

Figure 11 .

 11 Figure 11. Tradeoffs between the basal layer thickness and core size for the (a) value of the degree-two Love number, and (b) the Moment of Inertia (MoI) factor I assuming 11% of sulfur in the planet's core. Black curves delineate the acceptable solutions based on k2 and I estimates. (c) Core sulfur content associated with a fixed value of normalized MoI. It is assumed that the basal layer is entirely molten and that the crust and lithospheric thicknesses are fixed to recent estimates (Konopliv et al., 2016) within a 2-σ range: k2 = 0.169 ± 0.0012 and I/(M R 2 p ) = 0.36379 ± 0.0002, respectively. Mantle rheological parameters are η0 = 10 21 Pa s, E * = 300 kJ mol -1 , V * = 5 cm 3 /mol. (d) Solutions that simultaneously satisfy constraints on k2 and I estimates (in red) or satisfying only either k2 or I estimates (in green) assuming xS =11% of sulfur in the planet's core. (e-f) Same as (d), but for xS =14% and xS =17%, respectively. See text for further details.

Figure S1 .

 S1 Figure S1. Magma-ocean crystallization profiles, and composition of the deep layer. Hypothetical (i.e., no mixing) compositional and density profiles before (thin solid) and after (thick solid lines) overturn for D d =500 km (red) and D d =300 km (blue). The dashed lines display the compositional profile after final mixing during overturn with stable stratification in the basal layer, and a homogeneous overlying mantle (see text for details). (a) Fe#UL = 25. In this case, Fe# d is ∼82 and ∼90.9 for D d =500 km and D d =300 km, respectively. (b) Fe#UL = 20. In this case, Fe# d is ∼74.8 and ∼85.5 for D d =500 km and D d =300 km, respectively.

Figure S2 .

 S2 Figure S2. Fe# d as a function of D d for the range considered in this study. The iron content of the Martian mantle above the basal layer is (green) Fe#UL = 20 (blue) and Fe#UL = 25 (green).

  S3 we solve the following set of dimensionless conservation equations for mass, momentum, composition (C, which is a proxy for the chemical enrichment: C = 0 refers to the regular mantle, 1 ≥ C > 0 corresponds to enriched mantle) and energy, respectively under the Boussinesq approximation:∇ • u = 0,(S7)-∇p + ∇ • η(∇u + ∇u T ) + Ra(B C -T )e r = 0

  Figure S4. Influence of the mantle rheology (η0 and E * ) on the ability of the layer to convect. The fields displayed show the ratio of the thermal Rayleigh number in the basal layer, Ra d , (without compositional gradient within the layer) to the critical values Racr = 10 5 appropriate of variable viscosity in spherical geometry (Yanagisawa et al., 2016) for different values of layer thickness, D d . Warm colors are indicative of convective motions within the basal layer, while cold colors correspond to sub-critical Rayleigh numbers (no convection) in the basal layer. The black lines delineate the boundary between sub-and super-critical regions. The calculation assumes V * = 5 cm 3 /mol. See text for further details.

Figure S5 .

 S5 Figure S5. Melt-solid density difference (∆ρ) for a Martian mantle composition(Morgan & Anders, 1979) along a 1706.56 K isentrope. (a) ∆ρ for liquids with the same composition as (black), and FeO-enriched (colored) compared to the bulk solid composition ofMorgan and Anders (1979). Fe# of the liquid as labeled. The uncertainty bounds of the K term in the silicate liquid EOS are considered (K =6: solid; K =10.1 dashed)(Suzuki et al., 1998). (b) ∆ρ for a liquid basal layer with D d =500 km, and for an iron content of the mantle above the basal layer of Fe#UL = 25 (green) and Fe#UL = 20 (orange). The corresponding compositional profiles are displayed in FigureS1a) and in FigureS1b, respectively.

Figure S6 .

 S6 Figure S6. Effect of the mantle grain size on the viscoelastic response of Mars. For this figure we used a set of 1024 models with radius of the core RC = 1700 km, a stratified enriched layer of 300 km, a temperature of 1800 K, and an activation volume of 5 × 10 -6 m 3 /mol. The models have varying activation energy E * and reference viscosity η0. We calculated the values of k2 and Q considering two different grain sizes, namely 1 cm and 1 mm. The upper panels show the corresponding variation in k2, calculated as ∆k2 = 100 × (k2 mm -k2 cm)/k2 mm, as a function of E * and η0. Values are typically below 2%, and increase with increasing activation energy. The bottom panel shows the absolute variation in Q, calculated as ∆Q = Q2 cm -Q2 mm. It is clear that the grain size has a considerably larger effect on the dissipation.

Figure S7 .

 S7 Figure S7. Comparison between the present-day thermal states obtained with the dynamic modeling in cylindrical geometry displayed in Fig. 2 (solid curve) and the parameterized convection model (dashed curve) after 4.5 Gyr of evolution. See text for further details.

Figure S8 .

 S8 Figure S8. Influence of the basal layer thermal conductivity, k d , on the thermo-chemical evolution. Present-day areotherms and mantle melting curves (a) and crustal thicknesses (b) corresponding to different values of k d . All other model parameters are those corresponding to the case displayed in Figure 4. The thin horizontal dotted line in (a) marks the core-mantle boundary.

Figure S9 .

 S9 Figure S9. Influence of water, on the thermo-chemical evolution. (a) Time evolution of the uppermost mantle temperature (orange) and average planet temperature (black) for the same case displayed in Figure 4. A second case with the same characteristics is displayed (dashed curves) except that the influence of water is accounted for. (b) Time evolution of the corresponding crustal thicknesses without (plain) and with (dashed) the influence of water. A third case shown (thin dotted curve) has the same characteristics as the evolution influenced by the presence of water, except for an initial uppermost mantle temperature 50 K colder. See text for further details.

Table 1 .

 1 Symbols meaning, and values of the main physical parameters used for the modeling of the thermo-chemical evolution of Mars.

	Symbol	Meaning	Value	Unit
	T m0	Initial Mantle temperature	1800	K
	T d0	Initial denser layer temperature	1800	K
	T c0	Initial Core temperature	2100	K
	T s	Surface temperature	220	K
	ρ m	Mantle density	3500	kg m -3
	ρ c	Core density	7200	kg m -3
	ρ cr	Crust density	2900	kg m -3
	ρ d	Denser layer density	Variable (S1)	kg m -3
	g	Surface gravity	3.7	m s -2
	g c	Core surface gravity	3.1	m s -2
	g d	Denser layer surface gravity	Variable	m s -2
	C pm

manuscript submitted to Journal of Geophysical Research: Planets

Acknowledgments

We are grateful to Julian Lowman and two anonymous reviewers for their comments that helped improving the manuscript, and to the editor Laurent Montési and associate editor Scott King for their useful comments and for the handling of the manuscript. We thank Francis Nimmo for discussions and useful advice and ChEd Boukaré for discussions on solid-melt density contrasts. H.S. acknowledges the support from the Deutsche Forschugsgemeinschaft (grant number SA 2042/3-1). S.P. acknowledges funding from the DFG within the Research Unit FOR 2440 "Matter Under Planetary Interior Conditions". A.-C.P. gratefully acknowledges the financial support and endorsement from the DLR Management Board Young Research Group Leader Program and the Executive Board Member for Space Research and Technology. Numerical computations were partly performed on the S-CAPAD platform, IPGP, France. Main text figures were made using the Generic Mapping Tools (Wessel & Smith, 1995). This is InSight contribution number 130.

Derived data files used in this study (Samuel, 2020) are available using the DOI number: 10.5281/zenodo.4271582 . The numerical codes used to compute the results in this study are described in detail in Ballmer et al. (2017) (Pre-and post-magma ocean HPE and iron enrichment profiles); Hüttig et al. (2013) (stagnant-lid thermochemical convection); Samuel et al. (2019) (Mars parameterized convection); Padovan et al. (2014) (Love number and tidal dissipation); Crotwell et al. (1999) (ray paths and travel times).

Supporting information references

Supporting Information for "The thermo-chemical evolution of Mars with a strongly stratified mantle" Henri Samuel 1 , Maxim D Ballmer 2,3,4 , Sebastiano Padovan 5 , Nicola Tosi 5 , Attilio Rivoldini 6 , Ana-Catalina Plesa

S1 Computation of the basal layer HPE enrichment and iron content

For the layered cases, we consider fractional crystallization of the magma ocean (MO) (Elkins Tanton et al., 2003;Solomatov, 2015)). Accordingly, melt-solid segregation is dominated by crystal settling, and is therefore controlled by mineral-melt partitioning at small crystal fractions. As several chemical elements, such as iron, are preferentially partitioned into the liquid, the MO evolves compositionally towards a progressive enrichment. Thin lines in Figure S1 show the cumulate compositional profile after MO crystallization in terms of the molar ratio Fe# = 100×Fe/(Mg+Fe) for four example cases (with thicknesses of the basal layer of D d =500 km and D d =300 km, and for different bulk-Mars compositions).

To calculate these profiles, we assumed that the Martian mantle was initially entirely molten, and consider a liquid-solid distribution coefficient for Fe-rich and Mgrich end-members

whose value is set to be constant (a simplifying assumption) at 0.3 (Beattie et al., 1991). This value is between those relevant for olivine and pyroxene minerals for Mars-like bulk compositions (Filiberto & Dasgupta, 2011;Elkins Tanton et al., 2003).

We numerically compute the composition upon solidification by removing incrementally crystals at small crystal fractions. The crystal composition that is in equilibrium with the MO is:

Upon the incremental removal of crystals, the MO composition evolves accordingly (by mass balance), becoming gradually enriched in FeO. We consider that meltsolid segregation is efficient, and the MO is always internally well-mixed. The next incremental layer of crystals then also displays a slightly higher Fe# than the previous one. The final layer ultimately reaches a Fe# of 100 (Figure S1a) based on our assumption of constant K D . diffusive, similar to the stably stratified case displayed in Figure 2 (see supporting material Section S4).

S4 Likelihood of occurrence of convective vs. conductive basal layer

We assess the propensity of convective vs. conductive style of heat transfer within the basal layer in the case where the basal layer has a homogeneous composition. Evidently, the presence of a stable compositional stratification depends on poorly constrained aspects related to the origin of the compositional layering (e.g., magma ocean cooling dynamics (Maurice et al., 2017;Ballmer et al., 2017;Boukaré et al., 2018) and differentiation (Elkins Tanton et al., 2003, 2005;Zeff & Williams, 2019)). Therefore, determining robustly the likelihood of occurrence of a stable stratification may be questionable. Nevertheless, while a stably stratified basal layer can only transfer heat through conduction, a basal layer of homogeneous composition will not necessarily be convecting. Indeed, if the layer is too thin and/or too viscous, its associated Rayleigh number may be too small to allow for convection, regardless of any compositional stratification.

To estimate the ability of a homogeneous basal layer to convect, one can simply compare the value of its corresponding thermal Rayleigh number, Ra d , to the critical threshold for convection to occur, Ra cr . For simplicity we assumed a constant value Ra cr = 10 5 . The latter corresponds to a minimum bound in the case of variable viscosity in spherical geometry with modest values of activation energy, because Ra cr is expected to increase with increasing E * (Yanagisawa et al., 2016). Figure S4