Two-phase flow in Aragog

A magma ocean is rarely fully molten or fully solid. Across most of its evolution, the silicate mantle sits in a mushy regime where solid grains and silicate melt coexist over a finite depth range. Aragog represents this regime as a two-phase mixture of melt and solid that coevolves at every radial node, rather than as a single-phase fluid with a depth-cutoff between molten and solid layers. This page explains the conceptual difference between one-phase and two-phase formulations, lists the additional physics terms the two-phase representation enables, and points to the source files that implement them.

One-phase versus two-phase formulations

In a one-phase magma-ocean model, the mantle is treated as a single homogeneous fluid that is either entirely above the liquidus or entirely below the solidus. A solidification front (a depth or a single state variable) tracks where the transition lives, and the mantle on each side of the front is closed by its own equation of state. The bulk fluid evolves through a global energy balance, and the latent heat of fusion either appears as an effective heat-capacity spike at the phase boundary or as a global solidification term subtracted from the surface flux. Most of the boundary-layer-theory (BLT) lineage of magma-ocean codes (Elkins-Tanton (2008)⁴, Hamano et al. (2013)⁵, Schaefer et al. (2016)⁷) is one-phase by construction: the well-mixed assumption that defines BLT is a one-phase assumption.

In a two-phase model, both the melt and the solid coexist on the same spatial element, with a melt fraction \(\phi(r,t) \in [0,1]\) describing how much of each is present. Solid and melt have independent equations of state, blended by the lever rule into composite material properties. Solid grains and melt droplets can move with respect to each other in response to gravity (gravitational separation) and chemical diffusion (chemical mixing of melt fraction). Phase-change latent heat enters the energy budget as a continuous function of \(\phi\), not as a delta function at the solidus or liquidus. Aragog and SPIDER (Bower et al. (2018)²) are two-phase magma-ocean codes; this is the second defining feature, alongside mixing-length theory, that distinguishes them from the BLT lineage.

The trade-off is the usual one: the two-phase representation is necessary to capture the partial-melt physics that controls the rheological transition and the late-stage solidification window, but it costs an extra state variable per cell, two phase-boundary lookups per cell per call, and three additional flux contributions on the right-hand side.

How Aragog defines the two-phase state

The two-phase state at each radial node is defined by the local pressure \(P\) and specific entropy \(S\). The phase-boundary entropies \(S_\mathrm{sol}(P)\) (solidus) and \(S_\mathrm{liq}(P)\) (liquidus) are tabulated as functions of pressure in the EOS file. The melt fraction is computed from the un-truncated lever rule

The clipped \(\phi\) is the physical melt fraction; the un-truncated \(g\phi\) is retained for the smoothing functions that gate the two-phase fluxes (see below). This single ratio determines how all the other material properties are blended.

The single-pass evaluation lives in EntropyPhaseEvaluator._update_eos (file src/aragog/eos/entropy_phase.py), which mirrors SPIDER's EOSEval_Composite_TwoPhase (eos_composite.c:177-290) so the same blending rules are used in both codes. The phase-boundary lookups happen exactly once per call, and all derived quantities are computed from the cached intermediates rather than via repeated table calls.

The viscosity blend deserves special note. Costa et al. (2009) Figure 2 shows the relative viscosity of crystal-bearing magmas rising by \(\sim 10^{10}\) across the rheological transition near the critical melt fraction \(\phi_c\), while Aragog's solid/liquid end-member viscosities are \(\eta_\mathrm{sol} \sim 10^{21}\) Pa s and \(\eta_\mathrm{liq} \sim 10^{-1}\) Pa s by default, giving an end-member contrast of \(\sim 10^{22}\) across \(\phi\). This is what makes the partial-melt regime numerically stiff and motivates Aragog's smooth tanh blend over a hard if/else.

Two-phase flow contributions

In a one-phase formulation, the heat flux is conduction plus convection only. Two-phase flow adds two additional flux components and one additional source term in the partial-melt regime:

Gravitational separation of melt and solid

Solid and melt have different densities, so they separate vertically by gravity. The mass flux of melt (positive upward when melt is buoyant) is

The permeability factor \(K(\phi)\) spans three asymptotic regimes (Bower et al. (2018) §2.1):

• Stokes settling at high \(\phi\): melt is the matrix and solid grains settle as isolated spheres, derived from Stokes' law (Bower et al. (2018) Eq. 13a).
• Rumpf-Gupte at intermediate \(\phi\): power-law permeability fit to particle-bed measurements (Rumpf & Gupte, 1971, Chem. Ing. Tech. 43, 367)⁶.
• Blake-Kozeny-Carman at low \(\phi\): melt as a percolating fluid through a near-rigid solid matrix.

