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 is positive upward; Aragog treats melt as the lighter phase, so separation moves melt up and solid down:

The tabulated densities can cross so that the melt is locally the denser phase in a shallow region. The absolute value \(\lvert\rho_\mathrm{liq} - \rho_\mathrm{sol}\rvert\) keeps \(v_\mathrm{rel}\) non-negative there, so the separation direction stays fixed (melt up, solid down). Aragog floors the porosity denominator \(\rho_\mathrm{sol} - \rho_\mathrm{liq}\) at about 1 kg m⁻³ in that region, so the porosity there is \((\rho_\mathrm{sol} - \rho)/(1\ \mathrm{kg\,m^{-3}})\) clipped to \([0, 1]\): about 0 where the node density is at or above \(\rho_\mathrm{sol}\), and up to 1 below it. The flux \(j_\mathrm{grav}\) carries the factor \(\phi(1-\phi)\) of the melt fraction, so it vanishes in solid and fully molten nodes whatever the porosity.

\(\eta_\mathrm{mix}\) is the drag viscosity selected by separation_viscosity. The default, "melt", is the fixed liquid viscosity \(\eta_\mathrm{liq}\), matching SPIDER's GetGravitationalHeatFlux. The alternative, "mixture", is the rheological-transition-blended mixture viscosity: it ramps from the liquid value near and above phi_rheo to the solid value below it, so settling locks up through the same transition that sets the bulk rheology. In a coupled caps-off simulation, "melt" keeps the drag viscosity fixed at its low liquid value below phi_rheo, so gravitational separation keeps moving melt out of the deepest cells (melt rises, solid settles; the core-mantle boundary itself carries no mass flux) and the CMB temperature collapses; "mixture" ties the drag viscosity to the same solid-fraction rise that stiffens the bulk rheology, so separation locks up at the same melt fraction and the collapse does not occur. PROTEUS defaults separation_viscosity to "mixture" for this reason; Aragog's own default stays "melt" for SPIDER parity.

The mobility factor \(F(\phi)\) (the permeability \(K(\phi)\) divided by porosity) 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 two regime boundaries, \(\zeta_1 = 0.0769452\) and \(\zeta_2 = 0.771462\), are the porosities at which Aragog switches between adjacent permeability laws, placed at their crossings: Blake-Kozeny-Carman with Rumpf-Gupte at \(\zeta_1\), Rumpf-Gupte with Stokes settling at \(\zeta_2\) (Bower et al. (2018)² §2.1, Eqs. 13a to 13c). For the coded forms, \(\zeta_2\) is the exact crossing and \(\zeta_1\) sits \(1.7 \times 10^{-5}\) below it in porosity (0.02 percent; exact value 0.0769618), far inside the blend width of 0.02. The boundaries follow from the functional forms of the three laws. 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 general 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 fixed porosity thresholds, \(\zeta_1 = 0.0769452\) (BKC to RG) and \(\zeta_2 = 0.771462\) (RG to Stokes), close to the equal-density-ratio limits of these expressions (0.0769645 and 0.7714621), and blends through them with a tanh switch rather than recomputing per radial node. The resulting mobility-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.