# Solid-state mantle convection coupled with a crystallising basal magma ocean

Stéphane Labrosse, Adrien Morison, Paul James Tackley

► To cite this version:

Stéphane Labrosse, Adrien Morison, Paul James Tackley. Solid-state mantle convection coupled with a crystallising basal magma ocean. Comptes Rendus. Géoscience, 2024, 356 (S1), pp.1-17. 10.5802/cr-geos.275. hal-04795672

HAL Id: hal-04795672

https://hal.science/hal-04795672v1

Submitted on 21 Nov 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Comptes Rendus
## Géoscience
### Sciences de la Planète

Stéphane Labrosse, Adrien Morison and Paul James Tackley

**Solid-state mantle convection coupled with a crystallising basal magma ocean**

Published online: 12 November 2024

**Part of Special Issue:** Magma degassing and its impact on the Earth's atmosphere: from magma oceans to lava lakes

**Guest editors:** Manuel Moreira (Institut des Sciences de la Terre d'Orléans Université d'Orléans-CNRS-BRGM 1a rue de la Férollerie 45071 Orléans France), Bruno Scaillet (Institut des Sciences de la Terre d'Orléans Université d'Orléans-CNRS-BRGM 1a rue de la Férollerie 45071 Orléans France) and Clive Oppenheimer (Department of Geography, University of Cambridge, Downing Place, Cambridge CB2 3EN, UK)

https://doi.org/10.5802/crgeos.275

<figure><img src="image_2.png" /></figure>

This article is licensed under the

CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
<figure><img src="image_3.png" /></figure>

Research article

Magma degassing and its impact on the Earth's atmosphere: from magma oceans to lava lakes

# Solid-state mantle convection coupled with a crystallising basal magma ocean

Stéphane Labrosse<sup>a, *, a</sup>, Adrien Morison<sup>b, b</sup> and Paul James Tackley<sup>c, c</sup>

<sup>a</sup> LGLTPE, ENS de Lyon, Université de Lyon, 46 allée d'Italie, 69003 Lyon, France

<sup>b</sup> Univ Exeter, Phys & Astron, Exeter, Devon, England

<sup>c</sup> Department of Earth and Planetary Sciences, ETH Zürich, Sonneggstrasse 5, Zürich, 8092, Switzerland

*E-mails*: stephane.labrosse@ens-lyon.fr (S. Labrosse), A.Morison@exeter.ac.uk (A. Morison), paul.tackley@eaps.ethz.ch (P. J. Tackley)

**Abstract.** Fractional crystallisation of a basal magma ocean (BMO) has been proposed to explain the formation of large scale compositional variations in the mantle and the persistence of partially molten patches in the lowermost mantle. We present a complete set of equations for the thermal and compositional evolution of the BMO and show that it can be implemented in a mantle convection code to solve the long term mantle evolution problem. The presence of the BMO modifies the dynamics of the mantle in several ways. The phase equilibrium at the bottom of the solid mantle implies a change of mechanical boundary condition, which helps solid state convection. The net freezing of the BMO implies a change of computational domain, which is treated by mapping the radial coordinate on a constant thickness domain. Fractional melting and freezing at the boundary makes the composition of the BMO and the solid mantle evolve, which is treated using Lagrangian tracers. A sample calculation shows that the persistence of the BMO and its long term evolution drastically changes the dynamics of the solid mantle by promoting downwelling currents and large scale flow. The gradual increase of the FeO content in the BMO and in the solid that crystallises from it leads to the stabilisation of large scale thermo-compositional piles at the bottom of the mantle, possibly explaining the observations from seismology.

**Keywords.** Thermal evolution, Mantle convection, Core cooling, Basal magma ocean.

*Manuscript received 21 May 2024, revised 3 October 2024, accepted 7 October 2024.*

## 1. Introduction

The long term evolution of the Earth is paced by convection in the solid mantle, which is much slower than the dynamics of the underlying liquid core. The present day structure of the mantle is rather well constrained, thanks to the tremendous progress of seismic tomography over the last few decades [e.g.

Fichtner et al., 2024], and is generally well understood in the context of mantle convection in a plate-tectonic regime [e.g. Coltice et al., 2019]. Paleogeographic reconstructions have been used to constrain models of mantle convection going back 1 Gyr [Flament et al., 2022] although, of course, the uncertainties regarding paleogeography increase drastically with ages larger than about 200 Myr, the age of the oldest oceanic plate. These studies generally consider the presence of chemically denser

<footnote>* Corresponding author</footnote>
material in the form of thermochemical piles at the bottom of the mantle as a way to explain the seismically imaged large low shear velocity provinces (LLSVPs) [e.g. Hernlund and McNamara, 2015, for a review].

If we accept the chemical interpretation of LLSVPs, the question of their formation remains. Several scenarios have been proposed. The formation of the Earth could make the mantle initially chemically stratified in two layers and the slow erosion of the stratification by mantle convection could lead to the present state [e.g. Le Bars and Davaille, 2004]. Instead of gradually mixing an initial stratification, dense material produced by the extraction of the oceanic crust could accumulate at the bottom of the mantle after being subducted [Christensen and Hofmann, 1994, Nakagawa et al., 2010, Li and McNamara, 2013]. A last scenario, which is the topic of the present paper, proposes that the fractional crystallisation of a basal magma ocean (BMO) can lead to the stabilisation of chemically different piles at the bottom of the mantle [Labrosse et al., 2007].

In addition to LLSVPs, seismological studies of the deep mantle have uncovered regions of very reduced seismic velocity, on a much smaller scale than LLSVPs, termed ultra low velocity zones (ULVZs). These regions, that can be tens to hundreds of km wide, have seismic velocities about 30% (for S) and 10% (for P) smaller than the surrounding mantle, and a density a few percent larger [Rost et al., 2005], while typical variations of the seismic velocities on the large scale are of the order of a few percent. These extreme velocity reductions, and the fact that the velocity of S waves is more reduced than that of P waves, has been used to imply the presence of partial melt [Williams and Garnero, 1996]. If indeed the mantle in contact with the core is currently partially molten, more melt should have been present in the past when the core was hotter. Indeed, we know that the Earth has been cooling down [Jaupart et al., 2015] and the possibility of the core to be cooling faster than the mantle has been invoked [Driscoll and Bercovici, 2014, Labrosse, 2016, Patočka et al., 2020] to solve the long standing thermal catastrophe problem for the thermal evolution of the Earth [Christensen, 1985]. The rapid core cooling is also a consequence of the necessity to maintain a convective dynamo with a large thermal conductivity of the core [Labrosse, 2015, Patočka et al., 2020]. All these argu-

ments together led to the scenario of a basal magma ocean [Labrosse et al., 2007].

The aim of this paper is to present a first step toward including a basal magma ocean in a fully dynamical model of mantle evolution. Starting with a mantle convection model, in this case StagYY [Tackley, 2008], a first ingredient to add is the possibility of a solid-melt phase change at the boundary with the underlying magma ocean. Compositional changes associated with the phase change are also required. Both aspects have already been reported and the implications of such a boundary have been explored [Labrosse et al., 2018, Agrusta et al., 2019, Bolrão et al., 2021, Lebec et al., 2023, 2024]. The possibility of melting and freezing at one of the horizontal boundaries helps convection in the solid. This effect is included by applying a phase change boundary condition [Alboussière et al., 2010, Mizzon and Monnereau, 2013, Deguen et al., 2013, Deguen, 2013, Labrosse et al., 2018] controlled by a single dimensionless parameter, the phase change number Φ. For small values of this parameter (Φ ≲ 10), the phase change is fast and the critical Rayleigh number for the onset of convection is reduced compared to the situation without phase change or with a slow phase change (Φ ≳ 10³). The heat and mass transfer is also enhanced by the phase change at the boundary.

These previous studies on the effect of a solid-liquid phase change at the boundary did not include the net evolution of the planet, with the possibility of volume change of the basal magma ocean. Including this effect requires several important modifications of the model, which are presented below. Firstly, a numerical treatment of the moving boundary at the bottom of the solid mantle is necessary and this is the topic of Section 2.1. This requires knowing the moving rate of the boundary, which can be determined from the energy balance of the BMO. The relevant theory is presented in Section 2.2. An example calculation is then presented in Section 3 before discussing the implications and limitations of the model in Section 4.

## 2. Physical and numerical model

### 2.1. Convection in the solid mantle

We start with a mantle convection numerical code, StagYY, which can solve the equations for mass,
composition, energy and momentum balances in an infinite-Prandtl-number fluid, like planetary mantles, in various geometries and with many complexities [Tackley, 2008]. This code is widely used in the geodynamics community and needs not be detailed here; only features directly relevant to the present study are presented. The code uses a finite volume approach for the mass, energy and momentum balances and a particle-in-a-cell (PIC) approach for composition [Tackley and King, 2003, Gerya and Yuen, 2003, Ismail-Zadeh and Tackley, 2012]. In the present study, we use the spherical annulus geometry [Hernlund and Tackley, 2008] since the sphericity is important for a proper surface to volume scaling in an evolving planet. In order to keep the study as simple as possible, we consider an incompressible mantle using the Boussinesq approximation, with all physical quantities uniform. With these assumptions, the solid part of the model is controlled by three dimensionless numbers, the Rayleigh number Ra, the internal heating rate H and the buoyancy number B, their usual definition being,

