RESEARCH ARTICLE

History and Future of the Martian Dynamo and
Implications of a Hypothetical Solid Inner Core

Key Points:

Douglas J. Hemingway1

10.1029/2020JE006663• W
 e develop a new internal structure
and thermal evolution model tostudy dynamo regimes on rocky
planets, using Mars as a firstexample• An early, thermally driven Martian
dynamo could be followed by a later,
compositionally driven dynamo,
either in the past or the future• Detection of a solid inner core on
Mars could imply that the Martian
dynamo will reactivate at some pointin the future
Supporting Information:
Supporting Information may be foundin the online version of this article.
Correspondence to:
D. J. Hemingway,dhemingway@carnegiescience.edu
Citation:
Hemingway, D. J., & Driscoll, P. E.
(2021). History and future of theMartian dynamo and implications of ahypothetical solid inner core. Journalof Geophysical Research: Planets,
126, e2020JE006663. https://doi.
org/10.1029/2020JE006663
Received 21 AUG 2020
Accepted 2 MAR 2021
Author Contributions:Conceptualization: Douglas J.
Hemingway, Peter E. Driscoll
Formal analysis: Douglas J.
Hemingway
Investigation: Douglas J. Hemingway,
Peter E. Driscoll
Methodology: Douglas J. Hemingway,
Peter E. Driscoll
Software: Douglas J. Hemingway
Supervision: Peter E. Driscoll
Visualization: Douglas J. Hemingway
Writing – original draft: Douglas J.
Hemingway
Writing – review & editing: Peter E.
Driscoll© 2021. American Geophysical Union.
All Rights Reserved.

HEMINGWAY AND DRISCOLL

and Peter E. Driscoll1

1

Earth & Planets Laboratory, Carnegie Institution for Science, Washington, DC, USA

Abstract Although Mars does not possess a global magnetic field today, regions of its crust arestrongly magnetized, consistent with an early dynamo, likely powered by rapid heat flow from the core. If
the core is undergoing crystallization, the associated compositional changes would provide an additional
mechanism for driving convection—probably the dominant driver for Earth's dynamo today. This raises
the question: does the lack of a global dynamo field on Mars suggest the absence of a partially crystallized
core? More generally, what is the range of possibilities for the history and future of the Martian dynamo
and which scenarios could be ruled out by the presence or absence of a solid inner core? Here we develop
a new internal structure, thermal evolution, and buoyancy flux model to investigate the conditions
under which the Martian core could experience compositionally driven convection, either in the pastor the future. We show that the presence of a partially crystallized core is compatible with the lack of a
dynamo today but that such a scenario implies the Martian dynamo could reactivate at some point in the
future. We find that top-down core crystallization (iron snow) requires weak light element partitioning,
introduces limited buoyancy flux, and is unlikely to be effective at driving convection. Our model
demonstrates how sulfur content & partitioning and core conductivity & expansivity determine which
dynamo regimes are possible, can help in assessing implications of future observations relating to theMartian core, and forms the basis for further comparative study across rocky planets.
Plain Language Summary Although Mars does not possess a magnetic field today, parts of

its crust are strongly magnetized, suggesting an early field, likely powered by rapid heat flow from thecore. For Earth, an additional, and possibly dominant, mechanism for powering the magnetic field arises
from core crystallization. Does the lack of a magnetic field on Mars therefore imply its core is entirely
liquid and not undergoing freezing? More generally, what is the range of possibilities for the history andfuture of the Martian magnetic field and which scenarios would be ruled out by confirming the presenceor absence of a solid inner core? Here we develop a new model to investigate the conditions under which
the Martian magnetic field could be driven by core freezing, either in the past or the future. We showthat a partially crystallized core is compatible with the lack of a magnetic field today but that this would
imply the Martian magnetic field could reactivate in the future. Our model demonstrates how sulfur
content & partitioning and core conductivity & expansivity determine which scenarios are possible, canhelp in assessing the implications of observations relating to the Martian core, and forms the basis forcomparative study across rocky planets.

1. Introduction
Planetary magnetism is a fundamental but imperfectly understood phenomenon that has the power toinform our understanding of the deep interiors and thermal histories of planetary bodies of all kinds(e.g., Connerney, 2015; Stevenson, 2003; Stevenson et al., 1983). For terrestrial planets, convective motions in the electrically conductive liquid metallic core can drive a global scale self-sustaining magnetic
field, in a process known as a magnetohydrodynamic dynamo. This field shields the planet's atmosphereand surface from the impinging solar wind and imparts thermoremanent magnetization to rocks as theycool at or near the planet's surface, preserving a record of how the dynamo field may have changed overtime (e.g., Purucker & Whaler, 2015; Weiss & Tikoo, 2014). Paleomagnetic evidence has established that
the Earth's dynamo field has been active for at least the last ∼3.5 Gyrs, and possibly longer (e.g., Tarduno et al., 2010, 2015). By contrast, the case of Mars is interesting because it has no global dynamo field
operating today and yet magnetometer and electron reflectometer data from the Mars Global Surveyor

1 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

and MAVEN missions have established that large portions of the crust—particularly the ancient southern highlands—are extremely strongly and coherently magnetized (Acuña, 1999; Connerney et al., 2005;
Langlais et al., 2019; Stevenson, 2001), strongly suggesting that a dynamo was active early in Mars' history, but has since shut down.
The lack of magnetization in large impact basins like Hellas, Argyre, and Isidis, whose formation shouldhave produced significant local heating and melt, is generally interpreted as evidence that the early Martian dynamo had shut down by the time these basins formed, roughly 500 Myrs after the formation ofMars (Lillis et al., 2008, 2013; Vervelidou et al., 2017), or even that the formation of the basins themselves
could have initiated the shut down (Roberts et al., 2009). Several authors have previously studied thethermal evolution of Mars (e.g., Hauck & Phillips, 2002; Schubert et al., 1992; Zuber et al., 2000), withsome studies focusing on the evolution of the Martian dynamo and what controls its timing (e.g., Dietrich& Wicht, 2013; Hori & Wicht, 2013; Ke & Solomatov, 2009; Landeau & Aubert, 2011; Nimmo & Stevenson, 2000; O'Rourke & Shim, 2019; Stevenson et al., 1983; Williams & Nimmo, 2004). Williams and Nimmo (2004), for example, showed that if the initial core-mantle temperature difference was sufficiently
large, strong heat flow across the core-mantle boundary (CMB) could have driven thermal convection in
the liquid core for a limited period before the core became thermally stably stratified, naturally shutting
down convection and the dynamo. However, it has also been suggested that a weak dynamo persisted
or reactivated after the formation of the large basins (Lillis et al., 2006; Mittelholz et al., 2020; Schubert
et al., 2000), raising the question of what mechanisms and what conditions could permit such different
dynamo histories.
The primary mechanisms for driving core convection are generally thought to be thermal (i.e., due to rapid
core-mantle heat flow) and compositional (i.e., arising from compositional changes associated with core
crystallization), though other mechanisms could be important as well. If the core is initially much hotter
than the mantle, thermal convection will be effective at early times, but its effectiveness diminishes overtime as the core and mantle cool. Compositional convection becomes more important later when parts ofthe core begin to drop below the melting temperature. The Earth's core is believed to have been entirely
liquid for an extended period until its central temperature fell below the high pressure melting curve foriron, which may have occurred as recently as ∼1 Gyr ago (Davies et al., 2015; Driscoll & Bercovici, 2014;
Gomi et al., 2013; Labrosse et al., 2001). Since then, a solid inner core has been growing, currently having aradius about 35% of the core's total radius. As the core crystallizes, there is a tendency for lighter elements topartition preferentially into the liquid rather than being incorporated into the solid inner core. This meansthat light element rejection introduces a source of buoyancy at the inner core boundary which, for thepresent-day Earth, is likely an important—if not dominant—source of buoyancy driving convection in theliquid outer core.
At Mars, geodetic observations and interior modeling suggest that the core remains at least partially liquid
today (Rivoldini et al., 2011; Van Hoolst et al., 2003; Yoder et al., 2003), but it is not known whether any
portion of the Martian core is solid. The InSight lander, now operating on the Martian surface, includes
a seismometer (SEIS; Banerdt et al., 2020) and a radio science experiment (RISE; Folkner et al., 2018)
designed to study the deep interior, offering the possibility of confirming the presence or absence of asolid inner core within Mars' otherwise liquid core. Anticipation of such results is one of the primary
motivations for the present study. Similarly, analysis of seismic wave propagation could help to constrain
the outer core's density, which could in turn help to better constrain its sulfur content—a key parameter
whose importance will be discussed extensively throughout this work. Finally, a heat probe such as theone on InSight (HP3; Spohn et al., 2018) could potentially measure the conductive heat flow through theMartian crust. As we will see, surface heat flow, which is currently not well constrained, relates directly
to the temperature structure and thermal evolution of the Martian mantle, and may therefore indirectly
constrain core properties.
The mechanism of compositionally driven convection, which is thought to be so effective for Earth, thus
raises the question of whether or not a solid inner core on Mars is compatible with the absence of a dynamo
field at present—a key question in this study. In other words, under what conditions, if any, is it possible
to have a crystallizing core without driving convection? If it becomes possible to detect a solid inner coreand estimate its radius, what would this tell us about the history and future of the Martian dynamo? EvenHEMINGWAY AND DRISCOLL

2 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

if compositionally driven convection is not operating today, is it possible or likely that it has operated at
any point in the past or that it will operate at some point in the future, and what determines which of thesescenarios is most likely? It has also been suggested that, owing to the much lower pressures in the Martian
core, crystallization may occur from the top-down rather than the bottom-up—two regimes that could havevery different effects on convection and dynamo generation (Breuer et al., 2015; Davies & Pommier, 2018;
Hauck et al., 2006; Stewart et al., 2007). What factors determine the crystallization regime and when various
types of convection could operate? How different would the results be with a different internal structure, orfor a body that is larger or smaller than Mars?
To answer these questions, we have developed a new internal structure and thermal evolution modelthat tracks core and mantle heating and cooling resulting from radiogenic heating, heat flow across thecore-mantle boundary, conductive and eruptive heat loss to the surface, and the effects of gravitational energy and latent heat flow associated with core crystallization. We assume a well-mixed core with an adiabatic
temperature profile and we track the relationship between this core adiabat and the composition-dependent
melting curve to determine when parts of the core may crystallize, whether from the bottom-up (i.e., Earthlike solid inner core growth), or from the top-down (i.e., as “iron snow”). We use an iterative procedure
involving equations of state to compute a simplified, but physically reasonable radial density structure.
Though not the focus of this first study, this will allow us to investigate the effects of varying the body's size
and internal structure. Finally, we compute both thermal and compositional buoyancy fluxes to determine
how each may contribute to driving convection within the liquid part of the core. Apart from the specifics
of approximating the iron-sulfur melting curve appropriate for Mars core conditions, the model we describein Section 2 is not specific to Mars. We aim for a planet-independent model in part because we intend touse it in future work to study other bodies like the Moon, Mercury, Venus, or terrestrial-like exoplanets,
but also because this allows us to benchmark our model by verifying that it produces reasonable results forthe case of the Earth. In Section 3, we illustrate a series of possible Martian thermal evolution scenarios
and undertake a systematic exploration of the large parameter space, identifying what types of convection
could drive a dynamo, under what conditions, and when. In Section 4, we summarize the various possible
dynamo regimes, the factors that control their timing, and what can be concluded from future observations
such as confirming the presence or absence of a solid inner core. Finally, we offer brief concluding remarks
in Section 5.

2. Methods
This section details the various components of our model. Because we use thermal and compositional buoyancy flux as a proxy for whether or not convection is possible, we begin in Section 2.1 with a discussion of
buoyancy flux and how it is calculated. We show how the compositional buoyancy flux depends on thedensity anomalies arising from unequal partitioning of light elements into the liquid and solid phases upon
freezing, and we describe how those density anomalies are calculated. Section 2.2 describes our method
for using a Vinet equation of state and an iterative procedure to obtain a self-consistent internal structure
describing the radial profiles of density, pressure, bulk modulus, compression, and thermal expansivity. We
use this method to build interior models for both Earth and Mars, demonstrating that the model can be used
for an arbitrary planet, with a minimum of body-specific parameters, permitting future comparative studyof internal structure on convection and dynamo regimes. Section 2.3 describes the internal temperature
structure and our approach for approximating the adiabatic core temperature profile (Section 2.3.1) and
its relationship to the melting curve, which depends on the core's sulfur content (Section 2.3.2), to identify when and where parts of the core may begin to crystallize. Core cooling is largely controlled by CMB
heat flow, which is a function of the mantle temperature structure, described in Section 2.3.3. Section 2.3.3
includes our procedure for obtaining an internally consistent mantle temperature profile, consisting of aconvecting region between two conductive thermal boundary layers, whose thicknesses are functions of thetemperature-dependent mantle viscosity. Finally, Section 2.4 details the energy budget calculations, which
include radiogenic heating, CMB heat flow, gravitational and latent heat terms associated with core crystallization (Section 2.4.1), conductive heat flow across the thermal boundary layers at the base and top of the
mantle and possible eruptive heat loss to the surface (Section 2.4.2), and how this energy budget is used tocompute the body's thermal evolution (Section 2.4.3).
HEMINGWAY AND DRISCOLL

3 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

2.1. Buoyancy Flux
Dynamo generation depends on movement in the electrically conductive liquid portion of the core, which
we take to arise from some combination of thermal and compositional convection (for this study, we neglect
other possible sources, such as mechanical stirring, e.g., Dwyer et al., 2011; Le Bars et al., 2011). We quantify
the introduction of convection-inducing thermal or compositional density perturbations in terms of thebuoyancy flux, defined as (e.g., Olson et al., 2017)F 
g  u(1)cwhere g is the local gravitational acceleration, u is the advective velocity (i.e., carrying the density perturbation), ρc is the core fluid's background density, and where the density perturbation is, in general, given by    c  T     
(2)where α is the thermal expansivity and β is the compositional expansivity (T′ represents a temperature perturbation and χ′ represents a perturbation in chemical composition—e.g., due to an influx of fluid enrichedin light elements). Buoyancy flux describes the rate at which buoyancy is being introduced into the liquid
portion of the core (in terms of force per unit mass times a volume flux, yielding units of m2/s3). As shownby Lister (2003), although quite different in form and emphasizing different terms, this buoyancy flux approach is fundamentally equivalent to the approach of tracking entropy production, at least as it relates tothermal and compositional convection. Rather than relating the dynamo strength to the total ohmic dissipation in the core, however, we evaluate whether or not thermal and/or compositional density anomalies
introduce sufficient buoyancy to drive convection.
Thermal convection occurs when heat is removed from the core more rapidly than it can be conducted along
an adiabatic temperature profile (i.e., when the temperature profile is superadiabatic), leading parcels of core
fluid near the core-mantle boundary (CMB) to cool, become dense, and sink. We express this effect in terms

of the thermal buoyancy flux, FT, obtained by requiring the thermal energy advection (V c cc    T u  dV )to equal the superadiabatic heat flow, yielding (e.g., Driscoll & Bercovici, 2014, equation 36)

 g
FT  c cmb qT
(3)
c cc
where αc and ρc are the thermal expansivity and mean density of the fluid outer core, cc is the core's specificheat capacity, gcmb is the gravity at the CMB, and where qT is the superadiabatic heat flux, given byqT qcmb  qad(4)where qcmb is the heat flux across the CMB and qad is the adiabatic heat flux, given by dT 
qad   kc (5) dr adwhere kc is the core thermal conductivity and the adiabatic temperature gradient is given by (e.g., Stacey, 1992, p. 304)
 dT 
 c gcmb
Tcmb
(6)

  
dr
cc

ad

where Tcmb is the temperature at the CMB. The CMB temperature and heat flux are functions of time andare obtained by considering the energy budget as described in Section 2.4.
A positive thermal buoyancy flux indicates that heat is being removed across the CMB sufficiently quickly(i.e., superadiabatically) to introduce a source of density anomalies that can drive convection. A negative

HEMINGWAY AND DRISCOLL

4 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

thermal buoyancy flux indicates that the CMB heat flux is small enough that the core can cool conductively
and will therefore remain thermally stratified (i.e., the temperature profile remains subadiabatic and there
is no thermal convection).
If the (primarily iron) liquid core is cool enough that parts of it begin to crystallize, the tendency forthe lighter elements to partition preferentially into the liquid phase becomes a source of compositional buoyancy. If the inner core is growing from the bottom-up (i.e., radially outward from the center),preferentially crystallizing iron, the rejection of some of the remaining lighter elements introduces a source of buoyant (light element-enriched) liquid at the inner core boundary. In this case, thedensity perturbation in (1) is the density contrast between the outer core (ρoc) and the rejected lightelement-enriched liquid (ρχ), which we write as Δρχ = ρoc − ρχ. As the solid inner core grows, thedensity of the remaining liquid outer core decreases as it is gradually enriched in lighter elements.
The decrease in the outer core density is initially extremely gradual, however, and can be neglected aslong as the inner core radius remains small compared to the total core radius. Because it simplifies themodel significantly, we adopt the approximation that the outer core's sulfur content remains constant,and thus ρoc and Δρχ likewise remain constant. Hence, we do not distinguish between the outer coredensity (ρoc) and the core's overall bulk density (ρc). (We discuss the consequences of this approximation below.)
If the mass fraction of light elements in the core is χS, the core's bulk density (neglecting compression,
which is discussed in Section 2.2) is given by1
1  S S


