Iron Snow in the Martian Core?

The decline of Mars' global magnetic field some 3.8-4.1 billion years ago is thought to reflect the demise of the dynamo that operated in its liquid core. The dynamo was probably powered by planetary cooling and so its termination is intimately tied to the thermochemical evolution and present-day physical state of the Martian core. Bottom-up growth of a solid inner core, the crystallization regime for Earth's core, has been found to produce a long-lived dynamo leading to the suggestion that the Martian core remains entirely liquid to this day. Motivated by the experimentally-determined increase in the Fe-S liquidus temperature with decreasing pressure at Martian core conditions, we investigate whether Mars' core could crystallize from the top down. We focus on the"iron snow"regime, where newly-formed solid consists of pure Fe and is therefore heavier than the liquid. We derive global energy and entropy equations that describe the long-timescale thermal and magnetic history of the core from a general theory for two-phase, two-component liquid mixtures, assuming that the snow zone is in phase equilibrium and that all solid falls out of the layer and remelts at each timestep. Formation of snow zones occurs for a wide range of interior and thermal properties and depends critically on the initial sulfur concentration, x0. Release of gravitational energy and latent heat during growth of the snow zone do not generate sufficient entropy to restart the dynamo unless the snow zone occupies at least 400 km of the core. Snow zones can be 1.5-2 Gyrs old, though thermal stratification of the uppermost core, not included in our model, likely delays onset. Models that match the available magnetic and geodetic constraints have x0~10% and snow zones that occupy approximately the top 100 km of the present-day Martian core.


Introduction
Low-altitude vector magnetometer measurements from Mars Global Surveyor show that Mars presently lacks a global dipole field, but reveal large regions of strongly magnetized crust located mainly in the southern highlands (Acuña et al. 1998). The prevailing view is that this magnetization was acquired as the rock cooled in the presence of a global magnetic field (Stevenson, 2001;Breuer and Moore, 2015). The global field was likely produced in the liquid core by a dynamo process in which thermal (and possibly chemical) buoyancy forces drive convective motion (Stevenson, 2001). Inferences based on the age of impact craters (Acuña et al., 1998;Langlais et al., 2012) and Martian meteorites (Weiss et al., 2002) suggest that the global field decayed around 3.8-4.1 Ga. This event marks the demise of the Martian dynamo and may have been contemporaneous with changes in the planets' heat loss (Ruiz, 2014) and oxidation state (Tuff et al., 2013).
Explanations of Mars' magnetic history are intimately linked to the thermal evolution and crystallization regime of its metallic core. A thermal dynamo can operate in an entirely liquid core, provided that the ancient core-mantle boundary (CMB) heat flow exceeded the heat lost by conduction down the adiabatic temperature gradient (assuming no radiogenic heating).
In this scenario the core cooled, perhaps from an initially superheated state compared to the mantle (Williams and Nimmo 2004) or modulated by an early episode of plate tectonics (Nimmo and Stevenson, 2000), until fell below . Impact-induced thermal insulation of the core (Monteux et al., 2013;Arkani-Hamed and Olson, 2010) would produce a similar outcome. On the other hand, a thermochemical dynamo can operate with < . It has been suggested that rapid growth of an inner core early in Mars' history led to dynamo termination when the size of the liquid region fell below a critical threshold (Stevenson, 2001). This scenario has not been favored because inner core growth provides additional power sources that lead to a long-lived dynamo (Williams and Nimmo, 2004). In this paper we investigate a third scenario: the topdown crystallization of the Martian core.
A necessary condition for top-down core freezing is / < / for all pressures , where is the liquidus temperature of the core alloy and is the ambient temperature. / is positive and for an adiabatic core / ∝ ∝ , where is the CMB temperature.
Melting curves for iron-sulfur systems, the mixture used throughout this work, have been extensively studied (Kamada et al., (2012); Morard et al., (2011); Figure 2d). Of particular interest are the results of Stewart et al., (2007) who found that / <0 at the − conditions of Mars' core using Fe-S mixtures with 10.6 % and 14.2 % S, which assures top-down cooling over bottom-up cooling. However, application of top-down crystallization to Mars depends critically on whether its core has cooled sufficiently over the last 4.5 billion years for to fall below . A further issue is that additional power sources accompanying topdown crystallization could have provided sufficient entropy to restart the dynamo.
Here we build a parameterized model of top-down crystallization in the Martian core. We consider the so-called "iron snow" regime that arises when the bulk sulfur concentration is smaller than the sulfur concentration of the eutectic composition: solid produced on freezing is heavier than the residual liquid and iron "snows" down onto the underlying liquid (Hauck et al., 2006;Dumberry and Rivoldini, 2015;Rückriemen et al., 2015). We follow the premise of previous work (Hauck et al., 2006;Rückriemen et al., 2015) and assume that crystallization in the snow zone produces a slurry: solid particles are suspended in a liquid Fe-S mixture and the solid fraction remains low enough that the system behaves as a liquid. The fluid dynamical behavior of a binary slurry is fundamentally different from that of a binary liquid mixture and so the theory must be developed from scratch, starting with the fundamental conservation equations.
We derive energy and entropy equations from an established slurry theory (Loper and Roberts 1977) that does not appear to have been utilized in previous models of iron snow in planetary cores.
Our model assumes that the snow layer is always in phase equilibrium and that freezing produces iron solid that quickly falls to the deeper liquid core (see Sections 2 and 4 for detailed discussion of the modelling approximations). Starting from an equilibrium state with the entire region at the liquidus temperature, cooling reduces below leading to the formation of solid iron ( Figure   1). The local increase in releases latent heat and elevates the sulfur concentration in the coexisting liquid phase, which in turn depresses until it reaches ( < implies the layer is fully liquid). Assuming no net mass exchange between core and mantle on the timescales of interest, the light residual liquid rises, producing a stable chemical stratification across the snow zone (Dumberry and Rivoldini, 2015;Rückriemen et al., 2015). The heavy solid sinks out of the snow layer into the underlying liquid region where it remelts, absorbing latent heat and causing a decrease in sulfur concentration and an increase in . Gravitational energy is liberated in the snow zone due to iron sinking and also in the liquid due to stirring induced by dense iron remelting (Rückriemen et al., 2015; Figure 1). These additional heat sources, together with variations in composition and temperature across the snow zone induced by freezing, influence the core cooling rate and the power available to generate a magnetic field.
The complexity of the iron snow equations together with uncertainties in thermal and material properties of Fe-S alloys at high − conditions mean that we do not expect (or attempt) to obtain a definitive thermal history for Mars. Rather, we seek to understand the conditions that could have led to snow zone formation. Nevertheless, viable models must be compatible with the magnetic history of Mars and with geodetic observations, which suggest that at least part of the Martian core is liquid at the present day (Yoder et al., 2003).

