Main

Isotopic anomalies in short-lived radiogenic isotopic systems in mantle rocks that record magmatic differentiation processes occurring in the first 100 Myr show that the mantle of Earth has preserved chemical heterogeneities7,8,9,10,11 dating back to its infancy. These findings are corroborated by the noble gas geochemical record that argues for the preservation of these early-formed geochemical reservoirs12,13,14. The solidification of a deep primitive magma ocean alone can explain this early silicate differentiation event8,9. This solidification process can also explain the current seismic structure of the deep mantle, in which large low shear velocity provinces (LLSVPs) and ultralow velocity zones can be interpreted as residual products of primordial magma ocean solidification2,3,4,5. As remnants of magma ocean solidification, the two antipodal deep-rooted (above the core–mantle boundary) LLSVPs1 must play a leading part in global mantle and core dynamics15, plate tectonics and hot-spot magmatism16,17,18,19 during the entire history of the Earth. Therefore, understanding magma ocean solidification from a dynamical and petrological point of view is essential for our comprehension of the long-term evolution of the mantle of Earth and its present-day state.

These geochemical and seismological observations indicate that the last remnants of the terrestrial magma ocean were located deep in the mantle, above the core–mantle boundary, but this remains debated both dynamically and petrologically3,20,21,22,23,24. Classical magma ocean solidification models, similar to those developed for the Moon, stipulate that in a cooling magma ocean, the first solids appear at the bottom of the mantle because of the intersection of liquidus and adiabat at depths, pushing the residual melt upwards21,25. As crystallization proceeds, solidification is expected to occur from the bottom (core–mantle boundary) upwards. An alternative scenario, based on the fact that the solids are more buoyant than the melt in the deep mantle, argues that solidification takes place in the middle of the magma ocean, separating it into basal and shallow magma oceans3,20,23,26. The shallow magma ocean solidifies upwards quickly because of efficient cooling at the surface, whereas the basal magma ocean solidifies slowly, pushing the residual meltdown to the core–mantle boundary.

The issue of top-down versus bottom-up solidification is ultimately controlled by thermodynamic properties that determine (1) where solidification takes place, that is, the intersection of liquidus with temperature, and (2) where solids and liquids accumulate, that is, the density contrast between melt and solids. Moreover, the efficiency of solid–liquid phase separation plays a fundamental part because no fractionation can occur if melts cannot gravitationally segregate from solids, regardless of the petrological or geochemical nature of solids and melts. Depending on the efficiency of solid–liquid phase separation27, the magma ocean can freeze as a homogeneous silicate reservoir (batch crystallization scenario) or form strongly fractionated reservoirs of distinct compositions (fractional crystallization scenario)25.

Magma ocean dynamics and multiphase flow

Previous fluid dynamics studies mostly focused on the issue of solid–liquid phase separation in magma oceans at high melt fraction25,28,29,30,31,32. The question of whether crystals settle or are suspended by the flow has been investigated. If crystals can settle; they are deposited at the bottom of the magma ocean at a rate that is a function of the settling velocity corrected by a re-entrainment flux29,33. Alternatively, in radial thermal evolution models, crystal fraction can be determined assuming local thermodynamic equilibrium by comparing the temperature to the solidus and liquidus at a given depth; convective mixing can be captured by a Fick’s law while a prescribed phase segregation flux can capture the gravitational segregation of crystals from the magma31.

Except for a few works34, less attention has been given to the later stages of magma ocean differentiation, in which the mantle enters a solid-like dynamic regime because of high crystal fraction. Although previous fluid dynamics simulations have investigated the issue of phase separation efficiency, they did not account for the transport of solids and liquids of distinct composition, crystal production and remelting. These assumptions prevent the possibility of scenarios in which (1) a clear distinction can be made between the ___location of crystal formation and that of crystal accumulation and (2) radial magma ocean structure is controlled by the ___location where crystals settle, rather than the ___location where they form.