(7)

c

 Fe

S

where we use the subscript S because we assume sulfur to be the principal light element species in theMartian core (Franz et al., 2019; Gaillard & Scaillet, 2009), and where ρFe and ρS are the densities of the
major and minor species (assumed to be iron and sulfur), respectively. Neglecting the small density differences between the solid and liquid forms of both the major and minor species (i.e., assuming ρFe, sol = ρFe,3
3
liq = 6,980 kg/m and ρS,sol = ρS,liq = 1,819 kg/m ), we can approximate the light element-enriched liquiddensity using
1   S,liq  S,liq
1


(8)



 Fe

S

where, assuming small control volumes partition into equal parts solid and liquid by mass, it can be shownthat
2
S,liq S(9)1  DSis the mass fraction of light elements in the residual liquid, where DS is the solid-liquid phase partitioning
coefficient for the light elements (i.e., the mole fraction of light elements that partition into the solid phasevs. the liquid phase) (Zhang et al., 2020). Note that if sulfur partitions into the solid and liquid phases inequal proportions (i.e., if DS = 1), then ρχ remains the same as the bulk core density and Δρχ = 0; if sulfur
partitions preferentially into the liquid, then DS < 1 and ρχ < ρc so the residual liquid is positively buoyant.
Note that the above-described freezing behavior and the use of (9) should be limited to the subeutectic FeFeS system only, since our model does not currently accommodate sulfur content exceeding the eutectic(see Section 2.3.2).
For this inner core growth-driven compositional buoyancy flux, the rate at which the inner core radius is
growing, R ic, is the rate at which the density anomalies are being introduced in (1). Hence, this compositional buoyancy flux, scaled up to the CMB (r = Rc) to enable direct comparison with (3), is given by (e.g.,
Driscoll & Bercovici, 2014, equation 37)

HEMINGWAY AND DRISCOLL

5 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

2

Δ   Ric  (10)F  gicb Ric c  Rc 

where the gravity is taken to be linear across the core such thatR
gicb  gcmb ic
(11)
Rc

where the radius and growth rate of the inner core, Ric and R ic, given by (40) and (41), are functions of thetime varying temperature structure of the core and its relationship to the melting curve, as described in
Section 2.3.2.
As described above, our model neglects the outer core's gradually increasing sulfur concentration such thatρoc = ρc and Δρχ are treated as constants. Whereas this approximation is very good for the first few billion
years of inner core growth, its accuracy begins to decrease once the inner core radius exceeds roughly half
the total core radius (Figure S1). Specifically, our model does not account for the gradually decreasing
density of the outer core and the resulting increase in the density perturbation Δρχ/ρc (Figure S1a). As wediscuss in Section 2.3.2, the sulfur content also affects the melting temperature. Hence, another result ofneglecting the change in sulfur content is that our model can begin to overestimate the melting temperature
at the inner core boundary during the later stages of thermal evolution (Figure S1b). In other words, ourmodel does not account for the fact that the increasing sulfur content in the outer core further depresses
the melting temperature, slowing inner core growth at later times. Both of these effects influence (10) suchthat our results relating to compositional buoyancy flux should be interpreted with caution when the innercore radius exceeds approximately half of the full core radius. As we discuss in Section 4, however, the limitations of this approximation do not affect our overall conclusions.
As we discuss further in Section 2.3, it is also possible for the core to crystallize from the top-down rather than from the bottom-up, in what is known as the iron snow regime (e.g., Breuer et al., 2015; Hauck
et al., 2006). This occurs when the core temperature profile is steeper than the melting curve such thatthe central part of the core remains above the melting curve even as crystallization begins to occur in theregion above the radius at which the two curves meet, Rsnow, which we call the base of the snow zone. In
this case, as the core cools, iron-rich crystals forming above Rsnow are negatively buoyant compared to thesurrounding liquid (since the lighter elements partition preferentially into the liquid), and therefore sinkas iron snow. If the falling crystals are warmed enough to remelt, they mix into the liquid layer below,
introducing positive density anomalies and driving convection from above (e.g., Olson et al. 2017)—aprocess analogous to the buoyancy generated by exsolution of a minor species, leaving behind a negatively buoyant iron-rich residual liquid (e.g., Du et al., 2017; Mittal et al., 2020). The crystals will remelt at
a particular radius, Rremelt, where their temperature (assumed to follow the adiabat) reaches the melting
curve for the iron-rich crystals. In general, because the crystals are relatively depleted in light elements,
the relevant melting curve is elevated compared with that of the initial bulk liquid core. If the difference
between the two melting curves is neglected (e.g., Breuer et al., 2015; Rückriemen et al., 2015), thenRremelt = Rsnow, meaning the crystals remelt immediately below the base of the zone in which they are
forming. We instead compute a separate melting curve (see Section 2.3.2) for the crystals and find its
intersection with the adiabat to obtain Rremelt, which, in general, is smaller than Rsnow and may in fact bezero if the crystals never warm enough to remelt. Although the crystals may generate small eddies as theyfall, they are not expected to contribute to large-scale convection unless they can remelt and mix into theunderlying fluid (e.g., Olson et al., 2017).
If the falling crystals do remelt, they introduce density anomalies at Rremelt at the rate at which the snow
zone is growing,  R snow. The snow-driven buoyancy flux, once again scaled up to the CMB to enable directcomparison with (3) and (10), is thus
2Δsnow  Rremelt (12)Fs gremelt  R snowc  Rc 

HEMINGWAY AND DRISCOLL

6 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

whereR
gremelt  gcmb remelt
(13)
Rc

and where Δρsnow = ρsnow − ρc is the excess density of the iron snow compared with the surrounding liquid,
where we approximate the snow density using
1   S,sol S,sol
1


(14)

snow

Fe

S

the solid component counterpart to (8), where2
S,sol S(15)1  DS1

As with bottom-up inner core growth, the base of the snow zone and its rate of change are given by (40) and(41) and are functions of time as described in Section 2.3.2. Likewise, the location of remelting, Rremelt, isgiven by (40) but using a melting curve that corresponds to the iron-rich snow rather than to the bulk liquidcore. If the crystals do not remelt, then Fs = 0, and a solid inner core grows by accumulation of the falling
solids at the center of the core. Unlike the case of bottom-up inner core growth, in the case of falling andremelting iron snow, there is no net gradual change in the liquid core's sulfur content, and so the approximations we adopt in (10) do not apply to (12).
Even if the CMB heat flux is subadiabatic, such that the upper region of the core is thermally stably stratified, any ongoing core crystallization could continue to act as a source of buoyancy that could drive convection. However, the greater the thermal stratification at the top of the core, the smaller the compositionallyconvecting region will be, and the more difficult it will be to sustain a dynamo. Hence, as a simplification,and following Driscoll and Bercovici (2014), we will consider convection and dynamo action to be possible
only when the net buoyancy flux, resolved at the CMB, is positive. That is, wheneverFc  FT  F  Fs  0
(16)and not otherwise. We emphasize that this is a simplification and that even a positive net buoyancy flux
(Fc > 0) may not be a sufficient condition for dynamo generation. When Fc > 0, the resulting dipole moment
magnitude may be estimated from scaling laws such as (Olson & Christensen, 2006)
1/2

 
1/3
(17)
Mdipole  4 Rc3  c   Fc d  0

where μ0 = 4π × 10−7 Hm−1 is the vacuum permeability, γ ≈ 0.2 is the saturation constant for fast rotating
dipolar dynamos, and d is the thickness of the convecting layer (i.e., Rc − Ric). However, we will not focus
much on the strength of the convection-driven dynamo. Instead, we are concerned mainly with how thetiming of the rising and falling of Fc, and thus the timing of the onset and cessation of convection and dynamo action, is controlled by parameters such as the initial temperature structure, radiogenic heating, mantle
viscosity, thermal conductivity, light element partitioning, and core composition.
2.2. Density Structure
We model the planet's interior as a series of spherically symmetric concentric shells. The models consist ofan outer rocky mantle and an inner primarily iron core, which may be further partitioned into solid andliquid phases. Owing to compression at high pressures, we allow for radially variable density within eachlayer according to the Adams-Williamson equation

HEMINGWAY AND DRISCOLL

7 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

 2 r  g r 
d

(18)
dr
Ks  r 

where Ks is the isentropic bulk modulus, given by





K
KT 1    x    x  T  r (19)s r 

where KT is the isothermal bulk modulus, T is the temperature, α is the thermal expansivity, and γ is the
Grüneisen parameter—another thermodynamic property describing nonlinear aspects of the relationshipbetween temperature and crystal volume. KT is in turn related to the density by the third-order Vinet equation of state (Oganov, 2007; Vinet et al., 1989)   





KT K0 1  1   x1/3 1  x1/3 exp  1  x1/3 x 2 /3
(20)where the compression is represented by




x  0
(21)

(equivalently, x = V/V0, as it is expressed by Vinet et al., 1989 and Oganov, 2007) and where





3 


K0  1(22)2where K0 is the zero-pressure bulk modulus and K′0 is its derivative with respect to pressure. The thermal expansivity and the Grüneisen parameter are functions of the compression according to the parameterizations

  x    0 x T
(23)
and

  x    0 x 1
(24)based on high pressure experiments (Chopelas & Boehler, 1992; Merkel et al., 2000). These parameterizations were previously used by Driscoll and Olson (2011) but the sign on the exponents in our (23) and (24)are reversed from those in their equation 8 due to a typo in the text after their equation 6 in which their x is
inverted from the one previously defined by Oganov (2007) and in our (21) (see also Boujibar et al., 2020).
Whereas Driscoll and Olson (2011) and Boujibar et al. (2020) assumed the exponent in (23) to be δT = 3, weallow it to vary for different materials. The reference values, ρ0, K0, α0, and γ0, and the exponents, δT and γ1,
are material and temperature dependent; the values given in Table 2 have been selected to approximately
satisfy experimental data for the expected temperatures and compositions of the cores and mantles of Earthand Mars (e.g., Boehler et al., 1990; Chopelas & Boehler, 1992; Stacey & Davis, 2004; Vočadlo, 2015).
The temperature required by (19) is estimated by integrating r  g r  r dT T r 
(25)drKs  r 
subject to the boundary condition of an assumed temperature at the top of the mantle and the density andgravity profiles obtained as discussed immediately below. For purposes of modeling the planet's energy
budget and thermal evolution, however, we will define a distinct temperature structure as described in
Sections 2.3 and 2.4. Although the temperature profile obtained by integration of (25) differs from the temperature profile that emerges from the thermal evolution calculation (in part because the latter includes
a steep temperature change across the thermal boundary layers), by substituting an evolved temperature
profile into (19) instead of using (25), we verified that the effect on the computed radial density structure is
HEMINGWAY AND DRISCOLL

8 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

negligible (densities are everywhere within <1%). Hence, integration of (25) is sufficient for estimating thetemperature profile for purposes of computing the radial density structure prior to considering the planet's
thermal evolution, avoiding the computationally expensive process of having to continually update thedensity structure as the temperature profile evolves.
At a given radius, r, the enclosed mass, M  r , and the hydrostatic pressure, P  r , are given by integration ofdM   r  4 r 2
(26)
dr

and
dP
  r  g r 
(27)
dr

where the gravitational acceleration isGM  r (28)g r  r2

Although it follows that
2GM  r dg(29)
 4 G   r  
drr3

it is more straightforward to compute g  r  directly from (28) since M  r  is determined at every step of theintegration.
We obtain the interior structure by integrating Equations (18) and (25)–(27) inward from the outer edge ofthe body, in 5-m steps, subject to the boundary conditions of the known total mass and an assumed density
and temperature at the top of the mantle (e.g., Boujibar et al., 2020; Driscoll & Olson, 2011). At the transition between the mantle and the core, there is a discontinuity in density and elastic properties, which weestimate by requiring continuity in the pressure profile. Specifically, given the known value of Ks on the
mantle side of the boundary, we find the value of Ks for the core side that yields equal pressures on eitherside of the boundary, where the pressure is approximated as
P  K S  K0  / K0(30)and then we solve (19)–(21) to obtain the density at the top of the core before resuming the inward
integration.
In general, when the integration terminates at r = 0, the enclosed mass does not go precisely to zero, as isrequired for a self-consistent physical solution. To ensure that M  0   0, we carry out the above procedure
iteratively, making small changes to the assumed density at the top of the mantle until M  0   0 to withina factor of 10−5. Finally, we force the density to be constant over the innermost 2% of the body's radius andthen recompute the mass, pressure, and gravity profiles by integrating outward from r = 0, thereby ensuringthat M  0  and g  0  are precisely zero. As a result, the enclosed mass at r = R differs slightly from the true
total mass, but only by a factor of ∼10−5.

As an example, Figure 1 shows the simplified radial density structures we obtained for Earth and Mars
when we adopt the parameters given in Tables 1 and 2. The Martian core radius is uncertain but we adopt
the middle of the wide range obtained by Brennan et al. (2020)—although we note that at least one recent InSight-based study suggests the core may be larger (Deng & Levander, 2020). Similarly, analysis ofseismic wave propagation (e.g., via InSight) could aid in constraining the density structure, including theouter core density, which could in turn help to constrain its sulfur content. In principle, the Martian coreradius should be a function of its sulfur content because greater sulfur content implies a lower density core,

HEMINGWAY AND DRISCOLL

9 of 37

Journal of Geophysical Research: Planets
Table 1
Coarse Internal Structure of Earth and Mars
Parameter

Symbol (units)

Earth

Mars

Full radius

R (km)

6,371

3,389.5

Core radius

Rc (km)

3,480

1,745

Bulk density (kg/m3)

5,513.6

3,934.1

10.1029/2020JE006663

which in turn requires a slightly larger core radius in order to satisfy themass and moment of inertia constraints (e.g., Brennan et al., 2020). Suchchanges in the radial density structure do have some effects on relevant
model parameters, but they are small compared with uncertainties in
the reference values for those parameters. For example, we determinedthat higher sulfur content lowers the reference core density as well as
lowering core pressures such that compression, x, and therefore thermal
expansivity, α, are largely unaffected by changes in core sulfur content:
adjusting wt%S in the range 7 ± 7 affects α at the ∼10% level, an effectthat is overwhelmed by uncertainty in the reference value, α0.