Model and Methods
We generalize an existing 1D thermochemical evolution model (Davies, 2015) to study crystallization of an iron-sulfur alloy (Taylor, 2013) in the Martian core. A standard averaging procedure (Nimmo, 2015) is used to obtained the equations governing changes in the reference or equilibrium state, in which variables depend only on radius . In regions where there is no solid and outside thin boundary layers it is assumed that vigorous convection maintains a reference state where the pressure is determined by a hydrostatic balance, the sulfur concentration is uniform and the radial entropy gradient is zero (Braginsky and Roberts 1995).
These conditions imply that temperature follows an adiabatic profile.
The global energy budget determines the core evolution by balancing against the sum of the heat sources within the core as defined below. The energy balance does not contain information about the dynamo because magnetic energy is converted to heat within the core. The entropy balance contains the irreversible processes of thermal, chemical, mechanical and Ohmic dissipation. Together, these equations describe the thermal and magnetic history of the core.
The general slurry theory describes the time-dependence of particle composition and local departures from phase equilibrium (Loper and Roberts, 1977) and must be simplified for application to planetary cores. We therefore adopt the following two approximations espoused by Loper and Roberts (1977): 1) No light element partitions into the solid phase on freezing; 2) "fast melting", i.e. instantaneous relaxation to phase equilibrium. The first approximation is supported by experiments that reported very low sulfur concentrations in the solid phase (Kamada et al., 2012;Li et al., 2001) and has also been used to model iron snow in Ganymede's core (Rückriemen et al., 2015). The second approximation means that material in an infinitesimal volume of the continuum will melt/freeze instantly.
In an equilibrium snow zone, the entire system is on the liquidus and the liquidus is collinear with the adiabat. Heavy sulfur-depleted solid sinks and eventually falls out of the layer where it remelts because the temperature of the underlying liquid is above the liquidus. We assume, as in Rückriemen et al. (2015), that the timescale for sinking and remelting is much faster than the timescale for changes to the equilibrium state. At each timestep, all of the newly created solid sinks out of the layer and remelts, leaving the layer on the liquidus. We refer to this third approximation as "fast remelting".
The temperature profile may not be adiabatic throughout the Martian core because compositional and/or thermal stratification can develop below the CMB. Consider first the case where > , i.e. the temperature profile is everywhere unstably stratified. Subsequent growth of a snow zone will produce a stable compositional stratification below the CMB. In this case it can be shown using equation (8) below that the isentropic condition requires where is the specific heat capacity, ̅ > 0 is the heat of reaction coefficient defined below, and the fast remelting approximation and ≪ 1 have been used to remove the contribution from radial variation in solid fraction. The first term on the right-hand side is the usual definition of the adiabatic temperature in a homogeneous fluid. The second term increases | / | since / > 0 and shows that there must be a greater variation in temperature in the presence of a stabilizing compositional gradient in order to keep the layer isentropic. Accounting for this second term is complicated because / is determined from the liquidus, which is itself related to / (see below). Instead of undertaking a complex iterative procedure, which seems unnecessary in light of significant uncertainties in several of the model parameters, we ignore variations in in the adiabat. A posteriori estimates (Section 4) reveal that this is a good approximation.
The configuration of unstable thermal stratification and stable compositional stratification is potentially susceptible to oscillatory double-diffusive instabilities since heat and mass have different diffusion coefficients (Turner, 1973). In the standard doubly-diffusive configuration the horizontally averaged temperature is not expected to deviate significantly from the original adiabatic profile in this case (Buffett and Seagle, 2010) and any effect on terms in the global energy and entropy budgets should be minor. These results may not apply to an equilibrium slurry where Fickian diffusion no longer holds (Loper and Roberts, 1977); however, in the absence of theoretical or experimental evidence to the contrary we assume that doubly diffusive effects do not influence the adiabatic profile.
If < a region at the top of the core will become stable to thermal convection. The base of the thermally stable layer is located where the stabilizing thermal buoyancy balances the destabilizing buoyancy forces that drive convection in the deeper core (Lister and Buffett, 1998).
Departures from an adiabat are expected to be small because the thermal diffusion time, = 2 / where is the thickness of the thermally stable layer and ≈ 10 −6 is the thermal diffusivity, is around 10 7 yrs even for layers as thin as 10 km, and should not affect estimates of terms in the global equations significantly. Thermal history calculations for the Earth's core with and without a thermally stratified layer showed only minor differences to the global energy balance (Labrosse et al., 1997), which were caused primarily by the assumption that gravitational energy release occurs only in the unstably stratified region rather than by departures from adiabaticity. A more important effect arises because the inability of mantle convection to evacuate all of the heat brought to the CMB by core convection requires that the top of the core must heat up. Since snow zones form when < at the CMB, formation will be delayed when thermal stratification is present compared to when it is absent. Unfortunately, a thermally stable layer will, in general, not grow at the same rate as a snow layer; creating a parameterization for the dynamics and couplings between regions of different thermal and compositional stability significantly complicates the model and obscures the effects associated with the slurry that we wish to investigate. Our model considers an entirely adiabatic core and will therefore predict a lower and earlier snow zone nucleation than would be obtained if thermal stratification were taken into account; we return to consider the impact of thermal stratification when applying the results to Mars.
The main assumptions used to develop a quantitative model for the equilibrium evolution of the snow zone are: 1) All sulfur remains in the liquid phase on freezing.
3) Fast remelting of sinking solid, i.e. rapid sinking and remelting of solid iron.
4) An adiabatic temperature profile exists throughout the core.
Using approximations (1) and (2) the general thermal energy equation in a slurry can be written (Loper and Roberts, 1977) = −∇ ⋅ + ∇ ⋅ + , where the density , temperature , entropy , and chemical potential of the liquid are all functions of radius and / denotes the material derivative. The heat flux vector and mass flux vector are determined by constitutive relations. The total dissipation is assumed to arise solely from Ohmic heating, where is the electric current density and is the electrical conductivity, since the viscous dissipation is expected to be small in planetary cores (Nimmo, 2015). Radiogenic heating contributes little entropy (Williams and Nimmo, 2004) and is not considered here.
The global energy equation for a slurry is obtained by summing the internal, mechanical and magnetic energies and integrating over the volume of the slurry. These equations are supplemented by the equations describing conservation of total mass and light element, , which can be written (Loper and Roberts, 1977) = − ∇ ⋅ Changes in the total internal energy can be expressed using the equation With the approximations above the total mechanical energy budget for a slurry is where is the gravitational potential, which is calculated locally as described in Davies (2015).
The total magnetic energy budget is the same as for a two-component liquid: where is the Lorentz force. The rate of change of kinetic and magnetic energy, / and / respectively, are small and can be neglected (Nimmo, 2015) from equations (5) and (6).
Adding the integral of equation (1) and (4)- (6) gives Equation (3) has been used to obtain the second term on the left-hand side and is the (outwardpointing) area element on the surfaces that bound . The first term in (7) can be rewritten using the entropy differential [equation 5.9 of Loper and Roberts (1977)], where is the latent heat, = −1/ ( / ) is the thermal expansion coefficient and ̅ = −( / ) is the heat of reaction coefficient. Equation (8) is identical to the entropy differential for a binary liquid mixture except for the last term, which represents changes in entropy produced by latent heat release (absorption) when solid forms (melts).
The total energy budget for the whole core is obtained by applying equation (7) to the liquid and slurry regions and applying boundary conditions at the interface and CMB. We denote using superscripts and quantities on the snow and liquid side of the interface respectively. The constitutive relation for in a binary slurry is (Loper and Roberts, 1977;Loper and Roberts, 1980) where is the flux of solid particles and is the thermal conductivity. Note that = = 0 outside the slurry. At the CMB, we assume for simplicity that there is no net mass exchange; where the unit vector points radially outward. To determine the boundary condition at we follow standard pill-box arguments (Loper and Roberts, 1987), obtaining where is the velocity of the interface and 〈 〉 denotes the jump in the quantity across .
The terms on the right-hand side represent respectively the latent heat and heat of reaction in the shell of freezing material.
Since the core temperature is assumed adiabatic the cooling rate at radius can be related to the CMB cooling rate (Gubbins et al., 2003) = . (13) Assuming that the interface moves in the radial direction then ⋅ = − / . The latent heat of melting is defined as = Δ , where Δ is the entropy change on melting, and hence / can be related to the core cooling rate in a manner analogous to the situation of inner core growth: The gravitational energy released due to rearrangement of light material can be re-expressed using the identity ⋅ ∇ = ∇ ⋅ − ∇ ⋅ and taking the part of the density change due to composition. We separate the contributions to from the freezing out of solid in the snow zone, denoted , and remelting of solid in the liquid region, , as where = −1/ ( / ) is the compositional expansion coefficient for sulfur, assumed constant, and the second term on the right-hand side of (17) gives the contribution due to motion of the interface.
The sulfur concentration in the liquid region below the snow layer is obtained by applying equation (3) to the snow zone and liquid layer and adding: Applying a standard pill-box analysis (Loper and Roberts, 1987), the boundary condition at can be written where we have used the fact that the total mass of sulfur is conserved.
The time-averaging process removes the ⋅ ∇ part of the first two terms in (18). Furthermore, assuming that the liquid region is well-mixed allows / to be taken outside the integral. The second term then becomes where is the mass of the liquid region. Equation (18) can then be written The second term on the right-hand of (20) is very small because is continuous across the interface and so ≈ ( and are to be evaluated on either side of the interface). / is positive because , and hence , decrease with time: more light element is needed to keep the layer at the liquidus. Therefore, as expected, decreases with time as the liquid region becomes more enriched in iron.
We obtain / from the liquidus relation where Δ is the change in volume on freezing and ̅ = / . ̅ is calculated from ideal solution theory as ̅ = × × × 1000/ (J kg -1 ), where is Boltzmann's constant, is the electron volt, is Avagadro's number and is the atomic weight of S (Gubbins et al. 2004). Solid is formed rapidly and subsequent changes in occur due to rearrangement of the solid fraction. We therefore assume that changes in occur on a timescale that is long compared to changes in . Then = (1 − ) and (neglecting pressure changes) This equation resembles the relation = used by Rückriemen et al. (2015), who estimated from an empirical liquidus curve. Relations (20) and (22) determine and .
On the short timescale we neglect variations in . Then We must distinguish between the latent heat released on freezing of solid particles, , and the latent heat absorbed on remelting, . The total mass created, ∫ / , is equal to the mass destroyed; the only difference is that freezing occurs throughout the snow zone whereas remelting occur at . We therefore have = ∫ / , and = −∫ ( ) / .
Substituting equations (13)- (17) and (20)-(23) into (12) allows the global energy equation to be written symbolically as The additional energy sources that arise due to iron snow are the latent heat released due to formation of solid, , latent heat absorbed as falling snow remelts, , gravitational energy released due to mixing of the remelted iron in the liquid region, and gravitational energy released due to the sinking of iron particles in the snow zone. All terms are proportional to the CMB cooling rate as determined by equations (13), (15), (22) and (23).
The entropy equation is obtained from equation (1) in the usual way and is Viscous and chemical dissipations are thought to be much smaller than and so are neglected (Nimmo 2015). Equations (24) and (25) are evolved forward in time using a timestep of 1 Myr.
The initial time is 4.5 Ga and the final time is the present-day unless a snow zone occupies the whole core in which case the calculation is terminated at that point. The thermal and chemical evolution of the coupled snow-liquid system is calculated such that the base of the snow zone is at the liquidus temperature at each time step. This evolution repeats at each model iteration as the core cools.

