Two-phase flow of melt and solid

With cfg.twoPhaseFlow = true (default false) MILEVA treats partially molten mantle as two interpenetrating phases. The solid matrix creeps as before, the melt percolates through it by porous (Darcy) flow, and the matrix compacts or dilates to make room. The porosity $\varphi$, the volume fraction of melt, becomes a field on the grid that is carried by the solid, pushed by the melt flux and fed by melting. This replaces the prescribed extraction law of Melting and crust production, in which melt above a retained porosity leaves its marker over a relaxation time and appears at fixed depths of its own column. The implementation is 2-D, spherical annulus and Boussinesq; validate! refuses every other combination (see Configuration).

Physical background

Why melt moves relative to the rock

Partial melt in the mantle sits along grain edges, where it forms a connected network of tubes already at a fraction of a per cent of melt. Basaltic melt is some 300 to 600 kg m${}^{-3}$ less dense than the peridotite around it at upper-mantle pressures (Stolper et al. 1981), and its viscosity, of order 1 to 10 Pa s, is about twenty orders of magnitude below that of the matrix. The melt therefore flows through the pore network as a fluid through a porous medium, driven by its buoyancy and by pressure differences in the melt, and its flux relative to the solid follows Darcy's law (McKenzie 1984):

$$\mathbf{q}=\varphi ({\mathbf{v}}_f-{\mathbf{v}}_s)=-\frac{k(\varphi )}{\mu }(\mathrm{\nabla }{p}_f-{\rho }_f\mathbf{g}).$$

The permeability $k$ grows steeply with porosity, because the tubes widen and their number does not change. Measurements on texturally equilibrated aggregates are written $k=d^2{\varphi }^n/C$ with grain size $d$. For olivine-basalt aggregates, Stokes-flow simulations on X-ray microtomography images of the melt network give $n=2.6\pm 0.2$ and $C={58}_{-22}^{+36}$ for melt fractions of 0.02 to 0.20, with an assumed grain size of 35 μm (Miller et al. 2014). Direct measurements on quartzite with an aqueous fluid give $k\approx d^2{\varphi }^3/200$ (Wark & Watson 1998), the cubic law that the default exponent $n=3$ follows. For $d=1$ mm and $\varphi =0.01$ the two give $5\times {10}^{-15}$ to ${10}^{-13}$ m${}^2$, so melt rises at centimetres to metres per year, much faster than the mantle convects.

A: the default permeability $k_0{\varphi }^3$ (from darcy_mobility times $\mu$) beside the laws of Miller et al. (2014) and Wark & Watson (1998) for a 1 mm grain size. B: the speed of buoyant melt relative to the solid, $\mathrm{\Delta }\rho g{K}_D/\varphi$, for Venus ($\mathrm{\Delta }\rho$ = 500 kg m${}^{-3}$, $g$ = 8.87 m s${}^{-2}$, $\mu$ = 1 Pa s); with the default $k_0$ it is 1.4 cm yr${}^{-1}$ at 1 per cent porosity, inside the range of mantle flow, and a hundred times more with $k_0={10}^{-7}$ m${}^2$. Script: docs/scripts/fig_two_phase_sections.jl.

Compaction

The melt can only leave if the matrix closes behind it, and it can only enter dry rock if the matrix opens. The solid is therefore compressible as a two-phase aggregate even though each phase is incompressible. The volume change of the matrix, $\mathrm{\nabla }\cdot {\mathbf{v}}_s$, is resisted by a bulk (compaction) viscosity $\zeta$, and the stress that resists it is the compaction pressure, the difference between the solid and the melt pressure,

$$p_c=-\zeta \mathrm{\nabla }\cdot {\mathbf{v}}_s.$$

Theory and experiments place $\zeta$ at the shear viscosity divided by the porosity, $\zeta \sim \eta /\varphi$ (Simpson, Spiegelman & Weinstein 2010), so a nearly dry matrix is very hard to compact. Melt also weakens the matrix in shear. In experiments on olivine-basalt aggregates the creep viscosity falls as $\exp (-{\alpha }_{\varphi }\varphi )$ with ${\alpha }_{\varphi }$ between 25 and 30, i.e. about a factor of 1.3 per per cent of melt (Mei et al. 2002; Hirth & Kohlstedt 2003).

A laterally uniform column of the Venus-sized annulus (Nr = 96, Nt = 8, uniform 30.6 km cells, Boussinesq, MUMPS) with ${\eta }_s={10}^{21}$ Pa s, $k_0={10}^{-7}$ m${}^2$ and a band of 1 per cent porosity at 1500 km on a 0.2 per cent background, solved once by solve_stokes_darcy!; an operator illustration, not a run. A: the porosity. B: melt rises through the band at up to 1 cm yr${}^{-1}$ and the matrix moves down at the same rate, since the column is closed. C: the compaction pressure is negative above the band, where the matrix dilates to let melt in, and positive below it, where the matrix compacts behind the melt. Script: docs/scripts/fig_two_phase_sections.jl.

The compaction length

The balance of Darcy flow against matrix compaction defines a natural length, the compaction length (McKenzie 1984)

$$\delta =\sqrt{\frac{k}{\mu }(\zeta +\frac{4}{3}\eta )}.$$

Over distances shorter than $\delta$ the matrix behaves as rigid and the melt pressure gradient adjusts to keep the flow continuous. Over distances longer than $\delta$ the melt rises at its buoyancy-driven speed and the matrix compacts as needed. All the structure that two-phase flow adds to mantle convection, compaction boundary layers at the top of a melting region, porosity waves, melt channels and the focusing of melt towards ridges or the base of a lid, has the scale $\delta$ (Spiegelman 1993; Katz 2008). With $\eta ={10}^{19}$ to ${10}^{21}$ Pa s, $\varphi =0.01$ and $k={10}^{-15}$ m${}^2$ the compaction length is of order 1 to 10 km. A model resolves this physics only where its cells are finer than $\delta$; see Resolution.

compaction_length with the default bulk-viscosity law ($\zeta =\eta /\varphi$), without melt weakening, for $k_0={10}^{-9}$ m${}^2$ (default, solid) and ${10}^{-7}$ m${}^2$ (dashed). Since $\delta \propto \sqrt{k_0\eta }$ the dashed curve for ${10}^{19}$ Pa s coincides with the solid one for ${10}^{21}$ Pa s. Only curves above the dotted line are resolved by the 30.6 km cells of the Venus-sized annulus at Nr = 96. Script: docs/scripts/fig_two_phase_sections.jl.

Porosity waves

A column of uniform porosity in which a local excess of melt is placed does not simply let it rise. The higher permeability lets melt into the matrix above, which must dilate, while the matrix below compacts behind it, and the excess travels as a solitary wave that keeps its shape (Barcilon & Richter 1986). For $k\propto {\varphi }^3$ and constant $\zeta$, a wave of amplitude $A$ (peak porosity over background) moves at $2A+1$ times the background melt speed. This closed-form solution is the standard test of a magma-dynamics code (Simpson & Spiegelman 2011) and one of MILEVA's benchmarks below.

The finest solitary-wave run of benchmarks/two_phase_benchmarks.jl (amplitude $A=4$, 8 cells per compaction length, CFL 0.5), replotted from output/two_phase_benchmarks/solitary_profiles.csv; heights and porosity are scaled by the background compaction length ${\delta }_0$ and porosity ${\varphi }_0$, and the parts next to the closed walls are cut off. A: the wave rises without changing shape. B: its crest moves at $2A+1=9$ times the background melt speed. The comparison with the exact profile is under Verification. Script: docs/scripts/fig_two_phase_sections.jl.