To address these challenges, we use a new fluid geodynamical approach with the numerical code Bambari (Methods and Supplementary Information section 1). Bambari implements a multiphase flow mathematical formalism based on the averaging method27,35,36,37. We systematically explore a range of mantle differentiation scenarios through two independent mechanisms: petrological differentiation (chemical partitioning) and mechanical differentiation (solid–liquid phase separation). Our simulations account for a density crossover between melt and solid at mid-mantle depth22,23,38 (Supplementary Information section 1). We focus on magma ocean solidification dynamics when the mantle starts to behave as a solid.

We find that an iron-rich basal magma ocean (BMO) always forms above the core–mantle boundary (Fig. 1) and that it corresponds to the last residue of mantle solidification. This finding is irrespective of the depth at which liquidus and adiabat intersect20,21 (that is, where solidification is expected to occur). A BMO forms at the core–mantle boundary (CMB) through the accumulation of dense, FeO-rich solids and liquids. The FeO-rich solids, having low melting points, remelt, whereas the FeO-rich melt, being denser than the solids, migrate downwards because of negative buoyancy. Both processes contribute to the development of the BMO. We first consider the least favourable scenario for producing a BMO21, in which the adiabat intersects the liquidus at the base of the mantle. For this scenario, unidimensional radial models predict that solidification begins at the CMB and progresses upwards.

Fig. 1: Numerical simulations of the solidification of the mantle of Earth from a mushy magma ocean state.
figure 1

Links to the associated videos can be found in the Supplementary Information. a–c, Initial stage: after a rapid early stage of solidification, the averaged melt fraction is approximately at the rheological critical melt fraction, that is, 50%. d–f, Early stage: solidification of the mush occurs in concert with thermocompositional convection at the global mantle scale. Although crystals accumulate in the deep mantle, they are formed at the surface of the planet in cold downwelling plumes. g–i, Late stage: progressive melt extraction from the cumulates differentiates the mantle. At the end of the upper magma ocean solidification, the mantle is heterogeneous. j–l, Final stage: fusible iron-rich silicates progressively pile up in the deep mantle, forming a BMO by downwards melt extraction and remelting of iron-rich solids.

Source data

Full size image

Magma ocean evolution

The initial conditions of our simulations represent a compositionally homogeneous mantle with a uniformly distributed melt fraction of 50%. Large solid production occurs at the surface in the earliest stage of evolution (Fig. 1a,b), because the temperature at the surface of the planet drops below the solidus, regardless of where liquidus and adiabat curves intersect at depth. The crystal-rich shallow layer becomes rapidly gravitationally unstable both thermally—because it is colder than the underlying interior—and compositionally—because it is composed of negatively buoyant crystals (Supplementary Information section 2.1). This crystal-rich unstable layer then sinks as cold downwellings. However, at this point, the solids do not accumulate in the deep mantle because they progressively remelt during their descent (Fig. 1a–c).

With further cooling of the magma ocean (Fig. 1d–f), solids formed at shallow depths gradually accumulate in the lower mantle. Our simulations show that these solids, which are still fed from the surface by cold downwellings, do not remelt anymore (Fig. 1d–f). Notably, these deep-seated solids are not the result of the crystallization of deep mantle melt (where the liquidus intersects the adiabat). These solids of shallow origin carry a chemical signature produced by low-pressure fractionation. More precisely, these solids are FeO-rich and enriched in incompatible elements. This is at odds with what we would expect solely from the values of partitioning coefficients. However, crystals forming near the cold thermal boundary tend to settle, propelling FeO-rich melt upwards by mass conservation. This FeO-rich melt approaches the cold thermal boundary, cools down by diffusion and eventually solidifies. As a result, the cold downwellings are relatively enriched in FeO as well as in incompatible trace elements. For trace elements that are very incompatible at low pressures (for example, Sm, Nd, Lu, Hf and W), their relative abundances in the residual liquid remain unchanged until the system reaches a very low degree of melting (Supplementary Fig. 3). Consequently, these shallow-origin solids, formed by the rapid quenching of melt with unfractionated trace-element ratios, are also expected to exhibit compositions that remain unfractionated relative to the bulk silicate Earth. Later, we demonstrate that the accumulation of these solids partially diminishes the chemical imprint of high-pressure fractionation processes occurring in the deep mantle.