Model Parameters
Mantle convection sets while core convection sets the CMB temperature and so the evolution of the two systems should strictly be solved simultaneously. However, significant uncertainties in the parameterization of mantle convection, particularly the appropriate rheological law and the scaling of surface and CMB heat flow with temperature, mean that we do not expect to obtain a definitive thermal history for Mars but seek to understand whether the snow regime is potentially compatible with existing geodetic and magnetic observations.
Focusing on the core alone allows us to elucidate the individual effects of the various physical processes that arise from snow zone growth.
We consider two time-series of from previous studies that both match the inferred dynamo cessation time, but nevertheless exhibit significant differences that embody some of the uncertainties in the mantle problem. The time-series of Williams and Nimmo (2004)  series by three piecewise linear segments that represent the initial rapid decline, the near-constant variation in recent times, and the intermediate transition period (Figure 2b).
We vary core properties individually to elucidate their influence on the snow regime. This has the potential to produce an inconsistency since both W04 and L14 used a particular core model, which produced a particular time-series of that is compatible with their time-series of .
To mitigate against this effect we set the initial CMB temperature, , to be close the original values. For W04 = 2400 K, while L14 did not quote a value and so we vary around the W04 value. We find that deviations in after 4.5 billion years of evolution differ by < 50 K from the original values in the majority of models. The paucity of independent observational constraints leads to some interdependencies between estimates of interior structure properties (e.g. assuming a temperature profile in order to estimate the density profile), but in this initial exploration we vary each parameter independently.
Values of density , CMB radius and CMB pressure (Table A1) are selected from W04 and also from a recent detailed analysis of the Martian interior (Rivoldini et al. 2011) that produced a range of models constrained by moment-of-inertia and 2 Love number data. Rivoldini et al. (2011) find that varies by only 5 − 10% across the Martian core and so there is little error in taking constant. These values determine the structure of the Martian core, i.e. the radial profiles of gravity, gravitational potential and hydrostatic pressure (Figure 2). For each model the pressure scale constructed in this manner is used to establish the temperature profiles discussed below.
The Martian core is thought to be composed primarily of an iron-sulfur alloy (Dreibus and Wanke 1985) and this simple chemistry has been used in almost all thermal history models to date (Breuer and Moore 2015). A more complicated core chemistry might be expected since the high temperatures achieved during the early stage of core formation may have facilitated the incorporation of Si, O, Ni, P, O, and H into liquid iron (Tsuno et al. 2007). The Martian core is expected to contain a negligible amount of Si (Sanloup and Fei, 2004) and O (Tsuno et al. 2007;Rubie et al. 2011), while Ni has only a minor effect on phase equilibria of Fe-S (Stewart et al. 2007). The amount of phosphorus in the Martian core is thought to be ten times that of Earth's core (0.16 vs. 0.02 wt% P2O5) (Dreibus and Wanke, 1985;Hart and Zindler, 1986). Both its abundance and the magnetic transitions of P-bearing phases may influence density distribution in the core (Gu et al., 2014), but we do not consider this additional complexity here. The density and melting point may also be lowered by the presence of hydrogen, though its content in the core is poorly constrained. In the absence of sufficient constraints regarding core equilibrium chemistry or a suitable theory for the melting point depression in ternary (or higher order) mixtures we model the evolution of a Fe-S mixture. The initial sulfur concentration is varied between 10 and 15%, which is within the estimates of previous models of the composition of Mars (e.g. Dreibus and Wanke, 1985;Taylor, 2013).
The adiabatic temperature is parameterized by the equation = 0 (1 + 0.02 ), which fits the published profiles of Williams and Nimmo (2004) and Fei and Bertka (2005). The coefficient 0 varies as the core cools. Note that the adiabatic gradient / is proportional to and therefore decreases as the core cools (Figure 2c).