$$Ra = \frac{\alpha g \Delta T d^3}{\kappa \nu}, \quad H = \frac{\rho h d^2}{k \Delta T}, \quad B = \frac{\beta}{\alpha \Delta T}, \quad (1)$$

with α and β the thermal and compositional (FeO) expansion coefficients, respectively, g the gravity, κ the thermal diffusivity, v the kinematic viscosity, ΔT the temperature scale, d the thickness of the whole mantle (solid and BMO), ρ the density and h the time-evolving radioactive heating rate. Among the parameters entering the definition of the dimensionless numbers, ΔT and d deserve a special discussion since the actual temperature drop across the solid mantle and its thickness, which are normally used as scales to define the dimensionless numbers, vary with the progressive crystallisation of the BMO. We use fixed values (see Appendix) and take into account the variations of the actual values as time-varying factors in the balance equations emerging from our way of dealing with the moving boundary, which is the main modification of the numerical model as presented in previous studies. The solid mantle shell is bounded by two spheres, the basal ocean-mantle boundary (BOMB) with dimensionless radius r<sub>BOMB</sub>(t), and the Earth surface with dimensionless radius R<sub>E</sub>. We can map the time-varying

spherical shell to a constant one, between rescaled radii 1 and 2, by using as radial coordinate

$$z = 1 + \frac{r - r_{\text{BOMB}}}{R_{\text{E}} - r_{\text{BOMB}}}. \quad (2)$$

This change of radial coordinate implies a change of both radial and temporal derivatives, obtained using the chain rule. Considering a function f(r, t) = ſ̃(z, t), we get

$$\frac{\partial f}{\partial r} = \frac{\partial z}{\partial r} \frac{\partial \tilde{f}}{\partial z} = \frac{1}{R_{\text{E}} - r_{\text{BOMB}}} \frac{\partial \tilde{f}}{\partial z}, \quad (3)$$

and so on for higher order derivatives. Consider now the time derivative,

$$\frac{\partial f}{\partial t} = \frac{\partial \tilde{f}}{\partial t} + \frac{\partial z}{\partial t} \frac{\partial \tilde{f}}{\partial z} = \frac{\partial \tilde{f}}{\partial t} - \dot{r}_{\text{BOMB}} \frac{R_{\text{E}} - r}{(R_{\text{E}} - r_{\text{BOMB}})^2} \frac{\partial \tilde{f}}{\partial z}, \quad (4)$$

the overdot standing for time derivative. Therefore, the change of BMO radius implies an additional advection term in the equation of energy balance when mapping the domain to a constant one. This additional advection is also applied to tracers. Introducing the time-dependent scaled solid mantle thickness, Γ = R<sub>E</sub> - r<sub>BOMB</sub>, the dimensionless balance equations for mass, momentum, energy, FeO mass fraction ξ and concentration c<sub>HPE</sub> of heat producing elements (HPE) become

$$0 = \nabla \cdot \mathbf{u}, \quad (5)$$

$$0 = -\nabla p + \nabla^2 \mathbf{u} + Ra\Gamma^3 [(T - B\xi) - \langle T - B\xi \rangle] \mathbf{e}_z, \quad (6)$$

$$\Gamma^2 \frac{\partial T}{\partial t} = -\mathbf{u} \cdot \nabla T - \dot{\Gamma} \Gamma(2-z) \frac{\partial T}{\partial z} + \nabla^2 T + H c_{\text{HPE}} \Gamma^2 \exp \left( - \frac{t \ln 2}{\tau_{\text{HPE}}} \right), \quad (7)$$

$$\Gamma^2 \frac{\partial \xi}{\partial t} = -\mathbf{u} \cdot \nabla \xi, \quad (8)$$

$$\Gamma^2 \frac{\partial c_{\text{HPE}}}{\partial t} = -\mathbf{u} \cdot \nabla c_{\text{HPE}}, \quad (9)$$

with u the flow velocity, T the temperature, p the dynamical pressure. For simplicity, we consider only one heat producing element with a unique half life parameter τ<sub>HPE</sub> that represents a mean of the four main ones in the Earth, <sup>235</sup>U, <sup>238</sup>U, <sup>232</sup>Th and <sup>40</sup>K. The time derivative ˙Γ<sub>BOMB</sub> is computed from the balance equations for the BMO, which are explained in the next section. This calculation requires knowing the fluxes of heat and FeO, which are directly computed
in StagYY from the temperature field and tracers, respectively. The change with time of bottom temperature, $T_{BOMB}(t)$, does not appear directly in the equations above but is applied as boundary condition at the bottom of the solid mantle. Therefore, the effective Rayleigh and buoyancy numbers at any given time are

$$Ra_{\text{eff}} = Ra \Gamma^3 T_{\text{BOMB}}, \quad B_{\text{eff}} = B \frac{(\langle\xi_{\text{bot}}\rangle - \langle\xi_{\text{top}}\rangle)}{T_{\text{BOMB}}}, \quad (10)$$

where the actual FeO mass fraction difference across the layer is introduced in the definition of the effective buoyancy number.

An important feature of the model, already included in StagYY for a few previous studies [Agrusta et al., 2019, Bolrão et al., 2021, Lebec et al., 2023, 2024], is the solid-liquid phase change boundary condition at the bottom of the solid shell [e.g. Labrosse et al., 2018],

$$2 \frac{\partial u_r}{\partial r} - p - \Phi u_r = 0, \quad (11)$$

$u_r$ being the radial velocity, $p$ the dynamic pressure and $\Phi$ the phase change number. This dimensionless parameter is the ratio of two timescales, $\Phi = \tau_\phi/\tau_\eta$, with $\tau_\eta$ the timescale to create a topography from viscous stress in the mantle and $\tau_\phi$ the timescale to erase it by convection in the liquid layer. The boundary condition, initially derived for the dynamics of the inner core [Deguen et al., 2013], expresses the competition between the generation of a topography by viscous stress in the solid and its removal by melting and freezing. With this equation, the boundary condition can be the classical no-penetrative one for $\Phi \to \infty$, in which case the phase change is effectively prohibited ($u_r = 0$), or of the flow-through type for small values of $\Phi$. This boundary condition strongly affects convection: for small values of $\Phi$, convection is easier to start [i.e. the critical Rayleigh number for the onset of convection is reduced, see Deguen, 2013, Labrosse et al., 2018, Morison et al., 2024], heat and mass transfer are increased and the wavelength of convection is increased compared to the situation usually considered in mantle convection [Agrusta et al., 2019]. For a purely thermal problem [Deguen et al., 2013], the phase change number can be expressed as

$$\Phi = \frac{\rho_s (\rho_l - \rho_s) L d}{\rho_l^2 C_{pl} (m_{ad} - m_p) u_l \eta_s}, \quad (12)$$

with $\rho_s$ and $\rho_l$ the density of the solid and the liquid, respectively, $L$ the latent heat of fusion, $C_{pl}$ the heat capacity of the liquid, $u_l$ the typical flow velocity in the magma ocean, $\eta_s$ the viscosity of the solid and $(m_{ad} - m_p)$ the difference between the adiabatic gradient in the liquid and the Clapeyron slope. Several parameters are rather uncertain but reasonable estimates give $\Phi \sim 10^{-8}$ (see Appendix). With such a very small value, we are clearly in the flow-through regime. In that regime, down-welling currents reaching the boundary do not turn and, instead, melt to reach the BMO. For this to happen, the temperature has to reach the boundary temperature, which happens on a short length-scale. Resolving that thin boundary layer can be challenging at large Rayleigh number and we use here the technique introduced by Agrusta et al. [2019]: the fixed temperature boundary condition is replaced by a laterally varying Robin boundary condition representing the behaviour on the inner edge of the boundary layer:

$$\mathcal{H}_s(u_r + u_0)\theta + [1 - \mathcal{H}_s(u_r + u_0)] \frac{\partial \theta}{\partial r} = 0, \quad (13)$$

with $\mathcal{H}_s$ a smooth approximation of the Heaviside function, $\theta$ the temperature anomaly, and $u_0$ the reference velocity at which the boundary condition switches behaviour. For $u_r < -u_0$, in downwelling regions, the outflow of solid material toward the BMO makes $\mathcal{H}_s = 0$ and therefore leads to an imposed zero gradient. Conversely, where $u_r > -u_0$, $\mathcal{H}_s = 1$ and the value of the temperature is set to the freezing one. See Agrusta et al. [2019] for more details and benchmark comparisons.