The shallow fractionation depth of the mantle of Earth is at odds with all previous assessments based on one-dimensional solidification models, highlighting the importance of lateral variations in temperature and composition that cannot be accounted for in one-dimensional models. The reason why solids do not substantially remelt during their descent is that (1) the average temperature in the mantle (once the rheological transition is reached) is lower than the liquidus at all depths; (2) the slope of the liquidus favours crystallization at larger depths; and (3) the crystal-rich downwellings sink too fast to heat up and remelt by thermal diffusion (Supplementary Information section 2.3).

The final stage involves the accumulation of denser solids sinking as downwellings in the deep mantle (Fig. 1e–f). This forms a thermal lid above the core, leading to efficient reheating of the lowermost mantle. This in turn melts mantle fusible and FeO-rich components that are so negatively buoyant that they accumulate at the base of the mantle, forming a BMO (Fig. 1g–h). This continuous transport of solids from the surface to the lower mantle, in which they remelt and settle above the CMB, gradually enriches the BMO in FeO compared with the bulk silicate Earth (BSE).

Let us now discuss the case of a steeper liquidus curve20 intersecting the adiabat at mid-mantle depths, in which a BMO is the expected final stage of mantle solidification. We find that the dynamical evolution is essentially the same (Supplementary Fig. 18), which establishes that the crossing of the liquidus and adiabat does not play a main part in the style of solidification of the magma ocean. A liquidus curve that intersects the adiabat in the mid-mantle facilitates remelting in the lower mantle, promoting the formation of a BMO. Although the crystal fraction is higher in the mid-mantle, this sluggish shell does not isolate the upper mantle from the lower mantle. Downwellings that form at the surface do not accumulate at mid-mantle depth as they remain thermally negatively buoyant (Supplementary Information section 3.1 and Supplementary Fig. 18). These downwellings pursue their descent and remelt in the deep mantle.

We conclude that the formation of the basal magma ocean of Earth is primarily controlled by the relatively low melting temperature of FeO-rich silicates and the density contrast between FeO-rich liquid and coexisting solid silicates (that is, solids in thermodynamic equilibrium with that melt) and not by the relation between liquidus and adiabat. At the end of magma ocean solidification, the resulting iron-rich thermochemical structures in the lowermost mantle (Fig. 1) are geophysically and mineralogically consistent with the properties of LLSVPs (Fig. 1) and ultralow velocity zones2,3,4,5,39,40.

Role of solid–liquid phase separation

The solidification sequence described above requires the solid–liquid phase separation to be faster than mantle solidification and remixing by thermochemical convection. Phase separation efficiency is controlled by the dimensionless melt mobility number, δ (Supplementary Information section 2.2, Supplementary Figs. 5 and 11), which is mainly governed by crystal size and melt viscosity. In our simulations, although the viscosity of partially or essentially molten regions is unrealistically large, we preserve a realistic balance of the thermochemical convection velocity and phase separation velocity—thus, the extent of chemical differentiation—by using a relatively large melt mobility number (Supplementary Information section 2.2). To validate our approach and explore the competition between chemical differentiation and convective mixing, we explored numerically the Ra–δ parameter space, where Ra is the thermal Rayleigh number and δ is the melt mobility number. For each thermal Rayleigh number we investigated, we were able to identify a critical melt mobility number above which substantial chemical differentiation is observed (Supplementary Information section 2.2).

By extrapolating our regime diagram to magma ocean conditions, assuming a liquid silicate viscosity of 1 Pa s (refs. 41,42), melts will segregate from the solids if the crystal size is larger than 0.01 μm. Even assuming a four orders of magnitude higher viscosity (104 Pa s), the critical crystal size above which melt segregation can occur is 1 μm (Supplementary Information section 2.2), consistent with previous estimates25,43. As typical grain size estimates in a magma ocean are around 1 mm (ref. 25), phase separation is expected to be efficient and dominate over remixing under realistic conditions.