Results
Input parameters for all models are listed in Table A1 and diagnostics are presented in Table A2.
We focus first on a model that uses the default parameters of W04 except for the melting curve and initial S concentration, which are set using the S07-10.6 profile. In this model the dynamo cessation time = 459 Myrs (~4 Ga), defined by falling below zero, while the present-day CMB temperature = 1822 K (Table A2, highlighted in bold); both values are very close to the original solution obtained by W04. Figure 3 shows profiles of temperature, solid fraction and S concentration for this model. Approximately 3.2 billion years into the evolution, falls below at the CMB and an iron snow layer begins to form and grows to 146 km by the present day. The solid fraction remains below ≈ 0.2 %, consistent with the modeling assumptions and with a recent model of iron snow in Ganymede's core (Rückriemen, et al., 2015), though the profiles of do not exhibit the curvature obtained by Rückriemen et al., (2015) near the base of the layer, which appears to stem from the different methods used to estimate / . The sulfur concentration increases across the snow zone by almost a factor of 1.5 at the present day, which arises partly due to the decreased S concentration in the deep core as Fe remelts and partly because of the enrichment in S with radius required to keep the snow zone on the liquidus. Figure 4 shows the contributions of individual terms to the energy and entropy balances for the model in Figure 3. The latent heat terms and make an order of magnitude larger contribution to the energy budget than the gravitational energy terms and in agreement with the study of Rückriemen et al. (2015), while is negligible. and almost balance since the same amount of mass is produced and destroyed; the small difference arises because the latent heat coefficient varies with depth. Therefore, these terms have little impact on the cooling rate at the onset of snow formation. The smallness of and reflects the slowness of compositional changes because the cooling rate is low and is a weak function of at conditions corresponding to the upper region of the Martian core.
The high thermodynamic efficiency of and means that the corresponding entropies are comparable to and . Nevertheless, the overall entropy produced from the formation and remelting of iron snow is small (Figure 4) and the dynamo does not restart as long as the snow zone remains relatively thin. The dynamo only restarts when the entropy produced by gravitational energy release due to the remelting snow, , which is proportional to the snow zone volume and growth rate [equations (17) and (20)], overcomes the conduction entropy . Rückriemen et al. (2015) inferred that dynamo action arose in their models of Ganymede. This finding is not comparable to our results since they used scaling laws to assess the onset and maintenance of dynamo action, rather than the entropy formulation employed here. Table A2 shows that the dynamo restarts when − > 400 km in our suite of models. Figure 5 shows a solution obtained with the same parameters as the model in Figure 3, except with a lower value of (highlighted in italics in Table A2). Lowering reduces ̃ [equation (24)] and therefore leads to faster cooling at early times and an older snow zone. The effect of decreasing from 7000 kg m -3 to 6000 kg m -3 is significant, which might partly reflect the lack of feedback on due to changes in in our model; however, decreasing also decreases the difference in gravity, gravitational potential, pressure, and adiabatic temperature across the Mars' thermal history.
The constraint that the Martian dynamo cannot restart (negative at the present-day) places a nominal upper bound on the thickness of the present-day snow zone of ∼ 400 km based on the limited model set available (Table A2). Occasionally, models with thick snow zones can produce thin layers below the CMB where exceeds the eutectic composition of ≈ 16 wt % at ≈ 20 GPa (Stewart et al. , 2007). The dynamics of this scenario are not included in our model, but it would produce light solid that floats to the CMB, thus reducing the estimates of gravitational energy compared to our calculations. However, the fact that such layers are very thin suggests that the associated entropy reduction would not prevent the dynamo from restarting. Figure 6 demonstrates the influence of parameter variations on the snow layer for models. Here the label for each symbol denotes the single quantity that was changed compared to the default model, which used the parameters highlighted in Table A1. The inferred dynamo cessation time ( ) is relatively insensitive to changes in , , and , but is very sensitive to changes in ; in Figure 6 the dynamo fails too late with = 30 W m -1 K -1 and too early with > 50 W m -1 K -1 . Aside from these cases all models in Figure 6 match the inferred dynamo cessation time, produce present-day CMB temperatures well above the eutectic value (Table A2) of 1300 − 1500 K at Martian CMB pressures (Rivoldini et al., 2011), and produce thin snow zones consistent with geodetic observations that suggest a predominantly liquid present-day core (Yoder et al., 2003). The iron snow regime is less likely to emerge for larger and certain time-series, which both cause the core cooling rate to decrease, though we obtained snow zones with all and values tested. Our results predict a strong sensitivity to ; however, this may be an artifact of the model assumption that does not change with .
The crucial parameter is the initial S concentration 0 , which determines the melting temperature and therefore strongly influences the initial difference between adiabatic and liquidus temperatures at the CMB Finally we consider whether iron snow zones arise using the preferred interior structure model of Rivoldini et al. (2011) and other default parameter values in Table A1. These runs use 0 = 0.142 and the S07-14.2 melting curve (Table A2). The high values of 0 and do not favor iron snow formation, but we do find relatively thin present-day snow zones in models with ≈ 2250 K, approximately 150 K below the value used in W04. This value still suggests a core that was initially superheated with respect to the mantle, consistent with the original modeling assumptions, and we have not attempted to 'optimize' our solution through a systematic parameter search as the uncertainties in several key variables do not warrant such a procedure. This model has a relatively late dynamo cessation time of 3.6 Ga, but increasing to 50 W m -1 K -1 , which is well within uncertainty, provides an acceptable value of 4 Ga while leaving the snow zone evolution unchanged.