2.3. Temperature Structure
The model planet's temperature structure can be divided into four regions: a well-mixed core, a conductive
lower mantle thermal boundary layer, a well-mixed (convecting) region of the mantle, and a conductive
upper mantle thermal boundary layer. An example model temperature profile is illustrated in Figure 2a
for a case approximating the present-day Earth, and a few hypothetical Martian temperature profiles areillustrated in Figures 2b–2d.
2.3.1. Core Temperature Structure
The adiabatic temperature profile across the core can be written (e.g., Labrosse et al., 2001)
 r2 

Tcore  r  Tcen exp   2 
(31)
 D 
N 


where Tcen is the temperature at zero radius and DN is a length scale given by
3cc
DN 
(32)
2 c cG

where cc is the specific heat capacity, G is the gravitational constant, ρc is the mean core density, and αc is the
mass-averaged thermal expansivity of the core, given by
1 Rc

c
 0   x  dM
(33)
Mc

where Mc is the total core mass and dM is the mass of each infinitesimally thin shell, and where   x  is given by (23) with x being a function of radius emerging from the density profile via (21). It can be shown thatthe radial gradient of (31) at r = Rc is equivalent to (6). Note that Labrosse et al. (2001) use the zero-radius
density, ρcen, in (32), which gives a similar result, provided the density does not vary too greatly across thecore. However, when ρcen is used in place of ρc, the equivalence with (6) is not strictly preserved.
Because we often prefer to reference the CMB temperature, we can also write (31) as R2  r 2 Tcore  r   Tcmb exp  c 2 
(34)
 D

N



The mass-averaged core temperature is given by
1 Rc

Tc
 0 Tcore  r  dM  r 
(35)
Mc

Substituting from 34, we can write

HEMINGWAY AND DRISCOLL

10 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Table 2Equation-of-State Parameters Used in Radial Structure Computation, Assuming a Perovskite-Like Mantle and Liquid
Iron Core for Earth, and a Peridotite-Like Mantle and Liquid Iron-Sulfur Core for Mars
Parameter

Symbol (units)

Earth mantle

Earth core

Mars mantle

Mars core

Zero-pressure density3ρ0 (kg/m )

3,900

6,900

3,226

6,000

Zero-pressure bulk modulus

K0 (GPa)

200

125

128

125

Bulk modulus gradientK′0 (n/d)

4.0

5.5

4.2

5.5

Reference Grüneisen parameter

γ0 (n/d)

1.0

1.6

0.99

1.6

Grüneisen parameter exponent

γ1 (n/d)

1.4

0.92

2.1

0.92

Reference thermal expansivity

α0 (K−1)

4 × 10−5

5 × 10−5

4 × 10−5

4 × 10−5

Thermal expansivity exponent

δT (n/d)

4.5

3.0

4.9

4.5

Note: Values are similar to those of Driscoll and Olson (2011) and correspond approximately with various experimental
results (Boehler et al., 1990; Chopelas & Boehler, 1992; Stacey & Davis, 2004; Vočadlo, 2015).

Figure 1. Radial density structures (a), pressure profiles (b), gravity profiles (c), and thermal expansivity profiles (d) for Earth (blue) and Mars (red) obtained
as described in Section 2.2. The dashed lines in (a–c) show PREM (Dziewonski & Anderson, 1981) as a point of comparison for our simplified model Earth. The
density, pressure, and gravity profiles (a–c) are shown as a function of radius but the thermal expansivity profile (d) is shown as a function of pressure to enable
more direct comparison with experimental results (e.g., Boehler et al., 1990; Chopelas & Boehler, 1992; Vočadlo, 2015).

HEMINGWAY AND DRISCOLL

11 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 2. Example internal temperature structures for Earth (a) and Mars (b–d); the Martian core is emphasized in (c and d). Open circles and labels indicate
the CMB temperature (Tcmb), the midmantle temperature (Tm), and the surface temperature (Tsurf ). Pale red shading highlights the two conductive thermalboundary layers (TBLs) at the top and bottom of the mantle (UM, LM, respectively). Solid lines represent the temperature as a function of radius; dashed lines
represent the melting temperature across the core; dotted line represents the remelting temperature for falling iron snow (visible only in d—a case in which
DS is close to unity). Pale blue shading at the left in (a–c) indicates the presence of a solid inner core where the core temperature, (34), falls below the melting
curve, (37) (a), (b) or where falling iron snow does not remelt and accumulates at the center (c). Whereas the melting curve is steeper than the core adiabatin (a) and (b), such that the core crystallizes from the bottom-up, it is shallower in (c) and (d) such that crystallization occurs in the “snow zone” (pale green
shading) above the intersection of the melting curve and the core adiabat. In (d), the melting curve for the falling iron snow (dotted line) is very close to thebulk core melting curve such that the core adiabat is able to cross both melting curves: the intersection with the bulk core melting curve defines the base of thesnow zone (pale green shading) and the intersection with the snow remelting curve defines the radius at which the snow remelts, and below which snow-driven
convection may take place (pale purple shading). The snow remelting curve was also calculated for (c) but is not visible because it is far above the temperature
range shown.

 R2  1 R
 r2 
Tc
 0 c exp   2  dM  r 
exp  c2 
(36)


 D 
Tcmb
N 
 DN  Mcwhich we integrate numerically to obtain the relationship between Tc and Tcmb, which is time invariant as
we neglect the small change in the core density profile that occurs during inner core growth. For the models
Tdefined by Tables 1 and 2 and illustrated in Figure 2, c is ∼1.12 for Earth and ∼1.03 for Mars.
Tcmb
HEMINGWAY AND DRISCOLL

12 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

2.3.2. Melting Curve and Core Crystallization
Parts of the core will crystallize if part of the core temperature profile falls below the melting curve (Figure 2). The melting curve may be expressed in a number of ways (e.g., Labrosse, 2003; Rivoldini et al., 2011;
Williams & Nimmo, 2004). For the sake of computational efficiency, we prefer to use a form that, whenequated with 34, yields a closed form analytical expression for the radius at which the two curves intersect.
In particular, we will express the melting curve as r2 

Tmelt  r  Tmelt, cen exp   2 
(37)
 D 
m 


where Tmelt, cen is the melting temperature at the center of the planet and where Dm is a length scale that, like
DN, depends on the planet's density structure and the material properties of its core. For Earth, we choose
Tmelt, cen = 5,700 K and Dm = 5,000 km (Table 5) in order to obtain good agreement between the melting
curve given by (37) and those based on laboratory experiments and ab initio calculations (e.g., Anzellini
et al., 2013; Boehler et al., 1990; Fei et al., 2000; Vočadlo, 2015) as shown in Figure 3, where we have plotted
Tmelt against the pressure profile obtained by integration of (27).
For Mars, we assume an Fe-FeS core for which the melting curve depends sensitively on the amount ofsulfur in the core. We compute the Martian Tmelt, cen and Dm according to the polynomials
Tmelt, cen t0  t1 S
(38)

and
(39)
Dm d0  d1S  d2 S2  d3 S3where χS is the sulfur mass fraction of the core and where the constants (Table 3) were obtained by fitting, in a least squares sense, to the results of several laboratory experiments (Boehler et al., 1990; Feiet al., 2000; Morard et al., 2008; Stewart et al., 2007) as summarized by Rivoldini et al. (2011). In particular, we fit to the melting curves described by equation 2 by Rivoldini et al. (2011), noting the following
corrections: in their Table 3, the constants a1 and a2 are reversed and the constant c1 should be positive;in their Table 4, the reference temperature for ≥21 GPa should be 1,355 K not 1,255 K. The resulting
Martian core melting curves are shown in Figure 3 for several different assumed sulfur mass fractions.
The melting curves become increasingly flat with increasing sulfur content. If the sulfur content is sufficiently high (approaching the eutectic from the pure iron side), the melting curve may be shallower (i.e.,
increasing more gradually with depth) than the core temperature profile, given by (31), whose steepness
is largely controlled by the core material's thermal expansivity via (32). In such a case, the core will not
solidify from the bottom-up, but will instead solidify first near the CMB in what may be called the ironsnow regime (Breuer et al., 2015; Davies & Pommier, 2018; Gilfoy & Li, 2020; Hauck et al., 2006; Stewart
et al., 2007; Williams, 2009) (Figure 4a). Like Rivoldini et al. (2011), our current model cannot accommodate sulfur content exceeding the eutectic.
The simplified melting curves given by Stewart et al. (2007, their Figure 2) would imply that, for subeutectic
compositions, melting temperature always decreases with depth in the Martian core, even when sulfur content is as low as 10.6 wt%, implying core crystallization occurs only as iron snow. The Rivoldini et al. (2011)
phase curve model, however, produces curves in which the melting temperature generally increases withdepth, except for sulfur content above roughly 13 wt%. We adopt the Rivoldini et al. (2011) model because it
incorporates the Stewart et al. (2007) experimental results along with those of several other studies (Boehler
et al., 1990; Fei et al., 2000; Morard et al., 2008), providing a more complete picture of the behavior acrossthe range of pressures relevant for the Martian core. Davies and Pommier (2018) favor the iron snow regime
because they adapted the Williams and Nimmo (2004) melting curve to fit the Stewart et al. (2007) curves(see liquidus curve parameters in Table 1 by Davies and Pommier (2018)). When the original parameters ofWilliams and Nimmo (2004) are retained (their Table 1), the melting curve closely resembles ours. Further
T T0 1  0.02 P ,
more, the Davies and Pommier (2018) core temperature profile (which they describe aswith pressure in units of GPa) implies a ∼22% drop in temperature across the Martian core, whereas our
HEMINGWAY AND DRISCOLL

13 of 37

Journal of Geophysical Research: Planets
Table 3
Constants Required for Equations 38 and 39 for Martian Core Melting
Curve
Constant (units)Value

t0 (K)

2,547

t1 (K)−7,689d0 (m)

5.243 × 106

d1 (m)
d2 (m)
d3 (m)

10.1029/2020JE006663

adiabatic core temperature profiles, given by (31), are significantly flatter, with temperature drops closer to ∼7%, due in part to differences inassumptions about thermal expansivity. Both of these effects mean that,while the iron snow regime may occur in our model, it is not guaranteed, and there remains the possibility of core crystallization from thebottom-up, depending on the relative slopes of the core adiabat and themelting curve.

Whenever the temperature profile is shallower than the melting curve(i.e., when the temperature curve length scale is greater than the melting−2.183 × 108
curve length scale: DN > Dm), crystallization will occur first at the center
2.725 × 109(zero radius). Conversely, whenever DN < Dm, the core temperature profile is steeper than the melting curve and crystallization will occur first at
the CMB. The transition between bottom-up inner core growth and theiron snow regime is thus determined mainly by the thermal expansivity and sulfur content of the core sincethese determine DN and Dm, respectively (Figure 4b). The length scale DN is equally sensitive to the core's
specific heat capacity, but we focus on thermal expansivity in Figure 4b because it is not as well determined
and because a wide range of values have been adopted in previous studies. As we discussed in Section 2.2,
we have simplified our model by neglecting the effect of sulfur content on the radial density structure—in
principle, increased sulfur content in the core should alter the density structure in a way that leads to slightly larger core thermal expansivity, but the effect is small compared to uncertainty in the reference thermal
expansivity.
1.689 × 107

We obtain Rint, the intersection of the adiabatic core temperature profile and the melting curve, by equating(34) and (37) and solving for r to get





1/ 2

2 
 ln T
melt, cen / Tcmb   Rc / DN  (40)Rint  Dm2  DN2

Differentiating with respect to time yields
 Tcmb
1 
1R int  
 2

(41)

2

2 Rint  Dm  DN  Tcmb

In the case of bottom-up core crystallization, these define the radius and growth rate of the solid inner core,
Ric and R ic. These expressions are similar to but differ slightly from equations 29 and 30 by Driscoll andBercovici (2014) because we have expressed the melting curve using our (37) instead of their equation 28.
Though we do not use an equivalent equation, we note also that their equation 32 contains a typo: the first
Rc on the right-hand side should be Ric. In the iron snow regime, the intersection of the core adiabat and themelting curve instead defines the base of the snow zone, above which iron crystals are forming. Note, however, that the signs are different: when the melting curve is steeper than the temperature profile, DN > Dm,
the term in parentheses in (41) is positive, meaning that the inner core growth rate is positive if the core iscooling; when the melting curve is shallower than the temperature profile, DN < Dm, the term in parentheses is negative, meaning that the radius of the base of the snow zone decreases as the core cools (i.e., thesnow zone grows downward from the top of the core).
As we discussed in Section 2.1, we also require the radius at which the snow crystals remelt, Rremelt, which
is once again obtained by equating (34) and (37) but with the melting curve for the crystals corresponding
to their reduced sulfur content, given by (15). For very small partitioning coefficients (DS ≪ 1), the snow isalmost pure iron, with a correspondingly higher melting temperature, and will therefore not remelt unlessthe adiabat is very steep. If the crystals do not remelt, a solid inner core grows as the falling solids accumulate at the center of the core. We assume the solids accumulate without porosity and therefore that an innercore built in this way follows the density profile obtained as described in Section 2.2, just as in the case of
crystallization from the bottom-up.
HEMINGWAY AND DRISCOLL

14 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

2.3.3. Mantle Temperature Structure
We assume the mantle consists of conductive thermal boundary layers at the top and bottom, separated bya convecting region in which the temperature profile is described by gTmantle  r  Tm exp  m m  Rm  r  
(42)
 cm

Rc  R), αm is2
the midmantle thermal expansivity, cm is the mantle's specific heat capacity, and gm is the gravity at Rm. We
assume linear temperature profiles across the thermal boundary layers with the temperature drops being