A fundamental aspect is that the fluid dynamics within or in the vicinity of the top and bottom thermal boundary layers are at the heart of the chemical evolution of the magma ocean. The composition of solids formed at the surface of the planet is controlled by chemical fractionation between melt and low-pressure mineral phases. To transport a fractionated composition into the deep mantle, solids and liquids must segregate from each other in the vicinity of the cold thermal boundary (Supplementary Figs. 7 and 8). Otherwise, it is the bulk composition of the solid and liquid mixture that is brought to the deep mantle. Scaling analysis shows that the solid–liquid segregation in the top thermal boundary layer is expected to be faster than the growth rate of Rayleigh–Taylor instabilities (Supplementary Information section 2.3) in the magma ocean of Earth, allowing for chemical shallow differentiation to occur. Although present, scaling analysis indicates that this process is less pronounced in our fluid dynamics simulations than it would be under real conditions. The symmetric version of the processes described above is responsible for the formation of the BMO. FeO-rich melt migrates downwards within the hot thermal boundary layer before being entrained in upwelling currents (Supplementary Figs. 9 and 10).

Geochemical consequences

We explore the geochemical consequences of our dynamical model on the production and nature of primordial mantle heterogeneities, inherited from magma ocean solidification. We use experimentally determined partition coefficients of trace elements between melts and the liquidus phase at upper mantle (olivine)44 and lower mantle (bridgmanite)45 conditions to track a simple (no garnet and no ferropericlase) evolutionary model of key trace-elements elemental ratios (Sm/Nd, Lu/Hf and Hf/W) during solidification (Supplementary Information section 1.2). At the end of mantle solidification, solids in the upper mantle show superchondritic Lu/Hf ratios (Fig. 2a, yellow to red), owing to olivine crystallization during bottom-up crystallization of the upper mantle. Conversely, in the lower mantle above the BMO, we observe subchondritic Lu/Hf ratios (Fig. 2, blue to white) that indicate bridgmanite fractionation. Furthermore, we observe that solidification and remelting at different depths produce a complex (marble cake-like) geochemical structure in the solid mantle from core to crust (Fig. 2), far from a vision inspired by geochemical two-box modelling of single enriched (residual liquid) and depleted (precipitating solids) reservoirs45,46. Although the extent and magnitude of trace-element ratio distribution depend markedly on phase separation efficiency (Supplementary Figs. 6 and 7), planetary scale solidification systematically generates heterogeneities at all depths (Fig. 2).

Fig. 2: Geochemical signature of magma ocean solidification for an intermediate (δ = 10 × 10−3) phase separation efficiency.
figure 2

a–c, Snapshots at the final stage of the simulations showing the distribution of various geochemical ratios (normalized to their ratio in the BSE) across the mantle Lu/Hf (a), Hf/W (b) and Sm/Nd (c). Sm/Nd and Hf/W ratios are always larger than 1 in solids and lower than 1 in residual melts. By contrast, Lu/Hf behaves similarly at low pressure (olivine–melt fractionation) but becomes less than 1 in solids and greater than 1 in residual melts at high pressures (bridgmanite–melt fractionation). The geochemical signature of both low-pressure and high-pressure fractionation (Lu/Hf ratio) as well as enriched solids and depleted liquids (Sm/Nd and Hf/W ratios) are preserved in the early solid mantle after magma ocean solidification. These are then expected to be further stirred by ensuing thermochemical convection in the solid mantle throughout Earth’s subsequent history. See Supplementary Figs. 6 and 7 for more extreme cases.

Source data

Full size image

The Lu/Hf ratio is a suitable geochemical tracer to quantify the extent of low-pressure (that is, olivine) mantle solidification on a global magma ocean47. This is because Lu is incompatible (that is, enriched in the melt) both in olivine and bridgmanite, whereas Hf is incompatible in olivine but compatible in bridgmanite. In this simple model, solids with superchondritic (high) Lu/Hf ratios stem from olivine crystallization and denote a shallow origin, whereas solids with subchondritic (low) Lu/Hf ratios originate from bridgmanite crystallization and can have formed only at depth. The extent of mixing between these two components on the planetary scale can be seen on a Sm/Nd-Lu/Hf correlation map (Supplementary Fig. 8). We should note, however, that the extent and amplitude of fractionation shown in Fig. 2 varies with the efficiency of solid–liquid phase separation (see Supplementary Figs. 6 and 7 for more extreme cases). The amount of stirring obtained in these two-dimensional Cartesian simulations may also be quantitatively affected by three-dimensional effects and sphericity. Unlike the lunar magma ocean, the magma ocean of Earth does not evolve as a stack of immobile layers each retaining the geochemical characteristics of the depths at which they formed. Rather, primordial heterogeneities are scattered at all depths. This marks a shift in our understanding: the way lunar-like magma oceans evolve, based on the differentiation observed in magma chambers on Earth, does not necessarily apply to mantles as deep as the ones of the Earth or Mars. There, despite stirring during solidification caused by thermochemical convection, early-formed geochemical reservoirs are preserved but scattered throughout the mantle48,49.