Typical values

quantity	typical value	source
melt-solid density contrast, upper mantle	300 to 600 kg m${}^{-3}$	Stolper et al. (1981)
melt viscosity, basalt	1 to 10 Pa s	McKenzie (1984)
permeability law $k=d^2{\varphi }^n/C$, olivine-basalt	$n=2.6\pm 0.2$, $C={58}_{-22}^{+36}$ ($d=35$ μm assumed)	Miller et al. (2014)
permeability law $k=d^2{\varphi }^n/C$, quartzite with aqueous fluid	$n=3$, $C\approx 200$	Wark & Watson (1998)
melt weakening $\eta \propto \exp (-{\alpha }_{\varphi }\varphi )$	${\alpha }_{\varphi }\approx 25$ to 30	Mei et al. (2002); Hirth & Kohlstedt (2003)
bulk viscosity	$\zeta \sim \eta /\varphi$	Simpson, Spiegelman & Weinstein (2010)
compaction length, $\eta ={10}^{19}$ to ${10}^{21}$ Pa s, $\varphi =0.01$	1 to 10 km	this page, $k_0={10}^{-9}$ m${}^2$
melt fraction retained in the source	about 0.1 to 2 %	McKenzie (1984); Kohlstedt & Holtzman (2009)

Symbols used on this page

symbol	meaning	unit
$\varphi$	porosity, the volume fraction of melt	-
${\mathbf{v}}_s$, ${\mathbf{v}}_f$	solid and melt velocity	m s${}^{-1}$
$\mathbf{q}$	Darcy flux $\varphi ({\mathbf{v}}_f-{\mathbf{v}}_s)$	m s${}^{-1}$
$p_f$, $p_c$	melt (fluid) pressure and compaction pressure	Pa
$\tilde{p}$	the solved fluid pressure, $p_f$ minus the reference hydrostat absorbed by the buoyancy	Pa
${\eta }_s$, $\eta$	shear viscosity of the dry matrix (from the rheology) and after melt weakening	Pa s
$\zeta$	bulk (compaction) viscosity	Pa s
$k$, $k_0$, $n$	permeability, its prefactor and exponent, $k=k_0{\varphi }^n$	m${}^2$, m${}^2$, -
$\mu$	melt viscosity	Pa s
$K_D$	Darcy mobility $k/\mu$	m${}^2$ Pa${}^{-1}$ s${}^{-1}$
$\delta$	compaction length	m
$\mathrm{\Delta }\rho$	melt-solid density contrast ${\rho }_s-{\rho }_f$	kg m${}^{-3}$
$\mathrm{\Delta }{\rho }_{\text{seg}}$	segregation density, the density difference that drives the Darcy flux	kg m${}^{-3}$
${\alpha }_{\varphi }$	melt-weakening exponent	-
${\zeta }_0$, $m_{\zeta }$	bulk-viscosity factor and exponent, $\zeta ={\zeta }_0\eta {\varphi }^{-m_{\zeta }}$	-, -
$\mathrm{\Gamma }$	melting rate (negative for freezing)	kg m${}^{-3}$ s${}^{-1}$
$\mathit{\tau }$, $\dot{\mathit{\epsilon }}$	deviatoric stress and strain rate of the matrix	Pa, s${}^{-1}$
$\mathbf{b}$	solid buoyancy from compute_buoyancy, a lateral density anomaly times $g$	N m${}^{-3}$
${\rho }_0$, $\alpha$, $c_p$, $L$, $g$	reference density, thermal expansivity, specific heat, latent heat, gravity
$F$, $C_{\text{dep}}$, $F_{\text{max}}$, $\mathrm{\Delta }{T}_{\text{dep}}$	melt fraction, depletion, fusible fraction and solidus rise of Melting
$D$, $G$	discrete divergence (cells from faces) and gradient (faces from cells)
$V_{ij}$	cell volume $(r_{F,i+1}^3-r_{F,i}^3)\mathrm{\Delta }\theta /3$, the measure of the marker weights and the energy diagnostics	m${}^3$ rad${}^{-1}$

Equations

MILEVA solves the three-field form of the McKenzie (1984) equations (solid velocity, fluid pressure, compaction pressure) in the Boussinesq form used for magma-dynamics benchmarks (Keller, May & Kaus 2013). Both phases are incompressible and share one reference density ${\rho }_0$ in every mass equation; density differences enter only the body forces. Both phases share one temperature (local thermal equilibrium).

Material laws

$$\eta =\mathrm{clamp}({\eta }_s{e}^{-{\alpha }_{\varphi }\varphi },{\eta }_{\text{min}},{\eta }_{\text{max}}),\frac{1}{\zeta }=\frac{{\varphi }^{m_{\zeta }}}{{\zeta }_0\eta },K_D=\frac{k_0{\varphi }^n}{\mu }.$$

The bulk viscosity $\zeta$ itself is never formed. The compaction relation is written with $1/\zeta$, which goes smoothly to zero as $\varphi \to 0$, so a dry cell is the exact limit of a wet one rather than a division by zero.

The material laws with the default parameters. A: melt_weaken! (${\alpha }_{\varphi }$ = 27) and the bulk viscosity from inv_bulk_viscosity (${\zeta }_0$ = 1, $m_{\zeta }$ = 1), both relative to the dry viscosity; $\zeta$ grows without bound as $\varphi \to 0$ while $1/\zeta$ goes to zero. B: darcy_mobility for $n$ = 2 and the default 3 with $k_0={10}^{-9}$ m${}^2$ and $\mu$ = 1 Pa s; below cfg.twoPhasePorosityThreshold a cell is dry. Script: docs/scripts/fig_two_phase_sections.jl.

Field equations

The total stress of the mixture is $-(p_f+p_c)\mathbf{I}+\mathit{\tau }$, with the trace-free stress of a compressible matrix $\mathit{\tau }=2\eta (\dot{\mathit{\epsilon }}({\mathbf{v}}_s)-\frac{1}{3}(\mathrm{\nabla }\cdot {\mathbf{v}}_s)\mathbf{I})$. The mixture pressure $p_f+p_c$ is the pressure $P$ of the single-phase solver. Then

$$\begin{array}{rlrl}\mathrm{\nabla }\cdot \mathit{\tau }-\mathrm{\nabla }\tilde{p}-\mathrm{\nabla }{p}_c+\mathbf{b}+\varphi \mathrm{\Delta }\rho g{\mathbf{e}}_r& =0& & \text{momentum of the mixture}\\ \mathbf{q}& =-K_D(\mathrm{\nabla }\tilde{p}-\mathrm{\Delta }{\rho }_{\text{seg}}g{\mathbf{e}}_r)& & \text{Darcy}\\ \mathrm{\nabla }\cdot {\mathbf{v}}_s+\mathrm{\nabla }\cdot \mathbf{q}& =0& & \text{total volume}\\ \mathrm{\nabla }\cdot {\mathbf{v}}_s+p_c/\zeta & =0& & \text{compaction}\\ \frac{\partial \varphi }{\partial t}+\mathrm{\nabla }\cdot (\varphi {\mathbf{v}}_s)& =\mathrm{\nabla }\cdot {\mathbf{v}}_s+\frac{\mathrm{\Gamma }}{{\rho }_0}& & \text{porosity}\end{array}$$