where Tm is a representative midmantle temperature (i.e., the temperature at radius Rm 

ΔT
Tcmb  TLM
(43)
LM

and
ΔT
TUM  Tsurf
(44)
UM

where Tsurf is a prescribed fixed value for each planet and where TLM and TUM are computed by extrapolating, via (42), the mantle temperature profile to the base and top of the convecting part of the mantle, respectively. Because the thermal boundary layer thicknesses are not yet known at this point in the calculation, we
initially assume their thicknesses to be zero (this results in two small discontinuities in the mantle temperature profile) and then we iterate to obtain a self-consistent solution, as described at the end of this section.
The boundary layer thicknesses δLM and δUM are given by (e.g., Driscoll & Bercovici, 2014)
1/3

 Ra   
(45) LM   crit LM LM   m gcmbΔTLM 

and
1/3

 Ra  (46) UM   crit UM UM gΔT m surf UM 

where Racrit = 660 is a commonly assumed critical Rayleigh number (e.g., Driscoll & Bercovici, 2014). Thek, where thethermal diffusivities, κLM and κUM are computed from the thermal conductivities as   cmpressure-dependent and therefore depth-dependent thermal conductivities are given by (Driscoll & Olson, 2011; Van den Berg et al., 2002)K 

k k0  1  P  r  0 
(47)
K 0 
where the reference conductivity is k0 = 3.3 W/m/K, P  r  is the radially dependent pressure (Section 2.2), K0
is the zero-pressure bulk modulus of the mantle, and K′0 is its derivative with respect to pressure (Table 2).
We compute the temperature dependent viscosity for the thermal boundary layer at the top of the mantle asA  11  UM  ref exp  v 


(48)
 Rg  TUM Tref  where νref is a prescribed reference viscosity corresponding to a reference temperature Tref, Av is the activation energy, and Rg is the ideal gas constant. As a means of accounting for the increase in viscosity withpressure, we then compute the lower mantle viscosity as

 LM  f UM
(49)HEMINGWAY AND DRISCOLL

15 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

where f is a free parameter of order 1–10.
Finally, having determined the boundary layer thicknesses, we update our estimates of TLM and TUM byextrapolating, via (42), the mantle temperature profile to the base and top of the convecting part of the
mantle, respectively, and then we recalculate the temperature drops and boundary layer thicknesses via(43)–(49). This process can be carried out iteratively until the two discontinuities in the temperature profile are eliminated. In practice, two iterations are sufficient to make the discontinuities negligibly small(∼2 K).

2.4. Energy Budget
2.4.1. Core Energy Budget
The temperature evolution of the core is controlled by the balance of internal radiogenic heating, heat flowfrom the core to the mantle (i.e., across the CMB), and the latent heat and gravitational potential energy
associated with core crystallization (or melting) (Figure 5). Specifically (e.g., Nimmo, 2015),
Qsec, c Qrad, c  Qcmb  QL  QG
(50)where Qsec, c is the secular warming (or cooling, if negative) of the core, and where the first term on the righthand side is the heat released by the decay of radioactive isotopes, which we compute as
Qrad, c  t  Qrad, c, init exp  t /  c 
(51)where τc is the decay timescale for the relevant radioactive materials in the core, which we assume tobe primarily 40K, and Qrad, c, init is a prescribed radiogenic heat flow at time zero. Although the option to
include radiogenic heating in the core is available, this is not the focus of our study and so we generally
assume zero radiogenic heating in the Martian core. For our Earth benchmark case, we assume someradioactivity in the core, corresponding to roughly 300 ppm Potassium (Nimmo et al., 2004), as this yields
reasonable thermal evolution results—though we note there may be other ways to obtain similar results.
The second term on the right in (50) is the heat flow across the CMB (i.e., across the lower mantle thermalboundary layer) and is obtained via (60), as described in Section 2.4.2. The third and fourth terms on theright-hand side of (50) are nonzero only while parts of the core are freezing or melting. QL is the latent
heat related to the phase change
QL  M c, sol LFe
(52)where LFe is the latent heat of fusion for iron. Similarly, QG is related to the gravitational energy associatedwith the separation of the crystallizing iron from the melt containing lighter elements, and is approximated as
QG  M c, sol EG(53)Both of these terms require the rate of mass change of the solid portion of the core, M c, sol . In the bottom-up
crystallization regime (i.e., when DN > Dm), this is effectively the mass flow rate across the phase boundary(the inner core boundary) and is obtained by multiplying the rate of change of the inner core radius, R ic,
with the surface area at that radius, Ric, and the local density, ρicb, obtained from the radial density structureT
T(Section 2.2). Noting that it can be verified that c  cmb , and using (41) for R ic, the result is
Tc
Tcmb
2 Ric  icb TcM c, sol  (54)Dm2  DN2 Tc



Notice that if the core adiabat and the melting curve intersect at a steep angle (i.e., if the length scales
DN and Dm are very different), then a given increment of cooling leads to a small increase in the mass of
HEMINGWAY AND DRISCOLL

16 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

the solid portion of the core, Mc, sol. However, if DN and Dm are very similar, the core adiabat and melting
curves become nearly colinear such that even a small increment of cooling corresponds to a large increase in Mc, sol.
The secular warming/cooling term relates to the core's heat capacity and rate of temperature change as
Qsec, c  Mc ccTc
(55)allowing us to solve equation (50) for the rate of temperature change in the core, obtaining
Qrad, c  Qcmb
Tc 


2 Rint  int  LFe  EG  (56)Mc cc  Dm2  DN2 Tc







where the term in square brackets is zero when there is no net freezing or melting in the core. If the core isexperiencing net freezing or melting, then the rate of its temperature change will be reduced because a large
fraction of the energy flow goes into the phase change part of the budget (i.e., latent heat and gravitationalenergy) rather than to the warming or cooling of the core. Notably, if freezing or melting is occurring and
DN and Dm are very similar (i.e., the core adiabat and melting curves are nearly colinear), then the rate of
temperature change in the core will become extremely slow as the energy budget becomes dominated byfreezing or melting.
In the iron snow regime (DN < Dm), if the falling snow remelts, it can drive convection in the underlying fluid as described in Section 2.3.2, but the net production of solids is zero ( M c, sol  0), and therefore there is nocontribution to the energy budget. M c, sol is nonzero only during bottom-up growth of a solid inner core orwhen iron snow crystals form in the snow zone and sink to the center without remelting. In the latter case,M c, sol is a function of the snow zone's volume and the rate of crystallization per unit volume in the snow
zone, which is related to the rate at which the core temperature is decreasing. We do not, however, attempt
to model this scenario in detail because the falling (but not remelting) snow does not drive convection, andis therefore not relevant to the problem of identifying different modes of dynamo generation (see Text S2).
2.4.2. Mantle Energy Budget
The temperature evolution of the mantle is controlled by the balance of the heat flow across the CMB, the internal radiogenic heating, and the surface heat flows (Figure 5). Specifically (e.g., Driscoll & Bercovici, 2014)



 



Qsec, m  Qcmb  Qrad, m  Qsurf, cond  Qsurf, melt(57)where Qsec, m is the secular warming (or cooling, if negative) of the mantle, Qrad, m is the radiogenic heatingin the mantle, Qcmb is the heat flow from the core into the mantle (i.e., across the lower mantle thermalboundary layer), and Qsurf, cond and Qsurf, melt represent the heat lost out the top of the mantle (i.e., to the surface via conduction and erupting melt, respectively). There is no latent heat term in (57) because we assume
the mantle cooled quickly from the magma ocean phase and is already entirely solid by the beginning of
our simulations.
The radiogenic heating is given by
Qrad, m  t  Qrad, m, init exp   t /  m 
(58)where τm is an effective decay timescale for the relevant radioactive materials in the mantle (primarily 238U,
235
U, 232Th, and 40K), and Qrad, m, init is a prescribed radiogenic heat flow at time zero.
The conductive heat loss through the top of the mantle is a function of the thermal conductivity andthe temperature gradient across the boundary layer at the top of the mantle. Approximating the latter

HEMINGWAY AND DRISCOLL

17 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

as the temperature jump (ΔTUM) divided by the thickness of the thermal boundary layer (δUM), we can
write
ΔTUM
Qsurf, cond  4 R 2 kUM
(59)

 UM

where we are using the planet's full radius to approximate the boundary's surface area. Similarly, the heatflow across the CMB can be written2
(60)Qcmb  4 Rcmb
qcmb

where we have approximated the boundary's surface area using the CMB radius. The CMB heat flux, also
appearing in equation (4), is
ΔTLM
qcmb  kLM
(61)

 LM

The boundary layer thicknesses, temperature drops, and thermal conductivities in (59) and (61) are obtained as described in Section 2.3.3.
Heat is also lost from the mantle in the form of the production of melt, which is then erupted to the surface.
The rate of heat loss due to eruptions (i.e., the advective heat loss) may be taken as proportional to the conductive heat loss through the upper mantle boundary layer (i.e., the diffusive heat loss) with the constant ofproportionality being effectively a Peclet number for heat transport out the top of the mantle
Qsurf, melt
Pe 
(62)
Qsurf, cond

where a large value corresponds to very efficient transport of melt to the surface and where a tectonically/
volcanically quiescent planet would have a Peclet number approaching zero. Although one could estimate
Qsurf, melt directly by modeling the details of how much melt should be produced based on the temperature
structure near the surface and how efficient eruptions might be, given the type of tectonic activity (e.g.,
Driscoll & Bercovici, 2014), these details will vary so widely between Earth and Mars (and any other planets
we may wish to consider using our model) that instead, we will simply treat this effective Peclet number as
an additional free parameter in the model and allow for a wide range of values.
Finally, the rate of temperature change of the mantle, Tm, is given by
Qsec, m  Mm cmTm
(63)where Mm and cm are the mantle's mass and specific heat capacity, respectively. Hence,
Qcmb  Qrad, m  Qsurf, cond  Qsurf, melt
Tm 
(64)
M mcm

2.4.3. Thermal Evolution
The thermal history of the interior, including any possible core crystallization, can now be obtained bynumerically integrating (56) and (64) through time (we use the first order Euler method, with dt = 1 Myrs,
which is sufficient because of the small curvature in the temperature histories) subject to some assumed initial values for Tcmb, Tm, Qrad, c, and Qrad, m. At each time step, we explicitly solve the coupled set of equations
to obtain Tcmb, Qcmb, Ric, and R ic (or Rsnow and Rsnow) as functions of time. Finally, we use (3), (10), (12), and(16) to determine when thermally driven and/or compositionally driven dynamos may be active.

HEMINGWAY AND DRISCOLL

18 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

3. Results and Analysis
3.1. Thermal Histories
Before applying our model to Mars, we first benchmark its performance against the case of Earth, where
we have better constraints. Figure 6 shows one possible thermal history for Earth, computed as described in
Section 2.4, where we have adopted the fixed properties given in Table 4 and the initial conditions and variable properties given in Table 5. The parameters were chosen to be similar to those adopted in previous work
(e.g., Driscoll & Bercovici, 2014)—though our mantle viscosity is defined in a way that differs from theirs
and should be regarded as an effective viscosity, calibrated to yield a temperature structure and heat flowsthat are broadly consistent with modern estimates (e.g., Jaupart et al., 2015). Figure 6c illustrates that themodel delivers a present-day inner core radius (open circle in Figure 6c) consistent with the known value(Dziewonski & Anderson, 1981). Following Nimmo et al. (2004), we assume some core radioactivity, which
delays the onset of core crystallization, helping to achieve the correct present-day inner core radius. Prior toinner core nucleation at ∼3.6 Gyrs, the dynamo is driven entirely by thermal buoyancy (red dashed line in
Figure 6d). As the inner core grows, compositional buoyancy (green dash-dotted line) becomes increasingly
important. In the relatively near future (at ∼4.57 Gyrs), the outer core becomes thermally stratified. Convection may persist after this, but must be driven entirely by compositional buoyancy. Figure 8a illustrates apresent-day snapshot of the thermal history shown in Figure 6, illustrating the corresponding temperature
structure, energy budget, and buoyancy fluxes. Our model neglects heat production within the lithosphere,
but our energy budget is otherwise broadly consistent with previously estimated ranges (e.g., Figure 1 by
Lay et al., 2008; Table 12 by Jaupart et al., 2015).
Figure 7 shows one possible thermal history for Mars which follows from the fixed properties in Table 4
and the set of initial conditions and variable properties specified under “Mars 1” in Table 5. In this example
scenario, a thermal dynamo is possible for the first 500 Myrs but then the thermal buoyancy flux becomes
negative (the outer core becomes thermally stratified), causing convection to shut down. An inner corethen nucleates at around 4 Gyrs, introducing compositional buoyancy into the outer core. At this point,
the rate of temperature decrease in the core slows dramatically (Figure 7a) because the heat flow out of thecore becomes dominated by latent heat and because the similarity in the slopes of the core adiabat and themelting curve (Figure 8b) means that a given increment of cooling must be accompanied by a great deal of
crystallization (as we discussed in Section 2.4.1). Although a solid inner core is growing after ∼4 Gyrs, thesum of the thermal and compositional buoyancy fluxes remains negative until ∼4.76 Gyrs, meaning thatouter core convection and a dynamo could reactivate after that time. Figure 8b illustrates a present-day
snapshot of the thermal history shown in Figure 7, illustrating the corresponding temperature structure,energy budget, and buoyancy fluxes.
To illustrate the variety of possibilities, Figure 9 shows four additional example thermal histories (buoyancy
flux histories only; see Figures S2–S5 for full histories) for Mars, which follow from the Mars 2–5 columns
in Table 5. Compared with the Mars 1 scenario shown in Figure 7, Mars 2 represents a case with a reduced
mantle viscosity, permitting more rapid cooling at early times such that, by the time the inner core nucleates, the thermal stratification is so strong that compositional convection is not sufficient to reactivate thedynamo. Note that this more rapid cooling scenario also requires higher thermal conductivity and lower light element density contrast (i.e., larger partitioning coefficient, DS) to ensure that the early thermal
dynamo still operates for the first 500 Myrs but without compositional buoyancy being so strong that it
would drive convection even today. In the Mars 3 scenario, the lower sulfur content and therefore higher
melting temperature means that the inner core nucleates earlier, some ∼500 Myrs after the demise of theearly thermal dynamo, causing a brief period of compositionally driven convection starting about ∼1 Gyrs
after thermal convection stopped. Because of the rapid cooling and smaller light element density contrast,
compositional buoyancy can no longer overcome the thermal stratification after ∼3.7 Gyrs, and the dynamo shuts down again permanently. In the Mars 4 scenario, the higher sulfur content lowers the melting
temperature significantly but the starting temperature is also lower, allowing core crystallization to begin(in the iron snow regime) at ∼3.2 Gyrs. At this point, because most of the energy coming out of the core is
now latent heat, the temperature decline slows dramatically (kink in Figures 9c and S4a), but because thesnow does not remelt (owing to its much lower sulfur content and therefore much higher melting point),
there is no compositionally driven convection. Finally, in the Mars 5 scenario, core crystallization begins at
HEMINGWAY AND DRISCOLL

19 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Table 4
Fixed Parameters Adopted for Thermal Evolution Models
Parameter

Symbol (units)

Value

Reference

Activation energy

Av (kJ/mol)

300

Karato and Wu (1993, Table 1)

Critical Rayleigh number

Racrit (n/d)

660

Turcotte and Schubert (1982)

Core specific heat capacity

cc (J kg−1 K−1)

840

Stacey (1995) and Nimmo (2015)

Mantle specific heat capacity

cm (J kg−1 K−1)

1,260

Tsuchiya et al. (2005, Figure 4)

Latent heat of fusion for iron

LFe (kJ/kg)

750

Nimmo (2015)

Gravitational energy density

EG (kJ/kg)

300

Driscoll and Bercovici (2014)

Core radioactive material half life

thalf,c (Gyrs)

1.251

Driscoll and Bercovici (2014)

Mantle radioactive material half life

thalf,m (Gyrs)

2.940

Driscoll and Bercovici (2014)∼2.7 Gyrs and, assuming a partitioning coefficient very close to unity, the falling snow remelts and is able
to drive convection for a brief period. As the remelting radius drops quickly, the snow-driven convection
weakens until convection shuts down again permanently. Owing to the simplifications adopted in our model, the accuracy of our results begins to decrease during the later stages of inner core growth, once the innercore radius exceeds roughly half the full core radius. However, as we will discuss further in Section 4, weare more concerned with the broad differences between the five scenarios we have shown here than withthe precise input parameter values required to produce any particular results.
3.2. Parameter Space Exploration
Whereas Figures 7 and 9 illustrate just a few possible thermal histories for Mars, exploring a wide rangeof parameters will allow us to determine more generally which initial conditions and material properties
permit convection during which periods of Mars' history and future. Although the large number of parameters makes the parameter space too large to explore exhaustively, we break the problem down into parts,
exploring subsets of the parameter space systematically in a way that yields a series of useful insights. In
particular, we will first study the parameters that control the timing of the demise of the early thermal
dynamo (Section 3.2.2) and then use the approximately known thermal dynamo end time to constrain theparameter space significantly before examining the parameters affecting the possible subsequent compositional dynamo (Section 3.2.3). Throughout the parameter space exploration, the Mars 1–5 cases defined inTable 5 will serve as points of reference, wherever possible.
3.2.1. Parameter Ranges
First, we note that several of the variable parameter ranges can be bounded to some degree. Assuming the Martian core is at least partly liquid, the CMB temperature must be above the Fe-FeS eutectic, which is around 1,600 K at 20 GPa (Boehler, 1996), and so the initial CMB temperature, Tcmb, init,
must have been at least somewhat larger, perhaps 1,700 K. It is more difficult to establish an upper
limit for the initial CMB temperature, but we consider values up to 2,500 K, which is approximately
the mantle solidus at the CMB—because a molten mantle will cool extremely rapidly, we assume theMartian mantle is entirely solid by the beginning of our simulations. In general, we assume the initialmidmantle temperature was similar to or somewhat smaller than the initial CMB temperature, though
we consider a wide range of possibilities, from 1,500 to 2,300 K. As a point of reference, the Earth's
upper mantle viscosity is estimated to be roughly η ≈ 1020–1021 Pa s (or, in terms of kinematic viscosity, ν = η/ρ ≈ 1017 m2/s, with ρ ≈ 3,300 kg/m3) at around 1,500 K (Mitrovica & Forte, 2004). However,
we allow for a wide range of possible reference mantle viscosities, from 1014 to 1019 m2/s (roughly
1017–1023 Pa s), and emphasize that we are concerned more with the effect of varying mantle viscosity
than in the absolute numbers. We consider a wide range of thermal conductivities for the Martian core(see e.g., Hsieh et al., 2020; Suehiro et al., 2017), from as low as 30 W/m/K to as high as 120 W/m/K,
allowing for the extreme possibility of nearly pure iron.

HEMINGWAY AND DRISCOLL

20 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Table 5
Variable Parameters Adopted for Thermal Evolution Models Illustrated in Figures 6–9
Parameter

Symbol(units)Earth

Mars 1Mars 2Mars 3

Mars 4

Mars 5Initial CMB temperatureTcmb (K)

5,900

2,200

2,200

2,200

2,000

2,000

Initial midmantle temperature

Tm (K)

3,300

1,800

1,800

1,800

1,600

1,600

Initial mantle radiogenic heat

Qrad, m (TW)

35

5

5

5

5

5

Initial core radiogenic heat

Qrad, c (TW)

25

0

0

0

0

0

Core sulfur mass fraction

χS (wt %)

–

6

6

4.5

8

9

Light element partitioning
coefficientDS (n/d)

–

0.05

0.5

0.5

0.5

0.99

Core melting temperature at
r=0

Tmelt, cen (K)

5,700

2,086

2,086

2,201

1,932

1,855

Melting curve length scale

Dm (km)

5,000

6,059

6,059

5,809

6,592

6,981

Light element density contrast

3Δρχ (kg/m )

500

694

276

225

330

5.7

Core thermal conductivity

kc (W/m/K)

120

56

89

89

56

89

18

17

17

17

2

νref (m /s)

10

19

10

10

10

10

1016

LM/UM viscosity ratiof (n/d)

15

2

2

2

2

2

Surface heat flow Peclet number

Pe (n/d)

0.1

0.0

0.0

0.0

0.0

0.0

Surface temperature

Tsurf (K)

288

210

210

210

210

210

Present inner core radius

Ric (km)

Mantle reference viscosity at
1,500 K

Core crystallization style
Thermally driven convection
Compositionally driven
convection

1,215

414

825

1,071

614

689

Bottom-up

Bottom-up

Bottom-up

Bottom-up

Iron snow

Iron snow

Until 500 Myrs

Until 500 Myrs

∼1.60–3.75 Gyrs

Until 500 Myrs

Until 500 Myrs

Until 500 Myrs

Until ∼4.57 Gyrs

From ∼3.65 Gyrs From ∼4.76 Gyrs

Never

Never

∼2.69–
2.82 Gyrs

Note: For Mars, Tmelt, cen, Dm, and Δρχ, are not specified directly, they are computed from χS and DS using (8), (9), (38), and (39) and shown here to the nearest
K, km, or kg/m3; reference viscosities are computed as described in Section 3.2.2 in order to ensure the early thermal dynamo terminates at 500 Myrs, and areshown here to the nearest power of 10.

3.2.2. Early Dynamo Driven by Thermal Convection
Before discussing the effects of compositional buoyancy arising from a crystallizing core, we focus on an
early Martian dynamo driven by thermal convection. We focus on thermal buoyancy flux and the model
parameters that control it because we take a positive thermal buoyancy flux to be a proxy for the presence ofa dynamo driven by thermal convection. From equation (3), it is clear that a positive thermal buoyancy flux
requires that qT, given by (4), be positive—that is, the CMB heat flux must exceed the adiabatic heat flux.
The latter is given by (5) and depends on the CMB temperature and the core thermal conductivity, thermal
expansivity, and specific heat capacity. The CMB heat flux is given by (61) and depends on the lower mantle
thermal boundary layer thickness and thermal conductivity, and the temperature drop across it. In general,
the CMB heat flux is initially superadiabatic (i.e., qT is initially positive) and decreases over time as the coreand mantle temperatures gradually converge, decreasing the temperature jump at the CMB. If the core andmantle temperatures are initially similar, however, qT (and therefore the thermal buoyancy flux) initially increases before reaching a peak and then decreasing for the rest of the planet's history (Figure 10). Similarly,
although we do not show it here, the presence of radioactive elements in the core would tend to delay thepeak in qT (and thermal buoyancy flux). The timing of the peak is likewise delayed by higher mantle viscosities (Figures 10 and 11). Increasing the Peclet number means that relatively more heat escapes through the
top of the mantle, which in turn keeps CMB heat flow higher, thus decreasing the rate at which qT drops.
Hence, although we do not illustrate this effect here, if qT is initially positive, a higher Peclet number would
lead qT to remain positive for slightly longer (e.g., for the cases illustrated in Figures 7 and 9, increasing thePeclet number from 0.0 to 0.5 would delay the end of thermal convection from ∼500 to ∼750 Myrs).
HEMINGWAY AND DRISCOLL

21 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 3. Melting curves for iron-rich cores as a function of pressure. The dashed black line in (a) represents the melting curve for pure iron obtained byAnzellini et al. (2013). The solid lines represent melting curves described by (37) and corresponding approximately to the cores of Earth and Mars. The solid
blue line, which we take to be representative for the Earth's core, is obtained from (37) by assuming that Tmelt, cen = 5,700 K and Dm = 5,000 km. The solid lines
in the brown/orange/red color scheme represent melting curves for the Martian core assuming different amounts of sulfur. The thin rectangle in (a) highlights
the region relevant to the Martian core, which is shown in greater detail in (b). Additional Martian core melting curves, computed from (37), are shown as solidlines in (b) and are compared with the Rivoldini et al. (2011) model (dashed lines) and the results from pure iron experiments (Boehler et al. (1990); Fei et al.
(2000); squares) and Fe-FeS eutectic experiments (Morard et al., 2008; Stewart et al., 2007; triangles). Lines and symbols are color-coded by sulfur content(darkest red = 0 wt%S; brightest orange = 14 wt%S).

The inference from crustal magnetic anomalies that the Martian dynamo was active early and ceased∼500 Myrs after formation thus suggests a relatively large initial core-mantle temperature difference (aspreviously argued by Williams and Nimmo, 2004) and/or a relatively low mantle viscosity. On the other
hand, if the mantle viscosity is too low, the top of the core remains superadiabatic (i.e., qT remains positive)
for a much longer period, unless the core thermal conductivity is very high (Figures 10a and 12). Hence, anearly dynamo driven by thermal convection and ending ∼500 Myrs after formation is permitted only undera fairly narrow range of conditions—namely, when ΔTinit is sufficiently large (i.e., upper left parts of Fig-

Figure 4. (a) Martian core melting curves as a function of radius (solid lines, color coded by wt%S as in Figure 3) compared with adiabatic temperature profiles
across the core (dash-dotted lines), where we have set Tcen = Tmelt, cen and where αc = 2 × 10−5 K−1. Core temperature profiles are shallower than the melting
curve when sulfur content is ≲8%, but steeper otherwise. (b) Ratio between core temperature profile and melting curve gradients at r = 0 as a function of core
thermal expansivity and sulfur content. Depending on sulfur content and thermal expansivity, core crystallization occurs either from the bottom-up (orange
shaded region below the dashed black line) or as iron snow precipitating from near the CMB (blue shaded region above the dashed black line).

HEMINGWAY AND DRISCOLL

22 of 37

Journal of Geophysical Research: Planets
Core

QG
QL

10.1029/2020JE006663

MantleQsec,cQrad,c

Qcmb

Qsec,mQsurf,cond
Qsurf,melt

Qrad,mFigure 5. Schematic illustration of the global energy budget. Qsec, c and Qsec, m represent the secular warming (ifpositive) or cooling (if negative) terms for the core and mantle, respectively. The heat flows into and out of the coreand mantle are detailed in Sections 2.4.1 and 2.4.2 and include: gravitational potential energy (QG), latent heat (QL),
radiogenic heating (Qrad, c, Qrad, m), CMB heat flow (Qcmb), and surface heat flows from the mantle (Qsurf, cond, Qsurf, melt).
The crust is not included in our model.

Figure 6. Example time evolution of Earth's (a) temperature structure, (b) heat flows, (c) inner core radius, and (d) buoyancy fluxes, where CMB = coremantle boundary; cen = zero radius; LM = top of lower mantle thermal boundary layer; UM = bottom of upper mantle thermal boundary layer;
rad = radiogenic, surf = surface. The open circle in (c) represents the known present-day inner core radius. For the buoyancy fluxes in (d), the red dashed line
represents thermal buoyancy flux, the green dash-dotted line represents compositional buoyancy flux, and the solid blue line represents the sum of the two.
The red shaded region represents the thermally driven dynamo regime; purple shading represents the compositionally driven dynamo regime; yellow shading
represents the period where both thermal and compositional buoyancy contribute to driving convection and the dynamo. A thin dashed line marks the present
day at 4.54 Gyrs.

HEMINGWAY AND DRISCOLL

23 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 7. Example time evolution of the Martian (a) temperature structure, (b) heat flows, (c) inner core radius, and (d) buoyancy fluxes, adopting the Mars 1
parameters from Table 5, where CMB = core-mantle boundary; cen = zero radius; LM = top of lower mantle thermal boundary layer; UM = bottom of upper
mantle thermal boundary layer; rad = radiogenic, surf = surface. For the buoyancy fluxes in (d), the red dashed line represents thermal buoyancy flux, thegreen dash-dotted line represents compositional buoyancy flux, and the solid blue line represents the sum of the two. The red shaded region represents thethermally driven dynamo regime operating for the first 500 Myrs; purple shading represents the compositionally driven dynamo regime; the absence of shading(white region) represents the period where there is no convection and therefore no dynamo. A thin dashed line marks the present day at 4.5 Gyrs.

ure 11) and for certain combinations of reference mantle viscosity and core thermal conductivity (i.e., along
the dashed white lines in Figure 12). Figure 13 shows this early dynamo constraint in terms of the impliedreference mantle viscosity as a function of initial CMB and midmantle temperatures. This requirement that
permissible combinations of initial temperatures imply particular values of reference mantle viscosity andcore thermal conductivity, helps to reduce the range of possibilities we will need to explore below when
discussing compositionally driven convection.
3.2.3. Compositionally Driven Convection
Up to this point, we have focused on the early thermal dynamo, whose timing is independent of core composition—provided the core does not begin to crystallize before the thermal dynamo shuts down. The purpose of Section 3.2.2 was mainly to constrain the relationship between several of the model parametersin order to make the rest of the parameter space exploration more tractable once we begin to consider thevariable of core composition (i.e., sulfur content).

HEMINGWAY AND DRISCOLL

24 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 8. Temperature structure, heat flows, and buoyancy fluxes for Earth (a) and Mars (b). These are present-day snapshots of the thermal histories
illustrated in Figures 6 and 7, respectively, and are the same as the cases illustrated in Figures 2a and 2b but additionally showing heat flows (energy budget)and buoyancy fluxes computed as described in Sections 2.4 and 2.1. Pale red shading highlights the two conductive thermal boundary layers at the top andbottom of the mantle. Pale blue shading at the left indicates the presence of a solid inner core. Although a solid inner core is growing in both cases, the net
buoyancy flux in (b) is negative, meaning the compositional buoyancy flux is not sufficient to power a dynamo at the present day.

3.3. Bottom-Up Crystallization Regime
When the sulfur content and thermal expansivity are sufficiently small, the Martian core crystallizes fromthe bottom-up (Figure 4b). A solid inner core nucleates when the central temperature, Tcen, drops below
the central melting temperature, Tmelt, cen. The time at which this occurs is a function of the melting curve(which depends on the core composition—i.e., its sulfur content), the initial core temperature, and the rate
at which the core temperature declines. In the absence of radiogenic heating in the core, the latter is entirely
a function of the CMB heat flow and hence depends mainly on the mantle viscosity and the temperature
difference between the core and mantle, which in turn has a small dependence on the surface heat flow Peclet number. However, if we again require that the early thermal dynamo ends at 500 Myrs, then any given
combination of initial CMB and midmantle temperatures constrains the mantle reference viscosity (Figure 13). Hence, the only remaining relevant variable is the sulfur content, χS, which determines the melting
temperature, Tmelt, cen, which then determines when the inner core nucleates and the implied present-day
inner core radius (Figure 14). Although core thermal conductivity has no effect on the timing of inner core
nucleation, its influence on the timing of the demise of the thermal dynamo means that core thermal conductivity is relevant to the constraint on mantle viscosity (Figure 13) and is therefore indirectly relevant toinner core growth (Figure 14). Though we do not show it here, the surface heat flow Peclet number has anadditional small effect, with larger Peclet numbers likewise corresponding to more rapid cooling and earlier
inner core nucleation, and therefore larger present-day inner core radii. A present-day solid inner core ispossible under a wide range of conditions as long as the melting temperature is sufficiently high (i.e., thesulfur content is sufficiently low) and the initial CMB temperature is sufficiently low (Figure 14).
As the solid inner core is growing, however, compositional convection becomes increasingly important. As
a result, not all combinations of parameters in Figure 14 are compatible with net buoyancy flux being negative at present or with the early thermal dynamo shutting down at ∼500 Myrs. The compositional buoyancy
flux, given by (10), and thus the vigor of compositional convection, is a function of the inner core's mass
growth rate (i.e., its surface area and rate of change of its radius) and the density contrast, Δρχ, between
the liquid outer core and the light elements that are being rejected as the inner core grows—an additionalparameter that we have not yet discussed because it does not affect any of the results we have seen so far inthis section. The density contrast is, however, what accounts for the much stronger compositional buoyancy
flux in the Mars 1 case compared to the Mars 2 and 3 cases (compare Figure 7d against Figures 9a and 9b).
Because it is the mass growth rate that matters, compositional buoyancy flux is small at first due to theinner core's initially negligible surface area, even though the inner core radius may be increasing rapidly.
However, if the density contrast is large enough and the inner core begins growing early enough, the net
HEMINGWAY AND DRISCOLL

25 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 9. Example time evolution of Martian buoyancy fluxes for the cases labeled in Table 5 as Mars 2 (a), Mars 3 (b), Mars 4 (c), and Mars 5 (d). The
red dashed line represents thermal buoyancy flux, the green dash-dotted line represents buoyancy flux arising from inner core growth, the dotted gold line
represents the snow-driven buoyancy flux, and the solid blue line represents the sum of the three. The red shaded region at the left represents the early
thermally driven dynamo regime operating for the first 500 Myrs; purple shading represents the compositionally driven dynamo regime; the absence of shading(white region) represents the period where there is no convection and therefore no dynamo. A thin dashed line marks the present day at 4.5 Gyrs.

buoyancy flux can eventually become positive again (Mars 1 case, Figure 7d). Inner core nucleation time isa function of the initial core temperature relative to the melting temperature, and the rate at which the core
is cooling. The latter is indirectly a function of core thermal conductivity: higher core thermal conductivity(which implies lower mantle viscosity to satisfy the 500 Myr early thermal dynamo requirement) translates
to more rapid core cooling and earlier nucleation of the inner core. Higher surface heat flow Peclet numbers likewise drive more rapid cooling and earlier inner core nucleation. The subsequent evolution of net
buoyancy flux then depends primarily on the light element density contrast, with higher density contrasts(a consequence of a small light element partitioning coefficient, DS) yielding larger net buoyancy fluxes.
Figures 14–16 illustrate the effects of varying the core sulfur content, χS, as well as either the initial CMB
temperature (Figure 14), the core thermal conductivity, kc (Figure 15), or the light element partitioning
coefficient, DS (Figure 16). If the sulfur content is too small, the inner core would nucleate early and compositional convection would bolster early thermal convection, precluding dynamo cessation at ∼500 Myrs
(blank white regions). Similarly, even if sulfur content and/or the initial CMB temperature were sufficiently high to permit the early thermal dynamo to shut down at ∼500 Myrs, large portions of the parameter
space must be excluded because they would yield a positive net buoyancy flux today (pale white shading).
Conversely, if the sulfur content is relatively higher, depressing the melting temperature, and/or the initial
HEMINGWAY AND DRISCOLL

26 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 10. Superadiabatic heat flux, qT, given by (4), for the first 1.5 Gyrs after the formation of Mars assuming various initial CMB and midmantle
temperatures. Panel (a) shows a case of rapid cooling (with reference mantle viscosity νref = 1016 m2/s, and core thermal conductivity kc = 89 W/m/K). Panel(b) shows a case of slower cooling (with νref = 1018 m2/s and kc = 56 W/m/K). Initial CMB temperatures are represented by different colors (2000, 2200, or
2400K; see inset) and solid and dashed lines represent cases where the CMB to midmantle temperature drops are 200 K or 400 K, respectively. Higher values of
kc would shift the curves downward. The Peclet number is assumed to be zero for all curves; higher Peclet numbers would shift the curves slightly to the right.
The peak in qT occurs earlier for lower mantle viscosities, higher initial CMB temperatures, and larger core-mantle temperature differences. The black circle in(b) corresponds to the Mars 1 reference case and the red triangle in (a) corresponds to the Mars 5 reference case (Table 5); in both cases, the superadiabatic heat
flux begins at its peak and then declines, reaching zero at 500 Myrs.

temperatures are too high, the inner core would nucleate late and compositional convection would neverbe sufficient to reactivate the dynamo (dark shading). Under certain conditions, it is also possible thatcompositional buoyancy, though not sufficient to drive convection today, was sufficient in the past (paleblue shading) or will be sufficient in the future (pale orange shading). If DS is small and a solid inner coreis present today, convection may well reactivate at some point in the future (e.g., the Mars 1 case). It is also
possible that there was a limited period of compositionally driven convection in the past (e.g., Mars 3), but

Figure 11. Timing of peak thermal buoyancy flux assuming either rapid cooling (with νref = 1016 m2/s and kc = 89 W/m/K) (a) or slower cooling (withνref = 1018 m2/s and kc = 56 W/m/K) (b). Superadiabatic heat flux (Figure 10), and thus thermal buoyancy flux, peaks early (i.e., peak time is ≲100 Myrs) only
when the initial core-mantle temperature difference (ΔTinit = Tcmb, init − Tm, init) is sufficiently large. Gray dashed lines illustrate contours of constant ΔTinit. The
black circle in (b) corresponds to the Mars 1 reference case and the red triangle in (a) corresponds to the Mars 5 reference case (Table 5); in both cases, thethermal buoyancy flux begins at its peak and declines steadily thereafter (Figure 10).

HEMINGWAY AND DRISCOLL

27 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 12. Timing of the end of thermal convection (i.e., when thermal buoyancy flux becomes negative) as a function of core thermal conductivity andreference mantle viscosity, and assuming initial CMB temperatures of either 2,000 K (a) or 2,200 K (b) and ΔTinit = 400 K. The blank white regions represent
parts of the parameter space where the mantle viscosity is too high and/or the thermal conductivity is too large for the thermal buoyancy flux to ever be positive(i.e., the core is thermally stratified to begin with and never undergoes thermal convection)—the rough edges are an artifact of the finite grid resolution. The
white dashed lines represent the approximate time at which the early Martian dynamo is thought to have shut down. The numbered symbols indicate the Mars
1–5 reference cases defined in Table 5, whose time histories are illustrated in Figures 7 and 9; Mars 2 and Mars 3 occupy the same point in (b).

this may be difficult as it requires DS to be somewhat large (∼0.3 or more) in order to maintain a small lightelement density contrast, and it requires relatively rapid cooling (high thermal conductivity and low mantle
viscosity), allowing net buoyancy flux to fall over time after an initial peak driven by compositional buoyancy (Figure 9b). Once again, correcting for some of the model's simplifications (e.g., dynamic outer coresulfur concentration) could shift the positions of the contours in Figures 14–16, especially in regions where
the inner core radius is a large fraction of the total core radius. But what is important here is the generallandscape of possibilities rather than the exact positions of any of the contours.

Figure 13. Reference mantle viscosity required to yield an early thermal dynamo that ends 500 Myrs after the formation of Mars. We consider: (a) a range ofinitial CMB and midmantle temperatures and a fixed core thermal conductivity of kc = 56 W/m/K; or (b) a range of initial CMB temperatures and core thermal
conductivities and a fixed ΔTinit = 400 K, ensuring the thermal dynamo initiates early (see Figure 11). Gray dashed lines in (a) illustrate contours of constant
ΔTinit. The numbered symbols indicate the Mars 1–5 reference cases defined in Table 5 and illustrated in Figures 7 and 9; Mars 2 and Mars 3 occupy the same
point in (b).

HEMINGWAY AND DRISCOLL

28 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 14. Present-day solid inner core radius as a function of initial CMB temperature and sulfur content, assuming ΔTinit = 400 K and the reference mantle
viscosities shown in Figure 13 (ensuring the early thermal dynamo terminates at ∼500 Myrs) and either: (a) relatively slow cooling with kc = 56 W/m/K; or(b) more rapid cooling with kc = 89 W/m/K. The black circle in (a) corresponds to the Mars 1 reference case, in which the present inner core radius is 414 km;
the green square and the blue triangle in (b) correspond, respectively, to the Mars 2 and 3 reference cases, in which the present inner core radii are 825 and1,071 km. The blank white region at the lower left is excluded because the early thermal dynamo does not terminate at 500 Myrs. The pale white shading leftof the dash-dotted white line represents the part of parameter space in which the present net buoyancy flux is positive, which we take to be incompatible withthe lack of a present-day dynamo. In the dark region to the right of the black dashed line, net buoyancy flux is never positive again after the demise of the early
thermal dynamo. In the region between the white and black curves, compositional buoyancy flux is positive either in the future (pale orange shading in a) orin the past (pale blue shading in b). The blue dashed line at ∼7.7 wt%S marks the boundary between bottom-up and top-down (iron snow, gray shading) core
crystallization. The numbered symbols indicate the Mars 1–3 reference cases defined in Table 5, whose time histories are illustrated in Figures 7 and 9.

Figure 15. Present inner core radius as a function of core sulfur content and thermal conductivity, assuming Tcmb, init = 2,200 K, ΔTinit = 400 K, and partitioning
coefficients of either 0.05 (a) or 0.5 (b). The blank white region at the left is excluded because the early thermal dynamo does not terminate at 500 Myrs. The
pale white shading left of the dash-dotted white line represents the part of parameter space in which the present net buoyancy flux is positive, which we take to
be incompatible with the lack of a present-day dynamo. In the dark region to the right of the black dashed line, net buoyancy flux is never positive again after
the demise of the early thermal dynamo. In the region between the white and black curves, compositional buoyancy flux is positive either in the future (paleorange shading) or in the past (pale blue shading). The blue dashed line at ∼7.7 wt%S marks the boundary between bottom-up and top-down (iron snow, grayshading) core crystallization. The numbered symbols indicate the Mars 1–3 reference cases defined in Table 5, whose time histories are illustrated in Figures 7
and 9.

HEMINGWAY AND DRISCOLL

29 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Figure 16. Present inner core radius as a function of core sulfur content and the partitioning coefficient, assuming an initial CMB temperature of either 2,200K(a), (b) or 2,000 K (c), (d) and assuming either kc = 56 W/m/K (a), (c) or kc = 89 W/m/K (b), (d). The blank white region at the left is excluded because the early
thermal dynamo does not terminate at 500 Myrs. The pale white shading left of the dash-dotted white line represents the part of parameter space in which
the present net buoyancy flux is positive, which we take to be incompatible with the lack of a present-day dynamo. In the dark region to the right of the black
dashed line, net buoyancy flux is never positive again after the demise of the early thermal dynamo, except in the iron snow regime with remelting (extreme top
right of d). In the region between the white and black curves, compositional buoyancy flux is positive either in the future (pale orange shading) or in the past
(pale blue shading). The blue dashed line at ∼7.7 wt%S marks the boundary between bottom-up and top-down (iron snow, gray shading) core crystallization.
The numbered symbols indicate the Mars 1–5 reference cases defined in Table 5, whose time histories are illustrated in Figures 7 and 9.

3.4. Top-Down Crystallization (Iron Snow) Regime
When the sulfur content is very high, and thus the core melting temperature is very low, core crystallization
may still occur if the initial temperatures are sufficiently low. In this case, however, the iron crystals will
form at the top of the core and sink as iron snow. If the falling snow crystals warm enough to remelt, they
can drive convection in the underlying fluid. However, this requires the melting curve for the solid iron to
be very similar to that of the bulk liquid core, which is the case only if the light element partitioning coefficient, DS, is close to unity (see Mars 5 at the extreme top right of Figure 16d). Hence, although large portionsof the parameter space permit iron snow, the falling snow generally does not remelt, and therefore is notable to drive convection. Moreover, even if such a large partitioning coefficient is possible, the resulting density contrast between the iron snow and the bulk liquid core is consequently small (e.g., in the Mars 5 case,with DS = 0.99, Δρsnow ≈ 5.7 kg/m3), meaning that iron snow would not be a strong driver of compositional
buoyancy flux in any case.
HEMINGWAY AND DRISCOLL

30 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

4. Discussion
4.1. Controlling Parameters
Our analysis has built upon previous Martian thermal evolution studies by further clarifying the relationships between several of the most important unknown or poorly constrained parameters. We showed that
there are multiple ways of ensuring that a thermal dynamo shuts down ∼500 Myrs after formation, even
without a large initial core-mantle temperature difference (e.g., via very low mantle viscosity; Figure 13a),
but a large initial core-mantle temperature difference is required to ensure a strong thermal dynamo that is
active for a significant period before shutting down (Figure 11; Williams & Nimmo, 2004). We showed thatfixing the time of the end of thermal convection places a strong constraint on the trade-off between mantle
viscosity and core thermal conductivity: lower mantle viscosity corresponds to higher CMB heat flow, andtherefore stronger thermal buoyancy flux, which will prolong thermal convection unless the core thermal
conductivity is also higher, reducing the superadiabatic heat flow (Figure 12). Although the initial temperatures are not well constrained on their own, they are closely linked to mantle viscosity and core thermal
conductivity: lowering the assumed starting temperature by 200 K, for example, requires either lowering the
mantle viscosity by an order of magnitude or lowering the core thermal conductivity by roughly 30 W/m/K
in order to maintain the same thermal dynamo end time (Figure 13). Hence, any constraints that can be
independently placed on mantle viscosity and core thermal conductivity effectively constrain the initial
CMB temperature as well. For the bottom-up inner core growth regime, our results also clarified the tradeoff between initial CMB temperature and core sulfur content as it relates to the radius of the solid innercore (Figure 14). As we discussed in Sections 2.1 and 3.2.3, however, owing to our assumed fixed outer coresulfur concentration, the positioning of the contours in Figures 14–16 should be considered approximate,especially where the inner core radius exceeds about half the full core radius. Hence, these contours should
not be taken as a precise mapping from input parameters to particular outcomes, but rather as illustrativeof the range of possible distinct dynamo regimes, as we discuss below.

4.2. Possible Dynamo Histories
Our analysis has demonstrated that there is a range of possibilities for the history and future of the Martian
dynamo that are compatible with an early thermal dynamo operating until ∼500 Myr after formation andthe lack of a dynamo operating at present. Some, but not all of these scenarios, are also compatible with thepresence of a solid inner core today. Below, we summarize the possibilities:
1. No Crystallization. If the core is sufficiently hot or its sulfur content is sufficiently high (depressing themelting temperature), the core may be entirely liquid and could remain so for billions of years to come.
In this scenario, compositionally driven convection never takes place.
2. Bottom-Up Crystallization. If the core's sulfur content, thermal expansivity, and temperature are sufficiently low, the core can crystallize from the bottom-up in a fashion similar to Earth's ongoing inner core
growth. In this bottom-up inner core growth regime, ongoing crystallization and light element rejection
may or may not generate sufficient buoyancy to reactivate the dynamo after the demise of the early thermal dynamo. Possibilities include:(a) Inner Core Growth, but no Dynamo: If the core cools sufficiently rapidly at early times such thatcooling is now slow, buoyancy arising from a crystallizing solid inner core is insufficient to overcome thermal stratification, meaning that no compositionally driven dynamo ever operates.
(b) Inner Core Growth, Future Dynamo: If the light element density contrast (Δρχ) is sufficiently large,
and the inner core nucleated only recently, compositionally driven convection may be insufficient
to drive a dynamo at present, but could become sufficient to reactivate the dynamo in the future.
(c) Inner Core Growth, Past Dynamo: If the light element density contrast is lower and the core has
already cooled a great deal, a compositionally driven dynamo may have operated for a brief period
in the past before shutting down permanently at some point before the present, owing to low coreheat flow.
3. Top-Down Crystallization: If the core's sulfur content and/or thermal expansivity are sufficiently high(see Figure 4b), but the temperature is sufficiently low, the core can crystallize from the top-down in theiron snow regime. In such a top-down crystallization scenario:
HEMINGWAY AND DRISCOLL

31 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

(a) Iron Snow, Brief Dynamo: If the light element partitioning coefficient (DS) is sufficiently close tounity (i.e., the light elements partition almost as much into the solid as into the liquid), then if
crystallization begins (either in the past or the future), the falling iron snow could remelt and mix
into the layer below, driving convection from above. Such snow-driven convection would continue
until the core is too cool for the snow to remelt, at which point the falling snow would accumulate
at the center of the core, shutting down all convection.
(b) Iron Snow, but no Dynamo: If no part of the core is ever warm enough for the falling snow to remelt (i.e., if DS is small), the snow would accumulate at the center of the core and no snow-driven
dynamo would ever operate. In this scenario, iron snow may be forming even now, but will never
be able to power a dynamo.
4.3. Implications of a Solid Inner Core
In light of the above list of possibilities, confirmation of the presence or absence of a solid inner core could
have several implications. For example:• T he detection of a solid inner core today immediately constrains the relationship between the initial coretemperature and the sulfur content (with higher sulfur content requiring lower initial temperatures, as
illustrated in Figure 14)• A solid inner core today is compatible with the bottom-up core crystallization regime (if sulfur content
is low) with compositionally driven convection occurring either in the past (e.g., Mars 3; requires smallΔρχ [large DS] and rapid cooling), future (e.g., Mars 1; requires recent inner core nucleation and largerΔρχ [small DS]), or never (e.g., Mars 2; requires rapid cooling but late inner core nucleation, perhaps due
to high initial temperatures)• A solid inner core today is also compatible with a period of snow-driven convection (if sulfur content ishigh) operating in the past (e.g., Mars 5; requires DS approaching unity), but not the future• Conversely, the absence of a solid inner core today is compatible with a period of compositionally driven
convection (in either top-down or bottom-up crystallization regimes) in the future, but not the past• The absence of a solid inner core today is, of course, also compatible with the possibility that no partof the Martian core ever crystallizes (if the core is sufficiently hot or its sulfur content is sufficientlyhigh)
Of particular interest is the possibility of compositionally driven convection operating for a limited period
following the demise of the early thermal dynamo but shutting down again before the present (exemplified
by the Mars 3 case; pale blue shading in Figures 14–16). This scenario is interesting in that it could help
explain observations of magnetized units that apparently postdate the formation of the large impact basins(e.g., Mittelholz et al., 2020). On the other hand, this scenario requires a relatively large light element partitioning coefficient (DS > 0.3), which is not expected at Mars core pressures (Zhang et al., 2020). A smaller
partitioning coefficient would yield a larger light element density contrast, leading to stronger compositional buoyancy, that would remain effective even to the present day. If the partitioning coefficient is small,then the only possibilities for a compositionally driven dynamo is that it becomes active at some point in
the future (pale orange shading in Figures 14–16), or not at all. If the partitioning coefficient is small anda solid inner core is detected, then, for almost any combination of the other parameters, compositionallydriven convection will become active within the next ∼1 Gyrs. Conversely, if the Martian core is confirmed
to be entirely liquid, then it would not be clear whether or not compositionally driven convection will evertake place.

4.4. Caveats and Future WorkOur exploration of the large parameter space was made tractable in large part by adopting the assumptionthat the early thermal dynamo shut down 500 Myrs after the formation of Mars. Relaxing this assumption
would affect our subsequent results related to the timing of the possible onset and cessation of compositionally driven convection. Nevertheless, any early thermal dynamo would shut down eventually unless the
mantle viscosity were very low. And if the mantle viscosity were much lower, the resulting rapid cooling

HEMINGWAY AND DRISCOLL

32 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

would cause a solid inner core to nucleate early, leading to strong compositionally driven convection thatwould continue operating even today, unless the sulfur content were high enough for crystallization to
occur in the iron snow regime.
Another caveat is that our model involves a number of simplifications including: neglecting the changes in
sulfur concentration of the liquid part of the core over time; neglecting the radial variation in core thermal
conductivity and expansivity; restricting our analysis to the subeutectic Fe-FeS system; and defining thethreshold for convection and dynamo action based on the simple sum of the individual buoyancy fluxes in(16). The effects of relaxing these assumptions are not necessarily obvious but will be the subject of future
work. As we discussed in Sections 2.1, 3.2.3, and 4.1, although it does not affect our general conclusions,
our results could be made more precise by dynamically tracking the gradual increase in the outer core's
sulfur concentration during inner core growth. This could allow for a more precise mapping between input
parameters and particular outcomes such as the present inner core radius. The addition of a time-variable
Peclet number could also be valuable as a means of accounting for changes over time in the level of volcanic
activity. Whereas the efficiency of thermally driven convection decreases over time in any case, a decreasein volcanic activity could compound this effect. Because this affects the timing of the demise of the early
thermal dynamo, it may therefore be worth investigating how such a variable Peclet number would trade
off against other parameters such as mantle viscosity. Similarly, our simplified radial density and thermal
structure may lead us to underestimate or overestimate the conductive heat flux at the surface. This is nota concern for our current model because this quantity is not yet well known for Mars. However, improving
the fidelity of this aspect of the model may permit future heat probe measurements (Spohn et al., 2018)
to be used to help constrain mantle thermal evolution, therefore providing an indirect constraint on coreproperties. In computing the mantle viscosity with (48), we assumed a constant mantle composition overtime and did not account for the possibility of dehydration stiffening or other compositional changes which
could have affected the mantle's thermal evolution as volatiles degassed into the atmosphere (Fraeman &Korenaga, 2010; Grott et al., 2011). In principle, future observational constraints on eruption and degassing
rates could be used to further constrain the evolution of Mars' mantle viscosity. However, even then, there
is sufficient uncertainty in mantle viscosity that it may be difficult to meaningfully relate such effects to our
results.
Finally, we intend to use our model to study the Moon, Mercury, Venus, Io, Ganymede, terrestrial-like
exoplanets, and super Earths, as well as alternate internal structures for Mars (e.g., Brennan et al., 2020;
Deng & Levander, 2020), to characterize how planet scale and core size (see Boujibar et al., 2020) interact
with the various model parameters we have considered here to determine how and when thermally driven
and/or compositionally driven dynamos may operate in general. The Moon is a particularly interesting casebecause, like Mars, it has no dynamo operating today, but parts of its crust are nevertheless strongly magnetized. The combination of age dating and measurement of the magnetic properties of returned lunar samples
has permitted some reconstruction of the Moon's magnetic history (e.g., Garrick-Bethell et al., 2009; Tikoo
et al., 2017; Weiss & Tikoo, 2014). This has raised questions about how the Moon could have generated a
dynamo field strong enough to magnetize not just the collected samples, but also the even more strongly
magnetized larger scale crustal magnetic anomalies (e.g., Hood et al., 2001; Halekas et al., 2010; Hemingway & Tikoo, 2018), and how and when the dynamo shut down. Hence, we are motivated to apply ourmodel to the Moon and show how it relates to the results we have discussed here. The case of Mars has been
an appropriate starting point, however, as it may be a boundary case for which there is an especially wide
range of possibilities. For example, in larger bodies, the much higher core pressures make bottom-up core
crystallization more likely whereas the much lower pressures in smaller bodies like Ganymede and Mercury mean that top-down core crystallization is more likely. Our model for driving convection via remeltingof the falling iron snow showed that the details determining whether or not, and where, the crystals canremelt are clearly important and worthy of more detailed study. Although we showed that this mechanism
may not be efficient for Mars due to the similarity in slopes of the core adiabat and the melting curve, thesame mechanism could be more efficient at Ganymede and Mercury, where the much lower core pressurespermit much higher thermal expansivity and therefore a much steeper adiabat (see equation 6). In any case,a comprehensive exploration of the parameter space spanned by this wider range of bodies could further
illuminate the range of possibilities for terrestrial-like planet dynamo histories and the factors that control
them.
HEMINGWAY AND DRISCOLL