The PIC method consists of using Lagrangian tracers in the domain to carry various quantities as they are transported by the flow, and averaging the relevant quantities on the finite volume grid to compute the flow solution. In the context of the present study, we use a unique compositional information, representing the FeO content of the mineral assemblage. A complexity added to this approach by the phase change boundary condition and already treated by Bolrão et al. [2021] comes from the fact that the flow crossing the phase change boundary implies exchange of FeO with the BMO. In practice, when solid material crosses the boundary by melting, the associated tracers are removed while new tracers are added to regions where crystallisation occurs. These new tracers are given the information associated with
their FeO content according to a simple phase diagram model (Equation (19)) [Bolrão et al., 2021]. The net flux of FeO is computed by compiling the information carried by both removed and added tracers. The composition acts on the density through a linearised equation of state to give the last term on the right-hand-side of the momentum equation (6). The buoyancy number B measures the tendency of FeO entering from the BMO to stabilise against entrainment by thermal convection.

Note that this approach assumes a spherical boundary between the solid mantle and the BMO, i.e. a negligible topography of that boundary. The dynamic topography associated to convection in the solid mantle is in fact included in the theory that leads to the phase change boundary condition Equation (11) [see Labrosse et al., 2018, for a full development], which assumes it to be small. This assumption is consistent with the fact that it is limited by convection in the liquid ocean that tends to erase it by mixing laterally solute and energy.

In addition to the phase change at the bottom boundary of the solid mantle, partial melting could in principle occur anywhere in the bulk, depending on the local temperature and the solidus of the local composition. This possibility is already implemented in StagYY [e.g. Nakagawa and Tackley, 2012], the melt being then extracted to the surface to form a crust, but has not been used in the present paper to focus on the phase change at the bottom boundary.

## 2.2. Evolution of the BMO

The flow of heat and FeO across the boundary between the solid mantle and the BMO makes the boundary move by melting or freezing according to the phase diagram and the conservation of heat and solute in the BMO. The phase diagram prescribes the temperature of the liquidus, $T_L$, and the solidus, $T_s$, as function of composition and pressure. At the boundary between the solid and the liquid, the temperature equals the solidus of the solid composition and the liquidus of the liquid composition. This dependence on composition makes this evolution depend also on the net flux of FeO between the solid and the liquid, which is due to the partitioning between the two phases upon crystallisation and melting. The heat flow out of the BMO controls the rate of cooling and crystallisation of the BMO to form the

solid mantle and the rate of cooling and crystallisation of core. The evolution of the BMO follows from the global energy balance in much the same way as the evolution of the core. The equations are therefore similar and we follow the theory presented by Labrosse [2015].

### 2.2.1. Evolving reference state

We consider the BMO to be composed of an entirely liquid magma, whose composition can evolve between two end-members, an MgO-rich one and an FeO-rich one. The composition is quantified by the mass fraction of FeO, $\xi_l$, the corresponding one in the solid mantle being $\xi_s$.

Fractional crystallisation at the top of the BMO releases FeO at the top, which drives compositional convection. Cooling from the top and heating by the underlying core also favor convection. Since both temperature and composition promote convection, we assume that the BMO stays well mixed at all times, such that $\partial\xi_l/\partial r = 0$, and that it is also isentropic on average. Alternatively, the BMO could start stably stratified [Laneuville et al., 2018] but we neglect this possibility in the present study.

We therefore consider a well-mixed isentropic magma ocean whose reference profile of density, temperature, chemical potential, mass fraction of FeO, etc., are linked by the phase equilibrium occurring at the top of the layer. The system is characterised by three state variables, the specific entropy $s$, the mass fraction of FeO in the magma $\xi_l$ and pressure $P$. The reference state is well mixed, $\partial\xi_l/\partial r = 0$, isentropic, $\partial s/\partial r = 0$, and hydrostatic, $\partial P/\partial r = -\rho g$. For the sake of simplicity, the density $\rho$ is kept constant, and so is the gravitational acceleration $g$. In the reference state, the radial derivatives of the temperature $T$ and chemical potential $\mu$ come only from the pressure variation. The chemical potential can be obtained using thermodynamic identities by integration of [e.g. Braginsky and Roberts, 1995, Lister and Buffett, 1995, Labrosse, 2015]

$$ \frac{\partial \mu}{\partial r} = \left( \frac{\partial \mu}{\partial P} \right)_{s,\xi} \frac{dP}{dr} = -\beta g, \quad (14) $$

$\beta$ being the coefficient of chemical contraction

$$ \beta = -\frac{1}{\rho} \left( \frac{\partial \rho}{\partial \xi_l} \right)_{P,s} = \rho \left( \frac{\partial \mu}{\partial P} \right)_{s,\xi}. \quad (15) $$

Note that for FeO in a magma, $\beta$ here defined is negative. This is opposite to the case of light elements
in the core. Assuming, for simplicity, $\beta$ and $g$ to be constant, Equation (14) can readily be integrated to give

$$ \mu = \mu^{\text{BOMB}} + \beta g(r_{\text{BOMB}} - r) \equiv \mu^{\text{BOMB}} + \mu' \quad (16) $$

with $\mu^{\text{BOMB}}$ the chemical potential at the top of the BMO and $\mu'$ the deviation from that value in the BMO.

Similarly, for the temperature we get the classical isentropic gradient:

$$ \frac{\partial T}{\partial r} = \left(\frac{\partial T}{\partial P}\right)_{s,\xi} \frac{dP}{dr} = -\frac{\alpha g T}{C_p}, \quad (17) $$

with $\alpha$ the thermal expansion coefficient and $C_p$ the heat capacity at constant pressure. We know that $\alpha$ generally decreases with pressure and therefore depth [Anderson et al., 1992, Chopelas and Boehler, 1992, Duffy and Ahrens, 1993, Ricard, 2007, Ricard et al., 2022] but since we consider here a magma ocean whose total thickness is a few 100 km only, we consider the assumption of a constant $\alpha$ as sufficient. Therefore, Equation (17) can be integrated to give

$$ T = T^{\text{BOMB}} \exp \left[ \frac{\alpha g (r_{\text{BOMB}} - r)}{C_p} \right], \quad (18) $$

with $T^{\text{BOMB}}$ the temperature at the top of the BMO. This equation can be safely linearised if that simplifies the expressions. $T^{\text{BOMB}}$ is equal to the liquidus corresponding to the composition of the magma ocean, which evolves with time. The liquidus also depends on pressure and, for simplicity, we assume a linear dependence of the liquidus on both pressure and mass fraction of FeO, yielding

$$ T^{\text{BOMB}} = T_{\text{L}}(r_{\text{BOMB}}) = T_{\text{L}}(r_0) - \frac{\partial T_{\text{L}}}{\partial P} \rho g (r_{\text{BOMB}} - r_0) + \frac{\partial T_{\text{L}}}{\partial \xi} (\xi_1 - \xi_{l0}), \quad (19) $$

with $r_0$ the initial position of the freezing front, $\xi_{l0}$ the initial FeO mass fraction of the magma ocean and $T_{\text{L}}(r_0)$ the corresponding liquidus value. This equation can be time differentiated to give

$$ \frac{dT^{\text{BOMB}}}{dt} = -\frac{\partial T_{\text{L}}}{\partial P} \rho g \frac{dr_{\text{BOMB}}}{dt} + \frac{\partial T_{\text{L}}}{\partial \xi} \frac{d\xi_1}{dt}. \quad (20) $$

The two time derivatives on the right-hand-side can be related to each other using the equation for FeO conservation. This and other balance equations are derived in the next subsection.

The different profiles obtained above are only strictly valid in the well mixed isentropic bulk of the magma ocean and are complemented by boundary

layers on both top and bottom. However, as soon as the magma ocean is unstably stratified, which we assume here from the start, the Rayleigh number is enormous and convection is very efficient, which makes the super-isentropic temperature difference across the BMO very small [Labrosse et al., 2007, Ulvrová et al., 2012]. Therefore, we neglect the thickness of boundary layers and their associated temperature, composition and chemical potential jumps.

### 2.2.2. Balance equations

The mass fraction $\xi_1$ of FeO in the BMO is assumed to be uniform, owing to the high efficiency of convective stirring in the liquid. However, it evolves with time due to interaction with the solid mantle due to fractional crystallisation at the boundary, and possibly exchange by diffusion through the core mantle boundary. The equations describing the evolution of $\xi_1$ are derived for their introduction in StagYY. The integrated fluxes over the CMB and top surface of the BMO are denoted $I^{\text{CMB}}$ and $I^{\text{BOMB}}$, respectively, and are counted positive upward. The global balance equation for the FeO content is

$$ \frac{dM_{\text{BMO}}\xi_1}{dt} = I^{\text{CMB}} - I^{\text{BOMB}}, \quad (21) $$

