\makeatletter
\@ifundefined{HCode}
{\documentclass[screen,CRGEOS,Unicode,Thematic]{cedram}
\usepackage[square,comma,authoryear]{natbib}
\renewcommand\bibsection{\section*{\refname}}
\newenvironment{noXML}{}{}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tbody{}%\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.92\linewidth}{#1}}
\def\xxtabnote#1#2{\vskip4pt\parbox{#1\linewidth}{#2}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\ndash{\text{--}}
\def\og{\guillemotleft}
\def\fg{\guillemotright}
\RequirePackage{etoolbox}
\usepackage{amssymb}
\def\jobid{crgeos20240413}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\usepackage{multirow} 
\usepackage{wasysym} 
\def\nrow#1{\@tempcnta #1\relax%
\advance\@tempcnta by 1\relax%
\xdef\lenrow{\the\@tempcnta}}
\def\morerows#1#2{\nrow{#1}\multirow{\lenrow}{*}{#2}}
\def\xmorerows#1#2{#2}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\sout#1{}
\csdef{Seqnsplit}{\\}
\def\botline{\\\hline}
\def\back#1{}
\def\rmmu{\upmu}
\def\rmpi{\uppi}
\def\rmDelta{\Delta}
\def\rmOmega{\Omega}
\def\0{\phantom{0}}
\def\mn{\phantom{$-$}}
\DOI{10.5802/crgeos.275}
\datereceived{2024-05-21}
\daterevised{2024-10-03}
\dateaccepted{2024-10-07}
\ItHasTeXPublished
\def\mathbi#1{\text{\textit{\textbf{#1}}}}
\def\Lbreak{\newline}
\def\ubreak{\break}
\makeatletter
\g@addto@macro{\UrlBreaks}{\UrlOrds}
\gappto{\UrlBreaks}{\UrlOrds}
\usepackage{hyperref}
\makeatletter
\newcommand*{\hyperlinkcite}[1]{\hyper@link{cite}{cite.#1}}
\makeatother
}
{
\PassOptionsToPackage{authoryear}{natbib}
\documentclass[crgeos]{article}
\usepackage[T1]{fontenc}
\def\CDRdoi{10.5802/crgeos.275}
\def\hyperlink#1#2{#2}
\def\hypertarget#1#2{#2}
\def\xxtabnote#1#2{\tabnote{#2}}
\def\xmorerows#1#2{\morerows{#1}{#2}}
\let\ubreak\relax
\def\sout#1{}
\makeatletter
\def\CDRsupplementaryTwotypes#1#2{}
\def\textperthousand{\permil}
\def\hyperlinkcite#1#2{\citeyear{#1}}
}
\makeatother

\usepackage{upgreek}

\begin{document}

\begin{noXML}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Magma degassing and its impact on the Earth's atmosphere: from magma oceans to lava lakes}
%\def\cdr@specialtitle@french{Impact atmosph\'erique  du d\'egazage magmatique : des oc\'eans de magma aux  lacs de lave}
%\makeatother

\CDRsetmeta{articletype}{research-article}

\title{Solid-state mantle convection coupled with a crystallising basal
magma ocean}

\alttitle{Convection dans le manteau solide coupl\'ee \`a la
cristallisation d'un oc\'ean de magma basal}

\author{\firstname{St\'ephane} \lastname{Labrosse}\CDRorcid{0000-0002-7620-4363}\IsCorresp}
\address{LGLTPE, ENS de Lyon, Universit\'e de Lyon, 46 all\'ee d'Italie, 69003 Lyon, France}
\email[S. Labrosse]{stephane.labrosse@ens-lyon.fr}

\author{\firstname{Adrien} \lastname{Morison}\CDRorcid{0000-0002-8261-9143}}
\address{Univ Exeter, Phys \& Astron, Exeter, Devon, England}
\email[A. Morison]{A.Morison@exeter.ac.uk}

\author{\firstname{Paul James} \lastname{Tackley}\CDRorcid{0000-0003-4878-621X}}
\address{Department of Earth and Planetary Sciences, ETH Z\"urich, Sonneggstrasse 5, Z\"urich, 8092, Switzerland}
\email[P. J. Tackley]{paul.tackley@eaps.ethz.ch}

\keywords{\kwd{Thermal evolution}
\kwd{Mantle convection}
\kwd{Core cooling}
\kwd{Basal magma ocean}}

\altkeywords{\kwd{\'Evolution thermique}
\kwd{Convection mantellique}
\kwd{Refroidissement du noyau}
\kwd{Oc\'ean de magma basal}}

\shortrunauthors

\begin{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. 
\end{abstract} 

\begin{altabstract}
La cristallisation fractionn\'ee d'un oc\'ean de magma basal (OMB) a
\'et\'e propos\'ee pour expliquer la formation de variations de
composition \`a  grande \'echelle dans le manteau et la persistance de
zones partiellement fondues \`a la base du manteau inf\'erieur. Nous
pr\'esentons un ensemble complet d'\'equations pour l'\'evolution
thermique et compositionnelle de l'OMB et montrons qu'il peut \^etre
impl\'ement\'e dans un code de convection du manteau pour r\'esoudre le
probl\`eme de l'\'evolution du manteau \`a long terme. La pr\'esence de
l'OMB modifie la dynamique du manteau de plusieurs fa\c{c}ons.
L'\'equilibre de phase \`a la base du manteau solide implique un
changement des conditions limites m\'ecaniques, ce qui favorise la
convection \`a l'\'etat solide. La cristallisation nette de l'OMB
implique un changement de domaine de calcul, qui est trait\'e en
utilisant un changement de variable pour la coordonn\'ee radiale pour
garder une \'epaisseur constante. La fusion et la cristallisation
fractionn\'es \`a la fronti\`ere font \'evoluer la composition de l'OMB
et du manteau solide, ce qui est trait\'e \`a l'aide de traceurs
lagrangiens. Un exemple de calcul montre que la persistance de l'OMB
et son \'evolution \`a long terme modifient radicalement la dynamique du
manteau solide en favorisant les courants descendant focalis\'es et les
\'ecoulements \`a grande \'echelle. L'augmentation progressive de la
teneur en FeO dans l'OMB et dans le solide qui cristallise \`a partir
de celui-ci conduit \`a la stabilisation de piles
thermo-compositionnelles \`a grande \'echelle en base de manteau, ce qui
pourrait expliquer les observations sismologiques. 
\end{altabstract} 
 
\maketitle

\twocolumngrid

\end{noXML}

\newcommand{\bbR}{\mathbb{R}}
\newcommand{\bbC}{\mathbb{C}}
\newcommand{\bbZ}{\mathbb{Z}}
\newcommand{\avg}[1]{\overline{#1}}
\newcommand{\vect}[1]{\boldsymbol{#1}}
\newcommand{\vvect}[1]{\mathbi{#1}}
\newcommand\Ray{{\mathit{Ra}}} % Rayleigh number
\newcommand\Di{\mathit{Di}} % Dissipation number
\def\dpd#1#2{\displaystyle\frac{\partial #1}{\partial#2}}
\def\dod#1#2{\displaystyle\frac{\mathrm{d} #1}{\mathrm{d}#2}}
\def\dif#1{\mathrm{d}#1}

\section{Introduction}\label{sec: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 \mbox{seismic} tomography
over the last few decades \citep[e.g.][]{Fichtner_etal2024}, and is generally well understood in the
context of mantle convection in a plate-tectonic regime
\citep[e.g.][]{coltice_etal2019}. Paleogeographic reconstructions have
been used to constrain models of mantle convection going back 1 Gyr
\citep{Flament_etal2022} 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 \mbox{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)
\citep[e.g.][for\break a review]{Hernlund_McNamara2015}.  

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
\citep[e.g.][]{LeBars_Davaille04b}. 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
\citep{Christensen_Hofmann94,Nakagawa_etal2010,Li_McNamara2013}. 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
\citep{Labrosse_etal07}. 

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
\citep{Rost_etal05}, 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 \citep{Williams_Garnero1996}. 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 \citep{Jaupart_etal2015} and the
possibility of the core to be cooling faster than the mantle has been
invoked \citep{Driscoll_etal2014,Labrosse2016,Patocka2020} to solve the
long standing thermal catastrophe problem for the thermal evolution of
the Earth \citep{Christensen85a}. The rapid core cooling is also a
consequence of the necessity to maintain a convective dynamo with a
large thermal conductivity of the core
\citep{Labrosse2015,Patocka2020}. All these arguments together led to
the scenario of a basal magma ocean \citep{Labrosse_etal07}.

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
\citep{Tackley2008}, 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
\citep{Labrosse_etal2018,Agrusta_etal2019,Bolrao_etal2021,Lebec_etal2023,Lebec_etal2024}.
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
\citep{Alboussiere_etal2010,Mizzon_Monnereau2013,Deguen_etal2013,Deguen2013,Labrosse_etal2018}
controlled by a single dimensionless parameter, the phase change number
$\Phi$. For small values of this parameter ($\Phi\lesssim 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 ($\Phi\gtrsim10^3$). 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~\ref{sec:conv-solid-mantle}. 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~\ref{sec:evolution-bmo}. An
example calculation is then presented in Section~\ref{sec:example-dynam-evol}
before discussing the implications and limitations of the model in
Section~\ref{sec:limit-way-ahead}.}

\section{Physical and numerical model}
\label{sec:model}

\subsection{Convection in the solid mantle}
\label{sec:conv-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 \citep{Tackley2008}. 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 \citep{Tackley_King03,Gerya_Yuen03,IsmailZadeh_Tackley2010}.
In the present study, we use the spherical annulus geometry
\citep{Hernlund_Tackley2008} 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 $\Ray$, the internal heating rate $H$
and the buoyancy number $B$, their usual definition being,
{\begin{equation}
  \label{eq:dimensionless-numbers}
  \Ray = \frac{\alpha g \rmDelta T d^3}{\kappa\nu}, \quad
  H = \frac{\rho h d^2}{k \rmDelta T}, \quad
  B = \frac{\beta}{\alpha\rmDelta T},
\end{equation}}\unskip
with $\alpha$ and $\beta$ the thermal and compositional (FeO) expansion
coefficients, respectively, $g$ the
gravity, $\kappa$ the thermal diffusivity, $\nu$ the kinematic
viscosity, $\rmDelta T$ the temperature scale, $d$ the thickness of
the whole mantle (solid and BMO), $\rho$ the density and $h$ the time-evolving
radioactive heating rate. Among the parameters entering the definition
of the dimensionless numbers, $\rmDelta 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 \hyperlink{AppA}{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_{\mathrm{BOMB}}(t)$, and
the Earth surface with dimensionless radius $R_{\mathrm{E}}$. We can map the
time-varying spherical shell to a constant one, between rescaled
radii 1 and 2, by using as radial coordinate
{\begin{equation}
  \label{eq:scaled-rad}
  z = 1 + \frac{r - r_{\mathrm{BOMB}}}{R_{\mathrm{E}}-r_{\mathrm{BOMB}}}.
\end{equation}}\unskip
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)=\tilde{f}(z,t)$, we get
{\begin{equation}
  \label{eq:radial-der}
  \dpd{f}{r} =
  \dpd{z}{r}\dpd{\tilde{f}}{z}
  = \frac{1}{R_{\mathrm{E}}-r_{\mathrm{BOMB}}}\dpd{\tilde{f}}{z},
\end{equation}}\unskip
and so on for higher order derivatives. Consider now the time
derivative,
{\begin{equation}
  \label{eq:time-der}
  \dpd{f}{t} = \dpd{\tilde{f}}{t} + \dpd{z}{t}\dpd{\tilde{f}}{z} =
  \dpd{\tilde{f}}{t} -
  \dot{r}_{\mathrm{BOMB}}\frac{R_{\mathrm{E}}-r}{(R_{\mathrm{E}}-r_{\mathrm{BOMB}})^2}\dpd{\tilde{f}}{z},
\end{equation}}\unskip
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, $\Gamma = R_{\mathrm{E}} - r_{\mathrm{BOMB}}$,
the dimensionless balance equations for mass, momentum, energy, FeO
mass fraction $\xi$ and concentration $c_{\mathrm{HPE}}$ of heat producing
elements (HPE) become
{\begin{eqnarray}
  \label{eq:mass}
  0 & =& \vect{\nabla}\cdot\vvect{u},\Seqnsplit
  \label{eq:momentum}
  0 & =& -\vect{\nabla}p + \nabla^2 \vvect{u}\nonumber\\
&&  +\, \Ray\, \Gamma^3
      \left[\left(T - B\xi\right) -
      \left\langle T - B\xi\right\rangle\right]
      \vvect{e}_z,\Seqnsplit
  \Gamma^2\dpd{T}{t} & =& - \vvect{u}\cdot\vect{\nabla}T 
- \dot{\Gamma}\Gamma (2-z) \dpd{T}{z}
+ \nabla^2 T\nonumber\\
&&+\, H c_{\mathrm{HPE}}\Gamma^2
\exp\left(-\frac{t \ln{2}}{\tau_{\mathrm{HPE}}}\right),\Seqnsplit
\Gamma^2 \dpd{\xi}{t} & =&  - \vvect{u}\cdot\vect{\nabla}\xi ,\Seqnsplit
\Gamma^2 \dpd{c_{\mathrm{HPE}}}{t} & =& - \vvect{u}\cdot\vect{\nabla}c_{\mathrm{HPE}},
\end{eqnarray}}\unskip
with $\vvect{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 $\tau_{\mathrm{HPE}}$
that represents a mean of the four main ones in the Earth, $^{235}$U,
$^{238}$U, $^{232}$Th and $^{40}$K. The time derivative
$\dot{\Gamma}_{\mathrm{BOMB}}$ 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_{\mathrm{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
{\begin{equation}
  \label{eq:numbers_eff}
  \Ray_{\mathrm{eff}}=\Ray\,\Gamma^3 T_{\mathrm{BOMB}}, 
  \qquad B_{\mathrm{eff}} = B \frac{(\langle\xi_{\mathrm{bot}}
  \rangle-\langle\xi_{\mathrm{top}}\rangle)}{T_{\mathrm{BOMB}}},
\end{equation}}\unskip
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
\citep{Agrusta_etal2019,Bolrao_etal2021,Lebec_etal2023,
Lebec_etal2024}, is the solid--liquid phase change boundary condition
at the bottom of the solid shell \citep[e.g.][]{Labrosse_etal2018},
{\begin{equation}
  \label{eq:phase-chage-BC}
  2 \dpd{u_{\mathrm{r}}}{r} - p - \Phi u_{\mathrm{r}} = 0,
\end{equation}}\unskip
$u_{\mathrm{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 \citep{Deguen_etal2013},
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_{\mathrm{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 \citep[i.e.\ the critical Rayleigh number for
the onset of convection is reduced, see][]{Deguen2013,Labrosse_etal2018,Morison_etal2024}, heat and mass transfer are
increased and the wavelength of convection is increased compared to
the situation usually considered in mantle convection
\citep{Agrusta_etal2019}.
For a purely thermal problem \citep{Deguen_etal2013}, the
phase change number can be expressed as 
{\begin{equation}
  \label{eq:Phi}
  \Phi = \frac{\rho_{\mathrm{s}} (\rho_{\mathrm{l}}-\rho_{\mathrm{s}})
  L d}{\rho_{\mathrm{l}}^2C_{\mathrm{pl}} (m_{\mathrm{ad}}-m_p) u_{\mathrm{l}}
  \eta_{\mathrm{s}}},
\end{equation}}\unskip
with $\rho_{\mathrm{s}}$ and $\rho_{\mathrm{l}}$ the density of the solid and the
liquid, respectively, $L$ the latent heat of fusion, $C_{\mathrm{pl}}$ the
heat capacity of the liquid, $u_{\mathrm{l}}$ the typical flow velocity in the
magma ocean, $\eta_{\mathrm{s}}$ the viscosity of the solid and
$(m_{\mathrm{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\sim10^{-8}$ (see \hyperlink{AppA}{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
\citet{Agrusta_etal2019}: 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:
{\begin{equation}
  \label{eq:Robin}
  \mathcal{H}_{\mathrm{s}}(u_{\mathrm{r}}+u_0) \theta + [1 - \mathcal{H}_{\mathrm{s}}(u_{\mathrm{r}}+u_0)] 
  \dpd{\theta}{r} = 0,
\end{equation}}\unskip
with $\mathcal{H}_{\mathrm{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_{\mathrm{r}} <
-u_0$, in downwelling regions, the outflow of solid material toward
the BMO makes $\mathcal{H}_{\mathrm{s}}=0$ and therefore leads to an imposed zero
gradient. Conversely, where $u_{\mathrm{r}}>-u_0$, $\mathcal{H}_{\mathrm{s}}=1$ and the
value of the temperature is set to the freezing one. See
\citet{Agrusta_etal2019} 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 \citet{Bolrao_etal2021} 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~(\ref{eq:TL-P-xi-gen})) \citep{Bolrao_etal2021}. 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~(\ref{eq:momentum}). 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~(\ref{eq:phase-chage-BC}) 
\citep[see][for a full development]{Labrosse_etal2018}, 
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
\citep[e.g.][]{Nakagawa_Tackley2012}, 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.

\subsection{Evolution of the BMO}
\label{sec:evolution-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_{\mathrm{L}}$, and the solidus,
$T_{\mathrm{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 \citet{Labrosse2015}.


\subsubsection{Evolving reference state}
\label{sec:evolv-refer-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_{\mathrm{l}}$, the corresponding one in the solid mantle 
being $\xi_{\mathrm{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\overline{\xi}_{\mathrm{l}}/\partial r=0$, and that it is also
isentropic on average. Alternatively, the
BMO could start stably stratified \citep{Laneuville_etal2018} 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_{\mathrm{l}}$ and pressure $P$. The reference state is well mixed,
$\partial\xi_{\mathrm{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
\citep[e.g.][]{Braginsky_Roberts95,Lister_Buffett95,Labrosse2015}
{\begin{equation}
  \label{eq:diff-mu-a}
  \dpd{\mu}{r} = \left(\dpd{\mu}{P}\right)_{s,\xi}\dod{P}{r} = - \beta g,
\end{equation}}\unskip
$\beta$ being the coefficient of chemical contraction
{\begin{equation}
  \label{eq:beta}
  \beta\equiv-\frac{1}{\rho}\left(\dpd{\rho}{\xi_{\mathrm{l}}}\right)_{P,s}=\rho\left(\dpd{\mu}{P}\right)_{s,\xi}.
\end{equation}}\unskip
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~(\ref{eq:diff-mu-a}) can readily be integrated to give
{\begin{equation}
  \label{eq:mu-r}
  \mu = \mu^{\mathrm{BOMB}} + \beta g (r_{\mathrm{BOMB}} - r) \equiv \mu^{\mathrm{BOMB}} + \mu^\prime
\end{equation}}\unskip
with $\mu^{\mathrm{BOMB}}$ the chemical potential at the top of the BMO and
$\mu^\prime$ the deviation from that value in the BMO.

Similarly, for the temperature we get the classical isentropic gradient:
{\begin{equation}
  \label{eq:isentropic-grad}
  \dpd{T}{r} = \left(\dpd{T}{P}\right)_{s,\xi}\dod{P}{r} = -
  \frac{\alpha g T}{C_p},
\end{equation}}\unskip
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
\citep{Anderson_etal1992,Chopelas_Boehler1992,Duffy_Ahrens1993,Ricard2007,Ricard_etal2022} 
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~(\ref{eq:isentropic-grad}) can be
integrated to give
{\begin{equation}
  \label{eq:isentrope}
  T=T^{\mathrm{BOMB}}\exp\left[\frac{\alpha g (r_{\mathrm{BOMB}}-r)}{C_p}\right],
  \end{equation}}\unskip
with $T^{\mathrm{BOMB}}$ the temperature at the top of the BMO. This equation
can be safely linearised if that simplifies the expressions.
$T^{\mathrm{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
{\begin{eqnarray}
 T^{\mathrm{BOMB}}&=&T_{\mathrm{L}}(r_{\mathrm{BOMB}}) = T_{\mathrm{L}}(r_0)
  - \dpd{T_{\mathrm{L}}}{P} \rho g \left(r_{\mathrm{BOMB}} - r_0\right)\nonumber\\
&&  +\, \dpd{T_{\mathrm{L}}}{\xi} \left(\xi_{\mathrm{l}} - \xi_{l0}\right),
    \label{eq:TL-P-xi-gen}
\end{eqnarray}}\unskip
with $r_0$ the initial position of the freezing front, $\xi_{l0}$ the
initial FeO mass fraction of the magma ocean and $T_{\mathrm{L}}(r_0)$ the
corresponding liquidus value. This equation can be time differentiated
to give
{\begin{equation}
  \label{eq:dTBMO-dt}
  \dod{T^{\mathrm{BOMB}}}{t} =  - \dpd{T_{\mathrm{L}}}{P} \rho g \dod{r_{\mathrm{BOMB}}}{t}
  +\dpd{T_{\mathrm{L}}}{\xi} \dod{\xi_{\mathrm{l}}}{t}.
  \end{equation}}\unskip
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 \citep{Labrosse_etal07,Ulvrova_etal2012}. Therefore, we neglect
the thickness of boundary layers and their associated temperature,
composition and chemical potential jumps. 

\subsubsection{Balance equations}
\label{sec:balance-equations}
The mass fraction $\xi_{\mathrm{l}}$ 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_{\mathrm{l}}$ are derived for
their introduction in StagYY. The
integrated fluxes over the CMB and top surface of the BMO are denoted
$I^{\mathrm{CMB}}$ and $I^{\mathrm{BOMB}}$, respectively, and are counted positive
upward. The global balance equation
for the FeO content is
{\begin{equation}
  \label{eq:FeO-Conserv}
  \dod{M_{\mathrm{BMO}}\xi_{\mathrm{l}}}{t} = I^{\mathrm{CMB}} - I^{\mathrm{BOMB}},
\end{equation}}\unskip
$M_{\mathrm{BMO}}$ being the mass of the BMO in the $\xi_{\mathrm{l}}=0$ limit.
For the time being, we will neglect $I^{\mathrm{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~(\ref{eq:FeO-Conserv}) can be developed to give
{\begin{equation}
  \label{eq:FeO-Conserv2}
  M_{\mathrm{BMO}}\dod{\xi_{\mathrm{l}}}{t} + 4\rmpi r_{\mathrm{BOMB}}^2\rho_0 \xi_{\mathrm{l}} \dot{r}_{\mathrm{BOMB}} = - I^{\mathrm{BOMB}}.
\end{equation}}\unskip

The flux to the solid mantle, $I^{\mathrm{BOMB}}$, is due to the phase
change happening at the boundary, with a crystallisation mass rate $w$
associated with the solid radial velocity $u_{\mathrm{r}}$ at the
boundary as computed in StagYY, $u_{\mathrm{r}}$, by $w =
u_{\mathrm{r}} - \dot{r}_{\mathrm{BOMB}}$ since, in StagYY, the boundary
is kept fixed by a continuous adjustment of the computation domain
thickness (Section~\ref{sec:conv-solid-mantle}). We get
{\begin{equation}
  \label{eq:I_BMO}
  I^{\mathrm{BOMB}}=\int_{\rmOmega_{\mathrm{BOMB}}} \rho_0 \xi_{\mathrm{s}} \left(u_{\mathrm{r}} -
    \dot{r}_{\mathrm{BOMB}}\right) \dif{\rmOmega},
\end{equation}}\unskip
with $\rmOmega_{\mathrm{BOMB}}$ the surface of the boundary and
$\xi_{\mathrm{s}}$ the mass fraction of FeO on the solid side of the
boundary:
{\begin{equation}
  \label{eq:xi_s}
  \xi_{\mathrm{s}} =
  \begin{cases}
    D \xi_{\mathrm{l}}, & u_{\mathrm{r}} \ge 0 \\
    \xi_{\mathrm{s}}, & u_{\mathrm{r}}<0.
  \end{cases}
\end{equation}}\unskip

In regions of solid upwelling, fractional crystallisation occurs,
producing a solid with a mass fraction $\xi_{\mathrm{s}}= D \xi_{\mathrm{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 \citet{Bolrao_etal2021}. In
regions of downwelling, the solid with a mass fraction $\xi_{\mathrm{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 \citep{Labrosse_etal2018}.

Combining Equations~(\ref{eq:FeO-Conserv2}) and~(\ref{eq:I_BMO})
gives, after rearrangement:
{\begin{eqnarray}
&&  \left(4\rmpi r_{\mathrm{BOMB}}^2 \rho_0 \xi_{\mathrm{l}} -
    \int_{\rmOmega_{\mathrm{BOMB}}}\rho_0\xi_{\mathrm{s}}\dif{\rmOmega} \right)\dot{r}_{\mathrm{BOMB}}\nonumber\\
&&\qquad  +\, \frac{4\rmpi}{3}\rho_0\left(r_{\mathrm{BOMB}}^3-r_{\mathrm{CMB}}^3\right) \dot{\xi}_{\mathrm{l}}\nonumber\\
&&\quad  = -\int_{\rmOmega_{\mathrm{BOMB}}}\rho_0\xi_{\mathrm{s}}u_{\mathrm{r}}\dif{\rmOmega}.  \label{eq:var-xi}
\end{eqnarray}}\unskip
The two integrals in Equation~(\ref{eq:var-xi}) are computed in StagYY
at each time-step and that gives the relationship between $\dot{r}_{\mathrm{BOMB}}$ and
$\dot{\xi}_{\mathrm{l}}$. Equation~(\ref{eq:dTBMO-dt}) then allows us to express
the rate of change of $T^{\mathrm{BOMB}}$ as a function of $\dot{r}_{\mathrm{BOMB}}$
only. 

Note that in the case $u_{\mathrm{r}}=0$, which is obtained for $\Phi=\infty$,
$\xi_{\mathrm{s}}=D\xi_{\mathrm{l}}$ is uniform and the mass flux of FeO to the solid mantle is
$I^{\mathrm{BOMB}}=-4\rmpi r_{\mathrm{BOMB}}^2 D \xi_{\mathrm{l}}\rho_0\dot{r}_{\mathrm{BOMB}}$. 
This flux is positive for $\dot{r}_{\mathrm{BOMB}}<0$, i.e.\ when crystallisation
occurs. Equation~(\ref{eq:var-xi}) leads to
{\begin{equation}
  \label{eq:var-xi-ur0}
  \dot{\xi}_{\mathrm{l}} = - \frac{3 r_{\mathrm{BOMB}}^2
    \rmDelta\xi}{r_{\mathrm{BOMB}}^3-r_{\mathrm{CMB}}^3} \dot{r}_{\mathrm{BOMB}}
\end{equation}}\unskip
with $\rmDelta\xi=\xi_{\mathrm{l}}-\xi_{\mathrm{s}}$. This equation was used in the original
BMO paper \citep{Labrosse_etal07} but the general
Equation~(\ref{eq:var-xi}) accounting for flow through the boundary is
used in the present study.

The flux $\vvect{i}$ of FeO across both boundaries contributes to the
total heat flux density $\vvect{q}$ as
{\begin{equation}
  \label{eq:heat-flow-chempot}
  \vvect{q}=-k\vect{\nabla}T+\mu\vvect{i}
\end{equation}}\unskip
and we note the thermal part as
{\begin{equation}
  \label{eq:q-T}
  \vvect{q}_T\equiv -k\vect{\nabla}T.
\end{equation}}\unskip
Integrated over a surface $\rmOmega$ (BOMB or CMB), we get
{\begin{equation}
  \label{eq:Q-integrated}
  Q^\rmOmega=Q_T^\rmOmega+\mu^\rmOmega I^\rmOmega
\end{equation}}\unskip
with $\mu^\rmOmega$ the average value of $\mu$ over the
surface $\rmOmega$. 

The long term thermal evolution of the BMO is controlled by the
integrated energy balance equation, which is written as
\citep[e.g.][]{Braginsky_Roberts95,Lister_Buffett95}
{\begin{eqnarray}
&&  \int_{V_{\mathrm{BMO}}}\rho\left(T\dpd{s_{\mathrm{l}}}{t}+\mu\dpd{\xi_{\mathrm{l}}}{t}\right)\dif{V}\nonumber\\
&&\quad=   \int_{\rmOmega_{\mathrm{BOMB}}} \rho L (u_{\mathrm{r}} - \dot{r}_{\mathrm{BOMB}}) \dif{\rmOmega}\nonumber\\
&&\qquad  +\, Q_{\mathrm{R}}+Q^{\mathrm{CMB}}-Q^{\mathrm{BOMB}},  \label{eq:energy}
\end{eqnarray}}\unskip
with $Q_{\mathrm{R}}$ the radiogenic heat production in the BMO and $L$ the
latent heat of melting, which contains two contributions:
{\begin{equation}
  \label{eq:latent-heat}
  L = T_{\mathrm{L}}\rmDelta s + \mu \rmDelta\xi
\end{equation}}\unskip
$\rmDelta s = s_{\mathrm{l}}-s_{\mathrm{s}}$ being the entropy of melting.
The surface integral on the right-hand-side of
Equation~(\ref{eq:energy}) is the total latent released or consumed by
phase change at the boundary and is related to the crystallisation change rate
$u_{\mathrm{r}}-\dot{r}_{\mathrm{BOMB}}$. The entropy contribution to the latent heat,
$T_{\mathrm{L}}\rmDelta s$, is uniform along the phase change boundary and, since
$\rho u_{\mathrm{r}}$ averages to 0 at the interface, only the net
crystallisation rate $\dot{r}_{\mathrm{BOMB}}$ contributes to that part. 
Combining Equation~(\ref{eq:energy}) with the global balance in FeO
Equation~(\ref{eq:FeO-Conserv}) and the expressions for the BOMB and
CMB fluxes Equation~(\ref{eq:Q-integrated}) gives
{\begin{eqnarray}
    Q_T^{\mathrm{BOMB}} &= & - \int_{V_{\mathrm{BMO}}}\rho T \dpd{s}{t}\dif{V}\nonumber\\
&&-\, T_{\mathrm{L}}(r_{\mathrm{BOMB}})\rmDelta s 4\rmpi
    r_{\mathrm{BOMB}}^2\dot{r}_{\mathrm{BOMB}} +Q_{\mathrm{R}}+ Q_T^{\mathrm{CMB}}\nonumber\\
&&+\,\left(\mu^{\mathrm{CMB}}-\mu^{\mathrm{BOMB}}\right) I^{\mathrm{CMB}}
    -  \dod{\xi_{\mathrm{l}}}{t} \int_{V_{\mathrm{BMO}}}\rho\mu^\prime\dif{V}\nonumber\\
&& +\, \int_{\rmOmega_{\mathrm{BOMB}}}\rho\mu^{\mathrm{BOMB}}\left(\rmDelta\xi 
u_{\mathrm{r}}+\xi_{\mathrm{s}}\dot{r}_{\mathrm{BOMB}}\right)\dif{\rmOmega}.
  \label{eq:energy2}
\end{eqnarray}}\unskip
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}_{\mathrm{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^{\mathrm{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^{\mathrm{CMB}}$, for now. $Q_T^{\mathrm{BOMB}}$ is
provided at each time step by the convection model of the solid
mantle, StagYY. Computing $Q_{\mathrm{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^\prime$ in
Equation~(\ref{eq:energy2}), is easily related to
rate of change of $r_{\mathrm{BOMB}}$ using Equations~(\ref{eq:mu-r}):
{\begin{eqnarray}
  \label{eq:E-chi}
  E_{\chi 1} &\equiv& -  \dod{\xi_{\mathrm{l}}}{t} \int_{V_{\mathrm{BMO}}}\rho\mu^\prime\dif{V}
  =   - \rho_0\beta g \left[r_{\mathrm{BOMB}} V_{\mathrm{BMO}}
  \vphantom{r_{\mathrm{BOMB}}^4}
  \right.\nonumber\\ 
  &&\left.-\,\rmpi 
  \left(r_{\mathrm{BOMB}}^4-r_{\mathrm{CMB}}^4\right)\right] \dod{\xi_{\mathrm{l}}}{t}.
\end{eqnarray}}\unskip
Using Equation~(\ref{eq:var-xi}), this term can be related to
$\dot{r}_{\mathrm{BOMB}}$.

The second compositional term (last on the right-hand-side of
Equation~(\ref{eq:energy2})) can be simplified assuming $\mu^{\mathrm{BOMB}}$,
$\rho$ and $\xi_{\mathrm{l}}$ to be constant and using the fact that $u_{\mathrm{r}}$ is 0 on
average on the boundary. This gives
{\begin{equation}
  \label{eq:E-chi2}
  E_{\chi 2} = \rho_0 \mu^{\mathrm{BOMB}}\int_{\rmOmega_{\mathrm{BOMB}}} 
  \xi_{\mathrm{s}} \left(\dot{r}_{\mathrm{BOMB}}-u_{\mathrm{r}}\right)\dif{\rmOmega}.
\end{equation}}\unskip


The secular cooling term is expressed as
{\begin{equation}
  \label{eq:sec-cool-s}
  Q_C \equiv - \int_{V_{\mathrm{BMO}}}\rho T \dpd{s}{t}\dif{V} = -  
  \int_{V_{\mathrm{BMO}}}\rho C_p \dpd{T}{t}\dif{V},
\end{equation}}\unskip
which can be computed using the Equations~(\ref{eq:isentrope}) and 
(\ref{eq:dTBMO-dt}):
{\begin{eqnarray}
  Q_C &=&  - 4\rmpi\rho_0 C_P L_{\mathrm{M}}^3\left[\dod{T^{\mathrm{BOMB}}}{t} + \frac{\alpha
          g T^{\mathrm{BOMB}}}{C_p} \dot{r}_{\mathrm{BOMB}}\right]\nonumber\\
&&\times\, \int_{r_{\mathrm{CMB}}/L_{\mathrm{M}}}^{r_{\mathrm{BOMB}/L_{\mathrm{M}}}}x^2 
          \exp\left[\frac{\alpha g L_{\mathrm{M}}}{C_P}
          \left(\frac{r_{\mathrm{BOMB}}}{L_{\mathrm{M}}}-x\right)\right]\dif{x}.\nonumber\\
  \label{eq:sec-cool-T2}
\end{eqnarray}}\unskip
Using Equation~(\ref{eq:dTBMO-dt}) and~(\ref{eq:var-xi-ur0}), this
term can be expressed as an afine function of $\dot{r}_{\mathrm{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 \citep{Labrosse2015} 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
{\begin{equation}
  \label{eq:core-cooling}
  Q_T^{\mathrm{CMB}} = - \frac{4\rmpi}{3} \rho_N C_{pN} L_\rho^3\frac{f_C\left(\frac{r_{\mathrm{CMB}}}{L_\rho},
      \gamma\right)}{\left(1-\frac{r_{\mathrm{CMB}}^2}{L_\rho^2}-A_\rho \frac{r_{\mathrm{CMB}}^4}{L_\rho^4}\right)^\gamma} \dod{T^{\mathrm{CMB}}}{t},
\end{equation}}\unskip
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\"uneisen parameter of the core and
{\begin{eqnarray}
  f_C(r, \gamma ) &= &3
  \int_0^r x^2 \left( 1 - x^2 - A_\rho x^4 \right)^{1+\gamma}\,\dif{x}\nonumber\\
&  =& r^3 \left[1-\frac{3}{5} (\gamma + 1) r^2 \right.\nonumber\\
&& \left.-\frac{3}{14} (\gamma+1)
    (2 A_\rho - \gamma ) r^4 + O(r^6) \right].\label{eq:fC}
\end{eqnarray}}\unskip
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~(\ref{eq:isentrope}):
{\begin{equation}
  \label{eq:T-CMB}
  T^{\mathrm{CMB}} = T^{\mathrm{BOMB}}\exp\left[\frac{\alpha g}{C_p}\left(r_{\mathrm{BOMB}}-r_{\mathrm{CMB}}\right)\right].
\end{equation}}\unskip
Therefore, the rate of change of $T^{\mathrm{CMB}}$ is
{\begin{eqnarray}
  \label{eq:dTCMB}
  \dod{T^{\mathrm{CMB}}}{t} &=&  \left[\dod{T^{\mathrm{BOMB}}}{t} + \frac{\alpha g
      T^{\mathrm{BOMB}}}{C_p} \dot{r}_{\mathrm{BOMB}}  \right] \nonumber\\
&&\times\,  \exp\left[\frac{\alpha g}{C_p}\left(r_{\mathrm{BOMB}}-
r_{\mathrm{CMB}}\right)\right].
\end{eqnarray}}\unskip


The equations for the evolution of the BMO and the core presented in
this section have been made dimensionless (see \hyperlink{AppA}{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. 

\section{Example of dynamical evolution}
\label{sec:example-dynam-evol}

The model is controlled by many input parameters, the
classical ones for mantle convection using StagYY \citep{Tackley2008} 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_{\mathrm{FeO}}=0.3,\ K_{\mathrm{HPE}}=10^{-2}$, respectively and the phase
change parameter is $\Phi=3\times10^{-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 $\Ray=3\times 10^{7}$, $B=5$,
respectively. These are defined using, as a temperature scale, $\rmDelta
T_{\mathrm{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 $\rmDelta T/\rmDelta T_{\mathrm{M}}$ is of
order one. The chosen value is $\rmDelta T_{\mathrm{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
\hyperlink{AppA}{Appendix}. 

Figure~\ref{fig:evol} 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. 

\begin{figure*}
\includegraphics{fig01}
\caption{\label{fig:evol}
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. }
\vspace*{-2pt}
\end{figure*}

Convection first sets in at large scale (Figure~\ref{fig:evol}), which is
usual for convection with a phase change at the bottom boundary
\citep{Labrosse_etal2018,Morison_etal2019,Morison_etal2024}. 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 of \mbox{Figure~\ref{fig:evol}}, 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 \citep{Agrusta_etal2019,Lebec_etal2023}.  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 \citep{Labrosse_etal07}, 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~\ref{fig:profiles}. 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.

\begin{figure*}
\includegraphics{fig02}
\vspace*{-4pt}
\caption{\label{fig:profiles}
Temperature (left) and composition (right) mean profiles for
the time steps corresponding to those of Figure~\ref{fig:evol}. The
radial position reflects the crystallisation of the BMO with
time.}
\vspace*{-4pt}
\end{figure*}

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~\ref{fig:evol}. Figure~\ref{fig:C-time} shows the
evolution of the minimum FeO mass fraction at the \mbox{upper} \mbox{boundary},
$\min(\xi_{\mathrm{FeO}, \mathrm{top}})$, and the maximum at the bottom boundary,
$\max(\xi_{\mathrm{FeO}, \mathrm{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~(\ref{eq:numbers_eff}). The value obtained is also shown on
Figure~\ref{fig:C-time}. Comparison with the snapshots of Figure~\ref{fig:evol}
suggests that a critical value of $B_{\mathrm{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 \mbox{layered} initial
condition \citep{Tackley1998,Davaille1999,LeBars_Davaille04b}.

\begin{figure}
\includegraphics{fig03}
\caption{\label{fig:C-time}
Time evolution of the minimum FeO mass fraction at the upper
boundary, $\min(\xi_{\mathrm{FeO}, \mathrm{top}})$, the maximum at the bottom boundary,
$\max(\xi_{\mathrm{FeO}, \mathrm{bot}})$ (left axis) and the effective buoyancy number
$B_{\mathrm{eff}}$, right axis. Empty circles on the black line gives the
position of the snapshots presented on Figure~\ref{fig:evol}.}
\vspace*{-8pt}
\end{figure}

The bottom temperature, equal to the liquidus, also varies with time
because of the combined \mbox{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~\ref{fig:profiles} 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
\citep{Labrosse_etal2018,Agrusta_etal2019} 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\gtrsim10^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 thermo-chemical 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~\ref{fig:evol}, 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.

\section{Discussion and the way ahead}
\label{sec:limit-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~\ref{fig:evol}. 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
\citep{Nomura_etal2014}.

On the other hand, the thermal structure shown in
Figures~\ref{fig:evol} and~\ref{fig:profiles} 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
\citep{Labrosse_etal2018,Agrusta_etal2019,Lebec_etal2023}, 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 \citep{Davies1990,
  Sleep2000}. 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 \citep{Agrusta_etal2019} 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~\ref{fig:profiles}). 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
mecanism 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
\citep[e.g.][]{Hes_Par95,Elkins-Tanton_etal03,Elkins-Tanton2012}
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
\citep{Ballmer_etal2017,Maurice_etal2017,Boukare_etal2018,Morison_etal2019}. 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~\ref{sec:evolution-bmo} for the surficial magma ocean and including phase change
boundary condition. Some progress has been donc in that direction
\citep{Morison2019} but a systematic exploitation of this model is
still needed and will be the topic of a future study.


\section{Conclusions}
\label{sec: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 \citep{Tackley2008}, 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 \citep{Labrosse_etal2018}.

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 \citet{Labrosse_etal07}. 

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
\citep{Nomura_etal2014}. 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. 

\begin{table*}
\caption{\label{tab:adimParamsBMO}
Dimensionless parameters for the Basal Magma Ocean evolution}
\begin{tabular}{llll}
\thead
Symbol & Expression & Description & Run value\\
\endthead
$\Gamma^-_0$ &\raisebox{1.4pc}{} $\dfrac{R^--R_{\mathrm{c}}}{d}$\vspace*{5pt} & Initial BMO thickness & 0.3\\
$T_{\mathrm{BOMB}, 0}$ & $\dfrac{T_{\mathrm{BOMB}}(t=0)}{\rmDelta T}$\vspace*{5pt} & Initial BOMB temperature & 0.37\vspace*{2pt}\\
$\xi_0$ & $\xi(t=0)$ & Initial FeO content & 0.07\vspace*{2pt}\\
$c_{\mathrm{HPE}, 0}$ & $c_{\mathrm{HPE}}(t=0)$ & Initial HPE content & $5\times10^{-3}$\\
$B$ & $-\dfrac{\beta}{\alpha\rmDelta T}$\vspace*{5pt} & Compositional buoyancy & 5\\
$\Di$ & $\dfrac{\alpha g L_{\mathrm{M}}}{C_p}$\vspace*{5pt} & BMO dissipation number & 1.74\\
$T_{LP}$ & $\dfrac{\rho gL_{\mathrm{M}}}{\rmDelta T_{\mathrm{M}}}\dpd{T_{\mathrm{L}}}{P}$\vspace*{5pt} & 
\parbox[t]{5cm}{\raggedright Pressure dependence of melting temperature} &
\parbox[t]{5cm}{\raggedright Depth-dependent from \citet{Fiquet_etal2010}}\\
$T_{L\xi}$ & $\dfrac{1}{\rmDelta T_{\mathrm{M}}}\dpd{T_{\mathrm{L}}}{\xi}$\vspace*{5pt} & 
\parbox[t]{5cm}{\raggedright Composition dependence of melting temperature} & $-$0.125\\
$\mathrm{St}$ & $\dfrac{\rmDelta s}{C_p}$\vspace*{5pt} & Stefan number & 0.3\\
$\tau_{\mathrm{hpe}}$ & $\dfrac{\tau_{{\mathrm{hpe}},d}\kappa}{L_{\mathrm{M}}^2}$\vspace*{5pt} & HPE half-life & 0.02\\
$P_{\mathrm{c}}$ & $\dfrac{P_{\mathrm{c},d}}{\rho C_p L_{\mathrm{M}}^3}$ & Core cooling contribution & 3.67 \vspace*{2pt}
\botline
\end{tabular}
\vspace*{-2pt}
\end{table*}

\section*{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.

\section*{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\^ole Scientifique de
Mod\'elisation Num\'erique) of the ENS de Lyon.


\back{}

\appendix{}

\gdef\appendixlabel{Appendix}

\hypertarget{AppA}{}
\section{Dimensionless numbers and parameter values for the
example calculation} \label{sec:dimens-numb-param}

We use dimensionless variables in the whole model, both for convection
in the solid mantle and the \mbox{evolution} of the BMO and core. The solid
mantle part is explained in Section~\ref{sec:conv-solid-mantle}. The
value expected for the phase change number $\Phi$ can be estimated
using its expression in Equation~(\ref{eq:Phi}) using reasonable
values for the various parameters:
$\rho_{\mathrm{s}}\simeq5\times10^3~\mathrm{kg}{\cdot}\mathrm{m}^{-3}$,  
$\rho_{\mathrm{l}}\simeq5.5\times10^3~\mathrm{kg}{\cdot}\mathrm{m}^{-3}$,
$C_{\mathrm{pl}}\simeq10^3~\mathrm{J}{\cdot}\mathrm{K}^{-1}{\cdot}\mathrm{kg}^{-1}$,
$\eta_{\mathrm{s}}\simeq10^{22}~\mathrm{Pa}{\cdot}\mathrm{s}$,
$d\simeq3000$~km. For the latent heat, we use 
$L= T \rmDelta S$, with $\rmDelta S=652~\mathrm{J}{\cdot}\mathrm{kg}^{-1}
{\cdot}\mathrm{K}^{-1}$ the
entropy of melting \citep{Stixrude_etal09} and
$T\simeq4000$~K the melting temperature. The slopes of
the isentropic profile and melting temperatures are estimated from
\citet{Stixrude_etal09}:
$m_{\mathrm{ad}}\simeq 15~\mathrm{K}{\cdot}\mathrm{GPa}^{-1}$,
$m_p\simeq 4~\mathrm{K}{\cdot}\mathrm{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_{\mathrm{l}}\sim10^{-4}~\mathrm{m}{\cdot}\mathrm{s}^{-1}$.
This gives $\Phi\sim4\times10^{-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\lesssim10^{-1}$
\citep{Labrosse_etal2018,Agrusta_etal2019,Morison_etal2024}. 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\times10^{-2}$, which is
sufficiently low to get the proper behaviour. 

The magma ocean evolution equations are made dimensionless with the
following scales:
\begin{itemize}
    \item $d$, the total thickness of the mantle (solid and magma ocean),
        considered constant in this study for simplicity.  For the Earth,
        $d = 2900$~km.
    \item $\rmDelta 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 $\rmDelta T=4000$~K.
    \item $\rho$, the reference density of the considered materials.  For the
        Earth's mantle, $\rho=4000~\mathrm{kg}/\mathrm{m}^3$. 
\item $\kappa = k/(\rho C_p)$, the diffusion coefficient  $\kappa =
10^{-6}~\mathrm{m}^2/\mathrm{s}$, associated with the conductivity
$k=4$~W/(mK) and the heat capacity
$C_p=10^3~\mathrm{J}{\cdot}\mathrm{K}^{-1}{\cdot}\mathrm{kg}^{-1}$.
\end{itemize}

The dimensionless parameters arising from the non-dimensionalization
are shown in Table~\ref{tab:adimParamsBMO}. More detailed can be found in
\citet{Morison2019}.

%\nocite*{}

\bibliographystyle{crgeos}
\bibliography{crgeos20240413}
\refinput{crgeos20240413-reference.tex}

\end{document}