33 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

5. Conclusions
Our analysis has demonstrated that there is a range of possibilities for the history and future of the Martian
dynamo. Consistent with previous work, we have shown that a large initial core-mantle temperature difference can help to ensure that thermal convection is able to power a dynamo at early times, and that, given
an appropriate combination of mantle viscosity and core thermal conductivity, thermal stratification can bemade to shut down this early thermal dynamo ∼500 Myrs after the formation of Mars. We have additionallyshown that the presence of a solid inner core inside Mars' otherwise liquid core is compatible with the lackof an active dynamo at present and that its presence could imply that either: (a) a compositionally driven
dynamo will become active at some point in the next ∼1 Gyrs or, (b) if the light element partitioning islarge enough, such a dynamo may have been active for a limited period in the past. We have argued that, if
the light elements partition strongly into the liquid phase, it may be difficult to drive convection in the topdown crystallization (iron snow) regime because the falling snow crystals may not warm enough to remelt
and because, even if they do, their density may not contrast enough with that of the liquid they are mixing
into. Our analysis has clarified the effects of, and illustrated the importance of obtaining better constraints
on, many of the most important parameters, including: core composition (i.e., sulfur content), light element
partitioning, mantle viscosity, core thermal conductivity, and core thermal expansivity. Future experimental
and modeling efforts to better constrain these parameters could help to more definitively determine which
Martian dynamo scenarios are most likely. More generally, expanding our analysis to include a wider arrayof bodies, ranging from the Moon to super Earths, could help to further clarify the conditions that determine which types of dynamos may operate on which types of bodies, and when.