with ${\mathbf{e}}_r$ the outward unit vector. The porosity equation is the solid mass balance rewritten for $\varphi$, equivalently ${\mathrm{D}}_s\varphi /\mathrm{D}t=(1-\varphi )\mathrm{\nabla }\cdot {\mathbf{v}}_s+\mathrm{\Gamma }/{\rho }_0$, so melt accumulates where the matrix dilates and is squeezed out where it compacts. The total-volume equation carries no melting source $\mathrm{\Gamma }(1/{\rho }_f-1/{\rho }_s)$. The annulus is closed, so the integral of $\mathrm{\nabla }\cdot ({\mathbf{v}}_s+\mathbf{q})$ over it vanishes identically, and a volume source would make the equations inconsistent with the pressure gauge. Dropping it is the Boussinesq mass approximation of the benchmark literature.

A laterally uniform column of the Venus-sized annulus (Nr = 96, Nt = 8, uniform 30.6 km cells, Boussinesq, MUMPS) with ${\eta }_s={10}^{21}$ Pa s, $k_0={10}^{-7}$ m${}^2$ and a band of 1 per cent porosity at 1500 km on a 0.2 per cent background, solved once by solve_stokes_darcy!; an operator illustration, not a run. A: the divergence of the solid velocity and minus the divergence of the Darcy flux, per cell. B: the same divergence and $-p_c/\zeta$. Both equations hold to $6\times {10}^{-15}$ of the largest divergence. Script: docs/scripts/fig_two_phase_sections.jl.

Buoyancy and the segregation density

MILEVA's solid buoyancy $\mathbf{b}$ from compute_buoyancy is a lateral anomaly, in which the radial mean of each density term is removed and absorbed into the solved pressure. In single-phase flow that is harmless, because a purely radial body force is a pressure gradient. In two-phase flow the same pressure also drives the melt, so the absorbed hydrostat must be put back into the Darcy law. Writing $\tilde{p}=p_f-P_{\text{ref}}(r)$ with $\mathrm{d}{P}_{\text{ref}}/\mathrm{d}r=-{\rho }_{\text{ref}}(r)g$ gives the exact Darcy law above, with the segregation density

$$\mathrm{\Delta }{\rho }_{\text{seg}}(r,\theta )={\rho }_{\text{ref}}(r)-{\rho }_f(T)=\mathrm{\Delta }\rho +{\rho }_0\alpha (T-⟨T{⟩}_{\theta })+⟨\delta {\rho }_{\text{phase}}+\delta {\rho }_{\text{pile}}{⟩}_{\theta },$$

where $⟨\cdot {⟩}_{\theta }$ is the azimuthal mean of a row of cells. The last term restores the radial means of the phase and pile density terms, the only ones compute_buoyancy subtracts (it records them through its mean_profile keyword). Melt and solid share $\alpha$ at first order, so the thermal reference cancels except for its lateral anomaly. For the same reason the melt buoyancy $\varphi \mathrm{\Delta }\rho g$ enters the momentum equation without mean subtraction. Products of the porosity with the thermal and compositional anomalies, $O(\varphi \alpha \mathrm{\Delta }T)$, are dropped.

segregation_density on the Venus-sized annulus (Nr = 96, Nt = 384) for a temperature rising 0.3 K km${}^{-1}$ with depth plus a +200 K plume of the full depth, with the phase-density means recorded by compute_buoyancy through mean_profile (cfg.useClaperonBuoyancy = true). A: at 414 km the plume raises $\mathrm{\Delta }{\rho }_{\text{seg}}$ by its thermal anomaly ${\rho }_0\alpha (T-⟨T{⟩}_{\theta })$, and the ambient mantle sits slightly below the mean. B: with depth the restored radial means of the phase terms add to the 500 kg m${}^{-3}$ contrast, which is what the Darcy law sees. Script: docs/scripts/fig_two_phase_sections.jl.

Energy

The mixture energy equation gains the heat carried by the melt relative to the solid,

$${\rho }_0{c}_p(\frac{\partial T}{\partial t}+({\mathbf{v}}_s+\mathbf{q})\cdot \mathrm{\nabla }T)=\mathrm{\nabla }\cdot (k_T\mathrm{\nabla }T)-L\mathrm{\Gamma }+H,$$

so the extra source is $-{\rho }_0{c}_p\mathbf{q}\cdot \mathrm{\nabla }T$ (darcy_heat_source). It is added to the diffusion solve beside the other volumetric sources, and the marker hand-back delivers it like every other source (see Temperature transport and time stepping). The latent heat of melting and freezing goes through the heat-pipe channel H_magma of Melting.

A laterally uniform column of the Venus-sized annulus (Nr = 96, Nt = 8, uniform 30.6 km cells, Boussinesq, MUMPS) with ${\eta }_s={10}^{21}$ Pa s, $k_0={10}^{-7}$ m${}^2$ and a band of 1 per cent porosity at 1500 km on a 0.2 per cent background, solved once by solve_stokes_darcy!; an operator illustration, not a run, with a linear temperature from 740 K at the surface to 4025 K at the CMB. A: the Darcy flux. B: darcy_heat_source (line) and $-{\rho }_0{c}_p{q}_r\mathrm{d}T/\mathrm{d}r$ with the cell means of $q_r$ (dots); melt rising through a temperature that falls upward warms the rock it passes. Script: docs/scripts/fig_two_phase_sections.jl.

How MILEVA discretises it

Unknowns and operators

The two-phase system lives on the staggered mesh of the Stokes solver: $v_r$ on radial faces ($(N_r+1)\times N_{\theta }$), $v_{\theta }$ on azimuthal faces, and $\tilde{p}$, $p_c$ and $\varphi$ at cell centres. The unknown vector is

$$x=[v_r\mid v_{\theta }\mid \tilde{p}\mid p_c],$$

so the first three blocks are index-identical to the single-phase system and _assemble_stokes_core (solvers/stokes_assembly.jl) supplies them, with the stress made trace free (trace_free = true). The continuity divergence and the pressure gradient are the pair the single-phase operator already uses,

$$(Du{)}_{ij}=\frac{r_{F,i+1}^2{u}_{r,i+1,j}-r_{F,i}^2{u}_{r,i,j}}{r_{C,i}^2\mathrm{\Delta }{r}_i}+\frac{u_{\theta ,i,j+1}-u_{\theta ,i,j}}{r_{C,i}\mathrm{\Delta }\theta },(G_{r}p{)}_{i,j}=\frac{p_{i,j}-p_{i-1,j}}{\mathrm{\Delta }{r}_{F,i}},(G_{\theta }p{)}_{i,j}=\frac{p_{i,j}-p_{i,j-1}}{r_{C,i}\mathrm{\Delta }\theta },$$

with $r_F$, $r_C$ face and centre radii, $\mathrm{\Delta }{r}_i$ the cell width and $\mathrm{\Delta }{r}_{F,i}$ the centre spacing (stretched grids included). $G$ is minus the adjoint of $D$ in the volume-weighted inner product, so every flux difference below telescopes over the annulus and the Darcy operator $-D{K}_{D}G$ is symmetric positive semi-definite.