We analysed the statistical distribution of geochemical heterogeneities in the solid and molten (BMO) lower mantle, and more specifically the Lu/Hf population statistics (Fig. 3). If we account for only bridgmanite fractionation (Fig. 3a) and disregard olivine fractionation in the upper mantle, then the distribution of Lu/Hf heterogeneities is centred around 0 for the solid mantle, and the BMO composition is shifted towards higher Lu/Hf ratios (70% higher than the BSE) as we would expect from bridgmanite–melt fractionation. However, when both olivine–melt and bridgmanite–melt fractionations are accounted for (Fig. 3b), the distribution of Lu/Hf heterogeneities for the solid mantle and the BMO are both centred around 0, showing that the opposite effects of low- and high-pressure solid–melt fractionation cancel each other out.

Fig. 3: Statistical distribution of Lu/Hf heterogeneities.
figure 3

The statistical distribution is shown as the logarithm of the Lu/Hf normalized to BSE in the lower mantle (grey) and the BMO (red) at the end of magma ocean solidification. The histograms are normalized by the maximum number of occurrences. a, The bridgmanite–melt partition coefficients of trace elements are used at all mantle depths. b, The olivine–melt partition coefficients are used at low pressures and the bridgmanite–melt partition coefficients are used at high pressures. This second and more realistic case highlights the important contribution of low-pressure mineral–melt fractionation in the lower mantle and shows that bridgmanite–melt partitioning coefficients alone are not sufficient to predict the chemical evolution of the lower mantle during magma ocean solidification.

Source data

Full size image

The role of low-pressure chemical fractionation on deep mantle composition has a fundamental consequence on the quantification of the extent of magma ocean differentiation from the standpoint of trace-element ratios in mantle rocks. The primitive upper mantle (PUM) has unfractionated refractory lithophile trace-element ratios50, and because bridgmanite crystallization strongly fractionates some of these ratios45, this has been used to constrain the maximum amount of bridgmanite–melt differentiation that took place in the early Earth. On this basis, and to not disturb those ratios in the PUM (within their uncertainties), no more than 8% (refs. 45,51) bridgmanite can be crystallized alone in the magma ocean, which is inconsistent with experimental melt relations and the melting phase diagrams of pyrolite. Our work relaxes this constraint because bridgmanite–melt partitioning is not the only relevant process that describes the chemical differentiation between the PUM and the lower mantle during magma ocean solidification. Olivine–melt partitioning plays an important part as well, because cold downwellings form at the surface of the planet, constantly feeding the lower mantle with a low-pressure signal. These two effects have signatures that partially cancel each other out, allowing for a much larger extent of lower mantle silicate differentiation (that is, bridgmanite melt) to take place without affecting chondritic ratios in the PUM (Fig. 3).

Magma ocean outgassing

Finally, we conducted simulations to quantify the outgassing of volatile species during magma ocean solidification (see Supplementary Information section 3.2). We used a Lagrangian approach to track and estimate the amount of mantle materials that would rise to the exsolution depth. We assumed that if mantle parcels remain below the exsolution depth, they remain undegassed and retain their primordial noble gas signature and volatile content. In agreement with previous work52, our results indicate that a negligible fraction of the mantle can be expected to degas during magma ocean solidification beyond the rheological transition. This is consistent with the noble gas geochemical record that suggests the preservation of early-formed geochemical reservoirs12,13,14, although this topic remains highly debated53.

Conclusions

