Primordial and recycled material in the lower mantle (2-D)

The lower mantle may hold two kinds of chemically distinct material: primordial rock left over from the planet's early differentiation, and basaltic crust recycled from the surface. Gülcher, Ballmer & Tackley (2021) asked how the two coexist over the age of the Earth and what decides whether the primordial material survives. They started from a dense, viscous, bridgmanite-enriched layer that fills the lower half of the mantle, and followed it for 4.5 Gyr while convection, melting and crustal recycling went on above it. This example repeats their reference case M100dD in the spherical annulus.

Temperature (300 to 4000 K), the fraction of primordial material (0 to 1) and the viscosity (10¹⁹ to 10²⁵ Pa s, log scale) over 4.5 Gyr, Nr = 64 (64 × 512 cells).

Set-up

shell	$r$ = 3481 to 6371 km, free slip, net rotation removed
walls	300 K at the surface, 4000 K at the core-mantle boundary, no internal heating
reference state	anelastic liquid, Perple_X pyrolite table, adiabat from 1900 K
viscosity	Arrhenius diffusion creep with a pressure-dependent yield stress (below)
primordial layer	the lowest 1650 km, 5 % of it ambient mantle, density excess 63.3 kg m⁻³, 100 times more viscous
ambient mantle	15 % basalt, melting and crustal recycling on
initial state	adiabat with thermal boundary layers, a degree-6 perturbation of 20 K and 10 K noise
run	4.5 Gyr

The viscosity is Newtonian and depends on temperature and pressure,

$$\eta ={\eta }_0{\lambda }_C\exp [\frac{E_a+p{V}_a}{RT}-\frac{E_a}{R{T}_0}],{\eta }_{\mathrm{eff}}=min(\eta ,\frac{{\tau }_y}{2{\dot{\epsilon }}_{II}}),{\tau }_y=30\mathrm{MPa}+0.01p,$$

with ${\eta }_0=5\times {10}^{20}$ Pa s at $T_0$ = 1600 K, $E_a$ = 140 kJ mol⁻¹, $V_a$ = 1.8 cm³ mol⁻¹, and ${\lambda }_C$ = 100 for primordial material and 1 otherwise. The yield stress lets a cold lid break where it is stressed enough.

The density excess of the primordial material is set by the buoyancy ratio

$$B=\frac{\mathrm{\Delta }{\rho }_C}{\rho \alpha \mathrm{\Delta }{T}_{\mathrm{sa}}}=0.28,$$

evaluated at 1500 km depth on the reference adiabat ($\rho$ = 4913 kg m⁻³, $\alpha =1.77\times {10}^{-5}$ K⁻¹), with $\mathrm{\Delta }{T}_{\mathrm{sa}}$ = 2597 K the temperature drop across the shell minus the adiabatic increase. Primordial material that crosses the pyroxenite solidus, $T>920{}^{\circ }\mathrm{C}+130{}^{\circ }\mathrm{C}{\mathrm{GPa}}^{-1}p$, melts and becomes 40 % basalt and 60 % harzburgite, as in the paper. It then no longer counts as primordial. The fraction that survives is the preservation factor

$$\chi (t)=\frac{V_{\mathrm{prim}}(t)}{V_{\mathrm{prim}}(0)},$$

with the volume integrated using the shell's measure $r^{2}drd\theta$.

Stagnant lid, overturn and mixing

Temperature (top), primordial material (middle) and viscosity (bottom) at 1.5, 2.25, 2.6, 3.2 and 4.5 Gyr.

A: lid mobility (rms surface velocity over rms velocity of the whole mantle) and rms velocity. B: surface heat flow. C: primordial volume, whole mantle and below 660 km, as a percentage of its initial value. Computed by docs/scripts/fig_primordial.jl from diagnostics.csv and the 245 VTK frames.

The run has three phases.

Two layers under a stagnant lid, 0 to 2.1 Gyr. The upper mantle convects under a lid that does not move (mobility below 0.01), while the dense, viscous primordial layer barely deforms. The layer is heated from the core and loses little heat through its top, so it warms, and the surface heat flow decays from 60 to 22 mW m⁻².

Overturn, 2.1 to 2.6 Gyr. Once the primordial layer is warm enough, its thermal buoyancy overcomes part of its compositional excess. Broad hot domes rise from it, reach the base of the lid and break it. The lid mobility jumps to 2 to 3.5, the rms velocity rises tenfold, and the surface heat flow peaks at 135 mW m⁻² near 2.5 Gyr. The hot primordial material that reaches shallow depth melts, so the preserved volume begins to fall.