The trace-free stress is the compressible deviator that the anelastic (ALA/TALA) path also uses; there is one implementation, switched on by trace_free in _viscous_block_deviatoric. It adds $-\frac{2}{3}\eta \mathrm{tr}\dot{\mathit{\epsilon }}$ to each of the three normal stresses ${\tau }_{rr}$, ${\tau }_{\theta \theta }$ and ${\tau }_{\phi \phi }$, with the trace built from the same discrete strain rates as the stress itself, $\mathrm{tr}\dot{\mathit{\epsilon }}={\dot{\epsilon }}_{rr}+{\dot{\epsilon }}_{\theta \theta }+{\dot{\epsilon }}_{\phi \phi }$. The work of the stress is then the dissipation $2\eta (\dot{\mathit{\epsilon }}:\dot{\mathit{\epsilon }}-\frac{1}{3}{\mathrm{tr}}^2\dot{\mathit{\epsilon }})\ge 0$ cell by cell. This trace differs from the continuity divergence $Du$ at $O(h^2)$ on this grid, which matters only for the reduction to single phase below. With two-phase flow on, compute_strain_rate removes the same trace, so the strain rate the rheology reads is the deviatoric one.

Where the unknowns live, drawn at the true cell sizes of the top three rows of the Venus-sized annulus at Nr = 96, Nt = 384 ($\mathrm{\Delta }r$ = 30.6 km, $r_C\mathrm{\Delta }\theta$ = 99 km): pressures, porosity and temperature at the cell centres, the radial solid velocity and Darcy flux on the radial faces, their azimuthal components on the azimuthal faces. Script: docs/scripts/fig_two_phase_sections.jl.

The linear system

After the momentum rows are scaled as in single phase ($S_v$, the row scaling of the momentum equations), the system is

$$\left[\begin{array}{ccc}{A}^{\prime }& -S_{v}G& -S_{v}G\\ D& -D{K}_{D}G& 0\\ D& 0& Z\end{array}\right]\left[\begin{array}{c}u\\ \tilde{p}\\ p_c\end{array}\right]=\left[\begin{array}{c}-S_v(\mathbf{b}+\varphi \mathrm{\Delta }\rho g{\mathbf{e}}_r)\\ -D(K_D\mathrm{\Delta }{\rho }_{\text{seg}}g{\mathbf{e}}_r)\\ 0\end{array}\right],Z=\mathrm{diag}(\frac{{\eta }_{\text{ref}}}{{\zeta }_{ij}}),$$

assemble_stokes_darcy in solvers/two_phase_assembly.jl:

• Momentum. $A^{\prime }$ is the deviatoric operator with the trace-free stress; the compaction pressure enters with the same coefficients as $\tilde{p}$.

• Total volume. $Du-D{K}_{D}G\tilde{p}$, with the gravity term of the Darcy law on the right-hand side. $\mathrm{\Delta }{\rho }_{\text{seg}}$ is averaged onto each face.

• Compaction. $Du+({\eta }_{\text{ref}}/\zeta )p_c=0$ in wet cells ($\varphi >{\varphi }_{\text{th}}$ = cfg.twoPhasePorosityThreshold, default ${10}^{-5}$). In dry cells the row is the identity $p_c=0$, the exact $\varphi \to 0$ limit, which keeps the matrix non-singular.

• Scaling. The existing non-dimensionalisation (shell depth, diffusivity, median viscosity ${\eta }_{\text{ref}}$) is reused; then the same point pin and max-abs row equilibration as in single phase are applied.

The non-zero pattern of assemble_stokes_darcy on a 6 by 12 annulus with porosity in rows 3 and 4 (300 unknowns, 2340 non-zeros), blocks in the order $v_r$, $v_{\theta }$, $\tilde{p}$, $p_c$. The Darcy entries $-D{K}_{D}G$ appear in the $\tilde{p}$ diagonal block only where cells are wet, and the $p_c$ rows of dry cells are the identity $p_c=0$. A deliberately small grid, so the pattern can be read. Script: docs/scripts/fig_two_phase_sections.jl.

Face permeability

The Darcy flux lives where the velocity lives, $q_r$ on radial faces and $q_{\theta }$ on azimuthal faces. Its mobility on a face is chosen from the two cells that share it (face_mobility):

neighbours	face mobility	why
both wet	arithmetic mean of $K_D(\varphi )$	second order where the porosity is smooth
one wet, one dry	donor cell: $K_D$ of the upstream cell	melt can enter dry rock, a dry cell cannot lose melt it does not have
both dry	0
wall	0	no flux through the CMB or the surface

The flow direction at a wet/dry face is the one of the previous solve (upward before the first). solve_stokes_darcy! refreshes the directions from the new flux and solves again until no wet/dry face flips, at most 10 solves. A solve that reaches the cap keeps its last solution, in which some wet/dry faces carry a flux against the direction their mobility was taken from. It is counted in donor_capped, those faces are counted in donor_mismatch, both appear in the [2-phase] log line (capped, mismatch), and a warning names them. The mobilities the returned flux was computed with are kept with it, so the time step reads the melt speed that flux actually implies. Arithmetic averaging would let a dry cell below a melting zone drain into it and go to negative porosity on every step; harmonic averaging would forbid melt from ever entering a dry cell, so a melt front could not advance.

face_mobility on one column of eleven cells (default $k_0$ and $n$). A: the porosity. B: the face mobility with the donor direction set upward and downward, against the arithmetic and harmonic means. Inside the wet layer all choices agree; at the lower wet/dry face only downward flow carries a mobility, and at the upper one only upward flow, so melt can advance into dry rock and a dry cell cannot lose melt. Script: docs/scripts/fig_two_phase_sections.jl.

Boundary conditions and the pressure gauge

Free slip is unchanged: $v_r=0$ and zero shear traction on both walls, periodicity in $\theta$. The Darcy flux is zero on both walls, in the matrix and in the gravity right-hand side, so the buoyancy-driven part of $\mathbf{q}$ also vanishes there. A constant added to $\tilde{p}$ changes nothing, while a constant added to $p_c$ changes the compaction rows, so there is one pressure null mode, as in single phase. It is removed by one point pin on $\tilde{p}$, which replaces that cell's total-volume row. The pin is consistent because both $Du$ and $D{K}_{D}G\tilde{p}$ telescope to the walls, where $v_r$ and $q_r$ vanish, so the volume-weighted sum of the total-volume rows and of their right-hand side is identically zero. Rigid rotation keeps its pin on one $v_{\theta }$; it has zero strain rate, zero divergence and zero Darcy flux.

The column of the panels above with the band moved to 2700 km depth. A: $q_r$ and $-v_r$ over the whole depth. B: the faces closest to the CMB, where both vanish on the wall face; the volume-weighted sum of the total-volume rows over the annulus is zero to round-off, which is what makes the single pressure pin consistent. Script: docs/scripts/fig_two_phase_sections.jl.

Reduction to single-phase Stokes

With $\varphi =0$ everywhere every face is dry, so $D{K}_{D}G=0$ and its right-hand side vanishes, every compaction row is $p_c=0$ and the melt buoyancy is zero. The system differs from the single-phase one only by the trace term. In the continuum it is proportional to $\mathrm{\nabla }\cdot {\mathbf{v}}_s$, which vanishes, so the exact solutions coincide. Discretely the trace is $\mathrm{tr}\dot{\mathit{\epsilon }}$ while the total-volume rows hold $Du=0$, and the two differ at $O(h^2)$, so the computed solutions differ at that order: on the variable-viscosity test of the suite (contrast ${10}^6$) by $9\times {10}^{-9}$ in the velocity and $4.4\times {10}^{-5}$ in the pressure at Nr = 24, and $6\times {10}^{-10}$ and $1.05\times {10}^{-5}$ at Nr = 48; see Verification.