Our modelling was performed using the least favourable conditions to produce a BMO, from both a geochemical (full equilibrium fractionation, see Supplementary Information section 1.2) and a petrological (adiabat intersects the liquidus at the base of the mantle) standpoint; yet (1) it systematically produces a BMO in the final stage of evolution, which appears to be inevitable on Earth, and extensible to other large terrestrial planets with Earth-like composition; and (2) the geochemical imprint of this solidification on the solid mantle is far less marked than that predicted by two-box geochemical models47, owing to extensive vertical mixing during solidification. The composition of the PUM integrates the signature of a complex mixture of shallow and deep geochemical components. Concomitantly, the signature of olivine–melt trace element fractionation is present in the lower mantle at the end of magma ocean solidification and could be preserved over geological times.

These findings call for a re-interpretation of the available geochemical record and geophysical observations to better reconstruct the thermal and chemical history of Earth from its infancy to the present day, and more generally to better understand the diversity of terrestrial bodies.

Methods

Recent mathematical and numerical advances27,37,54,55 allow extending the multiphase flow formalism to geodynamic regimes relevant to magma ocean solidification dynamics with vigorous convection (expressed by the thermal Rayleigh number, which is set to 109 in this work) and high melt fraction (0–100%).

We model the solid–liquid multiphase physics using the numerical code Bambari27. Bambari implements a multiphase flow mathematical formalism based on the averaging method27,35,36,37,56,57,58,59,60,61 under infinite Prandtl number and Boussinesq approximations. The original implementation62,63 was considerably optimized in recent years, in particular with the use of the numerical stencils implemented for the momentum conservation (equation (6)) in the Finite-Volume code StreamV64,65. The latter allows handling reliably large and sharp viscosity contrasts with negligible spurious pressure effects66.

Bambari accounts for thermochemical convection in both liquid and solid states, solid–melt phase equilibrium, mineralogical phase change and solid–melt phase separation. Convective motion in a solidifying magma ocean is driven by three types of density difference. These originate from thermal expansivity (due to temperature differences), compositional differences (due to changes in iron content) and phase changes (due to the varying melt fraction)67. An additional mechanism that generates motion in a partially molten convective medium is the solid–liquid phase separation driven by shear deformation68, or density contrasts between the melt and the solid. Phase separation is limited by matrix deformation, that is, compaction, and viscous drag between the melt and the solid using Darcy’s law. We implemented a depth-dependent density contrast between liquids and solids22,23,38 and used a thermodynamically self-consistent compositional evolution of melts and solids based on recent experimental melting phase diagrams22,69 (Supplementary Information section 1 and Supplementary Figs. 1 and 2).

The model consists of two mechanical phases: a liquid phase and a solid phase. Each mechanical phase is composed of two main compositional phases corresponding to a FeO-rich end-member and a MgO-rich end-member. Trace elements are distributed among the mechanical phases according to their solid–liquid partition coefficients. Apart from the advection and diffusion of temperature, our model tracks the advection of 14 compositional fields: 7 liquids (active MgO-rich end-member, FeO-rich end-member and passive chemical species Hf, Sm, Nd, W and Lu) and the 7 solids counterparts.

In the following, we recall the governing equations solved by the code in their dimensionless form. We perform non-dimensionalization of lengths by the thickness of the mantle H, time by the thermal diffusion scale H2/κ (where κ is the average coefficient of thermal diffusivity), velocities by κ/H, temperatures by the total super-adiabatic temperature contrast across the upper and lower thermal boundary layers in the mantle ΔTm and pressure and deviatoric stress using viscous pressure scale ηsκ/H2, where ηs is the viscosity of the solid phase.

Mass conservation equations

The mass conservation equations for the four main chemical components are

$$\begin{array}{rcl} & & \frac{\partial {\phi }_{1}}{\partial t}+{\boldsymbol{\nabla }}\cdot {{\bf{v}}}_{{\rm{s}}}{\phi }_{1}=-\,{\varGamma }_{1}\\ & & \frac{\partial {\phi }_{2}}{\partial t}+{\boldsymbol{\nabla }}\cdot {{\bf{v}}}_{{\rm{s}}}{\phi }_{2}=-\,{\varGamma }_{2}\\ & & \frac{\partial {\phi }_{3}}{\partial t}+{\boldsymbol{\nabla }}\cdot {{\bf{v}}}_{{\rm{l}}}\,{\phi }_{3}=+\,{\varGamma }_{1}\\ & & \frac{\partial {\phi }_{4}}{\partial t}+{\boldsymbol{\nabla }}\cdot {{\bf{v}}}_{{\rm{l}}}\,{\phi }_{4}=+\,{\varGamma }_{2}\end{array}$$
(1)