Data Availability Statement
All model outputs are available as supplementary data files attached to the online version of this paper andat http://doi.org/10.5281/zenodo.4459482.
Acknowledgments
Douglas J. Hemingway and Peter E.
Driscoll were funded by the Carnegie
Institution for Science in Washington,
DC. We thank Yingwei Fei, Tushar Mittal, Francis Nimmo, Bruce Buffet, Tim
Jones, Nick Knezek, Naor Movshovitz,
Gary Glatzmaier, and Sabine Stanley for
helpful discussions. We also thank Rob
Lillis and an anonymous reviewer forthoughtful and constructive reviews.

HEMINGWAY AND DRISCOLL

References
Acuña, M. H. (1999). Global distribution of crustal magnetization discovered by the Mars global surveyor MAG/ER experiment. Science,
284(5415), 790–793. https://doi.org/10.1126/science.284.5415.790
Anzellini, S., Dewaele, A., Mezouar, M., Loubeyre, P., & Morard, G. (2013). Melting of iron at Earth's inner core boundary based on fast
X-ray diffraction. Science, 340(6131), 464–466. https://doi.org/10.1126/science.1233514
Banerdt, W. B., Smrekar, S. E., Banfield, D., Giardini, D., Golombek, M., Johnson, C. L., et al. (2020). Initial results from the InSight mission
on Mars. Nature Geoscience, 13(3), 183–189. https://doi.org/10.1038/s41561-020-0544-yBoehler, R. (1996). Fe-FeS eutectic temperatures to 620 kbar. Physics of the Earth and Planetary Interiors, 96(2–3 Spec. Iss.), 181–186.
https://doi.org/10.1016/0031-9201(96)03150-0
Boehler, R., Von Bargen, N., & Chopelas, A. (1990). Melting, thermal expansion, and phase transitions of iron at high pressures. Journal ofGeophysical Research, 95(B13), 21731–21736. https://doi.org/10.1029/jb095ib13p21731
Boujibar, A., Driscoll, P., & Fei, Y. (2020). Super-earth internal structures and initial thermal states. Journal of Geophysical Research: Planets, 125, e2019JE006124. https://doi.org/10.1029/2019JE006124
Brennan, M. C., Fischer, R. A., & Irving, J. C. E. (2020). Core formation and geophysical properties of Mars. Earth and Planetary Science
Letters, 530, 115923. https://doi.org/10.1016/j.epsl.2019.115923
Breuer, D., Rueckriemen, T., & Spohn, T. (2015). Iron snow, crystal floats, and inner-core growth: Modes of core solidification and implications for dynamos in terrestrial planets and moons. Progress in Earth and Planetary Science, 2(1), 39. https://doi.org/10.1186/
s40645-015-0069-y
Chopelas, A., & Boehler, R. (1992). Thermal expansivity in the lower mantle. Geophysical Research Letters, 19(19), 1983–1986. https://doi.
org/10.1029/92GL02144
Connerney, J. (2015). Planetary magnetism. In Treatise on geophysics (pp. 195–237). Amsterdam: Elsevier. https://doi.org/10.1016/
B978-0-444-53802-4.00171-8
Connerney, J. E. P., Acuña, M. H., Ness, N. F., Kletetschka, G., Mitchell, D. L., Lin, R. P., & Reme, H. (2005). Tectonic implications of Mars
crustal magnetism. Proceedings of the National Academy of Sciences of the United Stated of America, 102(42), 14970–14975. https://doi.
org/10.1073/pnas.0507469102
Davies, C. J., & Pommier, A. (2018). Iron snow in the Martian core? Earth and Planetary Science Letters, 481, 189–200. https://doi.
org/10.1016/j.epsl.2017.10.026
Davies, C., Pozzo, M., Gubbins, D., & Alfè, D. (2015). Constraints from material properties on the dynamics and evolution of Earth's core.
Nature Geoscience, 8(9), 678–685. https://doi.org/10.1038/ngeo2492
Deng, S., & Levander, A. (2020). Autocorrelation reflectivity of mars. Geophysical Research Letters, 47, e2020GL089630. https://doi.
org/10.1029/2020GL089630
Dietrich, W., & Wicht, J. (2013). A hemispherical dynamo model: Implications for the Martian crustal magnetization. Physics of the Earthand Planetary Interiors, 217, 10–21. https://doi.org/10.1016/j.pepi.2013.01.001
Driscoll, P., & Bercovici, D. (2014). On the thermal and magnetic histories of Earth and Venus: Influences of melting, radioactivity, andconductivity. Physics of the Earth and Planetary Interiors, 236, 36–51. https://doi.org/10.1016/j.pepi.2014.08.004

