Influence of Thermal Stratification on the Structure and Evolution of the Martian Core

The apparent end of the internally generated Martian magnetic field at 3.6–4.1 Ga is a key event in Martian history and has been linked to insufficient core cooling. We investigate the thermal and magnetic evolution of the Martian core and mantle using parameterized models and considered three improvements on previous studies. First, our models account for thermal stratification in the core. Second, the models are constrained by estimates for the present‐day areotherm. Third, we consider core thermal conductivity, kc , values in the range 5–40 W m−1 K−1 as suggested by recent experiments on iron alloys at Mars core conditions. The majority of our models indicate that the core of Mars is fully conductive at present with core temperatures greater than 1940 K. All of our models are consistent with the range of kc=16−35 W m−1 K−1 . Models with an activation volume of 6 (0) cm3 mol−1 require a mantle reference viscosity of 1019−1020(1020−1021) Pa s.

• Majority of our models have an entirely conductive core at present rather than the frequently assumed convecting core • Termination of internal magnetic field on Mars places limits on the thermal conductivity of the core • Pressure dependence of mantle viscosity and the abundance of U/K/ Th affect the early dynamo and so may be constrained by the cessation time

Supporting Information:
Supporting Information may be found in the online version of this article.
2 of 10 flow ( c > a ) to a sub-isentropic heat flow ( c < a ; Nimmo, 2015). As such, modeling the thermal evolution of both the core and mantle requires estimating c as a function of time, which in turn provides insight about the cooling history of Mars.
A variety of scenarios for the evolution of Mars have been proposed, involving different physical processes. Nimmo and Stevenson (2000) and Breuer and Spohn (2003) suggested that a brief period of plate tectonics was required to power a dynamo prior to ∼4 Ga, before a transition to stagnant lid tectonics. In contrast, based on the lack of evidence for plate tectonics on Mars, Williams and Nimmo (2004) showed that an initially super-heated core relative to the mantle could also provide the required short period of rapid cooling to power the dynamo. All previous core modeling studies generally agree that the Martian core should have undergone a rapid decrease in temperature early in its history (e.g., Breuer & Spohn, 2006), powering the dynamo, before the cooling rate of the planet became too low to sustain the dynamo after ∼3.6 − 4.1 Ga.
If the core has cooled down sufficiently to start crystallizing, then different scenarios can be considered. The growth of a large solid inner core (bottom-up crystallisation) would have likely restarted the dynamo (Williams & Nimmo, 2004), and is therefore not favored. Growth of a small inner core might not yet be able to restart the dynamo (Hemingway & Driscoll, 2021) but would be undetectable by current seismic data from InSight (Stähler et al., 2021). However, depending upon the relative slopes of the core temperature and its melting temperature, the Martian core may have started to solidify from the top-down rather than bottom-up (Stewart et al., 2007). In this regime, iron crystals nucleate at the top of the core and sink because of the density contrast with the residual liquid, forming an iron snow. Whilst this process provides a source of power for the dynamo, it is less efficient per degree of cooling than growth of an inner core and so may have formed at the top of the Martian core without restarting the dynamo (Davies & Pommier, 2018).
In this study, we focus on three factors that have not been previously considered by Martian core-mantle evolution studies. First, recent experimental studies on a variety of iron alloys at conditions relevant to the core of Mars have yielded values for the thermal conducitivity, c , from as low as 5 W m −1 K −1 to around 30 W m −1 K −1 , depending on the iron alloy composition in the Fe-S and Fe-S-O-Mg-Si systems (Pommier, 2018;Pommier et al., 2020). These experiments span a lower range of conductivities than previous models of Mars have considered, typically 30-120 W m −1 K −1 (Davies & Pommier, 2018;Hemingway & Driscoll, 2021;Nimmo & Stevenson, 2000;Williams & Nimmo, 2004). Given a is proportional to c , low c values will affect the ability of the core to generate a magnetic field , impacting which evolutionary histories are consistent with magnetic observations. Second, the previous parameterized models did not account for thermal stratification, which likely affects core temperature estimates. These previous models (Breuer & Spohn, 2006;Hemingway & Driscoll, 2021;Williams & Nimmo, 2004) predict that the Martian core was heavily sub-isentropic ( c < a ) for a significant part of its history, including at present. When sub-isentropic, the thermal state of the core deviates away from a convective state toward a conductive one, resulting in a stable, thermally stratified layer that grows from the top of the core downwards, as investigated for Earth's core (Greenwood et al., 2021;Labrosse et al., 1997;Nimmo, 2015). This approach assumes that chemical convection is not present or that it is insufficient to destabilize the thermal layer. Previous studies of Mars have assumed that even when c < a , the Martian core keeps convecting, which limits the predictions that can be made about the present-day temperature of the core. The effect of thermal stratification on core temperature also influences the onset of core crystallisation and the evolution of c . Accounting for thermal stratification as part of parameterized core models permits more precise estimates for the present-day core temperature, complementing other modeling techniques (e.g., Khan et al., 2018;Plesa et al., 2018;Rivoldini et al., 2011) and spacecraft observations (e.g., InSight mission) that are not as sensitive to the core temperature. Third, evolution models can now take advantage of recent improvements in geophysical inversions for the present-day areotherm (Khan et al., 2018). New information on the thermal state of the planet offers an additional constraint on evolution models that has not been previously utilized. In this paper, we performed parameterized modeling of Mars by accounting for the three improvements listed above, providing a revised view of the thermal history of Mars.
GREENWOOD ET AL.