Plate-like convection and mixing, 2.6 to 4.5 Gyr. The lid stays mobile (mobility 1.4 to 2.5), the mantle cools from 2180 to 1870 K, and the primordial material is drawn out into ridges on the core-mantle boundary, streaks and smaller blobs. By 4.5 Gyr 18 % of the primordial volume survives in the whole mantle and 13.5 % below 660 km.

The viscosity. During the first two billion years the primordial layer remains at 10²³ to 10²⁴ Pa s, roughly a hundred to a thousand times more viscous than the upper mantle above it. No low-viscosity layer forms above the core, as the hot boundary layer there spans only about one grid cell. After the overturn, primordial streaks retain their high viscosity throughout the mantle, including at shallow depths. By 4.5 Gyr the mantle has cooled and much of it has reached 10²⁴ Pa s or more.

Comparison with the paper

The run follows the same sequence as the paper: two-layer convection under a low-mobility lid, a compositional overturn that switches the lid to plate-like behaviour, and partial preservation of the primordial material afterwards. The timing differs. The overturn here occurs at 2.1 to 2.2 Gyr, about three times later than the 0.68 Gyr of case M100dD, though still within the range of 0.5 to 2.5 Gyr that the paper reports across its cases.

The later overturn leaves less time for mixing, and the mantle at 4.5 Gyr is cooler and slower than in the paper. Averaged over the last 250 Myr, the mean temperature is 1870 K against 2212 K in the paper, and the rms velocity is about 0.5 cm yr⁻¹ against roughly 2.3 cm yr⁻¹ (the paper's value is read from its figure over 4 to 5 Gyr). The lid is somewhat more mobile, with a mobility of 2.0 against 1.3. Slightly less primordial material survives: 18 % of the initial volume at 4.5 Gyr, compared with 23.7 % in the paper.

How the set-up relates to the paper

This example follows the set-up of case M100dD, and MILEVA is an independent code, so the two models share the configuration but not the implementation. The one process taken over directly is the treatment of primordial material: once it crosses the pyroxenite solidus, it is converted to 40 % basalt and 60 % harzburgite, as in the paper. The remaining processes use MILEVA's own formulations, and the main differences are summarised below.

Melting and density. Melting of the ambient mantle follows MILEVA's peridotite scheme rather than the composition-dependent parameterisation that the paper adopts from Nakagawa et al. (2010). Basalt and primordial material are given constant density excesses relative to a pyrolite reference, whereas the paper describes each with its own sequence of phase transitions. As a result, recycled crust in this run does not develop the buoyancy just below 660 km that it has in the paper.

Heat sources. Both models are purely bottom heated, with the core held at 4000 K, and both include adiabatic and viscous heating. This run does not include the latent heat of melting or of the solid phase transitions, nor the heat carried to the surface by erupted melt. MILEVA provides these terms, but they are switched off in this example.

Phase transitions and viscosity. The 660 km transition is represented without a Clapeyron slope, and post-perovskite is not included, so the thousandfold weakening the paper applies near the base of the mantle is absent here. The viscosity increase of primordial material is applied at all depths rather than only in the lower mantle. Where materials are mixed within a cell, viscosities are averaged arithmetically rather than geometrically, which makes partly mixed material noticeably stiffer, by about a factor of five for an even mixture.

Resolution. The radial grid is uniform, with 45 km cells at Nr = 64, while the paper refines its grid to about 10 km near the boundary layers and around 660 km.

The timing of the overturn. Taken together, these differences suggest that the primordial layer warms more slowly than in the paper. The coarser grid at the base and the absence of weak post-perovskite would both reduce the heat it receives from the core, and the stiffer mixing rule would slow the mixing that follows the overturn. These effects have not yet been tested individually. The header of the driver script lists all differences.

Running it

MILEVA_OUTDIR=output/primordial MILEVA_MAXTIME_MYR=4500 MILEVA_VTKEVERY=80 \
    julia --project=. -t 2 examples/annulus_primordial_recycled_2d.jl 64
julia --project=. docs/scripts/fig_primordial.jl output/primordial

Further arguments set $B$, ${\lambda }_C$ and the initial layer thickness in km. The run took 17 hours on 2 threads (19 557 steps).

References

Gülcher, A. J. P., Ballmer, M. D. & Tackley, P. J. (2021), Coupled dynamics and evolution of primordial and recycled heterogeneity in Earth's lower mantle, Solid Earth 12, 2087-2107, doi:10.5194/se-12-2087-2021.