Phase changes

MILEVA describes the solid-state transitions of the olivine system (the 410, 520 and 660 km transitions and post-perovskite), the basalt-to-eclogite transition, and the latent heat of these transitions and of melting. None of them is a marker state. Each is evaluated on the grid every step from the lithostatic pressure and the local temperature:

what	where it acts	function
integer phase index phase_id	gate for the smooth phase buoyancy; latent-heat counter when the smooth latent heat is off	melt_phase (physics/melt_phase.jl)
smooth phase functions ${\mathrm{\Gamma }}_k$	density in the momentum equation	compute_buoyancy (physics/compute_buoyancy.jl)
smooth phase fractions and their material derivative	latent heat in the energy equation	melt_phase
eclogite fraction $f_{\text{ecl}}$	density of basalt below the transition depth	compute_buoyancy

The pressure throughout is the lithostatic mantle_pressure(z, cfg) (physics/melt_phase.jl), the hydrostat of the reference density of the selected compressibility mode (${\rho }_{0}gz$ under Boussinesq). Transition pressures are evaluated with the same function at the transition depth, so a transition sits at its nominal depth on the reference temperature whatever the reference state.

Physical background

What a solid-state phase transition is

Mantle rock is a mixture of minerals, and each mineral is stable only over a limited range of pressure and temperature. When a parcel is carried deep enough, its atoms rearrange into a more compact crystal structure of the same chemical composition. This is a solid-state phase transition. Nothing melts, but the density jumps, heat is taken up or released, and the elastic and seismic properties change. The seismic discontinuities of Earth's mantle near 410 and 660 km depth are the clearest evidence for such transitions. Their depths, widths and density jumps fit transitions in the olivine component of a peridotite mantle (Ringwood 1991; Bina & Helffrich 1994). The Preliminary Reference Earth Model puts density jumps of about 180 kg m${}^{-3}$ at 400 km and about 390 kg m${}^{-3}$ at 670 km (Dziewonski & Anderson 1981), close to the defaults MILEVA uses.

Earth depth	reaction	Clapeyron slope	role
about 410 km	olivine to wadsleyite	positive (exothermic)	helps flow across it
about 520 km	wadsleyite to ringwoodite	positive, small density jump	weak effect
about 660 km	ringwoodite to bridgmanite plus ferropericlase	negative (endothermic)	resists flow across it
a few hundred km above the core	bridgmanite to post-perovskite	large, positive (Hirose 2006)	modifies the lowermost mantle