The three-regime \(\zeta_\mathrm{grav}(\phi)\) formulation that Aragog and SPIDER use is consolidated in Abe (1995)¹ and reviewed in Bower et al. (2018) §2.1 Eqs. 13a-c, where the regime transitions are density-ratio-dependent: \(\rho_\mathrm{liq}/(11.993\,\rho_\mathrm{sol} + \rho_\mathrm{liq})\) for BKC-to-RG and \(\rho_\mathrm{liq}/(0.29624\,\rho_\mathrm{sol} + \rho_\mathrm{liq})\) for RG-to-Stokes. Aragog hard-codes the equal-density-ratio limits (\(\rho_\mathrm{sol} = \rho_\mathrm{liq}\)) of these expressions, \(\zeta_1 = 0.0769452\) (BKC to RG) and \(\zeta_2 = 0.771462\) (RG to Stokes), and blends through them with a tanh switch rather than recomputing per radial node. The resulting permeability-vs-porosity curve (Figure 1 of Heat transport) is JAX-differentiable and shifts the transition location by a few percent in \(\phi\) relative to the density-ratio-dependent form. The corresponding heat flux is

with \(L(P)\) the EOS-tabulated, pressure-dependent latent heat of fusion.

bottom_up_grav_sep = true adds a SPIDER-parity gate that disables \(j_\mathrm{grav}\) below the rheological transition until the solid fraction is interconnected, reflecting the physical fact that melt cannot drain through a fully molten lower mantle.

Chemical mixing of melt fraction

Convective parcels carry melt fraction across the rheological transition. When parcels with different \(\phi\) mix into different ambient \(\phi\), latent heat is released or absorbed as the local lever rule is restored. Bower et al. (2018) Eq. 9 writes this as a Fick's-law flux on the melt fraction itself, \(F_\mathrm{mix} = -\rho T \Delta S_\mathrm{fus}\,\kappa_h\,\partial\phi/\partial r\). Aragog uses the SPIDER-parity bracket form, which vanishes at the solidus and liquidus by construction (no \(\partial\phi/\partial r\) singularity at phase boundaries) and reduces to the Bower (2018) Eq. 9 form in the strict mushy interior. The bracket form makes the smoothness of the RHS explicit at the cost of a more complex algebraic form:

The bracketed expression is the entropy-gradient excess relative to the linear (lever-rule) interpolation between the solidus and liquidus entropy gradients at the local pressure. Outside the mushy band \(\mathrm{smth}(\phi)\) vanishes and so does the flux. The chemical eddy diffusivity \(\kappa_c\) is the same MLT diffusivity as \(\kappa_h\) scaled by the configurable eddy_diffusivity_chemical.

Phase-boundary smoothing

The smoothing function \(\mathrm{smth}(\phi)\) is what keeps the two two-phase fluxes (\(F_\mathrm{grav}\) and \(F_\mathrm{mix}\)) differentiable at the solidus and liquidus. Aragog supports two forms:

Without smoothing, the raw \(\rho\phi(1-\phi)\,v_\mathrm{rel}\) flux drains the CMB cell off the EOS table edge in a single coupling step at first crystallisation; with smoothing, the flux ramps continuously and the JAX Jacobian remains finite across the transition.

Why this matters: where one-phase models lose physics

A one-phase magma-ocean model can be tuned to give a reasonable bulk cooling rate during the fully-molten phase, where the well-mixed assumption is well justified and a single Nu-Ra closure captures the surface flux. The lossy regime is the partial-melt phase: the time window between first crystallisation and final solidification, which for an Earth-mass mantle is several hundred kyr to tens of Myr. In that window:

• Gravitational separation moves melt up and solid down at speeds comparable to the convective velocity in the lowest-permeability regime, transporting both mass and latent heat that a one-phase model cannot represent.
• Chemical mixing carries entropy in the opposite direction to convection in parts of the mushy band, partially cancelling the convective flux on either side of the rheological transition; this is what produces the multi-component flux structure shown in Figure 2 of Heat transport.
• The relative viscosity rises by \(\sim 10^{10}\) across the rheological transition (Costa et al. (2009) Figure 2), well outside any single Ra-Nu fit, so the BLT closure that worked at fully-molten state is no longer accurate.
• The latent heat of solidification, encoded in Aragog as a continuous function of \(\phi\) via the lever rule, becomes a c\(_p\) delta function in a one-phase model, which is either ignored (incorrect) or imposed as a global term that misrepresents its spatial distribution.

In a coupled atmosphere-interior simulation, this is the regime where atmospheric H\(_2\)O outgassing, surface volatile budgets, and the timing of solidification are decided. The one-phase formulation can give the right qualitative answer but will miss the resolved time-dependence and the radial structure that controls the coupling to the atmosphere.

See also

For the algebraic forms of the four heat-flux components Aragog assembles, see Heat transport. For the closure that determines \(\kappa_h\) in the convective flux, see Mixing-length theory. For the SPIDER cross-check, see Aragog vs SPIDER.