34 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Driscoll, P., & Olson, P. (2011). Optimal dynamos in the cores of terrestrial exoplanets: Magnetic field generation and detectability. Icarus,
213(1), 12–23. https://doi.org/10.1016/j.icarus.2011.02.010
Du, Z., Jackson, C., Bennett, N., Driscoll, P., Deng, J., Lee, K. K. M., et al. (2017). Insufficient energy from MgO exsolution to power earlygeodynamo. Geophysical Research Letters, 44, 11376–11381. https://doi.org/10.1002/2017GL075283
Dwyer, C. A., Stevenson, D. J., & Nimmo, F. (2011). A long-lived lunar dynamo driven by continuous mechanical stirring. Nature,479(7372), 212–214. https://doi.org/10.1038/nature10564
Dziewonski, A. M., & Anderson, D. L. (1981). Preliminary reference Earth model. Physics of the Earth and Planetary Interiors, 25(4),297–356. https://doi.org/10.1016/0031-9201(81)90046-7
Fei, Y., Li, J., Bertka, C. M., & Prewitt, C. T. (2000). Structure type and bulk modulus of Fe3S, a new iron-sulfur compound. American
Mineralogist, 85(11–12), 1830–1833. https://doi.org/10.2138/am-2000-11-1229
Folkner, W. M., Dehant, V., Le Maistre, S., Yseboodt, M., Rivoldini, A., Van Hoolst, T., et al. (2018). The rotation and interior structure
experiment on the InSight mission to mars. Space Science Reviews, 214(5), 100. https://doi.org/10.1007/s11214-018-0530-5
Fraeman, A. A., & Korenaga, J. (2010). The influence of mantle melting on the evolution of Mars. Icarus, 210(1), 43–57. https://doi.
org/10.1016/j.icarus.2010.06.030
Franz, H. B., King, P. L., & Gaillard, F. (2019). Sulfur on Mars from the atmosphere to the core. Volatiles in the Martian Crust, 119–183.
https://doi.org/10.1016/B978-0-12-804191-8.00006-4
Gaillard, F., & Scaillet, B. (2009). The sulfur content of volcanic gases on Mars. Earth and Planetary Science Letters, 279(1–2), 34–43. https://doi.org/10.1016/j.epsl.2008.12.028
Garrick-Bethell, I., Weiss, B. P., Shuster, D. L., & Buz, J. (2009). Early lunar magnetism. Science (New York, N.Y.), 323(5912), 356–359.
https://doi.org/10.1126/science.1166804
Gilfoy, F., & Li, J. (2020). Thermal state and solidification regime of the Martian core: Insights from the melting behavior of FeNi-S at 20
GPa. Earth and Planetary Science Letters, 541, 116285. https://doi.org/10.1016/j.epsl.2020.116285
Gomi, H., Ohta, K., Hirose, K., Labrosse, S., Caracas, R., Verstraete, M. J., & Hernlund, J. W. (2013). The high conductivity of iron and thermal evolution of the Earth's core. Physics of the Earth and Planetary Interiors, 224, 88–103. https://doi.org/10.1016/j.pepi.2013.07.010
Grott, M., Morschhauser, A., Breuer, D., & Hauber, E. (2011). Volcanic outgassing of CO2 and H2O on Mars. Earth and Planetary Science
Letters, 308(3–4), 391–400. https://doi.org/10.1016/j.epsl.2011.06.014
Halekas, J. S., Lillis, R. J., Lin, R. P., Manga, M., Purucker, M. E., & Carley, R. A. (2010). How strong are lunar crustal magnetic fields at thesurface? Considerations from a reexamination of the electron reflectometry technique. Journal of Geophysical Research, 115, E07007.
https://doi.org/10.1029/2009je003516
Hauck, S. A., Aurnou, J. M., & Dombard, A. J. (2006). Sulfur's impact on core evolution and magnetic field generation on Ganymede.
Journal of Geophysical Research, 111, E09008. https://doi.org/10.1029/2005JE002557
Hauck, S. A., & Phillips, R. J. (2002). Thermal and crustal evolution of Mars. Journal of Geophysical Research, 107(E7). 5052. https://doi.
org/10.1029/2001JE001801
Hemingway, D. J., & Tikoo, S. M. (2018). Lunar swirl morphology constrains the geometry, magnetization, and origins of lunar magnetic
anomalies. Journal of Geophysical Research: Planets, 123, 2223–2241. https://doi.org/10.1029/2018je005604
Hood, L. L., Zakharian, A., Halekas, J., Mitchell, D. L., Lin, R. P., Acuña, M. H., & Binder, A. B. (2001). Initial mapping and interpretation oflunar crustal magnetic anomalies using Lunar Prospector magnetometer data. Journal of Geophysical Research, 106(E11), 27825–27839.
https://doi.org/10.1029/2000JE001366
Hori, K., & Wicht, J. (2013). Subcritical dynamos in the early Mars' core: Implications for cessation of the past Martian dynamo. Physics ofthe Earth and Planetary Interiors, 219, 21–33. https://doi.org/10.1016/j.pepi.2013.03.005
Hsieh, W. P., Goncharov, A. F., Labrosse, S., Holtgrewe, N., Lobanov, S. S., Chuvashova, I., et al. (2020). Low thermal conductivity ofiron-silicon alloys at Earth's core conditions with implications for the geodynamo. Nature Communications, 11(1), 1–7. https://doi.
org/10.1038/s41467-020-17106-7
Jaupart, C., Labrosse, S., Lucazeau, F., & Mareschal, J. C. (2015). Temperatures, heat, and energy in the mantle of the Earth (7). Amsterdam: Elsevier B.V. https://doi.org/10.1016/B978-0-444-53802-4.00126-3
Karato, S.-i., & Wu, P. (1993). Rheology of the upper mantle: A synthesis. Science, 260(5109), 771–778. https://doi.org/10.1126/
science.260.5109.771
Ke, Y., & Solomatov, V. S. (2009). Coupled core-mantle thermal evolution of early Mars. Journal of Geophysical Research, 114, E07004.
https://doi.org/10.1029/2008JE003291
Labrosse, S. (2003). Thermal and magnetic evolution of the Earth's core. Physics of the Earth and Planetary Interiors, 140(1–3), 127–143.
https://doi.org/10.1016/j.pepi.2003.07.006
Labrosse, S., Poirier, J.-P., & Le Mouël, J.-L. (2001). The age of the inner core. Earth and Planetary Science Letters, 190, 111–123. https://doi.
org/10.1016/s0012-821x(01)00387-9
Landeau, M., & Aubert, J. (2011). Equatorially asymmetric convection inducing a hemispherical magnetic field in rotating spheres andimplications for the past martian dynamo. Physics of the Earth and Planetary Interiors, 185(3–4), 61–73. https://doi.org/10.1016/j.
pepi.2011.01.004
Langlais, B., Thébault, E., Houliez, A., Purucker, M. E., & Lillis, R. J. (2019). A new model of the crustal magnetic field of mars using MGS
and MAVEN. Journal of Geophysical Research: Planets, 124, 1542–1569. https://doi.org/10.1029/2018JE005854
Lay, T., Hernlund, J., & Buffett, B. A. (2008). Core-mantle boundary heat flow. Nature Geoscience, 1(1), 25–32. https://doi.org/10.1038/
ngeo.2007.44
Le Bars, M., Wieczorek, M. A., Karatekin, Ö., Cébron, D., & Laneuville, M. (2011). An impact-driven dynamo for the early Moon. Nature,479(7372), 215–218. https://doi.org/10.1038/nature10565
Lillis, R. J., Frey, H. V., & Manga, M. (2008). Rapid decrease in Martian crustal magnetization in the Noachian era: Implications for the
dynamo and climate of early Mars. Geophysical Research Letters, 35, L14203. https://doi.org/10.1029/2008GL034338
Lillis, R. J., Manga, M., Mitchell, D. L., Lin, R. P., & Acuna, M. H. (2006). Unusual magnetic signature of the Hadriaca Patera Volcano:Implications for early Mars. Geophysical Research Letters, 33, L03202. https://doi.org/10.1029/2005GL024905
Lillis, R. J., Robbins, S., Manga, M., Halekas, J. S., & Frey, H. V. (2013). Time history of the Martian dynamo from crater magnetic field
analysis. Journal of Geophysical Research: Planets, 118, 1488–1511. https://doi.org/10.1002/jgre.20105
Lister, J. R. (2003). Expressions for the dissipation driven by convection in the Earth's core. Physics of the Earth and Planetary Interiors,
140(1–3), 145–158. https://doi.org/10.1016/j.pepi.2003.07.007
Merkel, S., Goncharov, A. F., Mao, H. K., Gillet, P., & Hemley, R. J. (2000). Raman Spectroscopy of Iron to 152 Gigapascals: Implications for
Earth's Inner Core. Science, 288(5471), 1626–1629. https://doi.org/10.1126/science.288.5471.1626