Methods
We investigate the thermal evolution of Mars by coupling parameterized convection models for the core and mantle. Our primary motivation is to investigate the evolution of the core and the growth of stable layers beneath the CMB. Therefore, we use a conventional mantle model as a reference point in order to highlight the impacts of our revised core evolution.
For the mantle, we adopt a well-established parameterization of convection in the stagnant lid regime (Breuer & Spohn, 2006;Nimmo & Stevenson, 2000;Thiriet et al., 2019) with a simplified model for the lithosphere. The isothermal convecting interior is heated/cooled by heat conducted away from the core and into the lithosphere within thermal boundary layers as well as internal heat generated by the decay of radiogenic isotopes of the heat producing elements K, U, and Th. Linear temperature profiles are assumed in the thermal boundary layers and across the stagnant lid. We follow Thiriet et al. (2019) by not attempting to model mantle melting and subsequent crustal growth and instead we assume a fixed crustal size. By neglecting these effects we avoid needing to specify several uncertain parameters, such as the latent heat of melting and mantle solidus (e.g., Morschhauser et al., 2011), in order to focus upon the effects of core stratification and c . We also assume a fixed stagnant lid thickness of 300 km as this makes little difference to the evolution of the convecting interior and core (see Figure S2 in Supporting Information S1 for a comparison of the mantle evolution assuming different stagnant lid thicknesses).
Temperature-and pressure-dependent mantle viscosity, η, in the upper/lower thermal boundary layers is evaluated at the average temperatures/pressures within the relevant boundary layer given by an Arrhenius equation, scaled by a reference viscosity, 0 , at 1600 K and atmospheric pressure where is the activation energy, taken to be 300 kJ mol −1 assuming diffusion creep (Karato & Wu, 1993), is the gas constant, and are the pressure and temperature, and is the activation volume. The core is assumed to be initially well mixed with an adiabatic temperature profile. When c becomes sub-isentropic, a stratified layer grows beneath the CMB with a conductive temperature profile, implemented following Greenwood et al. (2021). An entropy balance is then used to estimate the Ohmic dissipation produced by magnetic field generation, J (Greenwood et al., 2021;Gubbins et al., 2003;Williams & Nimmo, 2004): where k , s , and J refer to the changes in entropy arising from thermal conduction, secular cooling, and Ohmic dissipation within the core, respectively. Since k is proportional to c , a smaller c gives a smaller k and hence a larger J . Dynamo action is inferred when J becomes larger than 1 MW K −1 , which occurs close to the condition c > a (see Text S1 in Supporting Information S1 for further details on the entropy budget). Whilst our model for the core can account for solidification of either a solid inner core or iron snow (Davies & Pommier, 2018;Greenwood et al., 2021), we do not find scenarios where any core fluid freezes and so do not include their associated terms in Equation 2.
About 8 wt% nickel is expected in the Martian core (Wänke & Dreibus, 1994), along with 10-20 wt% sulfur (Khan et al., 2018;Rivoldini et al., 2011;Wänke & Dreibus, 1994) and as such we assumed an Fe-Ni-S core. We used the melting data of Gilfoy and Li (2020), who observed the melting point of Fe-Ni-S at ∼1500 K with 10 wt% S at CMB pressure. This low melting temperature suggests that the entire core is liquid at present and, as mentioned above, we find no scenarios where any solid is formed. One consequence of the absence of a solid phase in our calculations is that the composition of the liquid core is constant over time and radius since no partitioning of S between solid and liquid occurs.