$M_{\text{BMO}}$ being the mass of the BMO in the $\xi_1 = 0$ limit. For the time being, we will neglect $I^{\text{CMB}}$ since its treatment would require constraints on the partitioning between the liquid metal and the liquid silicate at the relevant pressure and temperature and a model for the vertical transfer of light elements at the top of the core. This could be a target for future developments.

With this assumption, Equation (21) can be developed to give

$$ M_{\text{BMO}} \frac{d\xi_1}{dt} + 4\pi r_{\text{BOMB}}^2 \rho_0 \xi_1 \dot{r}_{\text{BOMB}} = -I^{\text{BOMB}}. \quad (22) $$

The flux to the solid mantle, $I^{\text{BOMB}}$, is due to the phase change happening at the boundary, with a crystallisation mass rate $w$ associated with the solid radial velocity $u_r$ at the boundary as computed in StagYY, $u_r$, by $w = u_r - \dot{r}_{\text{BOMB}}$ since, in StagYY, the boundary is kept fixed by a continuous adjustment of the computation domain thickness (Section 2.1). We get

$$ I^{\text{BOMB}} = \int_{\Omega_{\text{BOMB}}} \rho_0 \xi_s (u_r - \dot{r}_{\text{BOMB}}) d\Omega, \quad (23) $$
with $\Omega_{BOMB}$ the surface of the boundary and $\xi_s$ the mass fraction of FeO on the solid side of the boundary:

$$ \xi_s = \begin{cases} D\xi_l, & u_r \ge 0 \\ \xi_s, & u_r < 0. \end{cases} \qquad (24) $$

In regions of solid upwelling, fractional crystallisation occurs, producing a solid with a mass fraction $\xi_s = D\xi_l$ with $D < 1$ the partition coefficient of FeO between the solid and the liquid. Note that $D$ depends on the composition and is related to the distribution coefficient $K$ of the phase change, as detailed by Bolrão et al. [2021]. In regions of downwelling, the solid with a mass fraction $\xi_s$ arrives in contact with the liquid and melts. Note that the composition of this solid may be different from that for equilibrium with the liquid at its liquidus temperature but melting can still proceed by pumping FeO from the liquid, which acts to reduce the mass fraction FeO in the liquid. This effect is balanced by the release of FeO in regions of crystallisation. The lateral transfer of FeO in the liquid occurs on a timescale, $\tau_\phi$, much shorter than the timescale for convection in the solid and the same timescale applies to the transfer of latent heat from regions of freezing to regions of melting. This is taken into account in the dimensionless number $\Phi = \tau_\phi/\tau_\eta$ which parameterises the boundary condition applied for convection in the solid, $\tau_\eta$ being the viscous timescale on which a topography is built as a result of viscous stress in the solid [Labrosse et al., 2018].

Combining Equations (22) and (23) gives, after rearrangement:

$$ \begin{aligned} & \left( 4\pi r_{\text{BOMB}}^2 \rho_0 \xi_l - \int_{\Omega_{\text{BOMB}}} \rho_0 \xi_s d\Omega \right) \dot{r}_{\text{BOMB}} \\ & \quad + \frac{4\pi}{3} \rho_0 (r_{\text{BOMB}}^3 - r_{\text{CMB}}^3) \dot{\xi}_l \\ & = - \int_{\Omega_{\text{BOMB}}} \rho_0 \xi_s u_r d\Omega. \end{aligned} \qquad (25) $$

The two integrals in Equation (25) are computed in StagYY at each time-step and that gives the relationship between $\dot{r}_{BOMB}$ and $\dot{\xi}_l$. Equation (20) then allows us to express the rate of change of $T^{BOMB}$ as a function of $\dot{r}_{BOMB}$ only.

Note that in the case $u_r = 0$, which is obtained for $\Phi = \infty$, $\xi_s = D\xi_l$ is uniform and the mass flux of FeO to the solid mantle is $I^{BOMB} = -4\pi r_{\text{BOMB}}^2 D\xi_l \rho_0 \dot{r}_{\text{BOMB}}$.

This flux is positive for $\dot{r}_{BOMB} < 0$, i.e. when crystallisation occurs. Equation (25) leads to

$$ \dot{\xi}_l = - \frac{3r_{\text{BOMB}}^2 \Delta\xi}{r_{\text{BOMB}}^3 - r_{\text{CMB}}^3} \dot{r}_{\text{BOMB}} \qquad (26) $$

with $\Delta\xi = \xi_l - \xi_s$. This equation was used in the original BMO paper [Labrosse et al., 2007] but the general Equation (25) accounting for flow through the boundary is used in the present study.

The flux $i$ of FeO across both boundaries contributes to the total heat flux density $q$ as

$$ q = -k\nabla T + \mu i \qquad (27) $$

and we note the thermal part as

$$ q_T \equiv -k\nabla T. \qquad (28) $$

Integrated over a surface $\Omega$ (BOMB or CMB), we get

$$ Q^{\Omega} = Q_{T}^{\Omega} + \mu^{\Omega} I^{\Omega} \qquad (29) $$

with $\mu^\Omega$ the average value of $\mu$ over the surface $\Omega$.

The long term thermal evolution of the BMO is controlled by the integrated energy balance equation, which is written as [e.g. Braginsky and Roberts, 1995, Lister and Buffett, 1995]

$$ \begin{aligned} & \int_{V_{\text{BMO}}} \rho \left( T \frac{\partial s_l}{\partial t} + \mu \frac{\partial \xi_l}{\partial t} \right) dV \\ & = \int_{\Omega_{\text{BOMB}}} \rho L (u_r - \dot{r}_{\text{BOMB}}) d\Omega \\ & \quad + Q_R + Q^{\text{CMB}} - Q^{\text{BOMB}}, \end{aligned} \qquad (30) $$

with $Q_R$ the radiogenic heat production in the BMO and $L$ the latent heat of melting, which contains two contributions:

$$ L = T_L \Delta s + \mu \Delta \xi \qquad (31) $$

$\Delta s = s_l - s_s$ being the entropy of melting. The surface integral on the right-hand-side of Equation (30) is the total latent released or consumed by phase change at the boundary and is related to the crystallisation change rate $u_r - \dot{r}_{BOMB}$. The entropy contribution to the latent heat, $T_L \Delta s$, is uniform along the phase change boundary and, since $\rho u_r$ averages to 0 at the interface, only the net crystallisation rate $\dot{r}_{BOMB}$ contributes to that part. Combining Equation (30) with the global balance in FeO Equation (21)
and the expressions for the BOMB and CMB fluxes Equation (29) gives

$$Q_T^{\text{BOMB}} = - \int_{V_{\text{BMO}}} \rho T \frac{\partial s}{\partial t} dV \\ - T_L(r_{\text{BOMB}}) \Delta s 4 \pi r_{\text{BOMB}}^2 \dot{r}_{\text{BOMB}} + Q_R + Q_T^{\text{CMB}} \\ + (\mu^{\text{CMB}} - \mu^{\text{BOMB}}) I^{\text{CMB}} - \frac{d\xi_1}{dt} \int_{V_{\text{BMO}}} \rho \mu' dV \\ + \int_{\Omega_{\text{BOMB}}} \rho \mu^{\text{BOMB}} (\Delta\xi u_r + \xi_s \dot{r}_{\text{BOMB}}) d\Omega. \quad (32)$$

This equation has a simple interpretation: the heat flow out of the BMO has contributions from secular cooling, latent heat of freezing (positive for $\dot{r}_{\text{BOMB}} < 0$), radiogenic heating, heat flow at the CMB, the flux of compositional energy due to the chemical flux at the CMB and the change of compositional energy content in the BMO. Note that the CMB flux term has a sign opposite to $I^{\text{CMB}}$, which is expected since mixing any dense component from the core would require energy.

This equation can be used to compute the thermal and compositional evolution of the BMO and the core for a given heat flow at the bottom of the solid mantle. Each term of the equation needs to be expressed as function of a minimum number of parameters, which is the next task here.

As stated above, we will neglect $I^{\text{CMB}}$, for now. $Q_T^{\text{BOMB}}$ is provided at each time step by the convection model of the solid mantle, StagYY. Computing $Q_R$ is straightforward for a given value of the partition coefficient of heat producing elements at the top of the BMO and their decay constants. The compositional energy term, the one involving $\mu'$ in Equation (32), is easily related to rate of change of $r_{\text{BOMB}}$ using Equations (16):

$$E_{\chi 1} \equiv -\frac{d\xi_1}{dt} \int_{V_{\text{BMO}}} \rho \mu' dV = -\rho_0 \beta g [r_{\text{BOMB}} V_{\text{BMO}} \\ -\pi (r_{\text{BOMB}}^4 - r_{\text{CMB}}^4)] \frac{d\xi_1}{dt}. \quad (33)$$

Using Equation (25), this term can be related to $\dot{r}_{\text{BOMB}}$.

The second compositional term (last on the right-hand-side of Equation (32)) can be simplified assuming $\mu^{\text{BOMB}}$, $\rho$ and $\xi_1$ to be constant and using the fact that $u_r$ is 0 on average on the boundary. This gives