HEMINGWAY AND DRISCOLL

35 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Mitrovica, J. X., & Forte, A. M. (2004). A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic
adjustment data. Earth and Planetary Science Letters, 225(1–2), 177–189. https://doi.org/10.1016/j.epsl.2004.06.005
Mittal, T., Knezek, N., Arveson, S. M., McGuire, C. P., Williams, C. D., Jones, T. D., & Li, J. (2020). Precipitation of multiple light elements
to power Earth's early dynamo. Earth and Planetary Science Letters, 532, 116030. https://doi.org/10.1016/j.epsl.2019.116030
Mittelholz, A., Johnson, C. L., Feinberg, J. M., Langlais, B., & Phillips, R. J. (2020). Timing of the Martian dynamo: New constraints for acore field 4.5 and 3.7 Ga ago. Science Advances, 6(18), eaba0513. https://doi.org/10.1126/sciadv.aba0513
Morard, G., Andrault, D., Guignot, N., Sanloup, C., Mezouar, M., Petitgirard, S., & Fiquet, G. (2008). In situ determination of Fe-Fe3S phase
diagram and liquid structural properties up to 65 GPa. Earth and Planetary Science Letters, 272(3–4), 620–626. https://doi.org/10.1016/j.
epsl.2008.05.028
Nimmo, F. (2015). Energetics of the core (8). Amsterdam: Elsevier B.V. https://doi.org/10.1016/B978-0-444-53802-4.00139-1
Nimmo, F., Price, G. D., Brodholt, J., & Gubbins, D. (2004). The influence of potassium on core and geodynamo evolution. GeophysicalJournal International, 156(2), 363–376. https://doi.org/10.1111/j.1365-246X.2003.02157.x
Nimmo, F., & Stevenson, D. J. (2000). Influence of early plate tectonics on the thermal evolution and magnetic field of Mars. Journal ofGeophysical Research, 105(E5), 11969–11979. https://doi.org/10.1029/1999JE001216
Oganov, A. R. (2007). Theory and practice—Thermodynamics, equations of state, elasticity, and phase transitions of minerals at high pressures and temperatures. Treatise on Geophysics, 2, 121–152. https://doi.org/10.1016/B978-044452748-6.00033-X
Olson, P., & Christensen, U. R. (2006). Dipole moment scaling for convection-driven planetary dynamos. Earth and Planetary Science
Letters, 250(3–4), 561–571. https://doi.org/10.1016/j.epsl.2006.08.008
Olson, P., Landeau, M., & Hirsh, B. H. (2017). Laboratory experiments on rain-driven convection: Implications for planetary dynamos.
Earth and Planetary Science Letters, 457, 403–411. https://doi.org/10.1016/j.epsl.2016.10.015
O'Rourke, J. G., & Shim, S. H. (2019). Hydrogenation of the Martian Core by Hydrated Mantle Minerals With Implications for the Early
Dynamo. Journal of Geophysical Research: Planets, 124, 3422–3441. https://doi.org/10.1029/2019JE005950
Purucker, M. E., & Whaler, K. A. (2015). Crustal magnetism (5). Amsterdam: Elsevier B.V. https://doi.org/10.1016/B978-0-444-53802-4.00100-7
Rivoldini, A., Van Hoolst, T., Verhoeven, O., Mocquet, A., & Dehant, V. (2011). Geodesy constraints on the interior structure and composition of Mars. Icarus, 213(2), 451–472. https://doi.org/10.1016/j.icarus.2011.03.024
Roberts, J. H., Lillis, R. J., & Manga, M. (2009). Giant impacts on early Mars and the cessation of the Martian dynamo. Journal of Geophysical Research, 114, E04009. https://doi.org/10.1029/2008JE003287
Rückriemen, T., Breuer, D., & Spohn, T. (2015). The Fe snow regime in Ganymede's core: A deep-seated dynamo below a stable snow zone.
Journal of Geophysical Research: Planets, 120, 1095–1118. https://doi.org/10.1002/2014JE004781
Schubert, G., Russell, C. T., & Moore, W. B. (2000). Timing of the Martian dynamo. Nature, 408(6813), 666–667. https://doi.
org/10.1038/35047163
Schubert, G., Solomon, S. C., Turcotte, D. L., Drake, M. J., & Sleep, N. H. (1992). Origin and thermal evolution of Mars. In H. H. Kieffer, B. M. Jakosky, C. W. Snyder, & M. S. Matthews (Eds.), Mars (pp. 147–183). Tucson, AZ: University of Arizona Press. https://doi.
org/10.1017/CBO9781107415324.004
Spohn, T., Grott, M., Smrekar, S. E., Knollenberg, J., Hudson, T. L., Krause, C., et al. (2018). The heat flow and physical properties package(HP3) for the insight mission. Space Science Reviews, 214, 96. https://doi.org/10.1007/s11214-018-0531-4
Stacey, F. D. (1992). Physics of the Earth. Brisbane: Brookfield Press.
Stacey, F. D. (1995). Theory of thermal and elastic properties of the lower mantle and core. Physics of the Earth and Planetary Interiors,89(3–4), 219–245. https://doi.org/10.1016/0031-9201(94)03005-4
Stacey, F. D., & Davis, P. M. (2004). High pressure equations of state with applications to the lower mantle and core. Physics of the Earthand Planetary Interiors, 142(3–4), 137–184. https://doi.org/10.1016/j.pepi.2004.02.003
Stevenson, D. J. (2001). Mars' core and magnetism, Nature. Nature, 412(6843), 214–219. https://doi.org/10.1038/35084155
Stevenson, D. J. (2003). Planetary magnetic fields. Earth and Planetary Science Letters, 208, 1–11. https://doi.org/10.1016/
S0012-821X(02)01126-3
Stevenson, D. J., Spohn, T., & Schubert, G. (1983). Magnetism and thermal evolution of the terrestrial planets. Icarus, 54(3), 466–489.
https://doi.org/10.1016/0019-1035(83)90241-5
Stewart, A. J., Schmidt, M. W., van Westrenen, W., & Liebske, C. (2007). Mars: A new core-crystallization regime. Science, 316, 1323–1325.
https://doi.org/10.1126/science.1140549
Suehiro, S., Ohta, K., Hirose, K., Morard, G., & Ohishi, Y. (2017). The influence of sulfur on the electrical resistivity of hcp iron: Implications for the core conductivity of Mars and Earth. Geophysical Research Letters, 44, 8254–8259. https://doi.org/10.1002/2017GL074021
Tarduno, J. A., Cottrell, R. D., Davis, W. J., Nimmo, F., & Bono, R. K. (2015). A Hadean to Paleoarchean geodynamo recorded by single
zircon crystals. Science, 349(6247), 521–524. https://doi.org/10.1126/science.aaa9114
Tarduno, J. A., Cottrell, R. D., Watkeys, M. K., Hofmann, A., Doubrovine, P. V., Mamajek, E. E., et al. (2010). Geodynamo, solar wind, andmagnetopause 3.4 to 3.45 billion years ago. Science, 327(5970), 1238–1240. https://doi.org/10.1126/science.1183445
Tikoo, S. M., Weiss, B. P., Shuster, D. L., Suavet, C., Wang, H., & Grove, T. L. (2017). A two-billion-year history for the lunar dynamo. Science
Advances, 3(8), e1700207. https://doi.org/10.1126/sciadv.1700207
Tsuchiya, J., Tsuchiya, T., & Wentzcovitch, R. M. (2005). Vibrational and thermodynamic properties of MgSiO3 postperovskite. Journal ofGeophysical Research, 110, B02204. https://doi.org/10.1029/2004JB003409
Turcotte, D. L., & Schubert, G. (1982). Geodynamics. Cambridge, UK: Cambridge University Press.
Van den Berg, A. P., Yuen, D. A., & Allwardt, J. R. (2002). Non-linear effects from variable thermal conductivity and mantle internal heating: Implications for massive melting and secular cooling of the mantle. Physics of the Earth and Planetary Interiors, 129(3–4), 359–375.
https://doi.org/10.1016/S0031-9201(01)00304-1
Van Hoolst, T., Dehant, V., Roosbeek, F., & Lognonné, P. (2003). Tidally induced surface displacements, external potential variations, andgravity variations on Mars. Icarus, 161(2), 281–296. https://doi.org/10.1016/S0019-1035(02)00045-3
Vervelidou, F., Lesur, V., Grott, M., Morschhauser, A., & Lillis, R. J. (2017). Constraining the date of the Martian dynamo shutdown by meansof crater magnetization signatures. Journal of Geophysical Research: Planets, 122, 2294–2311. https://doi.org/10.1002/2017JE005410
Vinet, P., Rose, J. H., Ferrante, J., & Smith, J. R. (1989). Universal features of the equation of state of solids. Journal of Physics: Condensed
Matter, 1(11), 1941–1963. https://doi.org/10.1088/0953-8984/1/11/002
Vočadlo, L. (2015). Earth's core: Iron and iron alloys. In Treatise on Geophysics (2, pp. 117–147). Oxford: Elsevier Ltd. https://doi.
org/10.1016/B978-0-444-53802-4.00032-4
Weiss, B. P., & Tikoo, S. M. (2014). The lunar dynamo. Science (New York, N.Y.), 346(6214), 1246753. https://doi.org/10.1126/science.1246753

HEMINGWAY AND DRISCOLL

36 of 37

Journal of Geophysical Research: Planets

10.1029/2020JE006663

Williams, J.-P., & Nimmo, F. (2004). Thermal evolution of the Martian core: Implications for an early dynamo. Geology, 32(2), 97–100.
https://doi.org/10.1130/G19975.1
Williams, Q. (2009). Bottom-up versus top-down solidification of the cores of small solar system bodies: Constraints on paradoxical cores.
Earth and Planetary Science Letters, 284(3–4), 564–569. https://doi.org/10.1016/j.epsl.2009.05.019
Yoder, C. F., Konopliv, a. S., Yuan, D. N., Standish, E. M., & Folkner, W. M. (2003). Fluid core size of Mars from detection of the solar tide.
Science (New York, N.Y.), 300(5617), 299–303. https://doi.org/10.1126/science.1079645
Zhang, Z., Csányi, G., & Alfè, D. (2020). Partitioning of sulfur between solid and liquid iron under Earth's core conditions: Constraints
from atomistic simulations with machine learning potentials. Geochimica et Cosmochimica Acta, 291, 5. https://doi.org/10.1016/j.
gca.2020.03.028
Zuber, M. T., Smith, D., Cheng, A., Garvin, J., Aharonson, O., Cole, T., et al. (2000). The shape of 433 eros from the NEAR-shoemaker laser
rangefinder. Science (New York, N.Y.), 289(5487), 2097–2101. https://doi.org/10.1126/science.289.5487.2097

HEMINGWAY AND DRISCOLL

37 of 37