Layering by an endothermic phase change (2-D)

A solid-state phase change in the mantle releases or absorbs latent heat and changes the density of the rock that crosses it. Whether the boundary helps or hinders flow depends on the sign of its Clapeyron slope $\gamma =dP/dT$. Cold downwellings reach an exothermic transition ($\gamma >0$) early and are pulled down by the dense phase, while at an endothermic transition ($\gamma <0$, the 660 km ringwoodite to bridgmanite reaction) they are held back. A strong enough endothermic transition stops material crossing it, and the mantle convects as two separate layers. Christensen & Yuen (1985) mapped this transition in terms of one number, the phase buoyancy parameter $P$. This example repeats their experiment in the spherical annulus at four values of $P$.

P = +0.4, exothermic

P = −0.2

P = −0.4

P = −0.8, layered

Temperature (300 to 1300 K) of the four runs, Ra = 10⁵, Nr = 64 (64 × 512 cells), to 1585 Myr (0.05 thermal diffusion times). The transition sits at 500 km depth, half way through the 1000 km shell.

The phase transition

The fraction of the high-pressure phase is a smooth step in pressure, shifted by temperature along the Clapeyron slope (Christensen & Yuen 1985):

$$\mathrm{\Gamma }=\frac{1}{2}[1+\tanh \frac{p-p_t-\gamma (T-T_{\mathrm{ref}})}{\delta p}],\rho ={\rho }_0[1-\alpha (T-T_0)]+\mathrm{\Gamma }\mathrm{\Delta }\rho .$$

Here $p_t$ is the pressure at 500 km, $T_{\mathrm{ref}}$ = 800 K is the mean of the wall temperatures, and the half-width $\delta p$ corresponds to 50 km. Where a cold downwelling enters a transition with $\gamma <0$, the boundary is deflected downwards. The downwelling is then surrounded by dense phase while it is still the light phase itself, and that buoyancy pushes it back up.

Two dimensionless numbers measure the effect. The phase Rayleigh number $R_b=\mathrm{\Delta }\rho g{D}^3/(\kappa \eta )$ compares the density jump with viscous resistance, the way $\mathrm{Ra}$ does for thermal buoyancy. The phase buoyancy parameter

$$P=\frac{R_b}{\mathrm{Ra}}\frac{\gamma \mathrm{\Delta }T}{{\rho }_{0}gD}$$

combines the density jump with the boundary deflection $\gamma \mathrm{\Delta }T/({\rho }_{0}g)$ in units of the depth. The driver fixes $R_b=2\mathrm{Ra}$, which gives $\mathrm{\Delta }\rho =2{\rho }_0\alpha \mathrm{\Delta }T$ = 198 kg m⁻³, and sets the slope from $P$:

$P$	+0.4	−0.2	−0.4	−0.8
$\gamma$ (MPa K⁻¹)	+6.47	−3.24	−6.47	−12.95

The latent heat of the transition is not included, so the experiment isolates its buoyancy.

Whole-mantle or layered

Averages over the last third of each run's frames, computed by docs/scripts/fig_phase_layering.jl. A: horizontally averaged temperature. B: rms radial velocity on each ring of grid nodes, divided by its largest value over depth. Dashed line: the transition depth.

$P$	radial flow at 500 km / maximum	${\mathrm{Nu}}_{\mathrm{top}}$	$V_{\mathrm{rms}}$ (cm yr⁻¹)	regime
+0.4	1.00	8.80	0.56	whole mantle
−0.2	0.92	8.75	0.27	whole mantle
−0.4	0.94	8.65	0.34	whole mantle
−0.8	0.26	4.71	0.19	layered

Nusselt number and rms velocity: means over the last third of the rows of diagnostics.csv.

With the exothermic transition the radial flow peaks right at the transition and the cells span the whole shell. A weak endothermic transition ($P$ = −0.2 and −0.4) leaves a small notch in the radial flow at 500 km, where some downwellings pause before they break through, but the mantle still convects as one layer and the heat flow is almost unchanged.

At $P$ = −0.8 the picture changes. The radial flow drops to about a tenth of its maximum just below the transition, so two stacks of cells form, one above the other. A thermal boundary layer grows on each side of the transition, and the mean temperature jumps by about 500 K across it. Heat must now cross the transition by conduction, an additional thermal resistance in series with the two outer boundary layers, so the surface Nusselt number falls to about half its whole-mantle value.

The runs are short (0.05 diffusion times), so these are snapshots of the regime, not converged statistics. Where exactly the change to layering happens also depends on $\mathrm{Ra}$ (Christensen & Yuen 1985), so the threshold between $P$ = −0.4 and −0.8 found here applies to Ra = 10⁵ in this shell.

Running it

for p in 0.4 -0.2 -0.4 -0.8; do
  MILEVA_OUTDIR=output/phase_P$p MILEVA_VTKEVERY=100 julia --project=. -t 2 \
      examples/annulus_phase_layering_2d.jl $p 1e5 64
done
julia --project=. docs/scripts/fig_phase_layering.jl output

Arguments are $P$, the Rayleigh number and Nr. Each run took 13 to 26 minutes on 2 threads.

References

Christensen, U. R. & Yuen, D. A. (1985), Layered convection induced by phase transitions, J. Geophys. Res. 90, 10291-10300, doi:10.1029/JB090iB12p10291.