$$E_{\chi 2} = \rho_0 \mu^{\text{BOMB}} \int_{\Omega_{\text{BOMB}}} \xi_s (\dot{r}_{\text{BOMB}} - u_r) d\Omega. \quad (34)$$

The secular cooling term is expressed as

$$Q_C = - \int_{V_{\text{BMO}}} \rho T \frac{\partial s}{\partial t} dV = - \int_{V_{\text{BMO}}} \rho C_p \frac{\partial T}{\partial t} dV, \quad (35)$$

which can be computed using the Equations (18) and (20):

$$Q_C = -4\pi\rho_0 C_p L_M^3 \left[ \frac{dT^{\text{BOMB}}}{dt} + \frac{\alpha g T^{\text{BOMB}}}{C_p} \dot{r}_{\text{BOMB}} \right] \\ \times \int_{r_{\text{CMB}}/L_M}^{r_{\text{BOMB}}/L_M} x^2 \exp \left[ \frac{\alpha g L_M}{C_p} \left( \frac{r_{\text{BOMB}}}{L_M} - x \right) \right] dx. \quad (36)$$

Using Equation (20) and (26), this term can be expressed as an affine function of $\dot{r}_{\text{BOMB}}$.

The last term to deal with is the CMB heat flow. The thermal evolution of the core can be parameterised using the CMB temperature [Labrosse, 2015] and we therefore have all the equations needed to solve the coupled evolution of the solid mantle, the BMO and the core. For a liquid core (only situation implemented for the moment), the core cooling term is related to the rate of change of the CMB temperature by

$$Q_T^{\text{CMB}} = -\frac{4\pi}{3} \rho_N C_{pN} L_{\rho}^3 \frac{f_C\left(\frac{r_{\text{CMB}}}{L_{\rho}}, \gamma\right)}{\left(1 - \frac{r_{\text{CMB}}^2}{L_{\rho}^2} - A_{\rho} \frac{r_{\text{CMB}}^4}{L_{\rho}^4}\right)^{\gamma}} \frac{dT^{\text{CMB}}}{dt}, \quad (37)$$

with $\rho_N$ the density at the center of the core, $C_{pN}$ the assumed constant heat capacity of the core, $L_\rho$ and $A_\rho$ the structure parameters describing the density variation with radius in the core, $\gamma$ the Grüneisen parameter of the core and

$$\begin{aligned} f_C(r, \gamma) &= 3 \int_0^r x^2 (1 - x^2 - A_\rho x^4)^{1+\gamma} dx \\ &= r^3 \left[ 1 - \frac{3}{5}(\gamma+1)r^2 - \frac{3}{14}(\gamma+1)(2A_\rho - \gamma)r^4 + O(r^6) \right]. \end{aligned} \quad (38)$$

The boundary layers on both sides of the CMB are tiny, with temperature differences across them that are negligible compared to the temperature differences across the whole layers due to compressibility effects. Therefore, the cooling of the core must follow that of the BMO. More precisely, the CMB temperature is related to that at the top of the BMO using Equation (18):

$$T^{\text{CMB}} = T^{\text{BOMB}} \exp \left[ \frac{\alpha g}{C_p} (r_{\text{BOMB}} - r_{\text{CMB}}) \right]. \quad (39)$$
Therefore, the rate of change of $T^{CMB}$ is

$$ \frac{dT^{CMB}}{dt} = \left[ \frac{dT^{BOMB}}{dt} + \frac{\alpha g T^{BOMB}}{C_p} \dot{r}_{BOMB} \right] \times \exp \left[ \frac{\alpha g}{C_p} (r_{BOMB} - r_{CMB}) \right]. \quad (40) $$

The equations for the evolution of the BMO and the core presented in this section have been made dimensionless (see Appendix) and implemented in StagYY. In practice, at each time step, the temperature and composition fields in the solid mantle are used to compute the body force responsible for convection. The velocity and pressure fields resulting from the momentum balance are computed, which permits to compute the evolution of the temperature and composition fields, using tracers for the latter. The heat flow at the bottom of the solid mantle is also used to compute the evolution of the BMO and core.

### 3. Example of dynamical evolution

The model is controlled by many input parameters, the classical ones for mantle convection using StagYY [Tackley, 2008] and the additional ones for the BMO evolution, so that the parameter space is effectively enormous. Although we have run many cases, the goal of the present paper is not to provide a comprehensive study of this complex system but to show with one example a possible evolution of the Earth with a basal magma ocean. As will appear clearly, the model results are in some ways encouraging, in the sense that we obtain some of the expected features, in particular thermo-chemical structures that might explain some of the current seismological observations. On the other hand, some of the outcomes point toward strong limitations of this model, at least with the choice of parameters of this specific calculation. Future developments to solve these issues will then be discussed. Even though all the complexities of mantle convection that are included in the numerical code StagYY are also accessible with the BMO model, we restrict ourselves here to the simplest case, notably with a constant viscosity and compositional variations only due to the fractional crystallisation of the BMO.

The parameters specific to the basal magma ocean and its evolution are as follows: the initial thickness is taken to be 30% of the total mantle, the partition coefficients for FeO and the heat producing element

are $K_{FeO} = 0.3$, $K_{HPE} = 10^{-2}$, respectively and the phase change parameter is $\Phi = 3 \times 10^{-2}$. The initial composition of the solid mantle and the BMO are uniform in both FeO and HPE and in equilibrium with each other according to these partition coefficients. The mean mass fraction of FeO is 0.1 while that for HPE makes the mean dimensionless internal heating rate equal to 5, which is somewhat smaller than the expected value for the bulk silicate Earth. The internal heating rate decays exponentially with a dimensionless half-life time of $10^{-2}$. The nominal Rayleigh and buoyancy numbers are $Ra = 3 \times 10^7$, $B = 5$, respectively. These are defined using, as a temperature scale, $\Delta T_m$, an estimate for the melting temperature difference between the top and the bottom of the mantle, so that the dimensionless temperature accross the solid $\Delta T/\Delta T_m$ is of order one. The chosen value is $\Delta T_m = 4000$ K. The compositional range implied by the definition of $B$ is 1, which corresponds to the density difference between the MgO and FeO end-members. Since the range of temperature and composition actually encountered in the model are different from these defining values, a rescaling needs to be performed after the calculations, as discussed below. Other parameters are detailed in the Appendix.

Figure 1 shows the evolution of several key fields with time, the temperature (top row), mass fraction of FeO (middle row) and HPE (bottom row). The time of each snapshot is written at the top of each column. The BMO is depicted by a uniform pink layer to visualise its shrinking by crystallisation. The total duration of this calculation, which lasts until nearly full crystallisation of the BMO, is about $t = 6.4 \times 10^{-3}$ (dimensionless), which is nearly half the dimensionless age of the Earth. This duration clearly depends on the value of several input parameters, and, in particular, is expected to decrease with an increase of the Rayleigh number, since this makes the heat flow increase, and decrease with the buoyancy number, since FeO buoyancy goes against thermal convection.

Convection first sets in at large scale (Figure 1), which is usual for convection with a phase change at the bottom boundary [Labrosse et al., 2018, Morrison et al., 2019, 2024]. Starting from a compositionally uniform mantle, variations of concentrations in FeO and HPEs soon develop from fractional crystallisation of the BMO. In the first three snapshots
<figure><img src="image_4.png" /><figcaption>**Figure 1.** Evolution of the temperature (top row), FeO content (middle row) and concentration in heat producing elements (bottom row) as function of time, as indicated on the top. Note that the colorbar is adapted to each snapshot for maximum lisibility.</figcaption></figure>

of Figure 1, variations in FeO content do not seem to affect the dynamics, which is controlled by cold plumes that reach the bottom boundary with the BMO and melt, as expected from previous studies of purely thermal convection with solid-liquid phase change [Agrusta et al., 2019, Lebec et al., 2023]. In the last two snapshots, FeO-enriched regions start forming dome-like structures from which hot plumes depart. This behaviour, that was predicted in the original BMO model [Labrosse et al., 2007], happens when the FeO concentration contrasts become sufficient to compete with the thermal buoyancy.

To better understand the evolution, it is useful to study the mean temperature and composition profiles, as shown on Figure 2. The radius is scaled with the total thickness of the mantle (solid and magma ocean) so that the figure shows the crystallisation of the magma ocean with time.

The FeO mass fraction profiles show a gradual increase with time and, more importantly for the dynamics, an increase of radial contrast. As the BMO concentration in FeO becomes larger, the solids that form from it become enriched and therefore denser, which leads to the partial stratification observed on Figure 1. Figure 3 shows the evolution of the minimum FeO mass fraction at the upper

