Thermal evolution of Mars and Mercury (2-D)

Mars and Mercury are one-plate planets: a thick, cold lid sits over a convecting mantle and the planet loses its heat by conduction through it. The two differ mainly in their mantles. Mars has a mantle 1550 km deep over a core about half the planet's radius; Mercury's mantle is only 420 km deep over a core that fills most of the planet. This example follows both over 4.5 Gyr in the spherical annulus, from a hot start, with decaying radiogenic heating and a core that cools or warms through the heat crossing the core-mantle boundary. Each planet is run with and without a 50 km crust that holds twice the bulk heat production.

Runs in progress

Mercury without crust has reached 4.5 Gyr; the other three runs are still computing. Their figures and numbers are from the state reached so far and are updated as the runs advance.

run	reached	status
Mars, no crust	1460 Myr	running
Mars, with crust	1090 Myr	running
Mercury, no crust	4500 Myr	finished
Mercury, with crust	120 Myr	running

Mars without crust at 60, 300, 690 and 1290 Myr, temperature between 1700 and 2300 K (the lid, colder than 1700 K, is the dark rim). Nr = 124 (124 × 1650 cells, 12.5 km radially). The mantle warms under its lid: the radiogenic heat is produced faster than the lid can conduct it away.

Mercury without crust at 30, 990, 2490 and 4500 Myr, temperature between 1000 and 2000 K. Nr = 84 (84 × 1533 cells, 5 km radially). The mantle is so thin that it holds a single row of small cells; they merge into fewer, wider ones as the lid thickens, and by about 3 Gyr the lid has taken the whole mantle and heat leaves by conduction alone.

Set-up

The physical set-up follows Fleury et al. (2024, their Table 1).

	Mars	Mercury
planet / core radius	3400 / 1850 km ($D=1550$ km)	2440 / 2020 km ($D=420$ km)
surface temperature	220 K	440 K
initial core temperature	2000 K	2000 K
gravity	3.7 m/s²	3.7 m/s²
thermal diffusivity	${10}^{-6}$ m²/s	$1.04\times {10}^{-6}$ m²/s
density core / mantle	6000 / 3500 kg/m³	6980 / 3380 kg/m³
heat capacity mantle / core	1142 / 850 J/kg/K	1142 / 850 J/kg/K
conductivity mantle / crust	4 / 3 W/m/K	4 / 3 W/m/K
present-day heat-producing elements	U 16 ppb, Th 58 ppb, K 309 ppm	U 7 ppb, Th 29 ppb, K 550 ppm
grid	124 × 1650	84 × 1533

Both planets share the rest:

energy	extended Boussinesq: incompressible, constant density, with adiabatic heating and viscous dissipation
viscosity	diffusion creep, $\eta ={\eta }_{\mathrm{ref}}\exp [\frac{E+pV}{RT}-\frac{E+p_{\mathrm{ref}}V}{R{T}_{\mathrm{ref}}}]$, ${\eta }_{\mathrm{ref}}={10}^{21}$ Pa s at 1600 K and 3 GPa, $E=300$ kJ/mol, $V=6$ cm³/mol, bounded to ${10}^{18}$ to ${10}^{26}$ Pa s
boundaries	free slip, fixed surface temperature
core	${\rho }_c{c}_c{V}_{c}d{T}_{\mathrm{CMB}}/dt=-q_c{A}_c$, no inner core
heating	U-238, U-235, Th-232 and K-40, each decaying with its own half-life from 4.5 Gyr before present
crust (crust runs)	a 50 km shell with twice the bulk heat production per kilogram and conductivity 3 W/m/K; the mantle below is depleted by mass balance
initial state	an adiabat from a potential temperature of 1745 K with 50 km boundary layers meeting the surface and core temperatures, plus a 20 K seed and 10 K of noise near the core

How MILEVA represents it: the extended Boussinesq equations are MILEVA's anelastic liquid approximation with a constant reference state. MILEVA's mantle heat capacity is fixed at 1250 J/kg/K, so the run scales density, gravity and heat production so that every product entering the equations (${\rho }_0{c}_p$, ${\rho }_{0}g$, $g/c_p$, the heat produced per unit volume and the conductivity) is the one in the table. The crust is a fixed shell defined by radius: it changes heat production and conductivity and nothing else. There is no melting, no compositional buoyancy and no free surface.

Running it

julia --project=. -t 2 examples/annulus_fleury_thermal_evolution_2d.jl mars 0 124 1650
julia --project=. -t 2 examples/annulus_fleury_thermal_evolution_2d.jl mercury 0 84 1533

Arguments are planet crust Nr Nt, with crust 0 or 1. Nr = 124 for Mars and 84 for Mercury put a cell face exactly at the base of the 50 km crust. Wall time on two threads with other runs sharing the machine: Mars reached 1460 Myr in about 73 h; Mercury reached 4.5 Gyr in about 45 h. A run interrupted by a reboot continues from its last checkpoint with MILEVA_RESUME=latest.

What to expect

Mars. The mantle starts at a volume-mean 1760 K and a stagnant lid grows to about 140 km within the first 400 Myr; surface velocities stay below ${10}^{-4}$ cm/yr, so lid_mobility reads ${10}^{-5}$ throughout. Radiogenic heating is strongest early, and the mantle warms for the first Gyr (1760 K to about 1910 K by 1300 Myr) while the surface heat flow falls from over 100 to about 34 mW/m². Radiogenic heating warms the mantle just above the core faster than the core itself, so heat flows into the core (q_cmb_mWm2 about −8 mW/m²) and the core warms, by about 90 K in the first 1300 Myr. The mantle convects in cells of order the mantle depth under the lid, at about 2 cm/yr.

Mercury. The thin mantle cools from the start and keeps cooling for the whole 4.5 Gyr. The core loses heat throughout (q_cmb_mWm2 between +5 and +10 mW/m² after the first 100 Myr) and cools from 2000 K to about 1710 K. The volume-mean temperature falls from about 1610 K to 1460 K at 1 Gyr and 1060 K at 4.5 Gyr, and the surface heat flow from about 100 to 13 mW/m². The lid grows from about 100 km at 300 Myr to 165 km at 2 Gyr. Convection under it weakens as the mantle cools: the root-mean-square velocity falls from about 0.5 cm/yr in the first few hundred Myr to 0.1 cm/yr at 2 Gyr and below 0.01 cm/yr after about 2.7 Gyr. From then on the mantle only conducts, and the lid (read from the 1300 K isotherm) reaches 270 km of the 420 km mantle by 4.5 Gyr.

The crust. The crust runs put a larger share of the heat production into the top 50 km, where it is conducted out directly instead of warming the convecting mantle. Their comparison with the runs without crust will be added when both reach 4.5 Gyr.

The energy budget in diagnostics.csv (budget_resid) stays within 0.3 % over the whole Mercury run. In the Mars runs single samples after the first Gyr scatter by up to 2 %, with both signs, while their means over 500 Myr windows stay within 0.2 %.