where t is the time, vl is the velocity of the liquid phase, vs is the velocity of the solid phase and Γi are associated with the rate of phase change. The subscripts 1 and 3 refer to the MgO end-members in the solid and liquid phases, respectively. The subscripts 2 and 4 refer to the FeO end-members in the solid and liquid phases, respectively. The densities of the four chemical end-members are

$$\begin{array}{rcl} & & {\rho }_{1}={\rho }_{0}(1-\alpha T)+\frac{1}{2}\Delta {\rho }_{1}-\frac{1}{2}\Delta {\rho }_{2}\\ & & {\rho }_{2}={\rho }_{0}(1-\alpha T)+\frac{1}{2}\Delta {\rho }_{1}+\frac{1}{2}\Delta {\rho }_{2}\\ & & {\rho }_{3}={\rho }_{0}(1-\alpha T)-\frac{1}{2}\Delta {\rho }_{1}-\frac{1}{2}\Delta {\rho }_{2}\\ & & {\rho }_{4}={\rho }_{0}(1-\alpha T)-\frac{1}{2}\Delta {\rho }_{1}+\frac{1}{2}\Delta {\rho }_{2},\end{array}$$
(2)

where ρ0 is the reference density, the thermal expansivity α is considered constant in the solid and liquid phases, and T is the temperature. The isochemical density contrast between the liquid and solid phases is equal to Δρ1 = ρ1 − ρ3 = ρ2 − ρ4. Similarly, the chemical density contrast between the dense and light end-members is Δρ2 = ρ2 − ρ1 = ρ4 − ρ3. We parameterized these density contrasts to fit a self-consistent thermodynamic model22 (Supplementary Fig. 1).

Energy conservation equation

We assume local thermal equilibrium between all the phases, that is, the temperatures at a given ___location in all the phases are the same (T1(x, z, t) = T2(x, z, t) = T3(x, z, t) = T4(x, z, t)). This leads to a single equation for the conservation of energy:

$${\partial }_{t}T+\overline{{\bf{v}}}\cdot {\boldsymbol{\nabla }}T={{\boldsymbol{\nabla }}}^{2}T-({\varGamma }_{1}+{\varGamma }_{2}){S}_{t},$$
(3)

where \(\bar{{\bf{v}}}\) is the velocity of the solid–liquid mixture, St is the dimensionless Stefan number, given by

$${S}_{t}=\frac{L}{{C}_{{\rm{p}}}\Delta {T}_{{\rm{m}}}},$$
(4)

where Cp is the thermal capacity at constant pressure for the mixture of the four components and L is the latent heat release on phase change.

The velocity of the solid–liquid mixture is related to the solid and liquid velocities using the following equation:

$$\overline{{\bf{v}}}=\phi {{\bf{v}}}_{{\rm{l}}}+(1-\phi ){{\bf{v}}}_{{\rm{s}}}\,.$$
(5)

Momentum conservation equations

As we account for two mechanical phases, two momentum conservation equations are required to describe the two-phase flow mechanical interactions. The following Stokes equation describes the momentum conservation for the averaged solid–liquid mixture:

$$-{\boldsymbol{\nabla }}\varPi +{\boldsymbol{\nabla }}\cdot \bar{\underline{{\boldsymbol{\tau }}}}-{\rm{Ra}}\,T\,{\bf{g}}-{\rm{Rp}}\,\phi \,{\bf{g}}+\frac{1}{2}\,{\rm{Rc}}\,({\phi }_{2}+{\phi }_{4}-{\phi }_{1}-{\phi }_{3}){\bf{g}}=0,$$
(6)

where ∇Π is the dynamic pressure gradient (∇Π = ∇P − (ρ0 + 1/2∆ρ1)g, where P is the total pressure, g is the gravity vector and ρ0 is the reference density) and \(\bar{\underline{{\boldsymbol{\tau }}}}\) is the deviatoric stress tensor for the solid–liquid mixture. The thermal Rayleigh number, Ra, and the phase and compositional Rayleigh numbers, Rp and Rc, respectively, are