For Mg${}_2$SiO${}_4$ the dissociation of spinel (ringwoodite) into perovskite and periclase was located experimentally at $P=27.6-0.0028T$ (GPa, with $T$ in °C, between 1000 and 1600 °C), a slope of $-2.8$ MPa K${}^{-1}$, and magnesian spinel was found to dissociate within less than 0.15 GPa at 1600 °C, a few kilometres of depth (Ito & Takahashi 1989). Post-perovskite was discovered in 2004 (Murakami et al.

1. and is thought to be stable in the D″ layer at the base of Earth's mantle

(Hirose 2006).

The density added by the solid-state transitions along the code's 1600 K adiabat, $\sum _k\mathrm{\Delta }{\rho }_k{\mathrm{\Gamma }}_k$, from phase_transitions and phase_fraction with the mesh-derived width of a 192-cell default shell and post-perovskite switched on. The steps of 180, 100, 400 and 50 kg/m³ sum to 730 kg/m³. The 410 km analogue sits slightly below its nominal 453 km because the adiabat there is hotter than its 1600 K pivot. In the buoyancy only the lateral departures from these steps act. Script: docs/scripts/fig_phase_changes_sections.jl.

Venus

The transitions are set by pressure, not by depth. A body with the mineralogy of Earth's mantle but a different gravity has the same transitions at different depths. Venus has a surface gravity of 8.87 m s${}^{-2}$ against Earth's 9.81, and a similar mantle density, so a given pressure lies deeper by roughly the ratio 9.81/8.87 = 1.106. This maps 410 km to 453 km and 660 km to 730 km, which are the default buoyancy depths in MILEVA. Whether the deepest Venus mantle reaches the post-perovskite field depends on its core size and interior temperature, both poorly known; post-perovskite is switched off by default (cfg.usePostPerovskite = false).

mantle_pressure for the default Earth shell (g = 9.81 m/s²) and a Venus shell (R_out = 6052 km, R_in = 3110 km, g = 8.87 m/s²), both under the default ALA reference state. The pressure of 13.9 GPa found at 410 km on Earth is reached at 454 km on Venus, and 22.9 GPa from 660 km at 730 km, which are the default transition depths of 453 and 730 km. Script: docs/scripts/fig_phase_changes_sections.jl.

The Clapeyron slope and boundary topography

Along a phase boundary the two phases have equal Gibbs free energy. Differentiating that equality gives the Clausius-Clapeyron relation,

$$\gamma =\frac{dP}{dT}=\frac{\mathrm{\Delta }s}{\mathrm{\Delta }v}=-\frac{L}{T\mathrm{\Delta }v}\simeq \frac{{\rho }^{2}L}{T\mathrm{\Delta }\rho },\mathrm{\Delta }v\simeq -\frac{\mathrm{\Delta }\rho }{{\rho }^2},$$

where $\mathrm{\Delta }s$ and $\mathrm{\Delta }v$ are the changes of specific entropy and specific volume from the low-pressure to the high-pressure phase, $\mathrm{\Delta }\rho >0$ is the density increase, and $L=-T\mathrm{\Delta }s$ is the latent heat released per unit mass on the downward crossing, the sign convention the code's $L_k$ use (Schubert, Turcotte & Olson 2001). The high-pressure phase is denser, so $\mathrm{\Delta }v<0$ on the way down. If it also has lower entropy, heat is released ($L>0$, exothermic) and $\gamma >0$. If it has higher entropy, heat is absorbed ($L<0$, endothermic) and $\gamma <0$.

The slope decides how the boundary responds to a temperature anomaly. Inside a cold downwelling a boundary with $\gamma >0$ is reached at lower pressure, so it rises; one with $\gamma <0$ is reached at higher pressure, so it is pushed down. The vertical shift is $\delta z\simeq \gamma \delta T/(\rho g)$, the formula in the next section. Seismology maps exactly this topography on the 410 and 660 km discontinuities, and their opposite signs are one of the strongest tests of the phase-change interpretation (Bina & Helffrich 1994).

(A) Depth of each boundary against the temperature anomaly from its pivot, found by inverting mantle_pressure for $P_k+{\gamma }_k\delta T$. The 410 and 520 km analogues move up in cold material and down in hot material, and the 660 km analogue does the opposite. (B) The deflection for $\delta T$ = +200 K against the slope, at the 730 km analogue (line) and for each default transition (dots). The 410 km analogue moves 14 km down and the 660 km analogue 13 km up, close to the estimate $\gamma \delta T/(\rho g)$. Default shell. Script: docs/scripts/fig_phase_changes_sections.jl.

How a phase boundary pushes on the flow

A boundary that is flat everywhere adds the same density at every point of a given radius. That changes the hydrostatic pressure but cannot drive flow, which responds only to lateral density differences (see Buoyancy and composition). A deflected boundary is different. In a cold slab sinking through the 410 km transition, the boundary has risen, so a slice of the slab is already in the dense phase where its surroundings are not. That extra mass pulls the slab down and adds to its negative thermal buoyancy. At 660 km the boundary inside the slab is pushed down, so the slab stays in the light phase while its surroundings at the same depth have already transformed. The slab is then buoyant relative to them and is held back. A hot plume rising from below is resisted in the same way. Inside it the 660 km boundary is raised, so the plume stays in the dense phase over a depth range where its surroundings are already in the light one. The endothermic transition is therefore a partial barrier to vertical flow, and the exothermic ones are accelerators.

The strength of the barrier is measured by the phase buoyancy parameter of Christensen & Yuen (1985). Compare the extra mass per unit area from the deflected boundary, $\mathrm{\Delta }\rho \delta z=\mathrm{\Delta }\rho \gamma \delta T/(\rho g)$, with the thermal mass deficit of a column of the same temperature anomaly over the mantle depth $h$, $\rho \alpha \delta Th$. Their ratio does not depend on $\delta T$:

$$\mathcal{P}=\frac{\gamma \mathrm{\Delta }\rho }{{\rho }^2\alpha gh}.$$

The phase buoyancy parameter $\mathcal{P}$ is positive for an exothermic and negative for an endothermic transition. With the default MILEVA values (${\rho }_0$ = 3300 kg m${}^{-3}$, $\alpha =3\times {10}^{-5}$ K${}^{-1}$, $g$ = 9.81 m s${}^{-2}$, $h$ = 2890 km) the 660 km analogue gives $\mathcal{P}\approx -0.11$ and the 410 km analogue $\mathcal{P}\approx +0.05$. Convection calculations find that a sufficiently negative $\mathcal{P}$ splits the mantle into separately convecting layers, and that intermediate values give intermittent layering, in which cold material accumulates above the 660 km boundary and then drains into the lower mantle in large avalanches (Machetel & Weber 1991; Tackley et al. 1993). The critical value depends on the Rayleigh number, the viscosity structure and the geometry (Christensen 1995; Faccenda & Dal Zilio 2017), which is why the 660 km slope is treated as a sensitivity parameter here.

Radial velocity on the axis of a column 300 K hotter (A) and one 300 K colder (B) than the 1600 K adiabat, from one direct solve of assemble_stokes with the body force of compute_buoyancy and a uniform ${10}^{21}$ Pa s. With thermal buoyancy alone the peak speed is 72.1 cm/yr. The 660 km analogue alone lowers it to 66.6 cm/yr in the plume and in the slab, and the profile bends where the column crosses 730 km. Adding the exothermic 410 km analogue gives part of it back (68.3 cm/yr). Default shell, 96 × 384 cells, an operator illustration rather than a convection result. Script: docs/scripts/fig_phase_changes_sections.jl.

Latent heat

The same entropy change that sets the Clapeyron slope releases or absorbs heat as material crosses the boundary. A parcel that transforms completely changes temperature by roughly $L/c_p$ if no heat is exchanged with its surroundings. Material sinking through the exothermic 410 km transition is warmed, which pushes its uplifted boundary back down and partly offsets the deflection. Material sinking through the endothermic 660 km transition is cooled, which reinforces the barrier. Melting is the same physics with a larger latent heat. Rising mantle that crosses its solidus cools as it melts, and melt that freezes heats its surroundings (Melting). Because $L$ and $\gamma$ are tied through $\mathrm{\Delta }v$ in the Clausius-Clapeyron relation, a transition with a steep slope and a large density jump also has a large latent heat.

$L/c_p$ with $c_p$ = CP_REF = 1250 J/(kg K). The transition values are measured from melt_phase by integrating its latent-heat rate over a column sinking through each boundary at 1 cm/yr, which recovers the 150, 100 and −250 kJ/kg fixed in the function. That is +120 K and +80 K of warming at the exothermic 410 and 520 km analogues, 200 K of cooling at the endothermic 660 km analogue. The melting bar is cfg.latentHeatMelt = 400 kJ/kg, 320 K of cooling for complete melting on the way up. Script: docs/scripts/fig_phase_changes_sections.jl.

Basalt to eclogite

Basaltic crust is not stable at depth either. With increasing pressure, plagioclase breaks down and garnet grows, and basalt or gabbro transforms through garnet granulite into eclogite, a garnet and clinopyroxene rock that is denser than the basalt it formed from and, through most of the upper mantle, denser than the surrounding peridotite (Ringwood & Green 1966). This is what lets oceanic crust sink with a subducting slab on Earth, and it is a candidate trigger for the delamination of thick crust on Venus, where there is no subduction to remove it. The reaction is slow in cold, dry rock, so basalt can survive metastably well below its equilibrium depth and eclogite can survive above it (Ahrens & Schubert 1975). MILEVA treats this transition as a function of depth by default and offers the second kind of metastability as an option, below.

Why a thickening crust can founder. (A) $\delta {\rho }_C$ of compute_buoyancy for a basaltic crust 40, 100 and 160 km thick, with an illustrative crust contrast cfg.crustDensityAnomaly = 300 kg/m³ (the default is zero). Basalt above the eclogite ramp is light and below it heavy. (B) The mass anomaly of the whole column, $\int \delta {\rho }_{C}dz$, against crust thickness. With the illustrative contrast the column is buoyant until the crust is about 170 km thick and heavy beyond. With the default zero contrast any crust thicker than about 60 km adds mass. 1500 radial cells on the default shell. Script: docs/scripts/fig_phase_changes_sections.jl.

Olivine-system transitions

A phase transition is evaluated in a fixed sequence of small steps. The code first finds the pressure at each cell from its depth, then moves the transition pressure up or down with the local temperature (the Clapeyron slope), and turns the distance from that boundary into a smooth fraction $\mathrm{\Gamma }$ of transformed material. Only the part of that fraction that differs from the horizontal average makes a density anomaly, and so a force on the flow. When material moves through the boundary, the same fraction, with the same depth, pivot and width, sets how much latent heat it takes up or gives off. The figure follows the 660 km analogue through these steps for one cold and one hot cell.

The 660 km analogue for one cell row, stage by stage, on the default shell with Nr = 192, a 1600 K potential temperature and the flags cfg.useClaperonBuoyancy, cfg.useMeltingPhase and cfg.usePhaseLatentHeat switched on. Right: the values for a cell at 730 km in a column 200 K colder than the adiabat sinking at 1 cm/yr (1673 K), and in a column 200 K hotter rising at 1 cm/yr (2073 K). The cold cell lies above its deflected boundary (745 km), is lighter than the row mean by 112 kg m${}^{-3}$ and is pushed up; the hot cell (boundary at 719 km) is denser by 116 kg m${}^{-3}$ and is pulled down. The integer phase index follows the same deflected boundary: 3 (ringwoodite) in the cold cell, 4 (bridgmanite) in the hot one. The latent heat uses the boundary, pivot temperature and mesh-derived width of the buoyancy and peaks at it (745 and 715 km); it cools the sinking cell and heats the rising one. The values equal what compute_buoyancy and melt_phase return for those cells. Script: docs/scripts/fig_phase_changes_pipeline.jl.

Phase boundaries

A transition $k$ is a straight line in pressure-temperature space through the reference point $(P_k,T_k)$ with Clapeyron slope ${\gamma }_k=dP/dT$ (Pa K${}^{-1}$):

$$P_{\text{b}}(T)=P_k+{\gamma }_k(T-T_k),P_k=P_{\text{litho}}(z_k).$$

A positive slope (exothermic, 410 and 520 km) puts the boundary deeper in hot material and shallower in cold material; a negative slope (endothermic, 660 km) does the opposite. For a lateral anomaly $\delta T$ the boundary moves by

$$\delta z\simeq \frac{{\gamma }_k\delta T}{\rho g},$$

about 15 km for $\gamma =2.5$ MPa K${}^{-1}$ and $\delta T=200$ K.

$P_{\text{b}}(T)$ for the default transition table of phase_transitions with post-perovskite on, each line passing through its pivot $(T_k,P_k)$ (dots), $P_k$ = mantle_pressure$(z_k)$. The dashed line is the code's adiabat at 1600 K potential temperature. Where it crosses a boundary line is where that transition sits in an adiabatic column. Default shell. Script: docs/scripts/fig_phase_changes_sections.jl.

Phase index

When cfg.useMeltingPhase = true, melt_phase assigns every cell an integer phase (1 olivine, 2 wadsleyite, 3 ringwoodite, 4 bridgmanite, 5 post-perovskite) by testing $P>P_{\text{b}}(T)$ for each transition in turn. The depths, slopes and pivots are those of the phase buoyancy below (phase_transitions, from mat and, for post-perovskite, cfg), so the phase index, the phase buoyancy and the phase latent heat sit on one boundary:

transition	$z_k$	${\gamma }_k$ (Pa K${}^{-1}$)	$T_k$
olivine to wadsleyite	453 km	$+2.5\times {10}^6$	1600 K
wadsleyite to ringwoodite	520 km	$+4.0\times {10}^6$	1600 K
ringwoodite to bridgmanite	730 km	$-2.5\times {10}^6$	1900 K
bridgmanite to post-perovskite	cfg.ppvDepth_km (2700 km)	cfg.ppvClapeyron (${10}^7$)	1600 K

cfg.useEarthPhaseBoundaries = false leaves every cell in phase 1 (for a body with a different mineralogy), and post-perovskite is assigned only with cfg.usePostPerovskite = true. The phase index does not enter the density: a step function of the cell index would put a whole density jump into a single cell. Its role is to switch on the smooth 520 km and post-perovskite terms below (they are added only when a phase field is passed to compute_buoyancy), and to drive the integer latent-heat counter when the smooth latent heat is off.

(A) phase_id returned by melt_phase for a section with a column 300 K hotter (left) and one 300 K colder (right) than the 1600 K adiabat, on the default shell with 192 × 768 cells. The boundaries are deflected exactly as the buoyancy's are, but in whole cells. (B) The index along the plume, the adiabat and the slab. The hot column enters wadsleyite and ringwoodite deeper and bridgmanite shallower than the adiabat, the cold column the reverse. Script: docs/scripts/fig_phase_changes_sections.jl.

Phase buoyancy

The density enters through smooth phase functions (Christensen & Yuen 1985),

$${\mathrm{\Gamma }}_k=\frac{1}{2}[1+\tanh (\frac{P-P_k-{\gamma }_k(T-T_k)}{\delta P})],\delta {\rho }_{\text{phase}}=\sum _k\mathrm{\Delta }{\rho }_k({\mathrm{\Gamma }}_k-⟨{\mathrm{\Gamma }}_k{⟩}_r),$$

with $⟨\cdot {⟩}_r$ the horizontal mean at each radius. A laterally uniform pressure gives ${\mathrm{\Gamma }}_k$ no lateral anomaly, so a transition drives flow only through the temperature deflection of its boundary. Behind an exothermic boundary a cold downwelling transforms early and is pulled down; behind the endothermic 660 km boundary it transforms late and is held back.

symbol	meaning	in the code
$P$	lithostatic pressure of the cell	mantle_pressure(depth, cfg)
$P_k$	transition pressure at the reference temperature, $P_{\text{litho}}(z_k)$	mantle_pressure(z_k, cfg), plus the optional basalt shift
$T$	cell temperature as stored on the grid	Tg
$T_k$	pivot temperature, at which the boundary sits exactly at $z_k$	mat.clap_T_ref_*
${\gamma }_k$	Clapeyron slope $dP/dT$	mat.gamma_*
$\delta P$	pressure half-width of the tanh	_clap_width
${\mathrm{\Gamma }}_k$	fraction of material in the high-pressure phase, 0 above, 1 below	phase_fraction
$\mathrm{\Delta }{\rho }_k$	density increase across the transition	mat.drho_*
$⟨{\mathrm{\Gamma }}_k{⟩}_r$	mean of ${\mathrm{\Gamma }}_k$ over the cells of the same radius	mean(Γ, dims=2) (2-D), dims=(2,3) (3-D)

The argument of the tanh is the pressure distance of the cell from its own deflected boundary, $P-P_{\text{b}}(T)$, in units of $\delta P$. A cell exactly on its boundary has ${\mathrm{\Gamma }}_k=1/2$; one more than about $2\delta P$ beyond it is fully transformed. Subtracting the horizontal mean removes the part of the jump that every column shares, which would only alter the hydrostatic pressure, and keeps the part that drives flow.

transition	depth $z_k$	${\gamma }_k$ (Pa K${}^{-1}$)	$T_k$	$\mathrm{\Delta }{\rho }_k$ (kg m${}^{-3}$)	fields
410 km analogue, exothermic	453 km	$+2.5\times {10}^6$	1600 K	180	mat.clap_depth_410, gamma_410, clap_T_ref_410, drho_410
520 km analogue	520 km	$+4.0\times {10}^6$	1600 K	100	mat.clap_depth_520, gamma_520, clap_T_ref_520, drho_520
660 km analogue, endothermic	730 km	$-2.5\times {10}^6$	1900 K	400	mat.clap_depth_660, gamma_660, clap_T_ref_660, drho_660
post-perovskite	cfg.ppvDepth_km (2700 km)	cfg.ppvClapeyron (${10}^7$)	1600 K	50	cfg.ppvDepth_km, cfg.ppvClapeyron, mat.clap_T_ref_ppv, mat.ppv_drho

The default depths 453 and 730 km are the Earth transition pressures mapped through Venus gravity (config/MaterialProperties.jl), and the density jumps follow Tian et al. (2023) Table A.1. The phase index and the latent heat use the same table. For an Earth model set mat.clap_depth_410 = 410e3 and mat.clap_depth_660 = 660e3. The 660 km slope is a sensitivity parameter: published values span roughly $-1$ to $-3$ MPa K${}^{-1}$.

Which terms are active.

flags	transitions in the buoyancy
cfg.useClaperonBuoyancy = true	410 and 660 analogues
cfg.usePhaseTransitionDensity = true (default) and a phase field, i.e. cfg.useMeltingPhase = true	520 analogue; the 410 and 660 analogues as well if the Clapeyron term is off; post-perovskite if cfg.usePostPerovskite

Each transition is counted once. useClaperonBuoyancy and useMeltingPhase both default to false, so a default run has no phase buoyancy. With cfg.compClapeyron_dP_410 and compClapeyron_dP_660 the 410 and 660 transition pressures of the Clapeyron term shift by that many pascals per unit basalt fraction (zero by default).

Width. The real transitions are a few kilometres wide, below any convection cell, so $\delta P$ is a numerical parameter. With mat.clap_width_Pa = 0 (default) it is derived from the mesh so that the 10-90 % width spans mat.clap_width_cells (default 3) radial cells at each transition's own depth, $\delta P_k=n_{\text{cells}}\mathrm{\Delta }r(d{P}_{\text{litho}}/dz)/2.2$ (_clap_width), where $\mathrm{\Delta }r$ is the radial cell width grid.drVec at $z_k$, $d{P}_{\text{litho}}/dz$ is the gradient of the same lithostatic model that gives $P_k$ (mantle_pressure), evaluated at $z_k$, and $2.2\approx 2\mathrm{artanh}(0.8)$ is the 10-90 % width of $\frac{1}{2}(1+\tanh x)$ in units of $x$; a positive mat.clap_width_Pa fixes it. The code warns once when the 660 km transition spans fewer than two cells.

Left: the phase index returned by melt_phase on a grid of depth and temperature (shading), with the ${\mathrm{\Gamma }}_k=1/2$ surfaces of the buoyancy phase functions for the default material fields (lines) and the code's adiabat at 1600 K potential temperature, on the default shell. Right: the phase part of compute_buoyancy (Clapeyron plus 520 km terms) for columns 200 K warmer and colder than an adiabatic column, relative to it, on a 300-cell radial grid. Dotted lines: the deflected boundaries. The exothermic 410 and 520 km analogues make the hot column lighter, the endothermic 660 km analogue makes it denser.

Basalt to eclogite

Basalt (the marker field basalt, projected to the cells) transforms to eclogite below cfg.eclogiteTransitionDepth (60 km) over a 20 km ramp:

$$f_{\text{ecl}}(z)=\mathrm{clamp}(\frac{z-z_{\text{ecl}}}{20\text{km}},0,1),\delta {\rho }_{\text{ecl}}=\mathrm{\Delta }{\rho }_{\text{ecl}}{C}_{\text{bas}}{f}_{\text{ecl}},$$

with $\mathrm{\Delta }{\rho }_{\text{ecl}}$ = cfg.basaltDensityAnomaly (200 kg m${}^{-3}$). The transition is a function of depth only, so it is instantaneous and reversible: basalt brought back above $z_{\text{ecl}}$ is neutral again. With cfg.useEclogiteKinetics = true (default false), eclogite above $z_{\text{ecl}}$ stays dense while it is cold and back-transforms smoothly around cfg.eclogiteBackT_K (1000 K) over cfg.eclogiteBackWidth_K (100 K):

$$f_{\text{ecl}}\leftarrow min\{1,f_{\text{ecl}}+{\mathbb{1}}_{z<z_{\text{ecl}}}\frac{1}{2}[1-\tanh (\frac{T-T_{\text{back}}}{W_{\text{back}}})]\}.$$

The same 20 km ramp defines the crust fraction used for crustal buoyancy and viscosity (Buoyancy and composition). The transition carries no latent heat. Its consequences for the melt cycle are on the Melting page.

(A) $f_{\text{ecl}}$ recovered from compute_buoyancy (compositional part divided by cfg.basaltDensityAnomaly) for a basalt fraction of 1, with cfg.useEclogiteKinetics = true. Below 80 km the basalt is fully eclogite, across the 60 to 80 km ramp (dotted) it follows the depth, and above 60 km it stays eclogite while colder than about 1000 K (dashed, cfg.eclogiteBackT_K). (B) The same at 43 km. With the option off (the default) basalt above the transition is never eclogite, and with it on the fraction falls from 1 to 0 over about ±200 K around 1000 K. 1500 radial cells. Script: docs/scripts/fig_phase_changes_sections.jl.

Latent heat

Both latent-heat sources are volumetric rates added to the energy equation; melt_phase returns their sum when cfg.useMeltingPhase = true. They are evaluated from the temperature at the start of this step and of the previous one, and the driver rescales the rate by $\mathrm{\Delta }{t}_{\text{prev}}/\mathrm{\Delta }t$ so that the energy applied over the step equals the energy booked.

Melting

$$Q_L=-{\rho }_{0}L\frac{DF}{Dt}=-{\rho }_{0}L(\frac{\partial F}{\partial t}+\mathbf{v}\cdot \mathrm{\nabla }F),$$

with $L$ = cfg.latentHeatMelt ($4\times {10}^5$ J kg${}^{-1}$) and $F$ the grid melt fraction of melt_phase (linear between the Hirschmann solidus and the Andrault liquidus, see Melting). Melting cools, freezing heats. $\partial F/\partial t$ differences the two start-of-step fields, and $\mathbf{v}\cdot \mathrm{\nabla }F$ is evaluated upwind; in a statistically steady state the Eulerian part alone is close to zero even where material melts on its way in. When the marker melting scheme owns the latent heat (cfg.trackBasalt and cfg.useHeatPipe, 2-D), it books $-{\rho }_{0}L$ per unit of extracted melt at the source (basalt_evolution!), and $Q_L$ is zeroed above cfg.meltExtractMaxDepth_km so the two do not count the same melt.

A column on the code's 1600 K adiabat rising at 1 cm/yr, with steady temperature so that only the advective part of $DF/Dt$ acts. (A) The adiabat crosses solidus_hirschmann near 55 km (liquidus_andrault in grey). (B) The grid melt fraction of melt_phase rises to 0.33 at the surface. (C) $Q_L$ returned by melt_phase is negative wherever $F$ grows upward, about −2.6 µW/m³, a cooling of about 20 K per Myr. The marker heat pipe is off here (cfg.trackBasalt = false). Default shell, 600 radial cells. Script: docs/scripts/fig_phase_changes_sections.jl.

Solid-state transitions

With cfg.usePhaseLatentHeat = true (default false) each transition contributes

$$Q_k={\rho }_0{L}_k\frac{D{\mathrm{\Gamma }}_k}{Dt},{\mathrm{\Gamma }}_k=\frac{1}{2}[1+\tanh (\frac{P-P_k-{\gamma }_k(T-T_k)}{\delta P_k})],$$

with ${\mathrm{\Gamma }}_k$ the phase function of the buoyancy (same depths, slopes, pivots and width $\delta P_k$), and $D{\mathrm{\Gamma }}_k/Dt$ the change between the two start-of-step temperature fields plus a centred radial $v_r\partial {\mathrm{\Gamma }}_k/\partial r$. The latent heats are fixed in the function:

transition	$L_k$ (J kg${}^{-1}$)	sign of $Q_k$ for downward crossing
410 km	$+1.5\times {10}^5$	heating (exothermic)
520 km	$+1.0\times {10}^5$	heating
660 km	$-2.5\times {10}^5$	cooling (endothermic)
post-perovskite (with cfg.usePostPerovskite)	$+0.5\times {10}^5$	heating

With the flag off, melt_phase instead counts integer phase-index changes between the two temperature fields and books $\pm {\rho }_0{L}_k/\mathrm{\Delta }t$ in a cell whose index changed. A boundary must then migrate about half a cell before any cell registers, so on mantle grids this counter is usually zero; switch the flag on whenever solid-state latent heat matters.

melt_phase for a column on the 1600 K adiabat sinking (A) and rising (B) at 1 cm/yr, with steady temperature, on the default shell with 192 radial cells. With cfg.usePhaseLatentHeat = true the sinking column is heated at the 410 and 520 km analogues and cooled at the 660 km analogue (about −5.5 µW/m³), the rising one the reverse; each peak has the mesh-derived width of the buoyancy. With the flag off the integer counter returns zero everywhere, because no cell changes its phase index while the temperature is steady. Script: docs/scripts/fig_phase_changes_sections.jl.

Assumptions and simplifications

• Univariant boundaries. Each transition is a single line in pressure-temperature space for one component. Real mantle transitions are divariant loops in a multicomponent rock, and near 660 km the garnet component transforms in parallel with ringwoodite. The tanh width stands for all of this, and it is set by the mesh rather than by the mineralogy.

• Pressure is lithostatic. The boundaries are located on the reference hydrostat, so dynamic pressure does not move them. The built-in reference densities that set the hydrostat carry no phase jumps; a Perple_X reference table (cfg.referenceTable) supplies its own pressure profile.

• Buoyancy only. The phase functions change the density in the momentum equation. They do not change the viscosity, thermal expansivity, heat capacity or conductivity. A viscosity increase into the lower mantle is a separate, depth-fixed option on the Rheology page (cfg.use660kmJump).

• One description of each boundary. The buoyancy, the phase index and the latent heat read the same depths, slopes and pivots (phase_transitions), and the buoyancy and the latent heat the same mesh-derived width. A wide width on a coarse grid therefore spreads the latent heat as far as the buoyancy, and transitions closer together than their widths (the 453 and 520 km analogues) overlap in both. The one exception is the optional basalt shift of the 410 and 660 transition pressures (cfg.compClapeyron_dP_410, _660, zero by default), which moves the buoyancy's boundary only.

• Latent heats are independent parameters. The values $L_k$ are fixed in melt_phase and are not derived from ${\gamma }_k$ and $\mathrm{\Delta }{\rho }_k$. With the default slopes and density jumps, ${\rho }_0$ = 3300 kg m${}^{-3}$ and 1600 K, the Clausius-Clapeyron relation would give about 66, 59 and $-147$ kJ kg${}^{-1}$ for the 410, 520 and 660 km transitions, against the 150, 100 and $-250$ kJ kg${}^{-1}$ in the code.

• Eclogite by depth. The basalt-to-eclogite transition depends on depth only, not on temperature or reaction time, apart from the optional cold metastable eclogite above the transition. Metastable basalt below it is not modelled.

(A) The latent heats fixed in melt_phase (measured as in the latent-heat figure above) against ${\gamma }_{k}T\mathrm{\Delta }{\rho }_k/{\rho }_0^2$ at 1600 K from the default slopes and density jumps. The code values are larger by factors of 2.3, 1.7 and 1.7. (B) $\mathrm{\Gamma }$ of the 660 km analogue along the adiabat on the default shell at Nr = 64, 128 and 256. The 10-90 per cent width is three cells in each case, 135, 68 and 34 km, so the transition sharpens as the mesh is refined and is never as narrow as the few kilometres of the real transition. Script: docs/scripts/fig_phase_changes_sections.jl.

2-D and 3-D

	2-D annulus	3-D Yin-Yang
phase index	yes	yes
phase buoyancy, mesh-derived width	yes, lateral mean over $\theta$	yes, lateral mean over $\theta$ and $\varphi$
eclogite gate and kinetics	yes	yes
crust term using the eclogite ramp	yes	no
advective part of $DF/Dt$	upwind, radial and azimuthal	upwind, radial only
advective part of $D{\mathrm{\Gamma }}_k/Dt$	centred, radial only	centred, radial only
latent heat of melting on the grid	zeroed above the extraction depth when the marker heat pipe is on	everywhere (no marker heat pipe in 3-D)

Under ALA and TALA the driver replaces the thermal buoyancy with its reference-state form and keeps the phase and compositional terms of compute_buoyancy unchanged (Buoyancy and composition).

The thermal and phase density of one column 200 K hotter than the 1600 K adiabat, from the 2-D dispatch of compute_buoyancy on an annulus of 432 columns and from the 3-D dispatch on a Yin patch of 12 × 36 columns, so that both lateral means are taken over 432 cells. Clapeyron buoyancy on, 64 radial cells on the default shell. The two profiles agree to the last bit. Script: docs/scripts/fig_phase_changes_sections.jl.

Configuration

field	default	meaning
cfg.useMeltingPhase	false	compute melt fraction, phase index and latent heat (melt_phase)
cfg.useClaperonBuoyancy	false	410 and 660 analogues in the buoyancy
cfg.usePhaseTransitionDensity	true	smooth 520 km (and post-perovskite) terms when a phase field exists
cfg.useEarthPhaseBoundaries	true	assign phases 2 to 4 in the phase index
cfg.usePostPerovskite	false	post-perovskite in the phase index, buoyancy and latent heat
cfg.ppvDepth_km, cfg.ppvClapeyron	2700 km, ${10}^7$ Pa K${}^{-1}$	post-perovskite reference depth and slope
cfg.usePhaseLatentHeat	false	smooth solid-state latent heat instead of the integer counter
cfg.latentHeatMelt	$4\times {10}^5$ J kg${}^{-1}$	latent heat of melting
cfg.compClapeyron_dP_410, compClapeyron_dP_660	0, 0 Pa	transition shift per unit basalt fraction
mat.clap_depth_410, clap_depth_520, clap_depth_660	453, 520, 730 km	transition depths
mat.gamma_410, gamma_520, gamma_660	$2.5$, $4.0$, $-2.5$ MPa K${}^{-1}$	Clapeyron slopes
mat.clap_T_ref_410, clap_T_ref_520, clap_T_ref_660, clap_T_ref_ppv	1600, 1600, 1900, 1600 K	pivot temperatures
mat.drho_410, drho_520, drho_660, ppv_drho	180, 100, 400, 50 kg m${}^{-3}$	density jumps
mat.clap_width_Pa, mat.clap_width_cells	0, 3	transition width, explicit or in cells
cfg.eclogiteTransitionDepth	60 km	top of the 20 km eclogite ramp
cfg.basaltDensityAnomaly	200 kg m${}^{-3}$	eclogite density excess
cfg.useEclogiteKinetics, eclogiteBackT_K, eclogiteBackWidth_K	false, 1000 K, 100 K	metastable eclogite above the transition

Where in the code

What	File, function
lithostatic pressure	physics/melt_phase.jl: mantle_pressure
phase index, latent heat of melting and of transitions	physics/melt_phase.jl: melt_phase, _phase_index, _radial_column, _advect_upwind, _advect_gamma, _apply_phase_latent!
transition table, phase function, phase buoyancy, width, eclogite gate	physics/compute_buoyancy.jl: phase_transitions, phase_fraction, compute_buoyancy, _phase_density, _clap_width, _warn_clap_resolution
marker latent heat of melting	physics/basalt_evolution.jl: basalt_evolution!
call order, rescaling by $\mathrm{\Delta }{t}_{\text{prev}}/\mathrm{\Delta }t$	main_SphCon.jl, steps 3 and 4 and the source assembly before diffuse_temperature
material defaults	config/MaterialProperties.jl

Tests

• tests/test_phase_transition_resolution.jl: the width is derived from the mesh and spans the requested cells at several resolutions, 2-D and 3-D agree, the buoyancy stays smooth under refinement and is a continuous function of temperature, and each transition is counted once.

• tests/test_config_timeloop.jl: a run with Clapeyron buoyancy, and the combined configuration, survive the time loop.

• tests/test_melt_frame.jl: melt_phase compares the solidus with the real temperature in the potential-temperature frame.

• tests/test_phase_transition_consistency.jl: every transition parameter is a mat field, a call without mat uses its defaults, the width follows the model's pressure and each transition's own depth, the latent heat and phase index sit on the buoyancy's boundary, and melt_phase measures depth from cfg.R_out.

Key references

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

• Ringwood, A. E. (1991). Phase transformations and their bearing on the constitution and dynamics of the mantle. Geochim. Cosmochim. Acta 55, 2083-2110. doi:10.1016/0016-7037(91)90090-R

• Murakami, M., Hirose, K., Kawamura, K., Sata, N. & Ohishi, Y. (2004). Post-perovskite phase transition in MgSiO3. Science 304, 855-858. doi:10.1126/science.1095932

• Tackley, P. J. (2008). Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid. PEPI 171, 7-18. doi:10.1016/j.pepi.2008.08.005

• Tian, J., Tackley, P. J. & Lourenço, D. L. (2023). The tectonics and volcanism of Venus: new modes facilitated by realistic crustal rheology and intrusive magmatism. Icarus 399, 115539. doi:10.1016/j.icarus.2023.115539

Further reading

• Ahrens, T. J. & Schubert, G. (1975). Gabbro-eclogite reaction rate and its geophysical significance. Rev. Geophys. 13, 383-400. doi:10.1029/RG013i002p00383

• Bina, C. R. & Helffrich, G. (1994). Phase transition Clapeyron slopes and transition zone seismic discontinuity topography. JGR 99, 15853-15860. doi:10.1029/94JB00462

• Christensen, U. (1995). Effects of phase transitions on mantle convection. Annu. Rev. Earth Planet. Sci. 23, 65-87. doi:10.1146/annurev.ea.23.050195.000433

• Dziewonski, A. M. & Anderson, D. L. (1981). Preliminary reference Earth model. Phys. Earth Planet. Inter. 25, 297-356. doi:10.1016/0031-9201(81)90046-7

• Faccenda, M. & Dal Zilio, L. (2017). The role of solid-solid phase transitions in mantle convection. Lithos 268-271, 198-224. doi:10.1016/j.lithos.2016.11.007

• Hirose, K. (2006). Postperovskite phase transition and its geophysical implications. Rev. Geophys. 44. doi:10.1029/2005RG000186

• Ito, E. & Takahashi, E. (1989). Postspinel transformations in the system Mg2SiO4-Fe2SiO4 and some geophysical implications. JGR 94, 10637-10646. doi:10.1029/JB094iB08p10637

• Machetel, P. & Weber, P. (1991). Intermittent layered convection in a model mantle with an endothermic phase change at 670 km. Nature 350, 55-57. doi:10.1038/350055a0

• Ringwood, A. E. & Green, D. H. (1966). An experimental investigation of the gabbro-eclogite transformation and some geophysical implications. Tectonophysics 3, 383-427. doi:10.1016/0040-1951(66)90009-6

• Schubert, G., Turcotte, D. L. & Olson, P. (2001). Mantle Convection in the Earth and Planets. Cambridge University Press. doi:10.1017/CBO9780511612879

• Tackley, P. J., Stevenson, D. J., Glatzmaier, G. A. & Schubert, G. (1993). Effects of an endothermic phase transition at 670 km depth in a spherical model of convection in the Earth's mantle. Nature 361, 699-704. doi:10.1038/361699a0