{'fraction_non_alphanumeric': 0.030890844873016764, 'fraction_numerical': 0.018248584403206817, 'mean_word_length': 4.842613822469702, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Can a fractionally crystallized magma ocean explain the thermo-chemical evolution of Mars?

The impact heat accumulated during the late stage of planetary accretion can melt a significant part or even the entire mantle of a terrestrial body, giving rise to a global magma ocean. [...] Assuming fractional crystallization of the magma ocean, dense cumulates are produced close to the surface, largely due to iron enrichment in the evolving magma ocean liquid (Elkins-Tanton et al., 2003). A gravitationally unstable mantle thus forms, which is prone to overturn. We investigate the cumulate overturn and its influence on the thermal evolution of Mars using mantle convection simulations in 2D cylindrical geometry. We present a suite of simulations using different initial conditions and a strongly temperature-dependent viscosity. We assume that all radiogenic heat sources have been enriched during the freezing-phase of the magma ocean in the uppermost 50 km and that the initial steam-atmosphere created by the degassing of the freezing magma ocean was rapidly lost, implying that the surface temperature is set to present-day values. In this case, a stagnant lid forms rapidly on top of the convective interior preventing the uppermost dense cumulates to sink, even when allowing for a plastic yielding mechanism. Below this dense stagnant lid, the mantle chemical gradient settles to a stable configuration. The convection pattern is dominated by small-scale structures, which are difficult to reconcile with the large-scale volcanic features observed over Mars' surface and partial melting ceases in less than 900 Ma. Assuming that the stagnant lid can break because of additional mechanisms and allowing the uppermost dense layer to overturn, a stable density gradient is obtained, with the densest material and the entire amount of heat sources lying above the CMB. This stratification leads to a strong overheating of the lowermost mantle [...]


Introduction
Geochemical analyses of the so-called SNC meteorites suggest the existence of three to four separate and isotopically distinct reservoirs in the martian mantle, which have been preserved over the entire planetary evolution (e.g., Jagoutz, 1991, Papike et al., 2009. Two of these reservoirs are depleted in incompatible elements and are most likely situated in the mantle, whereas the third one is enriched and could be located either in the crust or in the mantle (Foley et al., 2005). Until present, most of the dynamical simulations of Mars' mantle convection did not account for the formation or preservation of such reservoirs, and focused mainly on explaining, for instance, the crustal dichotomy or the Tharsis bulge (e.g., Breuer et al., 1998, Harder and Christensen, 1996, Schumacher and Breuer, 2006, Keller and Tackley, 2009,Šramek and Zhong, 2012. Scenarios that consider the early formation of chemical heterogeneities in the Martian mantle assume either a global magma ocean (Elkins-Tanton et al., 2003, 2005a, Debaille et al., 2009, or at least substantial partial melting in the earliest phase of the planet's evolution (Schott et al., 2001, Ogawa and Yanagisawa, 2011. It is suggested that early in the evolution of Mars, the large amount of primordial heat due to accretion, core formation, and possibly short lived radioactive elements can give rise to a magma ocean as a consequence of significant or perhaps even complete melting of the mantle (e.g., Breuer and Moore, 2007). This assumption is confirmed by studies on short-lived radionuclides such as 182 Hf suggesting that the separation between silicates and iron occurred in the first million years after accretion (e.g., Kleine et al., 2002). Furthermore, estimates of the timescale of the process of core formation suggest that, to achieve rapid separation, both silicates and metals need to be fluid, at least in the upper part of mantle (e.g., Stevenson, 1990).
In general, the magma ocean freezes from the bottom (i.e. from the CMB, if it comprises the entire mantle) to the surface because of the steeper slope of the mantle adiabat compared to the slope of the solidus (e.g., Solomatov, 2000). Solomatov (2000) discusses two main freezing mechanisms: equilibrium and fractional crystallization, depending on the size and the settling velocity of crystals that form upon solidification. In the case of equilibrium crystallization, the crystal size and the settling velocities are so small that freezing takes place before crystal-melt separation occurs. The magma ocean solidifies without differentiating and a chemically homogeneous mantle is formed. In the case of fractional crystallization, instead, denser cumulates are formed while the crystallization of the magma ocean proceeds from the bottom, mainly because of iron enrichment in the evolving magma ocean liquid (Elkins-Tanton et al., 2003). The consequence is the formation of an unstably stratified density profile. Whether equilibrium or fractional crystallization in a magma ocean is the dominant process also depends on the time scale of freezing and may change with time. A large degree of melting of the mantle implies a very low viscosity of the magma ocean. The viscosity of a melt/crystal mixture increases abruptly near a critical crystal fraction as suggested by theoretical and experimental studies of concentrated suspensions (e.g., Mooney, 1951, Roscoe, 1952, Campbell and Forgacs, 1990 and by experiments with partial melts (e.g., Arzi, 1978, Lejeune andRichet, 1995). This sudden increase in viscosity is defined by the so-called "rheological transition", which depends on the crystal size distribution, the crystal shape and orientation, and other factors. The critical crystal fraction for this tran-3 sition varies in the literature. Saar et al. (2001), using numerical models, found that for plagioclase crystals, the critical crystal fraction that allows a crystal network to form lies between 8 and 20 %. Other studies place this transition at 30% to 50% depending on the basalt composition (Philpotts et al., 1996). In Solomatov (2000), the rheological transition is assumed to take place at ∼ 60 %. Solomatov (2000) argues for instance that freezing of the Earth's lower mantle lasts only a few hundred years, until the crystal fraction becomes larger than the critical value marking the rheological transition. For larger crystal fractions, the time scale for further crystallization slows down by several orders of magnitudes (∼ 10 to 100 Ma) as the convection speed and cooling efficiency of the partially molten mantle is significantly reduced. Fractional crystallization may become dominant at this stage. If the crystals do not remain in suspension, a crystal network forms instead and efficiently compacts, ultimately resulting in a fractionation between crystals and melt. Therefore, depending on how the magma ocean freezes, i.e. via fractional or equilibrium crystallization, either a chemically homogeneous or stratified mantle is obtained and different evolutionary paths can be expected. So far, two end-member scenarios have been suggested in the literature for Mars: 1) In the first scenario, it is assumed that the early martian mantle is homogeneous apart from an upper depleted layer (Schott et al., 2001, Ogawa andYanagisawa, 2011). Such a structure can evolve if Mars did not experience a deep magma ocean but a large amount of melting causing secondary differentiation early in its evolution. Alternatively, if a deep magma ocean existed, the largest part of the mantle, from the CMB up to a certain depth, freezes rapidly by equilibrium crystallization and only melt from a remaining shallow layer rises via porous flow due to compaction of the silicate matrix or through cracks and channels.
This early melt produces a primordial crust and marks the start of the differentiation process.
As the mantle residue is in general less dense than the primordial mantle, chemically distinct reservoirs can form depending on the density contrast between depleted and primordial mantle . Typically, a depleted and compositionally buoyant upper mantle can develop, which overlies an undepleted and compositionally denser lower mantle (Schott et al., 2001, Ogawa andYanagisawa, 2011). These two layers convect separately with the lower one that tends to become thinner with time because of erosion, while its composition remains almost primordial. However, whether or not this scenario can actually explain the isotopic characteristics of the martian meteorites is unclear and, to the best of our knowledge, it is a problem that still needs to be investigated.
2) In the second scenario, the mantle freezes by fractional crystallization resulting in a stratification with Fe-rich cumulates close to the surface and Mg-rich cumulates at the CMB, a configuration that is gravitationally unstable. Being highly incompatible, radioactive elements are enriched in the evolving residual liquid of the magma ocean. Therefore, a substantial amount of radioactive heat-sources is concentrated in the uppermost layer of Mars' mantle (Elkins-Tanton et al., 2005b). The unstable layering results then in an global mantle overturn, which has been estimated to take place within ∼ 10 − 100 Ma (Elkins-Tanton et al., 2003, 2005b, Elkins-Tanton, 2008, Debaille et al., 2009. As a consequence of the overturn, the dense surface material sinks to the core-mantle boundary, while the light material from the CMB rises to the surface, thus leading to a stable, chemically layered mantle. It is further assumed that, during the overturn, rising upwellings can melt by adiabatic decompression and produce an early crust. The residual mantle of this process has been associated with the source region of the nakhlites (Debaille et al., 2009). The chemical layering following the overturn suppresses subsequent thermally driven convection and early-formed reservoirs are preserved throughout the planet's history. In this way, distinct reservoirs can be maintained stable over the rest of the thermal evolution of Mars.
An important feature of the mantle overturn is also the subduction of the dense surface material. The first models with which this process was investigated (Elkins-Tanton et al., 2003, 2005a) assumed a small viscosity contrast between the interior and the surface, implicitly allowing for the mobilization of the surface layers. A small viscosity contrast can be the consequence of high surface temperatures maintained during the solidification phase of the magma ocean. In fact, a large amount of volatiles is expected to be outgassed, thereby forming a steam atmosphere able to keep the surface temperature even above the solidus of silicates (Abe, 1993, Elkins-Tanton et al., 2005b, Elkins-Tanton, 2008. However, whether and how long a thick proto-atmosphere can be preserved is unclear since atmospheric loss processes are likely to be enhanced in the early stage of planetary evolution due to a stronger solar activity. Lebrun et al. (2013), using a convective-radiative atmospheric model in equilibrium with the magma ocean, showed that a thick proto-atmosphere could be lost in about 10 Ma. A recent study based on a 1D hydrodynamic upper atmosphere model taking into account the extreme XUV conditions of the young Sun suggests that a thick atmosphere produced by outgassing during magma ocean crystallization can be lost in 0.14 to 14 Ma (Erkaev et al., 2013). The time needed to overturn the densest cumulates has been computed assuming a constant viscosity of 10 20 Pa s and was found to be a fraction of a million year (Elkins-Tanton et al., 2003), thus comparable with the time predicted for the loss of the initial atmosphere. However, even assuming that a high surface temperature (900 • C according to Elkins-Tanton et al. (2005a)) can be maintained for longer than 10 Ma, when using an activation energy of 300 kJ/mol, as expected for an olivine-dominated rheology, the resulting surface viscosity would be about 8 orders of magnitude higher than the bulk viscosity of the interior, which still represent a major obstacle for the mobilization of the uppermost cumulates.
Another possibility to obtain surface mobilization, even in the presence of low surface temperatures, is by lithospheric failure. Only recently it has been shown that rising upwellings in the early stage of the overturn can localize surface deformation, inducing lithospheric-scale yielding and consequent surface mobilization (Debaille et al., 2009. These models, however, used for simplicity either a linear unstable density profile , or a density profile in which phase transitions were neglected (Debaille et al., 2009). Furthermore, an approximation of the viscosity law, the so-called Frank-Kamenetskii approximation (F-K), which can result in an underestimation of shallow viscosity gradients, is commonly used (Debaille et al., 2009).
In the present work, in contrast to previously published results that mainly focused on the overturn itself, we further investigate the consequences of a magma ocean cumulate overturn on the subsequent thermo-chemical evolution of Mars. We will show that considering a density profile that accounts for phase transitions, as well as a more appropriate strongly temperature-dependent rheology and a plastic deformation mechanism, the overturn tends 6 to be confined below a stagnant lid and does not lead to a mobilization of the surface layers.
Even if a whole mantle overturn takes place, assuming additional mechanisms discussed below, the density gradient obtained after the overturn suppresses convective heat transport.
In both cases, i.e. whole mantle overturn and overturn below the stagnant lid, the subsequent thermo-chemical evolution is hard to reconcile with observations. The findings of the present study demonstrates the difficulty to find plausible scenarios to fit results from geochemical and isotope analyses with those of geodynamical modelling.

Physical Model
We assume that the amount of heat generated by accretion and core formation processes has been large enough during the early planetary evolution to melt the whole mantle, creating a global magma ocean. Subsequent cooling and fractional crystallization of the magma ocean results in a stratified and gravitationally unstable mantle. We assume that the initial steam atmosphere created during the magma ocean crystallization phase is entirely lost at the time when the magma ocean has fully crystallized (Lebrun et al. (2013), Erkaev et al. (2013)), such that the surface temperature is set to a present day value of 250 K.
Radiogenic elements, being highly incompatible, will preferentially partition in the evolving liquid phase during the fractional crystallization of the magma ocean. Following the model by Elkins-Tanton et al. (2005a), we set the entire amount of heat producing elements in the uppermost 50 km.
The density profile established after freezing of the magma ocean is a crucial parameter, since it can dominate the entire dynamic and even suppress thermal buoyancy. However, such a profile is poorly constrained and varies in the literature from a linear function of the mantle depth (Zaranek and Parmentier, 2004) to a more complex one, strongly depending on the mantle mineralogy (Elkins-Tanton et al., 2003). In order to test the effects of the initial density profile and to compare our results with previous models, we compare in the following three different profiles ( Figure 1). We run one simulation assuming an unstable linear profile. For this we use the density contrast from Elkins-Tanton et al. (2005a) calculated at a reference temperature of 1 • C.
Therefore the linear profile has the same buoyancy ratio as the third profile presented below.
A second run is performed using the initial density profile from Elkins-Tanton et al. are also interested in the evolution after the magma ocean crystallization and overturn, we use densities at 1 atm and 1 • C rather than 1 atm and solidus temperature, which would be valid only in the overturn phase. Additionally, we account in our models for the transition from garnet and pyroxene to majorite and olivine to γ-olivine. As in the model of Elkins-Tanton et al. (2005a), we assume for simplicity that all phase transitions occur at a pressure of ∼ 13 GPa equivalent to the depth of ∼ 1000 km in Mars. We calculate the density difference of the mantle material according to the magnesium number and the volume percent of each phase (Figure 2b). Therefore, when material of the shallow regions of the mantle sinks to the core-mantle boundary, garnet, pyroxene and olivine increase their density as a consequence of the transition to majorite and γ-olivine. For material from the deep regions containing γ-olivine and majorite and rising to shallower depths, the density is reduced according to the phase transition. Although the transition from pyroxene and garnet to majorite is slightly endothermic, on average, the transition of all phases is exothermic (Elkins-Tanton et al. (2003)). The exothermic behavior is taken into account in such a way that the phase transition in cold downwellings takes place at a shallower depth compared to hot upwellings.

Mathematical Formulation
We model the thermo-chemical convection of Mars' mantle by solving the conservation equations of mass, momentum, energy and composition. These are scaled as usual using the thickness of the mantle D as length scale, the thermal diffusivity κ as time scale, the temperature drop ∆T and the maximum chemical density contrast ∆ρ across the mantle as temperature and compositional scales. Using the extended Boussinesq approximation (EBA), a Newtonian rheology and an infinite Prandtl number, as appropriate for highly viscous media with negligible inertia, the non-dimensional equations of thermo-chemical convection read: where u is the velocity vector, p is the dynamic pressure, T is the temperature, t is the time, η is the viscosity, u r is the radial velocity and e r is the radial unit vector. Ra and Ra C are the thermal and compositional Rayleigh number respectively, H is the internal heating rate, which is given by the ratio of Ra Q and Ra, where Ra Q is the Rayleigh number for internal heat sources (see Table 1). The three Rayleigh numbers read: In all simulations we use a fixed surface temperature, a cooling boundary condition at the core-mantle boundary, and decaying radioactive heat sources. Similarly to the model of Morschhauser et al. (2011), the concentration of long-lived elements is taken from Wänke and Dreibus (1994). The viscous dissipation is defined as follows: The dissipation number, Di, depends on the thermal expansivity α, the gravity acceleration g, the mantle thickness D and the heat capacity of the mantle c p : All variables and parameters used in equations (1)-(5) are listed in Table 1.
The viscosity is calculated according to the Arrhenius law for diffusion creep (Karato et al., 1986). The non-dimensional formulation of the Arrhenius viscosity law for temperature dependent viscosity reads (e.g. Roberts and Zhong, 2006): where E is the activation energy, T surf the surface temperature and T ref the reference temperature at which a reference viscosity is attained (see Table 1).
In some cases, we use a linearization of equation (8), the so-called Frank-Kamenetskii (F-K) approximation. For temperature dependent viscosity, this can be written as: where γ is the F-K parameter (Solomatov and Moresi, 1997), which in our simulations is set to 10 5 as in Debaille et al. (2009).
The reference viscosity is set to 10 22 Pa s at a reference temperature of 1600 K which results in a lower mantle viscosity of ∼ 10 18 − 10 19 Pa s. It should be noted that such reference viscosity, typical of dry materials, is likely too high for Mars' mantle, which is expected to be characterized by a wet rheology (Grott et al., 2013) and hence by a lower viscosity. Nevertheless, this choice only affects the timescale of mantle overturn (Sections 3.1 and 3.2) and does not impact our main conclusions. A mechanism of pseudo-plastic yielding is used to allow the stagnant lid to self-consistently fail if the convective stresses are high enough to overcome a depth dependent yield stress σ y . For this, the effective viscosity reads: where σ y = σ 0 + z∂σ/∂z, z is the depth, andε is the second invariant of the strain rate tensor. For the depth-dependent yield stress, we use a surface value σ 0 = 10 8 Pa and a gradient ∂σ/∂z = 160 Pa/m (Nakagawa and Tackley, 2005).

Numerical Model
To solve equations (1)-(4) we used the finite volume code Gaia (Hüttig and Stemmer, 2008). All simulations were run in a 2D-cylindrical geometry on a structured grid consisting of 170 equally spaced shells in the radial direction and 880 points in the lateral direction (corresponding to a lateral resolution of 10 km in the mid of the mantle) resulting in a total of ∼ 1.5 × 10 5 grid points. The cylindrical domain was also rescaled to match the inner-to outer-surface ratio of a three-dimensional spherical shell (van Keken, 2001).
This minimizes spurious effects due to numerical diffusion, which typically arise when using grid-based methods. In all tests, we assigned 20 particles to each cell, resulting in a total of ∼ 3 × 10 6 particles. Our PIC implementation has been tested using various benchmark runs of increasing complexity in both 2D cylindrical and 3D spherical geometry . In our simulations each particle has its own density (C) value. Furthermore, when accounting for an exothermic phase transition, the particles are additionally characterized by a phase parameter and a density difference (∆C). At each time-step and for each particle, the following phase function is computed (Christensen and Yuen, 1985): where z 0 (T ) = z 0 + γ 0 (T − T 0 ) is the temperature-dependent depth of the phase boundary, γ 0 , z 0 , T 0 , and w are the Clapeyron slope, reference depth, reference temperature, and width of the phase transition, respectively. When a particle crosses the phase boundary calculated from equation (11), it is assigned an appropriate density jump. Note that we neglected the release and consumption of the latent heat associated with the phase change. For an exothermic transition like the one assumed here, this results in an slight overestimation of its enhancing effect on mantle flow.

Results
We apply the thermo-chemical mantle convection model to simulate magma ocean cumulate overturn assuming three different density profiles (Figure 1). If not otherwise specified, the viscosity is calculated using an activation energy of 300 kJ/mol, as appropriate for dry olivine (Karato et al., 1986), and applying the Arrhenius law (equation (8)). Despite the presence of a high concentration of heat producing elements near the surface, this causes the formation of a stagnant lid, which can be mobilized when considering a plastic rheology (equation (10)). The initial temperature profile is the same in all tests and follows the mantle solidus, apart for a thin upper thermal boundary layer with a thickness of 20 km (black line in Figure 3a).

Effects of the Initial Chemical Density Profile
We investigate first the influence of the initial chemical density profile on the cumulate overturn of an unstably stratified mantle using a strongly temperature-dependent viscosity and taking into account plastic behaviour. As shown in Figure  A simple linear profile has been widely used in the literature (Hansen and Yuen, 2000, Zaranek and Parmentier, 2004. In contrast to previous studies in which simulations are started from a stable configuration, i.e. after the magma ocean overturn, here we start with an unstable profile with the densest material located at the top. Due to the temperature dependence of the viscosity, dense material lying initially just beneath the stagnant lid sinks into the deep mantle while lighter material rises to the surface. The convective stresses created by both the light material which rises to shallower depths and the return flow due to the sinking of heavy cumulates are high enough to overcome the imposed yield stress. The lithosphere experiences plastic yielding and the initially unstable chemical gradient evolves into a stable configuration. Figure 3a)  Upon overturn, the dense garnet layer located in the mid-mantle sinks to the core-mantle boundary.
In the third case, the density profile considering phase transitions is even more unfavorable for whole-mantle overturn than the previous one. Similarly to the second simulation, the overturn only involves the mantle below the stagnant lid. Upwelling material from the lower part of the mantle is blocked by lighter cumulates residing in the upper part. The garnet layer located in the upper mantle changes its density only slightly when undergoing the phase transition at greater depths. As a consequence, and in contrast to the results neglecting the phase transitions, the garnet layer does not sink to the core-mantle boundary but remains in the upper part of the lower mantle just below the phase transition. The smallest density of the lower mantle material, after this material rises and undergoes the phase transition, is about 3450 kg/m 3 , thus still higher than the upper mantle densities below the stagnant lid (see Figure 2). Some material from the lower part of the mantle moves above the garnet layer whose density is about 3660 kg/m 3 . However, only a small amount of material can be exchanged between the upper and lower mantle as the garnet layer, which remains at approximately the same depth, acts against advective transport. Exchange can only take place during the overturn and ceases as soon as a stable configuration is reached in the lower part of the mantle. The mantle then cools rapidly to a conductive state as shown in Figure 5 and melting in the upper mantle only occurs during the first 1000 Ma.

Arrhenius Law vs. Frank-Kamenetskii Approximation
In the following models we only use the third chemical density profile discussed in the previous section, and investigate the effects of the rheological law. It is important to note that the results presented in previous studies (e.g., Elkins-Tanton et al., 2005b, Debaille et al., 2009) are strongly influenced by the rheological parameters used. For example, the use of a low activation energy and high surface temperature, or of a small viscosity contrast all contribute to considerably reduce global viscosity contrasts. In Elkins-Tanton et al. (2005a), the use of an activation energy of only 100 kJ/mol and a surface temperature above 1000 K results in a viscosity contrast of about three orders of magnitude, which shifts the convection regime from stagnant-to sluggish-lid, in which also the surface layer becomes mobile.
In contrast to our findings, the densest material located at the surface overturns because of the mobilization of the lithosphere. More recently, Debaille et al. (2009) used a higher viscosity contrast (γ = ln(10 5 )) and showed that, when using a yielding mechanism, the lithosphere undergoes plastic deformation and whole-mantle overturn can take place as in Elkins-Tanton et al. (2005a). However, the use of the F-K approximation is highly debated, especially when considering plastic yielding as it can alter the viscosity profile in the lithosphere, leading to erroneous results (Noack and Breuer, 2013).
Indeed, when we use the Arrhenius law with an activation energy of 300 kJ/mol and a surface temperature of 250 K, the lithosphere remains rigid and does not experience plastic deformation. If instead we use the F-K approximation with a parameter γ = ln(10 5 ), the upper dense layer overturn. In this case, the less steep viscosity gradient implied by the F-K approximation allows for larger deformations at shallow depths than those observed using the Arrhenius law, with the consequence that the yield stress can be locally exceeded.
Nevertheless, it should be also noted that this behavior is valid for a particular choice of the yield stress. A smaller value of the latter would likely allow for the mobilization of the surface even when using the Arrhenius law . In Figure 4  The presented F-K case can however be used to study the consequences of the global mantle overturn with surface mobilization on the subsequent thermal evolution -in case 15 different mechanisms not considered here would lead to an early surface mobilization. During the overturn, the uppermost dense material sinks and accumulates at the CMB. The entire amount of heat sources, initially enriched in the uppermost 50 km, resides, after the overturn, above the CMB. Due to the stable density gradient established after the overturn, thermal convection is difficult to maintain. Although the upper mantle cools, the lower part of the mantle is strongly affected by the stable density gradient and by the enrichment of heat producing elements above the CMB. This results in a strong increase of the temperature to values that even exceed the liquidus (Figure 5b). Although our model does not include the effects of partial melting, we argue that even if melting were considered, the depth where this occurs would be below the density inversion depth (Ohtani et al., 1995) where the melt, due to its compressibility, is negatively buoyant. Therefore it would either remain in suspension with the silicate matrix or sink to the CMB.

Discussion and Conclusions
We have investigated the consequences of an unstable density profile resulting from fractional crystallization of a global magma ocean on the thermal evolution of Mars. We have found that considering a strongly temperature-dependent rheology with a high activation energy, plastic yielding with a surface yield stress of 10 8 Pa, and a surface temperature close to present-day value as suggested by recent modelling efforts (Lebrun et al., 2013, Erkaev et al., 2013, a global mantle overturn with surface mobilization is unlikely. Note, however, that even a high surface temperature of 900 • C, which can be maintained for longer Ga since the majority of the heat sources are also trapped in the stagnant lid. Although stable chemical reservoirs form early in the evolution -as required by the isotopic data of the martian meteorites -the described scenario has several drawbacks that appear to be inconsistent with our current knowledge of the planet. The strong chemical gradient implied by fractional crystallization (Elkins-Tanton et al., 2005a, Debaille et al., 2009) surppresses convection and partial melting. On the one hand, it is difficult to trace any reservoir in the lower part of the mantle. On the other hand, the rapid decline of convection and partial melting is at odds with the expected long-lived volcanic activity (e.g., Neukum et al., 2004).
Furthermore, a dense primordial lid is difficult to reconcile with the low crustal density estimated for the southern hemisphere (Pauer and Breuer, 2008). Such scenario is also not consistent with the formation of nakhlites as proposed in Debaille et al. (2009), since the garnet layer can not rise to shallow depths and the material most likely will not experience high melt fractions for garnet segregation to occur.
Although we argue that surface mobilization at the end of the magma ocean phase is unlikely, it remains a possibility that can not be ruled out. For instance, a dense atmosphere and a high surface temperature close to or higher than the solidus of silicates may allow surface material to be mobilized -this requires the overturn of the magma ocean within a time shorter than 10 Ma (Elkins-Tanton, 2008). Isotopic analyses of the SNC meteorites suggest, however, source reservoirs differentiation between 10 Ma (Shih et al., 1999, Foley et al., 2005 and 100 Ma (Blichert-Toft et al., 1999, Debaille et al., 2009. Similarly, external processes such as impacts may trigger the instability of the upper dense layer. In that case, the whole mantle overturn leads to a stable density profile characterized by a basal layer highly enriched in radiogenic heat sources. This results in a strong overheating of the lower mantle up to above the liquidus temperature. It is important to note that heating of this dense layer above the CMB does not destabilize it as suggested for example for the ilmenite-rich layer resulting from the solidification of the lunar magma ocean (e.g., Stegman et al., 2003). In the latter case, the difference in density between ilmenite-rich layer and the overlying mantle is assumed to be only 90 kg/m 3 (Stegman et al., 2003, and hence much lower than the overall density contrast we considered in our models (∼ 800 kg/m 3 following Elkins-Tanton et al. (2005a)).
The melt generated in this lower part of the mantle is highly enriched in iron, since the densest layer contains iron-rich cumulates. Moreover, this melt is located below the density inversion depth -the depth below which melt is denser than the surrounding silicate matrix.
For Mars, due to the high Fe-O content of the melts, this depth has been estimated to lie at 7 GPa, corresponding to a depth of about 540 km (Ohtani et al., 1995). Therefore, we speculate that this melt would remain trapped in the lower part of the mantle and, because of the slow cooling of this region, would probably persist until recent times (perhaps until present). The upper mantle instead rapidly cools to a conductive state due to the absence of heat sources. Volcanic activity is expected to end in the first few hundred million years, at odds with the martian volcanic history (e.g., Werner, 2009).
In both described scenarios, i.e. whole mantle overturn and overturn below the stagnant lid, we obtain the formation of stable chemical reservoirs. Nevertheless, because of the strong chemical gradient assumed by the fractional crystallization model (Elkins-Tanton et al., 2005a, Debaille et al., 2009, it is difficult to trace these reservoirs, particularly in the lower part of the mantle. Furthermore, our results can not be reconciled with the subsequent thermo-chemical evolution of Mars. Therefore, we conclude that a purely fractional crystallization most likely did not take place. If a global magma ocean existed, a combination of fractional and equilibrium freezing, allowing for a smaller density constrast, appears more plausible. In fact, in recent study, Tosi et al. (2013) shows that for a lower density gradient Another possibility would be that of a local or hemispherical magma ocean. This could be the result of a large impact in the early stage of the planetary evolution that has been for instance suggested to explain the crustal dichotomy (e.g. Golabek et al., 2011, Andrews-Hanna et al., 2008, Nimmo et al., 2008, Marinova et al., 2008. Such a scenario could possibly reduce the compositional gradient if only local parts of the mantle have been crystallized allowing for efficient mixing of mantle cumulates. In fact, isolated lithospheric oxygen reservoirs, inferred by the analysis of the recently recovered NWA7034 meteorite, also support the scenario of impact driven regional magmatic processes instead of a global magma ocean which would have homogenized these reservoirs (Agee et al., 2013). In addition, density variations in the mantle due to partial melting that are associated with density contrasts of about 60 kg/m 3 between primordial and depleted mantle can explain the early formation of long-lasting reservoirs (Schott et al., 2001, Ogawa and Yanagisawa, 2011. In fact, for a better understanding of the martian thermochemical evolution, both the likely density stratification after solidification of a magma ocean and by partial melting need to be considered.