The variable-viscosity case of tests/test_two_phase.jl (contrast ${10}^6$, Nt = 4 Nr), solved by solve_stokes_mumps and by solve_stokes_darcy! at $\varphi =0$, extended to Nr = 96. The pressure difference falls at second order, $5.5\times {10}^{-5}$, $1.4\times {10}^{-5}$ and $3.4\times {10}^{-6}$ at Nr = 24, 48 and 96; the velocity difference is $4.5\times {10}^{-9}$ at Nr = 24 and falls faster. Script: docs/scripts/fig_two_phase_sections.jl.

Solving

Every two-phase solve factorises afresh with MUMPS (solverType = "mumps") or UMFPACK ("direct"), twice per time step (see Porosity transport); a factor is valid only for the same viscosity, porosity and donor directions, which practically never recur between steps. Each solve is followed by one step of iterative refinement with the same factor. Where $\varphi$ is small and the mobility low, the compaction rows carry a small ${\eta }_{\text{ref}}/\zeta$ and are weak. Without the refinement, UMFPACK left a backward error of ${10}^{-9}$ to $2\times {10}^{-8}$ and a total-volume residual of up to $9\times {10}^{-5}$ of $max|D{v}_s|$ there (MUMPS up to $7\times {10}^{-11}$); with it both reach a backward error of about ${10}^{-14}$. The solve sits inside the Picard loop of the governing equations; melt weakening is applied to the viscosity before the loop, and the reported divergence is the total-volume residual $max|D{v}_s+Dq|$.

The field of the solver-accuracy test in tests/test_two_phase.jl (Nr = 32, Nt = 96, variable viscosity and porosity) on the Venus-sized annulus, assembled by assemble_stokes_darcy and solved by UMFPACK alone and by _solve_two_phase_linear with its refinement step. The backward error falls from $1.5\times {10}^{-9}$, $1.0\times {10}^{-8}$ and $8\times {10}^{-13}$ to below ${10}^{-14}$ for $k_0={10}^{-12}$, ${10}^{-9}$ and ${10}^{-7}$ m${}^2$. Script: docs/scripts/fig_two_phase_sections.jl.

Porosity transport

Porosity lives on the cells and is not carried by markers. Write the right-hand side of the porosity equation for a solid velocity ${\mathbf{v}}_s$ as

$$R(\varphi ,{\mathbf{v}}_s)=-[\mathrm{\nabla }\cdot (\varphi {\mathbf{v}}_s){]}_{\text{TVD}}+D{v}_s,$$

with the van Leer face reconstruction of the 2-D temperature advection in conservative form (zero-gradient wall ghost) and the solver's own divergence $D{v}_s$. update_porosity! advances it with Heun's method, the two-stage strong-stability-preserving Runge-Kutta scheme (Shu & Osher 1988):

$$\begin{array}{rl}{\varphi }^{\ast }& ={\varphi }^n+\mathrm{\Delta }tR({\varphi }^n,{\mathbf{v}}_s^n),\\ {\varphi }^{n+1}& ={\varphi }^n+\frac{1}{2}\mathrm{\Delta }t[R({\varphi }^n,{\mathbf{v}}_s^n)+R({\varphi }^{\ast },{\mathbf{v}}_s^{\ast })],\end{array}$$

where ${\mathbf{v}}_s^n$ is the solve of the step at ${\varphi }^n$ and ${\mathbf{v}}_s^{\ast }$ a second two-phase solve at ${\varphi }^{\ast }$, with the melt weakening moved to ${\varphi }^{\ast }$ and the buoyancy, segregation density and temperature of the step, put in the same no-net-rotation frame. Each stage is a forward-Euler step under the same CFL limit, so the positivity bound of the van Leer step carries over, and the scheme is second order in time. After the step the solver state (fluid and compaction pressure, Darcy flux) is again the one at ${\varphi }^n$, which the temperature and marker transport of the same step use; the residuals and donor counters cover both solves.

Both terms of $R$ telescope, so $\sum _{ij}{V}_{ij}{\varphi }_{ij}$ changes only through melting and freezing, to round-off and independently of the solver tolerance. The porosity measure is the exact cell volume $V_{ij}=(r_{F,i+1}^3-r_{F,i}^3)\mathrm{\Delta }\theta /3$: the transport increment is rescaled by $r_{C,i}^2\mathrm{\Delta }{r}_i\mathrm{\Delta }\theta /V_{ij}$, so the divergence telescopes in the same measure the marker weights carry. Negative porosity is never clipped. Values below $-{10}^{-12}$ are counted and their volume is reported, and every consumer treats $\varphi \le 0$ as dry.

The column of the panels above with a band of 1 per cent porosity at 2000 km on a dry background, advanced by 120 steps of update_porosity! at the step of two_phase_dt, with a fresh solve_stokes_darcy! before each. A: the band rises into the dry rock above it and steepens at its front, leaving a thin tail. B: $\sum V\varphi$ changes by less than $2\times {10}^{-15}$ of its value, and no cell goes below $-{10}^{-12}$. Script: docs/scripts/fig_two_phase_sections.jl.

Time step

On top of the mantle CFL limit, two_phase_dt imposes, with $C_{\varphi }$ = cfg.twoPhaseCFL (0.5, the van Leer positivity bound),

$$\mathrm{\Delta }t\le C_{\varphi }\frac{\mathrm{\Delta }x}{|v_s|+|q|/{\varphi }_{\text{face}}},\mathrm{\Delta }t\le C_{\varphi }\frac{\mathrm{\Delta }x}{|q|},\mathrm{\Delta }t\le \frac{0.5\varphi }{(1-\varphi )|D{v}_s|}(D{v}_s<0),$$

the melt CFL condition on every face that carries Darcy flux, the CFL condition of the heat carried by the melt, and a compaction limit that keeps a compacting cell from losing more than half its melt in one step. ${\varphi }_{\text{face}}$ is the porosity whose mobility the face used. The cost is set mainly by $k_0$. On Venus ($\mathrm{\Delta }\rho =500$ kg m${}^{-3}$, $g=8.87$ m s${}^{-2}$, $\mu =1$ Pa s) the melt speed $k_0{\varphi }^2\mathrm{\Delta }\rho g/\mu$ at $\varphi =0.01$ is 1.4 cm yr${}^{-1}$ with $k_0={10}^{-9}$ m${}^2$ and 1.4 m yr${}^{-1}$ with $k_0={10}^{-7}$ m${}^2$, a melt CFL step on 10 km cells of about 0.36 Myr and 3.6 kyr.

two_phase_dt after one solve_stokes_darcy! on the column of the panels above with a dry background, for peak porosities of 0.25 to 4 per cent and three values of $k_0$. The limit falls roughly as the inverse square of the porosity and in proportion to $1/k_0$ at small porosity, where the melt CFL condition sets it. Script: docs/scripts/fig_two_phase_sections.jl.

Melting and freezing

With cfg.trackBasalt on, melt_two_phase! (physics/two_phase.jl) replaces basalt_evolution! and returns the same quantities, so the time loop uses either. It melts and freezes to equilibrium every step, so the rates do not depend on the step length:

1. Melting, on markers shallower than cfg.meltExtractMaxDepth_km. The increment is the melt fraction that leaves the parcel on its depletion-raised solidus after paying its latent heat (Katz 2008),

$$\mathrm{\Delta }F=\frac{T-T_{\text{sol}}(P,C_{\text{dep}})}{\mathrm{\Delta }{T}_{\text{dep}}/F_{\text{max}}+L/c_p},$$

capped at the remaining fusible fraction $F_{\text{max}}-C_{\text{dep}}$, and $C_{\text{dep}}$ grows by $\mathrm{\Delta }F$. The solidus and its depletion rise are those of Melting.

2. Melt into the porosity. $\mathrm{\Delta }\varphi =\sum _m{w}_m\mathrm{\Delta }{F}_m(1-\varphi )/V$ over the cell's markers, with $w_m$ the marker weights.

3. Freezing, on cells that did not melt in the same call, hold melt and lie below their solidus (with the cell-mean depletion), again to equilibrium. The amount frozen is the volume whose latent heat brings the cell back to its solidus, $\mathrm{\Delta }{\varphi }_f=min(\varphi ,c_p(T_{\text{sol}}-T)/L)$. The fraction cfg.eruptionEfficiency of it erupts to the top populated cell of the column (stamping the surface age); the rest becomes basalt on the cell's own markers, with overflow through the cascade and per-column backlog of Melting. Intrusions therefore form where transported melt meets cold rock, not in a prescribed depth band.

4. Heat, returned as the temperature tendency H_magma. Melting removes ${\rho }_{0}L$ per unit melt volume produced, freezing in place releases it, and the erupted melt removes its sensible heat ${\rho }_0{c}_p(T-T_{\text{surf}})$ from its source (its latent heat was paid when it melted).

The ledger reports melt produced, erupted and frozen in place through the existing melt columns of the diagnostics. Without trackBasalt nothing melts and the porosity evolves from cfg.initialPorosity alone. Heat-producing elements stay on the markers; partitioning them into the melt is refused rather than dropped.

melt_two_phase! on the Venus-sized annulus (Nr = 96, Nt = 384) for eleven fertile melting cells at 150 km and nine lid cells at 107 km holding 1.5 per cent melt, each the only active cell of its column. A: the equilibrium increment follows $(T-T_{\text{sol}})/(\mathrm{\Delta }{T}_{\text{dep}}/F_{\text{max}}+L/c_p)$, one per cent per 8.2 K, up to the fusible cap; the 0.1 per cent per kelvin marker law of Melting is shown for comparison. B: freezing removes $c_p(T_{\text{sol}}-T)/L$, so 1.5 per cent of melt is used up by 4.8 K of undercooling. Script: docs/scripts/fig_two_phase_sections.jl.

How it is computed

One time step with cfg.twoPhaseFlow = true on a Venus-sized annulus (Nr = 96, Nt = 384, uniform 30.6 km cells, Boussinesq, MUMPS) with a prescribed viscosity of ${10}^{19}$ Pa s below a 125 km lid of ${10}^{23}$ Pa s. Right: the value each stage produces for a fertile, still dry cell at 138 km, 20 K above its solidus, and a lid cell at 107 km holding 1.5 % melt, 3 K below its solidus. The melting cell reaches 2.4 % porosity in one call and cools by 7.7 K; the freezing cell freezes 0.94 % of its volume and warms by 0.5 K. The compaction length is 1.8 km (0.06 cells) in the plume head and 52 km (1.7 cells) in the lid, so compaction is resolved only in the lid. The solved Darcy flux through the top face of the melting cell is 0.087 cm/yr against 0.198 cm/yr for the estimate $K_D\mathrm{\Delta }{\rho }_{\text{seg}}g$, which ignores the pressure the matrix builds as it opens above a single wet cell. In 20 kyr the cell passes 0.03 % of its volume to the cell above; a forward-Euler step passed 0.06 %, because it moved the melt into the dry lid cell at the rate of the dry cell, which carries no compaction resistance, while Heun's second solve sees that cell already wet and resisting dilation. The cell above, being the lid base 518 K colder, is heated by the arriving melt at $1.9\times {10}^{-6}$ W m${}^{-3}$. Every stage value is recomputed from the formula shown and checked against the code, the melt volume changes by ${10}^{-16}$ (relative) in the transport step, and the solver's total-volume residual is $4\times {10}^{-16}$ of $max|D{v}_s|$. Script: docs/scripts/fig_two_phase_pipeline.jl.

Verification

Three benchmarks with closed-form or exact discrete references are in benchmarks/two_phase_benchmarks.jl (CSVs in output/two_phase_benchmarks/, figures by docs/scripts/fig_two_phase.jl); tests/test_two_phase.jl guards them on small grids together with the properties below.

Compaction boundary layers, exact spherical solution

A laterally uniform column with constant ${\varphi }_0=0.01$, $\eta$, $\zeta$ and $K_D$, no thermal or compositional buoyancy and closed walls. With radial flow only, the annulus equations reduce exactly to

$${\delta }^2\frac{1}{r^2}\frac{\mathrm{d}}{\mathrm{d}r}(r^2\frac{\mathrm{d}{p}_c}{\mathrm{d}r})-p_c=-\frac{A}{r},p_c=\frac{A}{r}+\frac{B_1{e}^{-(r-R_{\text{in}})/\delta }+B_2{e}^{-(R_{\text{out}}-r)/\delta }}{r},$$

with $A=2\zeta K_D(1-{\varphi }_0)\mathrm{\Delta }\rho g$ and $B_1$, $B_2$ fixed by $q=0$ on both walls. The $A/r$ term is the spreading of a uniform upward melt flux through shells of area $r^2$, absent in a Cartesian column, so the test checks the spherical metric as well as the operator. Melt rises at its Darcy speed through the interior, piles up in a compaction layer of thickness $\delta$ below the surface and leaves one above the CMB.

Left: compaction pressure and Darcy flux at $\delta /D=0.05$, Nr = 64, against the exact solution. Right: maximum relative error against Nr for $\delta /D$ = 0.02, 0.05 and 0.2, uniform (solid) and stretched (dashed) radial grids.

The error falls at second order (observed orders 2.0 to 2.3 in $p_c$ and 1.7 to 2.0 in $q$ between successive grids) on both grids. At Nr = 256 on the uniform grid it is $6.6\times {10}^{-4}$, $1.0\times {10}^{-4}$ and $7.7\times {10}^{-6}$ in $p_c$ and $5.9\times {10}^{-4}$, $9.6\times {10}^{-5}$ and $8.4\times {10}^{-6}$ in $q$ for $\delta /D$ = 0.02, 0.05 and 0.2. The solid velocity matches $-q$ to the same error, the lateral velocity is ${10}^{-12}$ of the radial one or less, and the total-volume residual is at most $8\times {10}^{-12}$ of $max|D{v}_s|$.

Magmatic solitary wave

A wave of amplitude $A=4$ on a background ${\varphi }_0={10}^{-4}$, $n=3$, constant $\zeta$ ($m_{\zeta }=0$), no melt weakening, travelling about 100 compaction lengths up a radial column of height 250 ${\delta }_0$. The analytic profile of Barcilon & Richter (1986) in the form of Simpson & Spiegelman (2011) is first checked against the travelling-wave ODE it must satisfy (largest relative mismatch $3.5\times {10}^{-10}$), so a transcription error cannot pass. No exact solitary wave exists in a spherical column, because shells grow as $r^2$, so the column is placed at $R_{\text{in}}/{\delta }_0={10}^3$, ${10}^4$ and ${10}^5$, where the geometric departure falls as $1/R_{\text{in}}$.