of 10
A number of parameters for Mars are known well enough to be taken as constant in our study (Table S1 in Supporting Information S1). For example, the radius and density of the core have been constrained to 1830 ± 40 km and 6100-6500 kg m −3 (Khan et al., 2018;Stähler et al., 2021). Varying these parameters within their uncertainties does not significantly alter the evolution of the planet relative to some key unknowns which we will now discuss. Due to the thermostat effect, the mantle self regulates its temperature such that the initial core and mantle temperatures have little influence upon the present-day thermal state (Plesa et al., 2015). However, the initial temperatures have a strong influence on an early dynamo (Williams & Nimmo, 2004), in particular the initial super heat of the core, Δ = c − m , where c is the CMB temperature. We vary the initial mantle temperature, m,0 = 1900 − 2400 K and super heat of the core Δ = 0 − 400 K. The mantle viscosity imparts a strong control on both the thermal evolution of the mantle and the early dynamo by scaling the heat flow out of the core and into the base of the lithosphere. Estimates on 0 typically span the range 10 18 − 10 21 Pa s (Breuer & Spohn, 2006;Fraeman & Korenaga, 2010) and so we include 0 with this range in our parameter search. As discussed in the introduction, dynamo operation is heavily dependent upon the core thermal conductivity, and so, we consider a range from 5 to 40 W m −1 K −1 , based on laboratory-based electrical resistivity results for Fe-S and Fe-S-O-Mg-Si alloys (Pommier, 2018(Pommier, , 2020Pommier et al., 2020).
Estimates on the activation volume for a silicate mantle, , span a wide range from 0 to 20 cm 3 mol −1 (Hirth & Kohlstedt, 2003). Based on dynamical models, Plesa et al. (2018) suggested that 0 cm 3 mol −1 is required to explain some properties of Mars, such as the tidal love number 2 . Thermal history models constrained to InSight observations of the upper mantle/lithosphere statistically prefer ≤ 10 cm 3 mol −1 (Knapmeyer-Endrun et al., 2021; Khan et al., 2021). We consider 2 values of : 0 and 6 cm 3 mol −1 in order to produce sets of models with and without a pressure dependence on .
A final uncertain parameter we consider is the quantity of radiogenic heating in the mantle. Mantle melting can extract heat producing elements (HPEs) from the mantle and emplace them into the crust (Morschhauser et al., 2011), changing the abundance of HPEs in the interior. Observations using gamma spectroscopy suggest that the crust may be enriched by a factor of 10 (Taylor et al., 2006) relative to the bulk compositional model of Wänke and Dreibus (1994) (WD94). However, difficulties with the dating and analysis of Martian meteorites (Grott et al., 2013) make it difficult to assess whether they represent ancient or present Mars. In light of this, we consider two model configurations. In the first, we assume the abundance of HPEs available to the primitive mantle according to WD94. In the second, we assume that the HPEs available are reduced by 45% from that of WD94, equivalent to the reduction due to the instantaneous growth of a 50 km thick crust and an enrichment factor of 10. This should indicate the broad trends that changing HPEs abundance has on the simulations, and in particular dynamo generation, without introducing additional uncertain parameters required for modeling crustal growth.
We perform four sets of Monte-Carlo simulations each with 100,000 models in order to search the parameter space. Within each set, we drew a random value from within the previously stated ranges for m,0 , Δ , 0 , and c (uniform prior). 1. The "standard" configuration, where = 0 cm 3 mol −1 and HPE abundance of WD94 is used. 2. The "pressure dependent" set, with = 6 cm 3 mol −1 and the same HPE abundance as the standard configuration. 3. A "reduced HPE" configuration, where = 0 cm 3 mol −1 and the HPE abundance of WD94 is reduced by 45%. 4. A "combined" set, where both = 6 cm 3 mol −1 and the reduced WD94 abundance of HPEs by 45% are used.
Successful models satisfy two constraints. First, the dynamo cessation time, denoted Φ , must be consistent with the magnetic field history of Mars, where the dynamo shut off between 4.1 and 3.6 Ga (400-900 Myr after core formation) (Acuña et al., 1998;Langlais et al., 2012;Milbury et al., 2012;Mittelholz et al., 2020). Second, we choose models from our ensemble that are in agreement with estimates of the present-day temperature of Mars, using the recent estimates of the areotherm from Khan et al. (2018) based on inversions from geophysical data. Their model is not particularly sensitive to the CMB temperature because they cannot resolve the temperature change in the lower mantle thermal boundary layer, and so instead, we focus 5 of 10 on the temperature at the base of the stagnant lid. As a result, successful models correspond to those with m,u = 1650 − 1750 K. Figure 1 shows the time series of three cases that demonstrate the impact of varying thermal conductivity and mantle viscosity on the evolution of the planet. These results are important for interpreting the trends seen in the Monte-Carlo simulations. The reference case is a successful model taken from the "standard" Monte-Carlo simulation set ( m = 2327 K, Δ = 182 K, c = 24 W m −1 K −1 , 0 = 2.5 × 10 20 Pa s, = 0 cm 3 mol −1 , WD94 HPE abundance), whereas the other two cases are unsuccessful. The "low conductivity" case uses the same input parameters as the reference case, with the exception that c is lowered to an extreme of 5 W m −1 K −1 . Finally a "high viscosity" case uses the same input parameters as the reference case except that 0 is raised to 10 21 Pa s, the largest value we consider.

Results
In all cases there is an initially super-isentropic heat flow that drops rapidly, before flattening out (Figure 1a). Comparing the reference case to the "low conductivity" case shows the identical evolution of c until ∼500 Myr. After this time, the reference case becomes sub-isentropic ( c < a ) and a stable thermal layer begins to grow causing divergence between the two cases. The introduction of the stable thermal layer in the reference case relatively elevates the core temperature, leading to an increased temperature difference between the core and bulk mantle, and driving a larger heat flow at the CMB. A stable layer never grows in the "low conductivity" case due to the extremely low a values. Note that if a thermally stable layer had not been accounted for in the reference case, aside from the calculation of k and J , the evolution would have been identical to the "low conductivity" case. The presence of the thermal layer in the reference case elevates c and subsequently c relative to the "low conductivity" case. The relative ∼25% increase in c Figure 1. Time series for the three cases described in the text. Panel (a) shows the CMB heat flow, c (solid lines), and the heat flow down the isentropic temperature gradient at the CMB, a . Panel (b) shows the entropy due to Ohmic dissipation, J . Finally, panel (c) shows the growth of the thermally stratified layer by plotting the radius of the base of the layer, s , normalized to the core radius, c . The reference case is an example of a successful model, and the low c and high 0 cases are both considered unsuccessful. See Figure S4 in Supporting Information S1 for a present day temperature profile of the reference case.
6 of 10 enhances core cooling to offset much of the increase in c by present-day. Models that grow a stable layer have a present-day c ∼ 30 K hotter than models that do not account for a stable layer. This difference in temperature is small relative to uncertainties in geophysical estimations of the interior temperature. However, at the slow cooling rate of the core, this 30 K difference still represents ∼250 Myr worth of core cooling.
Comparison of the "high viscosity" case to the reference case indicates c is decreased because the highly viscous mantle produces thicker thermal boundary layers, through which less heat is conducted. Initially, c drops off more rapidly, before a short rise and a local maximum, followed by a steady decline, ending up comparable to the reference case value at the present day (Figure 1a). Figure 1c shows the growth of the stable layer through time. The lower heat flows in the "high viscosity" case result in the stable layer growing sooner and faster than in the reference case. The entire core is thermally stratified in both the reference and "high viscosity" cases by the present day. When the stable layer is thin, the expected growth rate is proportional to the square root of time. However, as the heat flow continues to drop and the layer grows, there is an acceleration in the growth rate when only a small proportion of the core is still convecting ( s∕ c < 0.4 ). The volume of the convecting region (proportional to 3 ) shrinks faster than the heat extracted from it (proportional to 2 ). This relation results in a decrease of the ratio of stored heat to heat extracted, leading to rapid cooling of the convecting region and subsequent rapid movement of s .
The effects of 0 and c on the dynamo entropy J are illustrated in Figure 1b. The "low conductivity" case predicts a dynamo operating at all times due to the low entropy associated with thermal conduction, k . This is inconsistent with magnetic field observations. The reference and "high viscosity" cases are both characterised by an early decrease in J , falling below our threshold for dynamo action of 1 MW K −1 (dashed line). However, the dynamo fails too early in the "high viscosity" case and only the reference case fits the time constraints on Φ . Note that J does not fall below 0 because accounting for a thermally stable layer when the core is sub-isentropic ensures entropy is correctly balanced (Text S1 in Supporting Information S1).
To describe the results from our ensemble of Monte-Carlo simulations we focus on the "standard" model configuration in Figure 2 (see Figure S1 in Supporting Information S1 for our other configurations). The left panel of Figure 2 illustrates the correlation between the dynamo cessation time, Φ , and the present-day upper mantle temperature m,u . Many models fit either the constraint upon m,u (n = 26,180) or Φ (n = 9040) and a small fraction of them (n = 2010) fall within the limits for both constraints. No correlation between m,u and Φ exists as Φ is primarily sensitive to heat flows rather than temperatures. The color scale indicates that there is also no preference for 0 based on models that fit the dynamo constraint ( Φ = 400 − 900 Myr) alone. However, there is a strong correlation between m,u and 0 , where m,u is proportional to log( 0) . Higher viscosities limit the heat release through the upper mantle thermal boundary layer, insulating the planet and producing a hotter present day mantle than an equivalent model with a lower viscosity. As such, the present day areotherm offers a complementary constraint to the dynamo cessation time since it limits the reference viscosity in our "standard" set to 0 ∼ 10 19.9 − 10 20.7 Pa s. The same relationships exist throughout our other sets of models although the viscosity is constrained to different ranges; "reduced HPE": 0 ∼ 10 20 − 10 21 Pa s, "pressure dependent": 0 ∼ 10 19.2 − 10 20 Pa s, "combined": 0 ∼ 10 19.6 − 10 20.7 Pa s.
The right panel of Figure 2 shows how all 2010 successful models of the "standard" set require c > 16 W m −1 K −1 . At c < 20 W m −1 K −1 the results fall onto two branches where the dynamo either fails early (<500 Myrs) or persists for too long (>1 Gyrs). On the lower branch, where the dynamo is short lived, models have the highest viscosities (high 0 and low m,0 ) and/or low initial superheats Δ . These conditions lead to low CMB heat flows that quickly become sub-isentropic and hence produce an early Φ . The absence of models in-between the branches arises from the behavior of c at higher mantle viscosities. As described for the "high viscosity" case on Figure 1, the local maximum in c at ∼ 1000 Myrs forces the dynamo to end ( c < a ) either early on the lower branch, or much later on the upper branch. For c > 20 W m −1 K −1 , a is sufficiently large that lower viscosities, which do not exhibit this local maximum in c , can provide the desired range of Φ . The color-scale on the right panel of Figure 2 indicates the proportion of the core that is convecting at present-day (see Figure S1 in Supporting Information S1 for the equivalent figure for the other configurations). We find that all successful models across all sets of simulations with c > 25 W m −1 K −1 are fully thermally stratified. For all sets of simulations except for the "combined" configuration, a fully stratified 7 of 10 core is the most common present-day state. When c > 20 W m −1 K −1 , there are fewer solutions we deem successful when = 6 cm 3 mol −1 ("pressure dependent" and "combined" sets; Figure S1 in Supporting Information S1). This trend arises from a balance of requiring a low enough lower mantle viscosity in order to sustain a dynamo for at least 400 Myrs, whilst maintaining a viscosity high enough in the upper mantle such that the mantle does not cool below our specified range for m,u . A viscosity contrast between upper and lower mantle makes this balance more likely for c < 20 W m −1 K −1 . In the "combined set," reducing the internal heating in the mantle makes this balance even more difficult since the mantle cools faster. Conversely, we would expect an increase in the abundance of HPE in the mantle to permit more successful solutions at c > 20 W m −1 K −1 when = 6 cm 3 mol −1 . Depending upon the specific configuration, a different range of values for c produce successful models fitting our constraints on both Φ and m,u . The "standard" and "reduced HPE" configurations both find a greater number of successful models with increasing c , requiring a minimum value of c of 16 or 7 W m −1 K −1 , respectively. The model sets including = 6 cm 3 mol −1 reverse this trend, instead producing more successful models when c < 20 W m −1 K −1 . We also note that for the "combined" set, only two successful models were found at c > 25 cm 3 mol −1 and that in general the inclusion of a pressure dependence on the viscosity makes it much harder to satisfy our two constraints on Φ and m,u (4632 successful models between the "standard" and "reduced HPE" configurations, vs. 213 successful models for the "pressure dependent" and "combined" configurations). Table S2 in Supporting Information S1 contains a summary of the ranges for c , 0 , s∕ c , and c attained from each configuration. In all of our models the core is far hotter than the liquidus temperature observed by Gilfoy and Li (2020), with the successful models giving c > 1940 K at the present-day. Our successful models are also above the liquidus curve given by Stewart et al. (2007) and those used by Rivoldini et al. (2011) and Hemingway and Driscoll (2021) for sulfur concentrations >10 wt% S. Note that we assumed an isothermal mantle except in 8 of 10 thermal boundary layers, and so, accounting for the roughly 100 K adiabatic increase in temperature across the convecting mantle would yield CMB temperatures roughly 100 K higher than we obtain. Figure 3 demonstrates the general impact of varying any one of the four variables 0 , m,u , Δ or c upon Φ and m,u relative to the reference case in Figure 1. The present-day temperature of the mantle is almost exclusively controlled by 0 which also exerts some control on the dynamo cessation time, particularly at higher values of 0 . The other inputs, m,u , c , and Δ , instead almost solely influence the cessation time, with the largest influence on Φ coming from c .

Discussion and Conclusions
Due to the large amount of FeO in the Martian mantle (Wänke & Dreibus, 1994), significant amounts of oxygen may have dissolved into the core (Tsuno et al., 2007). The presence of 0.5 wt% oxygen added to Fe-5 wt%S is expected to drastically reduce c to ≈18 W m −1 K −1 at 2000 K . Furthermore, nickel is present in metallic cores, and can significantly reduce c as well as the melting temperature (Gilfoy & Li, 2020). The addition of 10 wt% Ni to Fe-5wt%S would halve c to ≈20 W m −1 K −1 (Pommier, 2020). These experiments were conducted at 8-10 GPa, lower than the 20-40 GPa in the Martian core. Extrapolation of these results to 20-40 GPa suggests that c reduces by approximately 10% . Since >15 wt% of sulfur is predicted based on density estimates in the Martian core (Rivoldini et al., 2011;Khan et al., 2018), significantly lower values of c than 20 W m −1 K −1 may be expected. Our four different model configurations require differing ranges of c highlighting the need to further explore c as a function of composition and pressure to interpret the magnetic and thermal history of Mars.
Inference of the dynamo cessation time from our evolution models is dependent upon the scenario for cooling of the planet, where we have assumed stagnant lid convection in the mantle. Early plate tectonics has been proposed based upon observations of geological structures in the northern lowlands (Sleep, 1994) and magnetic anomalies in the southern highlands (Connerney et al., 1999) that were hypothesised to indicate ancient sea floor spreading. Plate tectonics would allow the mantle to cool more rapidly, increasing the heat flow from the core (Nimmo & Stevenson, 2000). Little evidence of plate tectonics has been subsequently discovered and any significant period of plate tectonics appears incompatible with the present-day crustal thickness (Breuer & Spohn, 2003). Given available information, assuming the stagnant lid regime for all time seems reasonable.
Water in olivine crystals can significantly impact (Mackwell et al., 1985), particularly in the Fe-rich Martian mantle (Kohlstedt & Mackwell, 2010). To some extent, this effect of H and Fe on viscosity is captured by our consideration of a wide range for 0 . However, we did not account for a viscosity change that may  Figure 1. Arrows and circular data points show the influence of the variable (indicated by colour) if only that variable is changed from the reference case. For example, the bottom right point marked 10 21 shows the results using all inputs the same as the reference case, except that 0 is changed to 10 21 Pa s. 9 of 10 arise from changing hydration levels of the mantle with time, yet in the absence of sufficient evidence of changing water content in the mantle we cannot include this effect at this stage.
The majority of our successful models are entirely thermally stratified and only a few (n = 325) models have a convecting region making up at least 50% of the core, which disproportionately originate from the "pressure dependent" and "combined" model configurations. Previous studies have assumed the core is isentropic (e.g., Khan et al., 2018) in order to construct an interior model for Mars. Our findings suggest that the temperature structure of the core is likely conductive and so an imposed conductive temperature profile should also be considered. Furthermore, the time evolution of c is significantly modified when thermal stratification is present in the core and should be taken into account by future thermal evolution models for Mars.
In summary, we have conducted a suite of models for the thermal evolution of Mars including thermal stratification in the core and considering a range of core thermal conductivities based on recent experimental data. In order to match estimates of the termination of the Martian dynamo, all of our model configurations agree upon the range c = 16 − 35 W m −1 K −1 . More successful models are found at progressively larger values for c if the activation volume is zero. Conversely, when the activation volume is 6 cm 3 mol −1 , this trend is reversed and more successful models have <20 W m −1 K −1 . Future studies employing parameterized thermal history models can therefore benefit from considering the cessation time to further constrain relevant parameters to the interior properties of Mars. Furthermore, we find it likely that the Martian core is entirely thermally stratified with a hot CMB temperature of >1940 K. Finally, within each of our different model configurations, the reference viscosity of the mantle was constrained to within 1 log unit, with all configurations spanning 0 = 10 19 − 10 21 Pa s.