boundary, $\min(\xi_{FeO,top})$, and the maximum at the bottom boundary, $\max(\xi_{FeO,bot})$. The minimum value is stable, expressing the fact that no diffusion is allowed in this model, while the maximum value, after a rapid increase in the early evolution, increases with an almost constant time rate. The Buoyancy number in this calculation is set to a nominal value of $B = 5$ but this corresponds to a change of composition between the MgO and FeO end-members and a dimensionless temperature change between 0 and 1. Since the model is evolving with time with fixed temperature and composition scales, the effective buoyancy number at each time is given by Equation (10). The value obtained is also shown on Figure 3. Comparison with the snapshots of Figure 1 suggests that a critical value of $B_{eff} \sim 0.5$ makes the convective regime transition between the well mixed situation to the stratified one. Determining the dependence of that number on other parameters of the problem, in particular the Rayleigh number, falls beyond the scope of the present paper. It is however grossly consistent with the value of the buoyancy number needed to get partial entrainment or a doming regime in convection with a compositionally layered initial condition [Tackley, 1998, Davaille, 1999, Le Bars and Davaille, 2004].
<figure><img src="image_5.png" /><figcaption>**Figure 2.** Temperature (left) and composition (right) mean profiles for the time steps corresponding to those of Figure 1. The radial position reflects the crystallisation of the BMO with time.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Time evolution of the minimum FeO mass fraction at the upper boundary, min($\xi_{\text{FeO,top}}$), the maximum at the bottom boundary, max($\xi_{\text{FeO,bot}}$) (left axis) and the effective buoyancy number $B_{\text{eff}}$, right axis. Empty circles on the black line give the position of the snapshots presented on Figure 1.</figcaption></figure>

The bottom temperature, equal to the liquidus, also varies with time because of the combined effect

of the pressure (depth) and composition of the magma. With the choice of parameters made for this calculation, the pressure effect dominates over that of FeO mass fraction, which implies an increase of the liquidus with time. At the end of the calculation, the mass fraction of FeO increases faster than pressure at the BMO boundary (BOMB) which makes the liquidus decrease with time. Other choice of parameters would change this evolution and a more realistic phase diagram could be implemented in the future.

A striking feature of the temperature profiles shown on Figure 2 is that, for the long early period where the solid stays well mixed, the mean temperature is nearly independent of the radius and varies significantly only in the upper boundary. Therefore, no boundary layer exists at the bottom and, consequently, no focused hot plumes. This behaviour has already been observed from both linear and non-linear models [Labrosse et al., 2018, Agrusta et al., 2019] in cartesian geometry and results from the phase change boundary condition. In this scenario, hot plumes can only develop once the basal magma ocean gets crystallised enough to get $\Phi \gtrsim 10^3$ or, as
in the present case, when a compositional stratification separates the phase change boundary from the bulk of the convecting mantle. It is worth recalling here that this model assumes an incompressible solid mantle within the Boussinesq approximation and the development with time of a significant temperature gradient at the bottom of the solid shell is due here solely to the stabilisation of dense thermochemical piles, not to any compressibility effect. To convert these temperature profiles to realistic ones for the Earth requires to add an isentropic temperature gradient.

The concentration in HPE, depicted on the bottom row of Figure 1, is somewhat similar to that of FeO since both come from fractional crystallisation of the BMO. However, because the partition coefficient is much smaller, the enrichment is slower and the mean concentration becomes significant only at the very end of the crystallisation process. Moreover, comparing the various FeO-enriched domes, we can see that their concentrations are not uniform: at the latest stages of their formation, small differences in crystallisation time result in large concentration variations, a phenomenon much less visible for the concentration in FeO which increases in a more gradual manner. Finally, the enrichment in HPEs of the domes makes them warm up and, in some calculations (not shown here), they become unstable at the end of the calculation.

## 4. Discussion and the way ahead

The model presented in the previous section is the first application of the complete set of equations coupling the evolution of a basal magma ocean to convection in the solid mantle. As such, it provides some interesting results but has some limitations when application to Earth are considered.

First of all, this model is a proof of concept for the possibility of forming large scale compositional variations in the lower mantle from fractional crystallisation of the BMO, as can be seen on Figure 1. The fact that this situation arises at the end of this calculation is linked to the shrinking of the BMO, which is necessary to make its FeO mass fraction increase sufficiently for solids that crystallise from it to be dense enough. Our calculation stops when the BMO crystallises almost entirely. In the present model, the BMO is imposed to fully freeze when a

minimum thickness is reached, while it is assumed to be fully liquid before. In reality, we expect that, at some point, a transition should occur toward a two-phase layer, the present interpretation of ULVZs as partially molten being the last remnant pockets of it. The present model does not include the physics necessary to deal with this transition. This two-phase situation could last quite a long time, the last melt getting enriched in all the incompatible elements, in particular volatiles, which significantly depress the freezing point of the magma [Nomura et al., 2014].

On the other hand, the thermal structure shown in Figures 1 and 2 may be difficult to reconcile with observations constraining the early Earth. Indeed, the absence of boundary layer at the bottom of the solid mantle implies that, except for thin downwelling plumes, the whole solid mantle in our model is at the same temperature as the bottom boundary. This situation, which results from the presence of phase change boundary condition at the bottom [Labrosse et al., 2018, Agrusta et al., 2019, Lebec et al., 2023], is difficult to sustain when applied to the mantle of the Earth since it would lead to very high temperature below the lithosphere, possibly making most of the upper mantle liquid. We expect that such a state would lead to massive eruptions and the associated rapid cooling would not allow it to persist for the length of time implied here [Davies, 1990, Sleep, 2000]. One way to avoid that problem is to increase the value of the phase change parameter, $\Phi$. Indeed, for $\Phi > 10^3$, the dynamics is similar to that obtained for a non-penetrative boundary condition [Agrusta et al., 2019] with a well-developed boundary layer on the solid side. However, the value of $\Phi$ is expected to be on the low side, similar to the value used here, although it is difficult to estimate it precisely.

In the present model, the situation with a high upper mantle temperature persists as long as the solid stays well mixed but once a compositional stratification develops, the bulk of the mantle can cool down compared to its bottom boundary (Figure 2). The main problem therefore comes from the initial condition of the calculation. The one shown here starts from a compositionally uniform solid mantle and the variations of composition only arise from fractional crystallisation of the BMO, since testing the feasibility of this mechanism is one of the goals of this paper. The freezing of the surficial magma ocean is, however, likely to produce a compositional stratification of the
solid that could be included as initial condition of our calculations. For example, models of fractional crystallisation of the lunar or martian magma oceans [e.g. Hess and Parmentier, 1995, Elkins-Tanton et al., 2003, Elkins-Tanton, 2012] suggest a gradual densification of the solid that make it prone to overturn and that should produce a solid mantle that is initially stably stratified. The degree of initial stratification can be rather extreme if the overturn happens after full crystallisation but several overturns are expected to occur during crystallisation [Ballmer et al., 2017, Maurice et al., 2017, Boukaré et al., 2018, Morison et al., 2019]. This initial stratification should therefore be considered as an input parameter whose effect should be quantified in future studies.

Alternatively, it could be desirable to run a model coupling the crystallisation of the surficial magma ocean to convection in the solid in order for the initial condition to result directly from this process. This means writing a set of equations similar to the one developed in Section 2.2 for the surficial magma ocean and including phase change boundary condition. Some progress has been donc in that direction [Morison, 2019] but a systematic exploitation of this model is still needed and will be the topic of a future study.

**5. Conclusions**

This paper presents, for the first time, a complete set of equations for the evolution of a basal magma ocean coupled to a convecting mantle and a cooling core. This theory has been used to build a model, starting from the mantle convection code StagYY [Tackley, 2008], for the evolution of the whole Earth that tracks the fractional crystallisation of the BMO and the implied compositional evolution of the mantle. The BMO is treated as a well-mixed layer whose evolution is controlled by general energy and composition balance equations. The solid mantle is treated with a full convection model that includes the possibility of melting and freezing at its bottom boundary as a phase change boundary condition [Labrosse et al., 2018].

The parameter space for this model is very large and, in this paper, we have only considered one typical case. It shows, for the first time in a self-consistent dynamical model, the feasibility of the scenario for the generation of large compositionally dense and

heterogenous anomalies proposed by Labrosse et al. [2007].

The phase change boundary condition at the bottom of the solid mantle has profound implications on its dynamics and thermal structure: as long as it is well mixed, its bulk temperature follows an isentropic upward extrapolation of the BMO liquidus. This would likely imply large amounts of melting in the upper mantle although, of course, a more thorough exploration of the parameter space is necessary to confirm that implication. In any case, that points toward the necessary discussion of the initial conditions, in particular concerning the compositional structure of the solid mantle, that results from the crystallisation of a surficial magma ocean. Models integrating that aspect in the evolution would help decipher this question.

Another interesting addition to this model would be to include a treatment of volatiles, and in particular water. Assuming water to be incompatible at high pressure, we can expect its concentration to increase in the basal magma ocean as it crystallises, which would help to maintain partial melt to the present day by decreasing the solidus [Nomura et al., 2014]. The solid forming from it would gradually see its H content increase, with important implication on its viscosity and possibly on the melting and degassing of the mantle near the surface.