Top: the wave at five times at $R_{\text{in}}={10}^5{\delta }_0$, $h={\delta }_0/8$, CFL 0.5, against the analytic profile. Bottom, left to right: crest position error against time; errors at $t=11$ and $R_{\text{in}}={10}^5{\delta }_0$ against the grid spacing at CFL 0.5 (filled) and Richardson-extrapolated to $\mathrm{\Delta }t\to 0$ from CFL 0.5 and 0.25 (open); errors at CFL 0.5, $h={\delta }_0/8$, against $R_{\text{in}}/{\delta }_0$.

The porosity step was forward Euler until 2026-10; it is now Heun's method (Porosity transport). Both at $R_{\text{in}}/{\delta }_0={10}^5$ and CFL 0.5, the time loop's Courant number:

$R_{\text{in}}/{\delta }_0={10}^5$, CFL 0.5	$h={\delta }_0/2$	$h={\delta }_0/4$	$h={\delta }_0/8$
amplitude error, forward Euler	+13.7 %	+6.7 %	+3.3 %
amplitude error, Heun	−0.31 %	−0.17 %	−0.13 %
speed error, forward Euler	+1.7 %	+0.71 %	+0.28 %
speed error, Heun	+0.25 %	−0.032 %	−0.10 %
shape error (L2), forward Euler	0.22	0.11	0.052
shape error (L2), Heun	0.0040	0.0020	0.0019

Forward Euler grew the crest by 3.3 % at $h={\delta }_0/8$: its first-order time error dominated and fell only as $\mathrm{\Delta }t\propto h$. With Heun the crest amplitude errs by −0.13 % and the speed by −0.10 %, and the shape error is 27 times smaller on the finest grid and 56 times on the coarsest. The time error is now negligible. Halving and quartering the step at $h={\delta }_0/4$ moves the amplitude error from −0.165 % to −0.175 % and −0.178 %, and the successive differences give an observed time order of 2.00 in the final profile (2.11 in the amplitude, 1.98 in the speed), as the scheme should. Extrapolated to $\mathrm{\Delta }t\to 0$ the $h={\delta }_0/8$ errors are −0.13 % (amplitude), −0.10 % (speed) and 0.0019 (shape), the same as at CFL 0.5. On the two finer grids the errors no longer fall with $h$: they sit on the geometric departure of the spherical column from the Cartesian wave, which at $h={\delta }_0/8$ and CFL 0.5 is −8.3 %, −1.0 % and −0.13 % in the amplitude and −6.9 %, −0.82 % and −0.10 % in the speed at $R_{\text{in}}/{\delta }_0={10}^3$, ${10}^4$ and ${10}^5$, a factor of about 8 per decade of radius. Melt volume is conserved to $2\times {10}^{-14}$ in every run.

Melt volume, positivity and infiltration

A Gaussian melt blob (peak 2 %, cut to zero below ${10}^{-4}$) in dry rock on Nr = 64 × Nt = 256, $k_0={10}^{-7}$ m${}^2$, $\eta ={10}^{21}$ Pa s, run for 200 steps (6.7 Myr) with the two-phase time-step limits, two solves per step.

Left: relative change of the melt volume $\sum V\varphi$ per step in the blob and in two solitary-wave runs. Right: the solver's total-volume residual in the blob run, relative to $max|D{v}_s|$.

The melt volume drifts by at most $7.3\times {10}^{-15}$ (relative), the lowest porosity reached is $-1.4\times {10}^{-14}$ (round-off; no cell is counted negative), 140 cells that were dry at the start become wet as melt infiltrates, the total-volume residual stays below $9\times {10}^{-12}$ of $max|D{v}_s|$, and no solve reaches the donor-iteration cap.

Further properties guarded by the test suite

property	measured
$\varphi =0$ reproduces single-phase Stokes at viscosity contrast ${10}^6$ (MUMPS and UMFPACK)	$v_r$, $v_{\theta }$ differ by $9\times {10}^{-9}$ (Nr = 24) and $6\times {10}^{-10}$ (Nr = 48) relative, $P$ by $4.4\times {10}^{-5}$ and $1.05\times {10}^{-5}$, the $O(h^2)$ of the trace; $p_c$ and $\mathbf{q}$ exactly zero
work equals dissipation: $\int {\mathbf{b}}_{\text{tot}}\cdot {\mathbf{v}}_s+\int \mathrm{\Delta }{\rho }_{\text{seg}}g{q}_r=\int \mathit{\tau }:\dot{\mathit{\epsilon }}+\int \zeta (\mathrm{\nabla }\cdot {\mathbf{v}}_s{)}^2+\int |\mathbf{q}{|}^2/K_D$	ratio 1.0041 and 1.0066 at Nr = 32, 1.0011 and 1.0024 at Nr = 64, for contrasts ${10}^3$ and ${10}^6$
solver accuracy on weak compaction rows ($k_0={10}^{-12}$, $\varphi \approx 0.01$)	backward error $1\times {10}^{-14}$, total-volume residual $2\times {10}^{-10}$ of $max|D{v}_s|$, both backends
melt_two_phase! ledger	porosity, basalt, depletion and heat close against the ledger
time loop with two-phase flow, checkpoint and resume	completes; melt volume conserved to $6\times {10}^{-16}$ without melting; the porosity after a resume is identical to a continuous run
flag off	the existing suite passes, and the single-phase matrix and right-hand side (Boussinesq, ALA and TALA), temperature advection, buoyancy, strain rate and a 20-step time loop hash identical to the code without two-phase flow

Resolution

Compaction-scale structure, the boundary layers, porosity waves and channels, is represented only where the compaction length spans several cells. The log line [2-phase] prints the largest compaction length over wet cells in units of the finest radial cell (delta_max/dr_min). With the default $k_0={10}^{-9}$ m${}^2$ and interior viscosities of ${10}^{19}$ to ${10}^{21}$ Pa s, $\delta$ is 1 to 10 km at $\varphi =0.01$, while production grids have radial cells of 10 to 45 km. Such a run is in the zero-compaction-length limit, where melt rises at its buoyancy-driven Darcy speed. Fluxes, porosity and the melt budget stay consistent, but channels and waves are not resolved.

Replotted from the benchmark CSVs of benchmarks/two_phase_benchmarks.jl (compaction_convergence.csv, solitary_runs.csv, uniform grids). A: the maximum error of $p_c$ against the exact spherical solution depends on $\delta /h$ alone, whatever $\delta /D$; it is 2 per cent at one cell per compaction length, ${10}^{-4}$ at about 13 cells and falls at second order. B: the solitary wave at 2, 4 and 8 cells per ${\delta }_0$; on the large annulus ($R_{\text{in}}/{\delta }_0={10}^5$) the speed and amplitude errors are already below $3\times {10}^{-3}$ at 2 cells, while on the small one ($R_{\text{in}}/{\delta }_0={10}^3$) the 7 to 8 per cent departure does not change with resolution, so it is the curvature of that shell rather than the grid. Script: docs/scripts/fig_two_phase_sections.jl.

Configuration