Application of the snow model to the Martian core
Relaxing the fast melting approximation (i.e. incorporating departures from phase equilibrium) introduces additional terms and equations into the slurry theory and drastically increases the complexity of the constitutive relations (Loper and Roberts, 1977). These additional effects require macroscopic parameterizations of microscopic processes (Loper, 1992) that are poorly understood and the resulting terms are hard to estimate for planetary cores. While the overall influence of fast melting is hard to quantify, we might expect that the effects may not be significant as long as the relaxation to phase equilibrium occurs on timescales that are much shorter than the long timescale of interest (Gyrs). Incorporating the effects of a multi-component solid phase also significantly complicates the theory by requiring that the history of individual particles is modeled. At present we believe both approximations are sensible compromises for modeling the long-term behavior of snow layers in planetary interiors. Rückriemen et al (2015) used scaling laws with simple assumed particle sizes and geometries to argue that the fast remelting approximation is appropriate for modeling iron snow in Ganymede's core. Solomatov and Stevenson (1993) analyzed the conditions required to perpetually suspend particles in a magma ocean, but our model does not predict quantities such as the convective velocity needed to apply their theory. If some of the solid particles remain suspended in the snow zone on long timescales the latent heat released on freezing, , will exceed the latent heat absorbed on remelting, . Since is released close to the CMB it has a low thermodynamic efficiency factor, suggesting a reduction in entropy available to power the dynamo compared to the fast remelting case considered in this paper. It therefore appears that relaxing the fast remelting assumption would not significantly change our results, though hydrodynamic simulations of slurry dynamics are needed to test the veracity of this statement.