**Declaration of interests**

The authors do not work for, advise, own shares in, or receive funds from any organization that could benefit from this article, and have declared no affiliations other than their research organizations.

**Acknowledgments**

We gratefully acknowledge two very constructive reviews that pushed us to improve and clarify the presentation of our paper. The calculation was performed on the cluster of PSMN (Pôle Scientifique de Modélisation Numérique) of the ENS de Lyon.

**Appendix A. Dimensionless numbers and parameter values for the example calculation**

We use dimensionless variables in the whole model, both for convection in the solid mantle and the
<table><caption>**Table 1.** Dimensionless parameters for the Basal Magma Ocean evolution</caption><thead><tr><th>Symbol</th><th>Expression</th><th>Description</th><th>Run value</th></tr></thead><tbody><tr><td>Γ<sub>0</sub><sup>-</sup></td><td><sup>R<sup>-</sup> - R<sub>c</sub></sup>&frasl;<sub>d</sub></td><td>Initial BMO thickness</td><td>0.3</td></tr><tr><td>T<sub>BOMB,0</sub></td><td><sup>T<sub>BOMB</sub>(t = 0)</sup>&frasl;<sub>ΔT</sub></td><td>Initial BOMB temperature</td><td>0.37</td></tr><tr><td>ξ<sub>0</sub></td><td>ξ(t = 0)</td><td>Initial FeO content</td><td>0.07</td></tr><tr><td>c<sub>HPE,0</sub></td><td>c<sub>HPE</sub>(t = 0)</td><td>Initial HPE content</td><td>5 × 10<sup>-3</sup></td></tr><tr><td>B</td><td>-<sup>β</sup>&frasl;<sub>αΔT</sub></td><td>Compositional buoyancy</td><td>5</td></tr><tr><td>Di</td><td><sup>αgL<sub>M</sub></sup>&frasl;<sub>C<sub>p</sub></sub></td><td>BMO dissipation number</td><td>1.74</td></tr><tr><td>T<sub>LP</sub></td><td><sup>ρgL<sub>M</sub></sup>&frasl;<sub>ΔT<sub>M</sub></sub>∂T<sub>L</sub></td><td>Pressure dependence of melting temperature</td><td>Depth-dependent from Fiquet et al. [2010]</td></tr><tr><td>T<sub>Lξ</sub></td><td><sup>1</sup>&frasl;<sub>ΔT<sub>M</sub></sub>∂T<sub>L</sub>/∂ξ</td><td>Composition dependence of melting temperature</td><td>-0.125</td></tr><tr><td>St</td><td><sup>Δs</sup>&frasl;<sub>C<sub>p</sub></sub></td><td>Stefan number</td><td>0.3</td></tr><tr><td>τ<sub>hpe</sub></td><td><sup>τ<sub>hpe,dK</sub></sup>&frasl;<sub>L<sub>M</sub><sup>2</sup></sub></td><td>HPE half-life</td><td>0.02</td></tr><tr><td>P<sub>c</sub></td><td><sup>P<sub>c,d</sub></sup>&frasl;<sub>ρC<sub>p</sub>L<sub>M</sub><sup>3</sup></sub></td><td>Core cooling contribution</td><td>3.67</td></tr></tbody></table>

evolution of the BMO and core. The solid mantle part is explained in Section 2.1. The value expected for the phase change number $\Phi$ can be estimated using its expression in Equation (12) using reasonable values for the various parameters: $\rho_s \approx 5 \times 10^3$ kg·m$^{-3}$, $\rho_l \approx 5.5 \times 10^3$ kg·m$^{-3}$, $C_{pl} \approx 10^3$ J·K$^{-1}$·kg$^{-1}$, $\eta_s \approx 10^{22}$ Pa·s, $d \approx 3000$ km. For the latent heat, we use $L = T\Delta S$, with $\Delta S = 652$ J·kg$^{-1}$·K$^{-1}$ the entropy of melting [Stixrude et al., 2009] and $T \approx 4000$ K the melting temperature. The slopes of the isentropic profile and melting temperatures are estimated from Stixrude et al. [2009]: $m_{ad} \approx 15$ K·GPa$^{-1}$, $m_p \approx 4$ K·GPa$^{-1}$. The least well constrained quantity is the flow velocity in the BMO and we tentatively set it to a value similar to that in the present-day core, $u_l \sim 10^{-4}$ m·s$^{-1}$. This gives $\Phi \sim 4 \times 10^{-8}$. For convection with a phase change boundary condition at one of the horizontal boundaries, the asymptotic limit of very small values of the phase change number is reached for $\Phi \lesssim 10^{-1}$ [Labrosse et al., 2018, Agrusta et al., 2019, Morison et al., 2024]. The value just estimated for the BOMB

puts it safely in that regime, irrespective of any possible overestimation of the flow velocity in the BMO. We therefore set this value to $3 \times 10^{-2}$, which is sufficiently low to get the proper behaviour.

The magma ocean evolution equations are made dimensionless with the following scales:

* $d$, the total thickness of the mantle (solid and magma ocean), considered constant in this study for simplicity. For the Earth, $d = 2900$ km.

* $\Delta T$, an arbitrary reference temperature difference. In practice it is taken as an estimate for the melting temperature difference between the top and the bottom of the mantle, so that the dimensionless temperature accross the solid is of order one. The chosen value is $\Delta T = 4000$ K.

* $\rho$, the reference density of the considered materials. For the Earth's mantle, $\rho = 4000$ kg/m$^3$.
• $\kappa = k/(\rho C_p)$, the diffusion coefficient $\kappa = 10^{-6}$ m$^2$/s, associated with the conductivity $k = 4$ W/(mK) and the heat capacity $C_p = 10^3$ J·K$^{-1}$·kg$^{-1}$.

The dimensionless parameters arising from the non-dimensionalization are shown in Table 1. More detailed can be found in Morison [2019].

## References

Agrusta, R., Morison, A., Labrosse, S., Deguen, R., Alboussière, T., Tackley, P. J., and Dubuffet, F. (2019). Mantle convection interacting with magma oceans. *Geophys. J. Int.*, 220, 1878–1892.

Alboussière, T., Deguen, R., and Melzani, M. (2010). Melting-induced stratification above the Earth’s inner core due to convective translation. *Nature*, 466, 744–747.

Anderson, O. L., Oda, H., and Isaak, D. (1992). A model for the computation of thermal expansivity at high compression and high temperatures: MgO as an example. *Geophys. Res. Lett.*, 19(19), 1987–1990.

Ballmer, M. D., Lourenço, D. L., Hirose, K., Caracas, R., and Nomura, R. (2017). Reconciling magma-ocean crystallization models with the present-day structure of the Earth’s mantle. *Geochem. Geophys. Geosyst.*, 18(7), 2785–2806.

Bolrão, D. P., Ballmer, M. D., Morison, A., Rozel, A. B., Sanan, P., Labrosse, S., and Tackley, P. J. (2021). Timescales of chemical equilibrium between the convecting solid mantle and over/underlying magma oceans. *Solid Earth*, 12, 421–437.

Boukaré, C.-E., Parmentier, E., and Parman, S. (2018). Timing of mantle overturn during magma ocean solidification. *Earth Planet. Sci. Lett.*, 491, 216–225.

Braginsky, S. I. and Roberts, P. H. (1995). Equations governing convection in Earth’s core and the geodynamo. *Geophys. Astrophys. Fluid Dyn.*, 79, 1–97.

Chopelas, A. and Boehler, R. (1992). Thermal expansivity in the lower mantle. *Geophys. Res. Lett.*, 19(19), 1983–1986.

Christensen, U. R. (1985). Thermal evolution models for the Earth. *J. Geophys. Res.*, 90, 2995–3007.

Christensen, U. R. and Hofmann, A. W. (1994). Segregation of subducted oceanic crust in the convecting mantle. *J. Geophys. Res.*, 99(B10), 19867–19884.

Coltice, N., Husson, L., Faccenna, C., and Arnould, M. (2019). What drives tectonic plates? *Sci. Adv.*, 5(10), article no. eaax4295.

Davaille, A. (1999). Simultaneous generation of hotspots and superswells by convection in a heterogeneous planetary mantle. *Nature*, 402, 756–760.

Davies, G. F. (1990). Heat and mass transport in the early earth. In Newsome, H. E. and Jones, J. H., editors, *Origin of the Earth*, pages 175–194. Oxford University Press, New York.

Deguen, R. (2013). Thermal convection in a spherical shell with melting/freezing at either or both of its boundaries. *J. Earth Sci.*, 24, 669–682.

Deguen, R., Alboussière, T., and Cardin, P. (2013). Thermal convection in Earth’s inner core with phase change at its boundary. *Geophys. J. Int.*, 194, 1310–1334.