field	default	unit	meaning
cfg.twoPhaseFlow	false		master switch
cfg.meltDensityContrast	500	kg m${}^{-3}$	$\mathrm{\Delta }\rho ={\rho }_s-{\rho }_f$
cfg.meltViscosity	1	Pa s	$\mu$
cfg.permeabilityK0	${10}^{-9}$	m${}^2$	$k_0$ in $k=k_0{\varphi }^n$; sets the cost
cfg.permeabilityExponent	3		$n$
cfg.meltWeakeningAlpha	27		${\alpha }_{\varphi }$
cfg.bulkViscosityFactor	1		${\zeta }_0$
cfg.bulkViscosityExponent	1		$m_{\zeta }$; with the defaults $\zeta =\eta /\varphi$
cfg.twoPhasePorosityThreshold	${10}^{-5}$		${\varphi }_{\text{th}}$, cells at or below carry $p_c=0$
cfg.twoPhaseCFL	0.5		$C_{\varphi }$
cfg.initialPorosity	0		uniform porosity at a fresh start

validate! refuses twoPhaseFlow with a 3-D grid, a geometry other than "spherical_annulus", a solver other than "direct" or "mumps", a compressibility mode other than "Boussinesq", radioHeatPartitionD > 0, meltSegregationMovesMarkers, useStickyAir, useDeviatoricStress = false, useVelocityLimit = true (the limiter rescales the solid velocity but not the Darcy flux or the compaction pressure), pPen other than 0 (the penalty would sit on the total-volume rows), useAntisqueeze (the matrix compacts physically, so the squeeze trigger would fire on every compacting cell), and, with trackBasalt, useHeatPipe = false or useMagmaOceanFreezing = true. Each refusal names its reason. meltExtractionTime_yr and intrusionDepth_km are not read in two-phase mode.

Output. VTK files gain Porosity, FluidPressure, CompactionPressure and DarcyFlux; diagnostics.csv gains melt_volume_m3rad, phi_max, neg_phi_cells, neg_phi_volume_m3rad, q_max_cmyr, mass_resid, compaction_resid and donor_capped, only when the flag is on, so the CSV of a single-phase run is unchanged. Checkpoints carry the porosity and the donor directions; a resume with twoPhaseFlow from a checkpoint without them is refused.

Assumptions and simplifications

• Boussinesq, one temperature. Both phases are incompressible with one reference density, the melt-solid density contrast is constant, so melt is buoyant wherever it exists (melt is produced only above meltExtractMaxDepth_km), and melt and matrix are in thermal equilibrium.

• The rest of the step stays first order. The porosity is advanced to second order in time, but the temperature, the markers and the melting coupling use the solid velocity and Darcy flux of the solve at ${\varphi }^n$, and the second solve reuses the viscosity of the step with only the melt weakening moved to ${\varphi }^{\ast }$.

• Heat-producing elements stay on the markers, and the melt carries no trace elements.

• Solid markers keep fixed weights while the matrix compacts, an error of second order in the porosity.

What using the solve at ${\varphi }^n$ for temperature, markers and melting misses. The column of the discretisation figures (band of 1 per cent porosity at 2000 km on a dry background, ${\eta }_s={10}^{21}$ Pa s, $k_0={10}^{-7}$ m${}^2$) advanced by one update_porosity! step of a given fraction of two_phase_dt and solved again: the Darcy flux changes by 17 per cent of its maximum over a full step and by 0.6 per cent over 1/32 of it, in proportion to the step, so the lag is first order and bounded by the two-phase time-step limit. Script: docs/scripts/fig_two_phase_sections.jl.

Scope

The implementation is 2-D, spherical annulus only, and is verified against 1-D solutions (compaction layers and solitary waves) and 2-D conservation properties. A 2-D porosity-wave or melt-channel benchmark is not part of it.

Tests

• tests/test_two_phase.jl: off by default and refusals; zero porosity against single-phase Stokes with both solvers, converging at $O(h^2)$; solver accuracy on weak compaction rows; the extended work-dissipation identity; compaction boundary layers on stretched grids; porosity conservation, positivity and infiltration; the melt_two_phase! ledger; the time loop, checkpoint and resume.

• benchmarks/two_phase_benchmarks.jl: the full compaction, solitary-wave and conservation series behind the figures (about three hours).

Further reading

• Barcilon, V. & Richter, F. M. (1986). Nonlinear waves in compacting media. Journal of Fluid Mechanics 164, 429-448. doi:10.1017/S0022112086002628

• Hernlund, J. W. & Tackley, P. J. (2008). Modeling mantle convection in the spherical annulus. PEPI 171, 48-54. doi:10.1016/j.pepi.2008.07.037

• Hirth, G. & Kohlstedt, D. L. (2003). Rheology of the upper mantle and the mantle wedge: a view from the experimentalists. In Inside the Subduction Factory, AGU Geophysical Monograph 138, 83-105. doi:10.1029/138GM06

• Katz, R. F. (2008). Magma dynamics with the enthalpy method: benchmark solutions and magmatic focusing at mid-ocean ridges. Journal of Petrology 49, 2099-2121. doi:10.1093/petrology/egn058

• Keller, T., May, D. A. & Kaus, B. J. P. (2013). Numerical modelling of magma dynamics coupled to tectonic deformation of lithosphere and crust. Geophysical Journal International 195, 1406-1442. doi:10.1093/gji/ggt306

• Kohlstedt, D. L. & Holtzman, B. K. (2009). Shearing melt out of the Earth: an experimentalist's perspective on the influence of deformation on melt extraction. Annual Review of Earth and Planetary Sciences 37, 561-593. doi:10.1146/annurev.earth.031208.100104

• McKenzie, D. (1984). The generation and compaction of partially molten rock. Journal of Petrology 25, 713-765. doi:10.1093/petrology/25.3.713

• Mei, S., Bai, W., Hiraga, T. & Kohlstedt, D. L. (2002). Influence of melt on the creep behavior of olivine-basalt aggregates under hydrous conditions. EPSL 201, 491-507. doi:10.1016/S0012-821X(02)00745-8

• Miller, K. J., Zhu, W., Montési, L. G. J. & Gaetani, G. A. (2014). Experimental quantification of permeability of partially molten mantle rock. EPSL 388, 273-282. doi:10.1016/j.epsl.2013.12.003

• Shu, C.-W. & Osher, S. (1988). Efficient implementation of essentially non-oscillatory shock-capturing schemes. Journal of Computational Physics 77, 439-471. doi:10.1016/0021-9991(88)90177-5

• Simpson, G. & Spiegelman, M. (2011). Solitary wave benchmarks in magma dynamics. Journal of Scientific Computing 49, 268-290. doi:10.1007/s10915-011-9461-y

• Simpson, G., Spiegelman, M. & Weinstein, M. I. (2010). A multiscale model of partial melts: 2. Numerical results. JGR Solid Earth 115, B04411. doi:10.1029/2009JB006376

• Spiegelman, M. (1993). Flow in deformable porous media. Part 1. Simple analysis. Journal of Fluid Mechanics 247, 17-38. doi:10.1017/S0022112093000369

• Stolper, E., Walker, D., Hager, B. H. & Hays, J. F. (1981). Melt segregation from partially molten source regions: the importance of melt density and source region size. JGR 86, 6261-6271. doi:10.1029/JB086iB07p06261

• 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

• Wark, D. A. & Watson, E. B. (1998). Grain-scale permeabilities of texturally equilibrated, monomineralic rocks. EPSL 164, 591-605. doi:10.1016/S0012-821X(98)00252-0