Thermal stratification at the CMB
The demise of the Martian dynamo is signified by the Ohmic dissipation dropping below zero.
However, ≥ 0 by definition and so negative values indicate an inconsistency in the modeling assumptions. The fact that < for most of the evolution suggests that thermal stratification ensues and the temperature profile deviates from the assumed adiabat profile near the top of the core in order to balance the entropy budget with = 0 after the dynamo fails. If = 0 prior to snow zone formation, the gravitational entropy terms , > 0 (Figures 4 and   5) that arise during snow zone growth would make > 0 and potentially restart the dynamo.
Strictly, must exceed some minimum value, denoted , for dynamo action to occur. is hard to estimate because it depends on the magnetic field morphology in the core, including the small-scale fields and the field components that remain inside the core, neither of which can be observed. Using just the observable field at the CMB gives ~10 6 MW/K (Gubbins, 1975;Backus et al., 1996), similar to the values of and in our models (Figures 4 and 5). The real value of is likely to be higher than this (Nimmo, 2015), suggesting that snow zone growth would not restart the dynamo. Since the dominant contributions to come from small-scale magnetic fields inside the core it may be possible that some field generation accompanies snow zone formation but produces an extremely weak signal at the planet's surface.
As discussed in Section 2, the thermally stable layer receives more heat through its base than can be removed at the CMB. Thus, the layer must heat up and the CMB temperature should be higher than predicted by our model, raising the question of whether the snow zone would still form. To address this issue we must first determine the relevant equations governing temperature variations in the conducting region. The temperature equation in a slurry ignoring pressure, radiogenic, and dissipative effects and assuming constant material properties is (Loper and Roberts, 1987) = ∇ 2 + ∇ ⋅ + where < 0 is the downward flux of solid material. The last two terms represent the total rate of change of solid mass per unit volume. A detailed analysis is complicated because depends on the size and distribution of solid particles. However, we observe that the fast melting approximation requires that solid freezes out quickly ( > 0) while fast remelting requires that solid falls out of the layer quickly (∇ ⋅ < 0) compared to the long timescale over which the temperature is changing. Since all solid leaves the layer after each timestep, on this long timescale the last two terms are expected to cancel out, leaving a standard diffusion equation for the temperature.
To estimate the temperature difference between an adiabatic and thermally stratified region , we first consider an infinite half-space with prescribed time-independent subadiabatic heat-flux at the boundary and zero initial temperature (corresponding to no departure from an initial adiabatic profile). In this case, the solution to (27) without solid ( = = 0) gives a boundary temperature 0 that can be written (Carslaw and Jaeger, 1959) where = / is the thermal diffusivity. In Figure 4, a thermally stable layer starts to grow at time = 400 Myrs into the calculation and a snow zone formed at approximately = 3.2 billion years, giving = − = 2.8 Gyrs. At time , = 0.257 TW and = 0.875 TW, corresponding to the strongest subadiabatic conditions (Figure 4).
With these values equation (27) shows that thermal conduction increases the CMB temperature by Δ = 0 ( ) − 0 ( ) ≈ 100K over 2.8 Gyrs above the value predicted from cooling on an adiabat. Using values for other models that match the cessation time for the Martian dynamo inferred from magnetic observations (Table A2) gives Δ = 100 − 250 K.
The analytical expression (27) ignores the effects of spherical geometry, finite stable layer thickness and temporal changes in CMB heat flow. We account for these effects by numerically solving the 1D conduction equation / = −2 / ( 2 / ) in a spherical shell of thickness . Using the time-series of in Figure 4 and an initial adiabatic from this run at time we obtain Δ = 70 − 140 K for 100 ≤ ≤ 700 km. The cooling rate in our models is 70-150 K Gyrs -1 at the time of snow zone formation and the snow zones form 0.5-1.5 Gyrs before present, suggesting that snow zone formation would be delayed until the recent past.
However, all of these Δ values are over-estimates because they ignore movement of the stable layer interface and omit the reduction in core cooling rate induced by stratification and by entrainment due to the underlying convection. We conclude that thermal stratification could delay, but not prevent, the onset of snow formation, though a more complete model of these effects is clearly needed.
The adiabatic temperature profile used in our calculations ignores the effect of a stabilizing compositional gradient, the second term in the relation = − − ̅ . The first term, calculated directly from the models, is ≈ 0.6 − 1 K km -1 at the CMB.  (Table A1) and / ≈ 4 × 10 −7 m -1 from Figure 2 means that the second term is ≈ 0.2 K km -1 . This is an overestimate since / depends on / in the model as discussed above, suggesting that the 'dry' adiabat assumed in the modelling is a good approximation to the 'wet' adiabat that includes compositional variations.
Departures from an adiabatic temperature profile affect the energy budget mainly through the term since this involves an integral over . To quantify the effect, we consider for simplicity a linear subadiabatic profile in the top 100 km of the present-day core with the CMB temperature 140 K above an adiabat, corresponding to the most extreme estimates above. The resulting 5% decrease in produces a change in the cooling rate of ~1 K Gyr -1 . Changes in sulfur concentration and solid fraction will also decrease in the presence of thermal stratification as the terms are proportional to −1 (equations (22) and (23)), but the overall effect on core cooling rate, and hence snow zone growth rate, is very small compared to other uncertainties in the calculation. The term in the entropy balance is reduced by a greater amount that , but this effect is countered by a reduction in since the conduction profile is shallower and hotter than an adiabat, resulting in a small change to the predicted dynamo entropy. Even weaker effects are predicted at earlier times or for younger stable layers. These simple calculations suggest that the assumption of neglecting departures from the adiabat in the energy-entropy balance is justified.
Finally, we expect that the presence of thermal stratification would reduce our estimates of the gravitational energy generated by migration of solid within the snow zone, though our calculations suggest that this term makes a negligible contribution to the overall energy and entropy budgets when the snow zone is only a few hundred km ( Figure 4). Nevertheless, the current parameterization of is simple at best and would benefit greatly from new experimental and/or numerical studies.