$$\begin{array}{rcl} & & {\rm{Ra}}=\frac{{\rho }_{0}\alpha \Delta {T}_{{\rm{m}}}\,g{H}^{3}}{\kappa {\eta }_{{\rm{s}}}},\\ & & {\rm{Rp}}=\frac{\Delta {\rho }_{1}g{H}^{3}}{\kappa {\eta }_{{\rm{s}}}},\\ & & {\rm{Rc}}=\frac{\Delta {\rho }_{2}g{H}^{3}}{\kappa {\eta }_{{\rm{s}}}},\end{array}$$
(7)

where g is the constant acceleration of gravity (the magnitude of g). We define the stress tensor of the solid–liquid mixture as

$$\bar{\underline{{\boldsymbol{\tau }}}}=(1-\phi ){\underline{{\boldsymbol{\tau }}}}_{{\rm{s}}}+\phi {\underline{{\boldsymbol{\tau }}}}_{{\rm{l}}}=\eta (\phi )({\boldsymbol{\nabla }}\overline{{\bf{v}}}+{[{\boldsymbol{\nabla }}\overline{{\bf{v}}}]}^{t}),$$
(8)

where η(ϕ) is the viscosity of fluid mixture that varies from η(ϕ = 0) = ηs to η(ϕ = 1) = ηl. We use

$$\eta (\phi )={\eta }_{{\rm{s}}}\times 1{0}^{\left(1+\tanh \left(\frac{\phi -0.5}{0.1}\right)\right)\frac{{\eta }_{{\rm{l}}}}{2{\eta }_{{\rm{s}}}}}.$$
(9)

For the second momentum conservation equation, we use a retro-action relationship27,36 that describes Δv = vs − vl, the velocity difference between the two phases:

$$\phi \Delta {\bf{v}}={\phi }^{2}(1-\phi ){\delta }^{2}\left(\zeta {\boldsymbol{\nabla }}\left[\frac{(1-\phi )}{\phi }{\boldsymbol{\nabla }}\cdot \phi \Delta {\bf{v}}\right]+X{\bf{g}}+{\boldsymbol{\nabla }}\cdot \bar{\underline{{\boldsymbol{\tau }}}}\right),$$
(10)

where

$$X={\rm{Rp}}+\frac{1}{2}\left(\frac{{\phi }_{3}-{\phi }_{4}}{\phi }+\frac{{\phi }_{2}-{\phi }_{1}}{1-\phi }\right){\rm{Rc}},$$
(11)

where δ is the melt mobility number and ζ is the dimensionless compaction viscosity. The melt mobility number and dimensionless compaction viscosity write

$$\begin{array}{r}\delta =\sqrt{\frac{{r}_{0}^{2}}{{C}_{0}{H}^{2}}\frac{{\eta }_{{\rm{s}}}}{{\eta }_{{\rm{l}}}}},\\ \zeta =\frac{{\eta }_{{\rm{c}}}}{{\eta }_{{\rm{s}}}},\end{array}$$
(12)

with r0 is the crystal size, C0 is a constant whose value is given in the Supplementary Information and ηc is a constant associated with the compaction viscosity, sometimes referred to as bulk viscosity, which describes the viscous resistance of a solid matrix to compact (that is, close the porous space)35,36,70,71. Here we follow ref. 36 in which the compaction viscosity, ηbulk, is defined as ηbulk = ηc/ϕ (the compaction viscosity must be infinite when the melt fraction is 0).

The above equations are discretized using conservative finite-difference schemes. The corresponding set of equations (6) and (10) involves the inversion of two sparse matrices, which are the most time-consuming operations. In Bambari, these matrix inversions are performed using the direct PARDISO library that is parallelized with OpenMP directives72,73. The numerical implementation of the two-phase flow physics has been benchmarked against one-dimensional analytical solutions (for example, section 3.4 in ref. 63). The new additions to the numerical implementation were successfully benchmarked against numerical and analytical solutions65,66.

Phase change and chemical partitioning

Methods related to the modelling of phase change and chemical partitioning can be found in Supplementary Information.