Driscoll, P. and Bercovici, D. (2014). On the thermal and magnetic histories of Earth and Venus: Influences of melting, radioactivity, and conductivity. *Phys. Earth Planet. Inter.*, 236, 36–51.

Duffy, T. S. and Ahrens, T. J. (1993). Thermal expansion of mantle and core materials at very high pressures. *Geophys. Res. Lett.*, 20(11), 1103–1106.

Elkins-Tanton, L. T. (2012). Magma oceans in the inner solar system. *Ann. Rev. Earth Planet. Sci.*, 40, 113–139.

Elkins-Tanton, L. T., Parmentier, E. M., and Hess, P. C. (2003). Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: Implications for Mars. *Meteorit. Planet. Sci.*, 38(12), 1753–1771.

Fichtner, A., Kennett, B. L. N., Tsai, V. C., et al. (2024). Seismic tomography 2024. *Bull. Seism. Soc. Am.*, 114, 1185–1213.

Fiquet, G., Auzende, A. L., Siebert, J., Corgne, A., Bureau, H., Ozawa, H., and Garbarino, G. (2010). Melting of peridotite to 140 GigaPascals. *Science*, 329, 1516–1518.

Flament, N., Bodur, Ö. F., Williams, S. E., and Merdith, A. S. (2022). Assembly of the basal mantle structure beneath Africa. *Nature*, 603(7903), 846–851.

Gerya, T. V. and Yuen, D. A. (2003). Characteristics-based marker-in-cell method with conservative finite-differences schemes for modeling geological flows with strongly variable transport properties. *Phys. Earth Planet. Inter.*, 140, 293–318.

Hernlund, J. and McNamara, A. (2015). 7.11 – the
core-mantle boundary region. In Schubert, G., editor, *Treatise on Geophysics. Mantle Dynamics*, pages 461–519. Elsevier, Oxford, 2nd edition.

Hernlund, J. W. and Tackley, P. J. (2008). Modeling mantle convection in the spherical annulus. *Phys. Earth Planet. Inter.*, 171(1–4), 48–54.

Hess, P. C. and Parmentier, E. M. (1995). A model for the thermal and chemical evolution of the Moon's interior: implications for the onset of mare volcanism. *Earth Planet. Sci. Lett.*, 134, 501–514.

Ismail-Zadeh, A. and Tackley, P. J. (2012). *Computational Methods for Geodynamics*. Cambridge University Press, Cambridge.

Jaupart, C., Labrosse, S., Lucazeau, F., and Mareschal, J.-C. (2015). 7.06–temperatures, heat, and energy in the mantle of the Earth. In Schubert, G. and Bercovici, D., editors, *Treatise on Geophysics. Mantle Dynamics*, pages 223–270. Elsevier, Oxford, 2nd edition.

Labrosse, S. (2015). Thermal evolution of the core with a high thermal conductivity. *Phys. Earth Planet. Inter.*, 247, 36–55.

Labrosse, S. (2016). Thermal state and evolution of the Earth core and deep mantle. In Terasaki, H. and Fischer, R., editors, *Deep Earth: Physics and Chemistry of the Lower Mantle and Core*, AGU Geophysical Monograph Series, pages 43–56. American Geophysical Union, Washington DC.

Labrosse, S., Hernlund, J. W., and Coltice, N. (2007). A crystallizing dense magma ocean at the base of Earth's mantle. *Nature*, 450, 866–869.

Labrosse, S., Morison, A., Deguen, R., and Albousière, T. (2018). Rayleigh-Bénard convection in a creeping solid with a phase change at either or both horizontal boundaries. *J. Fluid Mech.*, 846, 5–36.

Laneuville, M., Hernlund, J., Labrosse, S., and Gutenberg, N. (2018). Crystallization of a compositionally stratified basal magma ocean. *Phys. Earth Planet. Inter.*, 276, 86–92.

Le Bars, M. and Davaille, A. (2004). Whole layer convection in a heterogeneous planetary mantle. *J. Geophys. Res.*, 109, article no. B03403.

Lebec, L., Labrosse, S., Morison, A., and Tackley, P. J. (2023). Scaling of convection in high-pressure ice layers of large icy moons and implications for habitability. *Icarus*, 396, article no. 115494.

Lebec, L., Labrosse, S., Morison, A., and Tackley, P. J. (2024). Effects of salts on the exchanges through

high-pressure ice layers of large ocean worlds. *Icarus*, 412, article no. 115966.

Li, M. and McNamara, A. K. (2013). The difficulty for subducted oceanic crust to accumulate at the Earth's core-mantle boundary. *J. Geophys. Res.*, 118(4), 1807–1816.

Lister, J. R. and Buffett, B. A. (1995). The strength and efficiency of the thermal and compositional convection in the geodynamo. *Phys. Earth Planet. Inter.*, 91, 17–30.

Maurice, M., Tosi, N., Samuel, H., Plesa, A.-C., Hüttig, C., and Breuer, D. (2017). Onset of solid-state mantle convection and mixing during magma ocean solidification. *J. Geophys. Res.: Planets*, 122(3), 577–598.

Mizzon, H. and Monnereau, M. (2013). Implication of the lopsided growth for the viscosity of Earth's inner core. *Earth Planet. Sci. Lett.*, 361, 391–401.

Morison, A. (2019). *Convection dans le manteau primitif en interaction avec des océans de magma globaux*. PhD thesis, Université de Lyon, ENS de Lyon.

Morison, A., Labrosse, S., Deguen, R., and Alboussière, T. (2019). Timescale of overturn in a magma ocean cumulate. *Earth Planet. Sci. Lett.*, 516, 25–36.

Morison, A., Labrosse, S., Deguen, R., and Alboussière, T. (2024). Onset of thermal convection in a solid spherical shell with melting at either or both boundaries. *Geophys. J. Int.*, 238, 1121–1136.

Nakagawa, T. and Tackley, P. J. (2012). Influence of magmatism on mantle cooling, surface heat flow and urey ratio. *Earth Planet. Sci. Lett.*, 329–330, 1–10.

Nakagawa, T., Tackley, P. J., Deschamps, F., and Conolly, J. A. D. (2010). The influence of MORB and harzburgite composition on thermo-chemical mantle convection in a 3-D spherical shell with self-consistently calculated mineral physics. *Earth Planet. Sci. Lett.*, 296(3–4), 403–412.

Nomura, R., Hirose, K., Uesugi, K., Ohishi, Y., Tsuchiyama, A., Miyake, A., and Ueno, Y. (2014). Low core-mantle boundary temperature inferred from the solidus of pyrolite. *Science*, 343(6170), 522–525.

Patočka, V., Šrámek, O., and Tosi, N. (2020). Minimum heat flow from the core and thermal evolution of the Earth. *Phys. Earth Planet. Inter.*, 305, article no. 106457.
Ricard, Y. (2007). Physics of mantle convection. In Bercovici, D. and Schubert, G., editors, *Treatise of Geophysics. Mantle Dynamics*, volume 7, pages 253–303. Elsevier, Oxford.

Ricard, Y., Alboussière, T., Labrosse, S., Curbelo, J., and Dubuffet, F. (2022). Fully compressible convection for planetary mantles. *Geophys. J. Int.*, 230(2), 932–956.

Rost, S., Garnero, E. J., Williams, Q., and Manga, M. (2005). Seismological constraints on a possible plume root at the core-mantle boundary. *Nature*, 435, 666–669.

Sleep, N. H. (2000). Evolution of the mode of convection within terrestrial planets. *J. Geophys. Res.*, 105, 17563–17578.

Stixrude, L., de Koker, N., Sun, N., Mookherjee, M., and Karki, B. B. (2009). Thermodynamics of silicate liquids in the deep Earth. *Earth Planet. Sci. Lett.*, 278, 226–232.

Tackley, P. J. (1998). Three-dimensional simulations of mantle convection with a thermochemical CMB boundary layer: D''? In Gurnis, M., Wyession, M. E.,

Knittle, E., and Buffett, B. A., editors, *The Core-Mantle Boundary Region*, pages 231–253. American Geophysical Union, Washington DC.

Tackley, P. J. (2008). Modelling compressible mantle convection with large viscosity contrasts in a three-dimensional spherical shell using the yin-yang grid. *Phys. Earth Planet. Inter.*, 171(1–4), 7–18.

Tackley, P. J. and King, S. D. (2003). Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations. *Geochem. Geophys. Geosyst.*, 4(4), article no. 8302.

Ulvrová, M., Labrosse, S., Coltice, N., Raback, P., and Tackley, P. J. (2012). Numerical modeling of convection interacting with a melting and solidification front: application to the thermal evolution of the basal magma ocean. *Phys. Earth Planet. Inter.*, 206–207, 51–66.

Williams, Q. and Garnero, E. J. (1996). Seismic evidence for partial melt at the base of earth's mantle. *Science*, 273(5281), 1528–1530.