Conclusions
The presence of an approximately 100-km-thick snow layer at the top of the Martian core is consistent with the planets' magnetic history and available observational constraints on its core structure, temperature and composition. Snow zone nucleation is favored for lower initial sulfur concentrations and core temperatures and for smaller core sizes. Snow zones that grow to ≈ 400 km produce enough gravitational energy to restart the dynamo, suggesting that this is an upper limit on the layer depth in the Martian core.
Future work simulating slurry dynamics with and without thermal stratification should enable improved parameterizations of the thermal and compositional profiles in the snow zone and the gravitational energy terms in the energy balance and therefore mitigate the relaxation of some assumptions invoked in this study. Future parameterized models could also include coupled coremantle evolution. Considering the core in isolation has allowed us to focus on snow zone dynamics, divorced from the complexities and uncertainties in mantle evolution modeling, but at the expense of being restricted to a narrow range of CMB heat flow time-series and initial core temperatures. In particular, if solutions satisfying the available constraints can be obtained with lower initial CMB temperatures it will be possible to obtain thicker present-day snow zones than we have found.
Snow layers would not support seismic shear waves owing to the spatially dispersed nature of the solid phase, but could affect the core density. If these differences can be detected by observations from future spacecraft missions, it will provide profound inside into the thermochemical evolution of the Martian interior. Table A1. Input parameters used in the thermal history model. Gravity , gravitational potential and pressure are derived from the density assuming hydrostatic balance. Density is assumed depth-independent as interior structure models predict only 5 − 10% variation across the Martian core (Rivoldini et al. 2011); the constant value 7211 kg m -3 was sometimes used instead of the Williams and Nimmo value of 7011 kg m -3 as this accounts for the increase of with depth in their model. The thermal expansion coefficient , heat of reaction coefficient ̅ and volume change on melting Δ that appear in the governing equations are not included because the terms in the energy balance in which they appear are small enough to neglect. The latent heat is = Δ . Bold indicates the default value when multiple values have been used. Here W04 refers to Williams and Nimmo (2004), F05 is Fei and Bertka (2005) and S07 is Stewart et al. (2007).