# Birth and decline of magma oceans.
Part 1: erosion and deposition of crystal layers in evolving magmatic
reservoirs.

###### Abstract

This paper is the first of a two companion papers presenting a theoretical and
experimental study of the evolution of crystallizing magma oceans in
planetesimals. We aim to understand the behavior of crystals formed in a
convective magma and the implications of crystal segregation for the reservoir
thermal and structural evolution. In particular, the goal is to constrain the
possibility to form and preserve cumulates and/or flotation crusts by
sedimentation/flotation of crystals. We first use lab-scale analog experiments
to study the stability and the erosion of a floating lid composed of plastics
beads over a convective viscous fluid volumetrically heated by microwave
absorption. We propose an erosion law that depends only on two dimensionless
numbers which govern these phenomena: (i) the Rayleigh-Roberts number,
characterizing the strength of convection and (ii) the Shields number, that
encompasses the physics of the flow-particle interaction. We further consider
the formation of a cumulate at the base of the convective layer due to
sedimentation of beads that are denser than the fluid. We find that particles
deposition occurs at a velocity that scales with the Stokes velocity, a result
consistent with previous experimental studies. The theoretical framework built
on these experimental results is applied in a second paper on the evolution of
magma oceans in planetesimals and the fate of particles in this convective
environment.

JGR: Planets

Université de Paris, Institut de Physique du Globe de Paris, CNRS, F-75005,
France

Cyril<EMAIL_ADDRESS>
We study experimentally the transient behavior of a convective fluid bearing
particles that can float/sediment to form deposits.

The erosion process is described by a modified Shields number that compares
the convective shear to the beads buoyancy. The deposition rate of crystals
scales with the Stokes velocity of beads in absence of convection.

This theoretical framework yields scaling laws that can be applied to various
magmatic systems based on the dimensionless numbers that characterize them.

## Plain Language Summary

At early time in planetary formation, the heat supplied by $\rm{{}^{26}Al}$
was large enough to generate massive melting episodes in planetesimals and to
produce magma oceans. Following melting and iron core differentiation, the
proto-mantle of these 10-100’s km radius rocky bodies was a magma ocean that
behaved like a convective fluid bearing solid crystals. The goal of this study
is to understand the crystal-melt segregation in a magma ocean. In order to
characterize the behavior of crystals in a convective flow we build a
theoretical model that is based on the results of lab-scale experiments. We
use a viscous fluid heated by microwave absorption in order to simulate the
convection regime expected in magma oceans. We use plastic beads to simulate
crystals in the experiments and study the formation of either a flotation lid
or a sedimented basal cumulate. We demonstrate that both deposits form only if
the the buoyancy of the fluid balances the shear imposed by convection, and we
predict the transient evolution observed in the experiments.

##  1 Introduction

Crystal segregation in magmatic reservoirs, including magma chambers [Sparks .
(2019)], lava lakes [Helz . (1989)], and planetary-scale systems such as magma
oceans [Solomatov (2000)] is a key process in their thermal evolution. Their
dynamics can be described as the one of a convective fluid bearing silicate
crystals. Because the composition of crystals differs from that of the magma,
a density difference always exists between the two phases and induces a
buoyancy force whose strength and sign depend on the type of considered
crystal. For instance, olivine crystals are usually negatively buoyant and may
form a cumulate in magma chambers [Jellinek Kerr (2001)], whereas plagioclase
crystals are positively buoyant and may form a flotation mush at the surface
[Namur . (2011)]. At planetary scale, this segregation process may also occur
and heavy crystals can form cumulates of dunites and/or harzburgites at the
core-mantle boundary (CMB) by sedimentation [Righter Drake (1997), Mandler
Elkins-Tanton (2013)], whereas light crystals accumulate and form a crust at
the surface, like plagioclases that compose the anorthosite crust of the Moon
[Wood (1970), Wood . (1970), Warren (1985)].
The fate of a suspension in magma oceans has an overriding influence on its
geochemical and geophysical evolution. Studies that deal with the formation of
deposits do not usually consider the interaction between the sedimentation and
the convection and make the hypothesis of either instantaneous crystal
segregation [Hamano . (2013), Maurice . (2020)] or stable suspension of
crystals in the liquid [Elkins-Tanton . (2011), Bryson . (2019)]. However, the
dynamics of such deposition/flotation processes in a convective system are not
well constrained. Suckale12a,Suckale12b proposed a model that takes crystal-
fluid interaction into account by developing a numerical study describing the
suspension and applied their method to constrain the formation of the lunar
crust. However, as this model is based on the direct numerical resolution of
the conservation of mass, momentum and energy, it does not provide a physical
criterion that can be used in independent studies to determine the possibility
of forming a deposit or not. Hence, we aim to provide a general criterion that
predicts the ability of forming a crust and a cumulate.
A stability criterion for a suspension has been proposed by Solomatov93b, but
gives only the maximal concentration of particles that can be sustained by
convection at steady state. Sturtz21 revisited this criteria by adapting the
physical framework describing the stability of grains in a bed-load [Shields
(1936), Charru . (2004)] to the stability of a bead layer sheared by thermal
convection in an internally heated system. They used experiments to establish
that a modified Shields number xan be used to determine if at steady state a
particles layer is stable or not, depending on the Rayleigh-Roberts number
characterizing the convective vigor and on particle properties. However,
Sturtz21 did not describe the dynamics of erosion and deposition, and in
particular the characteristic timescales that are involved. The present paper
proposes this improvement.
In the following, we model the dynamics of particle erosion and deposition in
order to study the time evolution of a suspension and associated particle
layers. Based on an experimental study, we analyse the thermal evolution of a
convective system that is insulated by an erodible lid of particles. We show
how the framework developed by Sturtz21 at steady state can be adapted to
describe the thermal evolution of the system and the lid thickness dynamics.
These results will be used in a companion paper to constrain the evolution of
magma oceans and the formation of layered mantles in planetesimals.

##  2 Experimental approach

###  2.1 Experimental set-up

We study the thermal and mechanical evolution of a convecting viscous fluid
under an erodible lid using the method described in Sturtz21. A $300\times
300\times 50\,\rm{mm^{3}}$ tank is placed in a modified microwave oven
[Surducan . (2014)]. The top surface of the tank is composed of an aluminium
plate whose temperature is fixed and monitored. The lateral and basal walls of
the tank are made of plastic (poly-methyl metacrylate, PMMA), which insulate
thermally the interior of the reservoir. The fluid inside the tank is heated
volumetrically by microwaves absorption, that generates convection at high
Rayleigh numbers [Limare . (2013)].
To image the convective fluid, a laser sheet scans the tank, and the signal is
monitored by two CCD cameras equipped with filters sensitive to different
spectral ranges. The temperature field is measured using a two-dye laser
induced fluorescence method, and the velocity field is obtained by particle
image velocimetry (PIV). The spatial resolution for each field is 0.2 and 0.8
mm respectively. For a complete description of the experimental methods, the
readers can refer to Fourel17.
The particles used are PMMA beads, and the fluid is a mixture of 44 wt$\%$
glycerol and $56$ wt$\%$ ethylene-glycol. The beads and the fluid have
different thermal expansion coefficients (see Table 2, A). As shown in Figure
1, at low temperatures beads are lighter than the fluid and can float. Above
the inversion temperature $T_{\rm{inv}}=37.4\rm{{}^{o}C}$, beads become
negatively buoyant and can sink. We kept the surface temperature below
$37.4\rm{{}^{o}C}$ so that a floating lid can form at the surface of the
convective layer. Meanwhile, the bulk is heated by microwaves and its
temperature can eventually exceed the inversion temperature; in these
conditions a basal cumulate can form. The experimental system can thus be used
to study the formation and stability of both a flotation crust and/or a basal
cumulate.

Figure 1: Variation of the respective density of the beads and of the experimental fluid as a function of the temperature, taken from Sturtz21. The thermal expansion coefficient is $\alpha_{p}=3.17\,10^{-4}\ \rm{K^{-1}}$ and $\alpha_{f}=5.5\,10^{-4}\ \rm{K^{-1}}$ for the beads and the fluid respectively. Beads are positively buoyant for $T<37.4\rm{{}^{o}C}$ and negatively buoyant for $T>37.4\rm{{}^{o}C}$. Name | Beads radius ($\mu$m) |  $\delta_{0}$ (mm) |  $T_{s}$ (oC) |  $T_{bulk}$ (oC) |  $Ra_{H}$ ($10^{7}$) | Erosion? | Cumulate? | Heat bump?
---|---|---|---|---|---|---|---|---
IHB04 | $290$ | 5.3 | 21.8 | 44.5 | 4.3 | Partial | yes | no
IHB05 | $290$ | 5.3 | 21.8 | 47.7 | 8.9 | Partial | yes | no
IHB07 | $290$ | 1.9 | 22.3 | 35.9 | 2.8 | Partial | no | no
IHB08 | $290$ | 1.9 | 14 | 32.8 | 2.4 | Partial | no | no
IHB09 | $290$ | 1.9 | 34.1 | 43.4 | 4.1 | Partial | yes | no
IHB11 | $290$ | 4.7 | 22.8 | 43.9 | 4.2 | Partial | yes | no
IHB12 | $290$ | 4.7 | 22.8 | 35.3 | 1.4 | Partial | no | no
IHB13 | $290$ | 4.7 | 23 | 48.3 | 8.0 | Partial | yes | no
IHB14 | $290$ | 4.7 | 8.9 | 42.0 | 5.9 | Partial | no | yes
IHB16 | $290$ | 3.8 | 22.3 | 53.7 | 14.4 | Partial | yes | no
IHB17 | $290$ | 3.8 | 22.2 | 38.7 | 2.8 | Partial | no | yes
IHB19 | $145$ | 4.7 | 22.9 | 42.1 | 6.8 | Total | no | yes
IHB20 | $145$ | 4.7 | 29.2 | 45.7 | 8.1 | Total | no | yes
IHB21 | $145$ | 4.7 | 23.1 | 30.6 | 1.5 | Total | no | yes
IHB22 | $145$ | 4.7 | 30.5 | 45.3 | 10 | Total | yes | yes
IHB23 | $145$ | 4.7 | 23.2 | 33.1 | 2.6 | Total | no | yes
IHB24 | $145$ | 4.7 | 23.1 | 26.3 | 0.27 | Partial | no | yes
IHB25 | $145$ | 4.7 | 23.2 | 28.5 | 0.67 | Total | no | yes
IHB26 | $145$ | 4.7 | 33.2 | 36.3 | 1.0 | Total | no | yes
IHB27 | $145$ | 4.7 | 34.4 | 43.7 | 4.4 | Total | no | yes
IHB33 | $145$ | 6.0 | 22.7 | 36.8 | 3.1 | Partial | no | yes
Table 1: Experimental conditions: the beads average radius $r$ (note that two
families of beads are investigated), the initial floating bed thickness
$\delta_{0}$, the imposed surface temperature $T_{s}$, the mean bulk
temperature at steady state $T_{\mathrm{bulk}}$, the Rayleigh-Roberts number
calculated at steady state $Ra_{H}$. The next two columns inform about the
erosion of the floating lid (whether partial or total), and the formation of a
cumulate at the bottom of the tank. The last column deals with the shape of
the thermal evolution and designates whether the system experienced a heat
bump or not.

###  2.2 Typical experiments - qualitative behavior of the floating lid

Figure 2: Snapshots from experiment IHB33 ($Ra_{H}=3.10^{6}$,
$T_{s}=23\rm{{}^{o}C}$) showing the destabilization of the floating lid. Dune-
shape cusps form once beads have begun to move. Erosion is then triggered and
the lid begins to thin. The complete evolution of the lid along the experiment
is shown in the movie available in supplementary materials.

During experiments (summarized in Table 1), we record the evolution of the
floating lid thickness in order to study the mechanism of destabilization.
Snapshots of the lid are displayed in Figure 2 and the complete movie is
available in supplementary materials. Convection occurs quasi instantaneously
in all experiments. A bed-load at the base of the floating lid forms after a
delay of a few to several tens of minutes. At that stage, particles move at
the surface of the bed and form dunes. This mechanism facilitates beads re-
suspension, following the mechanism described by Solomatov93a. Once erosion
has been triggered, the lid thins until it reaches an equilibrium thickness at
steady state. The final thickness is zero in case of total erosion. We observe
that the destabilization of the lid obeys an erosion mechanism, rather than a
mechanism involving oscillation of domes and basins as it can be the case for
convection in stratified fluids [Jaupart . (2007), Limare . (2019)]. This
observation prompts to the modeling of the lid evolution by an erosion
mechanism such as the one described in the next section.

##  3 Thermal convection under an erodible lid

###  3.1 Convection driven by internal heating

In a fluid internally heated, the temperature scale $\Delta T_{H}$ is [Roberts
(1967)]:

$\Delta T_{H}=\frac{Hh^{2}}{\lambda_{\rm{f}}},$ (1)

with $h$ the reservoir thickness, $H$ the rate of internal heating and
$\lambda_{\rm{f}}$ the fluid thermal conductivity (table 2). Convection is
fully characterized by two dimensionless numbers: the Rayleigh-Roberts number
$Ra_{H}$ and the Prandtl number $Pr$ defined as follows:

$\displaystyle Ra_{H}$ $\displaystyle=$
$\displaystyle\frac{\alpha_{\rm{f}}\rho_{0,\rm{f}}gHh^{5}}{\eta_{\rm{f}}\kappa_{\rm{f}}\lambda_{\rm{f}}},$
(2) $\displaystyle Pr$ $\displaystyle=$
$\displaystyle\frac{\eta_{\rm{f}}}{\rho_{0,\rm{f}}\kappa_{\rm{f}}},$ (3)

where $\alpha_{\rm{f}}$ is the fluid thermal expansion coefficient,
$\rho_{0,\rm{f}}$ its density at a reference temperature, $g$ the acceleration
of gravity, $\eta_{\rm{f}}$ the fluid dynamical viscosity, and
$\kappa_{\rm{f}}$ the fluid thermal diffusivity. In an internally heated
convective system, there is a single thermal boundary layer (TBL) at the top
boundary of the convective layer. In the experiments performed in this study,
$Ra_{H}$ is high enough for the average bulk temperature to be quasi-
isothermal, except in the upper thermal boundary layer (Figure 3). In the
high-$Pr$ limit, its thickness of the TBL $\delta_{TBL}$ and the drop of
temperature across it $\Delta T_{TBL}$ at steady state scale with $Ra_{H}$ as
follows [Limare . (2015)]:

$\displaystyle\delta_{TBL}$ $\displaystyle=$ $\displaystyle
C_{\delta}\,h\,Ra_{H}^{-1/4},$ (4) $\displaystyle\Delta T_{TBL}$
$\displaystyle=$ $\displaystyle C_{T}\,\Delta T_{H}\,Ra_{H}^{-1/4},$ (5)

with $C_{T}=3.41$ [Vilella . (2018)] and $C_{\delta}=7.36$ [Limare . (2015)]
constants that only depend on the mechanical boundary condition at the top.
If a lid of thickness $\delta(t)$ is present at the top boundary, convection
occurs in a the fluid layer of thickness $h-\delta(t)$ and the conservation of
energy is written [Kerr . (19901), Kerr . (19902)]:

$\rho_{\rm{f}}c_{p,\rm{f}}\frac{\partial T_{\rm{bulk}}}{\partial
t}=H-\frac{Q_{s,\rm{conv}}}{h-\delta},$ (6)

with $T_{\rm{bulk}}$ the bulk volume averaged temperature, $c_{p,\rm{f}}$ is
the fluid specific heat, and $Q_{s,\rm{conv}}$ is the heat flux out of the
convective layer.
At steady state, $Q_{s,\rm{conv}}$ equals $(h-\delta)H$, a relationship that
does not hold if the temperature of the system or the thickness of the lid are
not constant. To determine a general expression for $Q_{s,\rm{conv}}$, we
first introduce the modified rate of internal heating $H^{*}$, including both
internal heating and secular cooling effects [Vilella Kaminski (2017),
Kaminski . (2020), Limare . (2021)]:

$H^{*}=H-\rho_{f}c_{p,f}\frac{\partial T_{\rm{bulk}}}{\partial t},$ (7)

and we furthermore consider that the temperature difference between the bulk
and the base of the lid, $T_{\rm{bulk}}-T_{\rm{lid}}$, scales as the
temperature drop through the TBL given by (5):

$T_{\rm{bulk}}-T_{\rm{lid}}=C_{T}^{*}\Delta T_{H}^{*}\,(Ra_{H}^{*})^{-1/4},$
(8)

with $C_{T}^{*}$ a constant that can depend on the mechanical boundary
conditions, and $Ra_{H}^{*}$ the modified Rayleigh-Roberts number obtained by
substituting $H$ by $H^{*}$ in (2). After some algebra, we get the following
equation for the time-dependent evolution of the bulk fluid temperature:

$\rho_{\rm{f}}c_{p,\rm{f}}\frac{\partial T_{\rm{bulk}}}{\partial
t}=H-\frac{1}{(C_{T}^{*})^{4/3}}\frac{\lambda_{\rm{f}}}{h-\delta}\left(\frac{\alpha_{\rm{f}}\rho_{0,\rm{f}}g}{\kappa_{\rm{f}}\eta_{\rm{f}}}\right)^{1/3}\,(T_{\rm{bulk}}-T_{\rm{lid}})^{4/3}.$
(9)

Comparing (6) and (9), we further obtain a the scaling law for
$Q_{s,\rm{conv}}$:

$Q_{s,\rm{conv}}=\frac{\lambda_{\rm{f}}}{(C_{T}^{*})^{4/3}}\,\left(\frac{\alpha_{\rm{f}}\rho_{0,\rm{f}}g}{\kappa_{\rm{f}}\eta_{\rm{f}}}\right)^{1/3}\,(T_{\rm{bulk}}-T_{\rm{lid}})^{4/3},$
(10)

with $C_{T}^{*}=3.59\pm 0.15$ according to Limare21.

Figure 3: Schematic view of the evolving system. The red solid line represents
the horizontally averaged temperature profile, $T_{\rm{lid}}$ is the basal
temperature of the lid, $T_{s}$ is the surface temperature, $Q_{s}$ is the
surface heat flux, $Q_{s,\rm{lid}}$ is the heat flux that enters the
conductive lid, $Q_{s,\rm{conv}}$ is the heat flux out of the convective bulk
and $\delta$ the thickness of the floating lid.

###  3.2 Evolution of the erodible lid

We now describe the evolution of the floating lid, whose thickness and
temperature are controlled by the coupling between the lid and the convective
fluid in two ways: (i) thermally, by the heat flux that is extracted from the
convective bulk and that imposes a thermal boundary condition at the base of
the lid, (ii) mechanically, by the erosion process occurring at the interface.

####  3.2.1 Thermal state of the floating lid

To determine the thermal state of the lid, we solve the Cauchy problem set as
follows. The temperature field in the lid $T_{b}(z,t)$ follows the heat
equation for conduction without internal heating source, noting that beads do
not absorb microwaves:

$\frac{\partial T_{b}(z,t)}{\partial
t}=\kappa_{b}\,\frac{\partial^{2}T_{b}(z,t)}{\partial z^{2}},$ (11)

where $\kappa_{b}$ is the lid thermal diffusivity. This equation is subject to
the initial and surface conditions:

$\displaystyle T_{b}(z,t=0)$ $\displaystyle=$ $\displaystyle T_{s}.$ (12)
$\displaystyle T_{b}(z=0,t)$ $\displaystyle=$ $\displaystyle T_{s},$ (13)

with $T_{s}$ the surface temperature. The basal condition is:

$\lambda_{\rm{b}}\,\frac{\partial T_{b}(z,t=0)}{\partial z}=Q_{s,\rm{lid}}$
(14)

with $\lambda_{b}$ the lid thermal conductivity, and $Q_{s,\rm{lid}}$ the heat
flux at the base of the lid. The global transient thermal state of the system
is shown in Figure 3. One should note that at this stage we do not write a
priori that the basal heat flux is equal to the heat flux out of the
convective bulk $Q_{s,\rm{conv}}$, as the thermal boundary layer at the base
of the lid and the top of the convective layer has to be taken into account.
Following Kerr90a,Kerr90b, heat conservation at the moving interface yields:

$Q_{s,\rm{lid}}(t)=\rho
c_{p}(T_{\rm{bulk}}-T_{\rm{lid}})\,\frac{\rm{d}\delta}{\rm{d}t}+Q_{s,\rm{conv}}(t).$
(15)

The left side term stands for the variation of internal energy of the small
moving thermal boundary layer due to change of temperature. In our
experiments, the first term in the righthand side is negligible compared to
the flux out of the convective bulk (see B) and we take
$Q_{s,\rm{lid}}(t)=Q_{s,\rm{conv}}(t)$.

####  3.2.2 Erosion/deposition model

To describe and quantify the erosional thinning of the lid in our experiments,
we use the formalism of Shields36. Within this framework, we consider that
particles are subject to two forces: (i) the convective shear stress that
tends to entrain particles, (ii) the friction forces, proportional to the
bead’s buoyancy, that tends to maintain the particles at the surface. The
ratio between these two forces defines the Shields number at the base of the
lid:

$\zeta_{\rm{lid}}=\frac{\eta_{\rm{f}}\dot{\gamma}}{\Delta\rho(T_{\rm{lid}})gr},$
(16)

where $\dot{\gamma}$ is the strain rate due to convection, and $r$ is the bead
radius. Erosion starts when the Shields number is larger than a critical value
$\zeta_{c}$.
Sturtz21 showed experimentally that the volume average root mean square (RMS)
strain rate scales with $\kappa/h^{2}\,Ra_{H}^{*,3/8}$ so that the Shields
number can be expressed as a function of the characteristics of convection:

$\zeta_{\rm{lid}}=\frac{\eta_{\rm{f}}\kappa_{\rm{f}}}{\Delta\rho(T_{\rm{lid}})grh^{2}}\,Ra_{H}^{*,3/8},$
(17)

and further obtained experimentally that for internally heated convective
system $\zeta_{c}=0.29\pm 0.17$.
To obtain an evolution equation for the lid thickness $-$ which will be used
later to describe the full evolution of the system $-$ we carry out the mass
balance of particles at the base of the lid following the erosion-deposition
model [Charru . (2004), Lajeunesse . (2010)]. The evolution of the lid
thickness is given by the balance between the flux of particles that settle at
the lid per unit area $\phi_{d}$ and the flux of eroded particles per unit
area $\phi_{e}$, i.e.:

$\frac{\rm{d}\delta}{\rm{d}t}=\phi_{d}-\phi_{e}.$ (18)

As the lid interface is located in the thermal boundary layer where the
convective movements are cold downwellings, we assume that if particles are
eroded, they are not able to re-settle as they are entrained by cold
downwellings. Hence we take $\phi_{d}\approx 0$. Besides, the flux of eroded
particles $\phi_{e}$ is given by the product of the volume flux of fluid
drained by the downwellings $\phi_{v}$ and the volume fraction of particles in
the TBL $c_{\rm{TBL}}$:

$\phi_{e}\sim\phi_{v}\,c_{\rm{TBL}}.$ (19)

The volume that is drained by the downwellings from the TBL depends on their
characteristics. If $N_{i}$ is the number of downwellings per unit area,
$A_{i}$ their cross section and $W_{i}$ their vertical velocity, then the
volume flux of fluid that is drained per unit area is given by $\phi_{v}\sim
N_{i}A_{i}W_{i}$. Using scaling laws proposed by Vilella18 : $N_{i}\sim
h^{-2}\,Ra_{H}^{1/4}$, $A_{i}\sim h^{2}\,Ra_{H}^{-3/8}$ and
$W_{i}\sim\kappa/h\,Ra_{H}^{3/8}$, we get:

$\phi_{v}\sim\frac{\kappa}{h}\,Ra_{H}^{*,1/4}.$ (20)

Besides, the volume fraction of particles in the TBL $c_{\rm{TBL}}$ is given
by:

$c_{TBL}\sim\frac{v_{1}\,n}{\delta_{TBL}},$ (21)

with $n$ the volume of eroded particles per unit area, and $v_{1}$ the volume
of one bead. Following Lajeunesse10, we assume that the number of particles
per unit area $n$ that are eroded scales with the difference between
$\zeta_{\rm{lid}}$ and the critical value:

$n\sim\frac{1}{r^{2}}\,(\zeta_{\rm{lid}}-\zeta_{c}).$ (22)

Combining (18)-(22), we get the following lid thickness evolution:

$\frac{\rm{d}\delta}{\rm{d}t}=-c_{e}\,\frac{\kappa
r}{h^{2}}\,Ra_{H}^{*,1/2}\,(\zeta_{\rm{lid}}-\zeta_{c}),$ (23)

with $c_{e}$ a constant to be determined experimentally.

Figure 4: Comparison of experimental data and model prediction for two
representative experiments (IHB05 top, and IHB33 bottom). The best fit
parameter $c_{e}$ is $0.1$ and $1$ for IHB05 and IHB33 respectively.

##  4 Comparison between model predictions and experimental data

###  4.1 Bulk temperature record

Figure 4 shows the two families of thermal evolution that are observed in
experiments. The first evolution is described by a monotonic increase of the
temperature towards a steady state thermal plateau reached the steady state.
The second one is characterized by an episode of transient thermal maximum
followed by a monotonic decrease towards a steady state. To confirm the
ability of the model to reproduce these two kinds of evolution, we solve
numerically the evolution equations using as input parameters: the surface
temperature $T_{s}$, the rate of internal heating $H$, the initial thickness
of the floating lid $\delta_{0}$, and the final thickness $\delta_{\rm{th}}$.
Note that the final thickness is not a direct input parameter but rather
deduced from the steady state, assuming conductive heat transfer in the
floating lid (see Sturtz21):

$\displaystyle Hh$ $\displaystyle=$
$\displaystyle\lambda_{b}\,\frac{T_{\rm{lid}}-T_{s}}{\delta_{\rm{th}}},$ (24)
$\displaystyle\frac{\delta_{\rm{th}}}{h}$ $\displaystyle=$
$\displaystyle\frac{\lambda_{b}}{\lambda_{f}}\
\left(\frac{T_{\mathrm{bulk}}-T_{\mathrm{s}}}{\Delta
T_{H}}-C_{T}Ra_{H}^{-1/4}\right),$ (25)

where $T_{\rm{lid}}$ is given by (8). We solve (6) and (23) using an explicit
finite differences scheme and we solve (11)-(14) using an implicit finite
difference scheme. At each step, we calculate the Shields number using the
scaling law (17) to obtain the erosion rate. More details on the numerical
scheme are given in C.
By construction our model gives the coupled evolution of $T_{\rm{lid}}$ and
$\delta$. The comparison between the evolution of $\delta$ during an
experiment and the model prediction thus provides an independent test of the
model’s consistency. In order to measure the lid thickness over time, we
analyse the light scattered by the lid. As illustrated in Figure 2, when the
beads are compacted and settled, they scatter more light than when they are in
suspension. Thus we use the number of pixels of high intensity at the top of
the frames as a proxy of the total intensity scattered by the lid. The thicker
the bed, the higher the intensity scattered by the floating lid and the higher
the number of pixels counted. To calibrate the link between the number of
pixels of high intensity $\mathcal{I}(t)$ and the thickness of the bed
$\delta(t)$, we measured the number of pixels of high intensity at the
beginning of the experiment $\mathcal{I}_{0}$ and at the end
$\mathcal{I}_{f}$. Assuming a linear dependence between the scattered
intensity and the thickness of the bed, the thickness of the floating lid is
obtained through:

$\delta(t)=\delta_{0}+(\delta_{th}-\delta_{0}).\frac{\mathcal{I}_{0}-\mathcal{I}(t)}{\mathcal{I}_{0}-\mathcal{I}_{f}}$
(26)

We then compare the numerical simulations of the floating lid thickness
evolution with the experimental results, as illustrated in Figure 4, right
panels. We note that contrary to the thermal case, there is only one kind of
evolution of the lid thickness: a transient period during which the lid
remains constant, followed by a delayed monotonic erosion of the lid until it
reaches the thickness at steady state.
The good agreement between experimental data and model prediction shown in
Figure 4 is verified in all experiments if $c_{e}$ is allowed to vary between
$0.06$ and $2$, leading to an average value $c_{e}=1.0\pm 0.8$, which confirms
that our model encompasses the main physical characteristics of the system.
Data do not allow us to identify a clear second order dependence of $c_{e}$,
for instance with respect to $Ra_{H}$. We thus conclude that the large
variability of $c_{e}$ mainly reflects the sensitivity of the model on the
determination of $T_{\rm{lid}}$ and $\delta_{\rm{th}}$. We further note that
here we describe the erosion of the bed thanks to a simplified 1D model,
whereas erosion induces a topography at the lid interface. Hence the variation
of $c_{e}$ could reflect this local bias in the estimate of $\delta$ for
example. Moreover, other parameters usually affects the erosion mechanism,
such as the packing of the bed [Agudo Wierschem (2012)] and we further measure
a dispersion in beads radius (see Sturtz21, supplementary materials) that can
also induce a variability in the determination of the Shields number
$\zeta_{\rm{lid}}$.
In both cases illustrated in Figure 4, we note a small increase of the lid
thickness before the erosion begins. This transient thickening can be
explained by several factors. At the onset of convection, under the action of
shear stress, if the bed is compacted enough, beads are not free to move one
around the others. Consequently, the first stage of evolution is a dilatation
(i.e. lid thickneing) that enables the motion of beads [Reynolds (1885),
Behringer Chakraborty (2018)]. Another explanation can come from the mechanism
of entrainment. As described by Solomatov93a, before being entrained, beads
are moving through the lid interface to form dune-shape cusps. This can lead
beads to reach the wall and accumulate in the upper corners of the tank. The
formation of this small meniscus of beads at the beginning of our experiment
could explain the apparent thickening of the lid which might be an artefact in
this case.

##  5 Deposition dynamics

Figure 5: Deposition dynamics in experiment IHB16. (a) Snapshot of the final
state of the experiment, showing a stable floating lid and a cumulate. The red
vertical line shows the one-pixel slice that is followed as a function of time
(b). The cumulate thickness is fitted by an exponential function (dashed
line).

In some experiments we performed, not only is the lid eroded, but also a
cumulate forms by sedimentation of the beads from the convective suspension
(see Table 1). Indeed the formation of a cumulate is expected when the bulk
temperature is high enough for the beads to become negatively buoyant and are
able to settle at the base of the tank. The formation of the cumulate always
happens in a late stage of the experiments, once the lid has reached its
steady state thickness. The evolution of the lid and of the cumulate can then
be modeled independently. However, as highlighted by Sturtz21, the criterion
for the cumulate formation is still related to the Shields number, that must
be subcritical in order to allow beads settle.
To model the formation of the cumulate, we consider a convective suspension
with a solid fraction $\phi(t)$ that can only decrease by particle deposition,
a valid hypothesis in our case as the floating lid has reached its steady
state thickness. We express the mass conservation of beads in the bulk as:

$\frac{\rm{d}V_{d}}{\rm{d}t}=-\frac{\rm{d}V_{\rm{sus}}}{\rm{d}t},$ (27)

where $V_{d}\sim\phi_{RLP}\delta_{d}$ is the volume of particles in the
cumulate of packing $\phi_{\rm{RLP}}$ and $V_{\rm{sus}}\sim\overline{\phi}h$
is the volume of particles in suspension with $\overline{\phi}$ the volume
averaged concentration of particles in the bulk.From the expressions of
$V_{\rm{sus}}$ and $V_{d}$ we get:

$\frac{\rm{d}\delta_{d}}{\rm{d}t}\sim-\frac{h}{\phi_{\rm{RLP}}}\,\frac{\rm{d}\overline{\phi}}{\rm{d}t},$
(28)

To go further, we use the model of Lavorel09 to determine $\phi(t)$. We
consider that beads can leave the bulk only if the Shields number is smaller
than its critical value, and do so at a settling velocity $v_{s}$ to form the
deposit. As the bulk Shields number is sub-critical, the cumulate is stable
and we consider that no re-entrainment occurs from the cumulate formed at the
bottom of the tank. All these assumptions put together leads to the evolution
equation of the volume averaged concentration of particles
$\overline{\phi}(t)$:

$\frac{\rm{d}\overline{\phi}}{\rm{d}t}=-\frac{v_{s}}{h}\overline{\phi}.$ (29)

Finally, using (27) and (29) we get the deposition law:

$\frac{\rm{d}\delta_{d}}{\rm{d}t}=c_{d}v_{s}\frac{\overline{\phi}}{\phi_{RLP}}.$
(30)

where $c_{d}$ depends on bead shape.
To further propose a first order solution of (30), we consider an initial
concentration of beads in suspension $\phi_{0}$ before deposition, which
yields:

$\phi=\phi_{0}\,\frac{h}{h-\delta_{d}}-\phi_{RLP}\frac{\delta_{d}}{h}.$ (31)

As the floating lid is supposed to be at steady state before cumulate starts
to form, $\phi_{0}$ is linked to the difference between the initial and the
final lid thicknesses. The solution of (30) is an exponential:
$\delta_{d}\sim\rm{exp}(-t/\tau_{s})$ where the deposition timescale is:

$\tau_{s}=\frac{1}{c_{d}}\frac{h}{v_{s}}.$ (32)

The problem thus boils down to the determination of the settling velocity
$v_{s}$. To that aim, we analyse the experiments where a stable deposit
formed, as illustrated in Figure 5. We choose a vertical slice of the frames
at position $x$ (vertical blue line in (a)), and we represent the evolution of
this slice over time (b) as a proxy of the evolution of the cumulate. As
illustrated in Figure 5 (b), the cumulate growth can be fitted by the
exponential law proposed above. Following several authors [Martin Nokes
(1988), Martin Nokes (1989), Lavorel Le Bars (2009)], we assume that the
settling velocity scales with the Stokes velocity: $v_{s}\sim\Delta\rho
gr^{2}/\eta_{f}$. Using the four experiments with a sufficient cumulate
thickness (the 3 other experiments produced a too thin deposit), we obtain
$c_{d}=0.24\pm 0.14$. The value of $c_{d}$ is close to the coefficient
$2/9\approx 0.22$ that appears in the expression of the Stokes velocity which
means that we do not have to introduce correcting a factor on the viscosity
due to convection. The dispersion is due to the inhomogeneity of the cumulate
thickness. To reduce the uncertainty, we take 3 values per experiments at
different position $x$. The beads radius distribution contributes also to the
dispersion (see Sturtz21, supplementary materials). Again, the dispersion can
be due to the distribution in bead radius that affects the sedimentation
velocity.

##  6 Conclusion

We develop an erosion/deposition model for particles sheared by convection by
adapting the Shields’ formalism to convective systems. We model the erosion
process in granular beds that lie in a thermal boundary layer. The erosion law
is a function of the Rayleigh-Roberts number that characterizes convection,
and the Shields number that quantifies the ability of beads to be freed from
frictional interactions. Within the same framework, we further established and
described the subsequent formation of basal cumulates.
The model has been validated using lab-scale experiments in a fluid bearing
plastic beads and heated from within. The theoretical framework and the
scaling laws proposed here can be applied to various geophysical systems, such
as magma oceans. Such applications are the subject of a companion paper.

## Declaration of interest

The authors report no conflicts of interest.

###### Acknowledgements.

This paper is part of Cyril Sturtz’s PhD thesis (Université de Paris, Institut
de Physique du Globe de Paris). The authors would like to thank Éric
Lajeunesse and Claude Jaupart for fruitful discussions that have made this
article possible. This study contributes to the IdEx Université de Paris
ANR-18-IDEX-0001. This work was supported by the Programme National de
Planétologie (PNP) of CNRS/INSU, co-funded by CNES.

##  Appendix A Experimental conditions

Properties | Symbol | Value | Unit
---|---|---|---
Fluid density ($20^{o}\rm{C}$) | $\rho_{0,f}$ | $1192$ | $\mathrm{kg\ m^{-3}}$
Beads density ($20^{o}\rm{C}$) | $\rho_{0,p}$ | $1187$ | $\mathrm{kg\ m^{-3}}$
Fluid thermal expansion | $\alpha_{f}$ | $5.5\ 10^{-4}$ | $\mathrm{K^{-1}}$
Beads thermal expansion | $\alpha_{p}$ | $3.2\ 10^{-4}$ | $\mathrm{K^{-1}}$
Fluid viscosity ($20^{o}\rm{C}$) | $\eta_{f}$ | $0.151$ | Pa s
Activation energy | $E_{a}$ | $41.7$ | $\mathrm{kJ\ mol^{-1}}$
Fluid thermal diffusivity | $\kappa_{f}$ | $9.1\ 10^{-8}$ | $\mathrm{m^{2}\ s^{-1}}$
Beads thermal diffusivity (*) | $\kappa_{p}$ | $1.\ 10^{-7}$ | $\mathrm{m^{2}\ s^{-1}}$
Fluid thermal conductivity | $\lambda_{f}$ | $0.276$ | $\mathrm{W\ m^{-1}\ K^{-1}}$
Beads thermal conductivity (*) | $\lambda_{p}$ | $0.21$ | $\mathrm{W\ m^{-1}\ K^{-1}}$
Table 2: Main physical properties of the fluid and beads. The activation
energy is obtained from the viscosity fit with an Arrhenius law:
$\eta(T)=\eta_{f}\exp\left[\dfrac{E_{a}}{R}\left(\dfrac{1}{T}-\dfrac{1}{T_{0}}\right)\right]$,
with $T_{0}=20^{o}\rm{C}$. Properties are all measured in the lab, except
those marked with (*) which are taken from [Mark (2007)]. See Supplementary
Materials for further information on the way properties measurements have been
carried on.

##  Appendix B Continuity of the flux at the base of the floating lid

We compare here that the flux out of the convective bulk to the . Following
Kerr90a,Kerr90b, the difference between the two fluxes is given by the
conservation of heat at the moving interface:

$Q_{s,\rm{lid}}(t)=\rho
c_{p}(T_{\rm{bulk}}-T_{\rm{lid}})\,\frac{\rm{d}\delta}{\rm{d}t}+Q_{s,\rm{conv}}(t).$
(33)

We take the dimensionless form of this equation by using the characteristic
flux $Hh$. Using the scaling law for convection
$T_{\rm{bulk}}-T_{\rm{lid}}\sim\Delta T_{H}Ra_{H}^{*,-1/4}$ and the evolution
equation of the lid thickness given by (23):

$Q_{s,\rm{lid}}^{*}(t)=\frac{r}{h}(\zeta-\zeta_{c})\,Ra_{H}^{*,1/4}+Q_{s,\rm{conv}}^{*}(t),$
(34)

where $Q_{s,i}^{*}$ stands for dimensionless fluxes. The first term in the
right-side of this equation equals $0.01-0.09$ in our experiments, and that is
why we assume that $Q_{s,\rm{lid}}(t)\approx Q_{s,\rm{conv}}(t)$.

##  Appendix C Numerical scheme for the experimental set-up

To simulate numerically what happens experimentally, we discretized the system
of equations over $N_{t}=6000$ points in time. Thus, the time at the time-step
$i$ is written $t_{i}=i\times\rm{d}t$. Similarly, all parameter $\mathcal{P}$
at $t_{i}$ is call $\mathcal{P}_{i}$. At each step, we deduce the evolution of
the system by the following numerical scheme, assuming that all parameters are
determined at $t_{i-1}$:

  1. 1.

Mechanical evolution of the bed :
if $\zeta_{i}>\zeta_{c}$, the bed is eroded, then we calculate the rate of
erosion by using the following explicit scheme of the equation (23):

$\delta_{i+1}=\delta_{i}+\frac{\rm{d}\delta}{\rm{d}t}\times\rm{d}t,$ (35)

with $\rm{d}\delta/\rm{d}t$ given by (23). For the sake of convenience, we can
re-express the (23) as a function of the lid thickness only. To do so, we re-
express $\zeta_{\rm{lid}}$ and $\zeta_{c}$ as a function of the $\delta$. In
this way, we assess that the scaling law determined by Sturtz21 for
$\zeta_{\rm{lid}}$ remains true in transient state, replacing $Ra_{H}$ by
$Ra_{H}^{*}$. Besides, as the bed is thin compared to the tank thickness, the
temperature profile in the bed can be assessed to be linear at any time, so
that: $T_{\rm{lid}}=T_{s}+Q_{s}(t)\delta(t)/\lambda_{b}$. This assumption is
justified by comparing conductive timescale
$\tau_{\rm{cond}}\sim\delta^{2}/\kappa_{b}$ to the convective one
$\tau_{\rm{conv}}\sim h^{2}/\kappa_{f}\,Ra_{H}^{-1/4}$ [Limare . (2019)]. We
get $\tau_{\rm{cond}}/\tau_{\rm{conv}}\ll 1$, which corroborates that the lid
is at thermal equilibrium at any time. Moreover,
$T_{\rm{lid}}=T_{bulk}-C_{T}^{*}\,\Delta T_{H}^{*}\,Ra_{H}^{*,-1/4}$ and
$\lambda_{b}$ the floating lid thermal conductivity. Consequentely, (17)
becomes:

$\zeta_{\rm{lid}}=\frac{\zeta_{s}}{1-\frac{\Delta(\rho_{0}\alpha)}{\Delta\rho_{0}}\frac{\delta(t)Q_{s}(t)}{\lambda_{b}}}.$
(36)

Similarly, $\zeta_{c}$ is expressed as a function of $\delta_{\rm{th}}$, the
thickness of the bed at steady state:

$\zeta_{c}=\frac{\zeta_{s}}{1-\frac{\Delta(\rho_{0}\alpha)}{\Delta\rho_{0}}\frac{\delta_{\rm{th}}Hh}{\lambda_{b}}}.$
(37)

The erosion of the bed (23) is therefore expressed as a function of the
thickness of the bed at steady state $\delta_{\rm{th}}$, which can be
calculated from the thermal steady state of each experiments by combining (5)
and assuming that the heat flux through the floating lid is given by
$Hh=\lambda_{b}(T_{\rm{lid}}-T_{s})/\delta_{\rm{th}}$, [Sturtz . (2021)]:

$\frac{\delta_{\rm{th}}}{h}=\frac{\lambda}{\lambda_{f}}\,\left(\frac{T_{\rm{bulk}}-T_{s}}{\Delta
T_{H}}-C_{T}\,Ra_{H}^{-1/4}\right),$ (38)

where $T_{\rm{bulk}}$ is the volume average temperature of the convective bulk
at steady state.

If $\zeta_{i}<\zeta_{c}$, deposition is possible.
As the lid does not absorbe microwave, we calculate the real rate of internal
heating $H_{i+1}=H_{0}\times h/\delta_{i+1}$.

  2. 2.

Basal temperature of the floating lid:
we use an explicit scheme of (6):

$T_{\rm{bulk,i+1}}=T_{\rm{bulk},i}+\frac{1}{\rho_{\rm{f}}c_{p,\rm{f}}}\left(H_{i}-\frac{Q_{s,\rm{conv},i}}{h-\delta_{i}}\right)\times\rm{d}t.$
(39)

  3. 3.

Thermal state of the bed :

     * •

If the lid thickness does not change: we compute equations (11)-(12) using a 1
dimension implicit finite difference scheme. The lid is discretized in
$N_{z}=500$ points, and the position in the lid is called
$z_{j}=j\times\rm{d}z$. The temperature field in the bed at time-step $i$ is
represented by the $N_{z}$-long vector $\mathcal{T}_{i}$. The temperature at
position $z_{j}$ and time $t_{i}$ is noted $T_{b,i}^{j}$. Thus, the
temperature in the floating lid at the time-step $t_{i+1}$ can be calculated
as follow :

$\mathcal{T}_{i+1}=\mathcal{A}^{-1}_{i}\cdot\mathcal{T}_{i},$ (40)

where

$\mathcal{A}_{i}=\begin{pmatrix}\epsilon&0&\cdots&&&0\\\
\hline\cr-\epsilon&C&-\epsilon&0&\cdots&0\\\
\vdots&\ddots&\ddots&\ddots&&\vdots\\\ \vdots&&\ddots&\ddots&\ddots&\vdots\\\
0&\cdots&0&-\epsilon&C&-\epsilon\\\ \hline\cr
0&\cdots&0&-\epsilon&-4\epsilon&-3\epsilon\\\ \end{pmatrix}\hskip
14.22636pt\rm{and}\hskip
14.22636pt\mathcal{T}_{i}=\begin{pmatrix}T_{s}\times\epsilon\\\ \hline\cr
T_{b,i}^{2}\\\ \vdots\\\ T_{b,i}^{j}\\\ \vdots\\\ T_{b,i}^{N_{z}-1}\\\
\hline\cr 2\rm{d}z\,\frac{Q_{s,\rm{conv},i}}{\lambda_{b}}\,\epsilon\\\
\end{pmatrix}$ (41)

where $\epsilon=\kappa_{b}\rm{d}t/\rm{d}z$, $C=1+2\epsilon$ and $T_{b,i}^{j}$
is the temperature field in the floating lid at the previous time-step. We
included the boundary conditions in the matrix. The Dirichlet condition at
$z_{j}=0$ is a fixed surface temperature : $T_{b,i}^{j=0}=T_{s}$. The Neumann
condition at the base of the floating lid is due to the fixed flux imposed by
the convection, and is calculated by using an order 2 approximation for the
temperature gradient at $z_{j=N_{z}}=\delta_{i+1}$ :

$\left[\frac{\rm{d}T_{b}}{\rm{d}z}\right]_{N_{z}}\approx\frac{T_{b,i}^{Nz-2}-4T_{b,i}^{Nz-1}+3T_{b,i}^{Nz}}{\rm{d}z}=\frac{Q_{s,\rm{conv},i}}{\lambda_{f}}.$
(42)

  4. 4.

If the lid thickness has evolved: we adapted the grid discretizing the lid in
order to always be regular and have $N_{z}=500$ points.
If the lid is eroded, we use only the temperature field at $t_{i}$.
Necessarily, the thickness $\delta_{i+1}$ is bracketed between
$z_{j_{max}}<\delta_{i+1}$ and $z_{j_{max}}>\delta_{i+1}$. Thus, we can
estimated the temperature at $\delta_{i+1}$ by an interpolation of the
temperature field between $z_{j_{max}}$ and $z_{j_{max}}$.
If the lid thickens, we add one row at $T^{i}$, corresponding to the position
$\delta_{i+1}$ and the corresponding temperature is estimated by an
extrapolation of $T_{b,Nz}^{i}$ using the basal heat flux. Then, this field is
resample on a regular grid of size $N_{z}$.
Once the lid evolves thermally, we deduce the basal temperature of the lid
$T_{lid,i+1}=T_{b,i+1}^{N_{z}}$. If the lid is totally eroded,
$T_{lid,i+1}=T_{s}$.

  5. 5.

Calculation of the convective heat flux : we use the scaling law (10) :

$Q_{s,\rm{conv},i+1}=\frac{\lambda_{f}}{C^{*,4/3}_{T}}\,\left(\frac{\alpha_{f}\rho_{f}g}{\kappa_{f}\eta_{f}}\right)^{1/3}\left(T_{\rm{bulk},i+1}-T_{\rm{lid},i+1}\right)^{4/3}$
(43)

  6. 6.

Calculation of all parameters :
We calculate all dimensionless numbers that can be important : the Rayleigh-
Roberts number $Ra_{H,i+1}^{*}$, and the Shields number $\zeta_{i+1}$.

## References

  * Agudo Wierschem (2012) Agudo12Agudo, JR. Wierschem, A.  2012\.  Incipient motion of a single particle on regular substrates in laminar shear flow Incipient motion of a single particle on regular substrates in laminar shear flow.  Physics of Fluids249093302.  https://doi.org/10.1063/1.4753941 10.1063/1.4753941
  * Behringer Chakraborty (2018) Behringer18Behringer, RP. Chakraborty, B.  2018nov.  The physics of jamming for granular materials: a review The physics of jamming for granular materials: a review.  Reports on Progress in Physics821012601.  https://doi.org/10.1088/1361-6633/aadc3c 10.1088/1361-6633/aadc3c
  * Bryson . (2019) Bryson19Bryson, J., Neufeld, J. Nimmo, F.  2019\.  Constraints on asteroid magnetic field evolution and the radii of meteorite parent bodies from thermal modelling Constraints on asteroid magnetic field evolution and the radii of meteorite parent bodies from thermal modelling.  Earth and Planetary Science Letters52168-78.
  * Charru . (2004) Charru04Charru, F., Mouilleron, H. Eiff, O.  2004\.  Erosion and deposition of particles on a bed sheared by a viscous flow Erosion and deposition of particles on a bed sheared by a viscous flow.  Journal of Fluid Mechanics51955-80.  10.1017/S0022112004001028
  * Elkins-Tanton . (2011) ElkinsTanton11Elkins-Tanton, L., Burgess, S. Yin, QZ.  2011April.  The lunar magma ocean: Reconciling the solidification process with lunar petrology and geochronology The lunar magma ocean: Reconciling the solidification process with lunar petrology and geochronology.  Earth and Planetary Science Letters3043326-336.
  * Fourel . (2017) Fourel17Fourel, L., Limare, A., Jaupart, C., Surducan, E., Farnetani, CG., Kaminski, EC.Surducan, V.  2017\.  The Earth’s mantle in a microwave oven: thermal convection driven by a heterogeneous distribution of heat sources The earth’s mantle in a microwave oven: thermal convection driven by a heterogeneous distribution of heat sources.  Experiments in Fluids588.  https://doi.org/10.1007/s00348-017-2381-3 10.1007/s00348-017-2381-3
  * Hamano . (2013) Hamano13Hamano, K., Abe, Y. Genda, H.  2013\.  Emergence of two types of terrestrial planet on solidification of magma ocean Emergence of two types of terrestrial planet on solidification of magma ocean.  Nature4977451607-610.
  * Helz . (1989) Helz89Helz, R., Kirschenbaum, H. Marinenko, J.  198904\.  Diapiric transfer of melt in Kilauea Iki lava lake, Hawaii: A quick, efficient process of igneous differentiation Diapiric transfer of melt in kilauea iki lava lake, hawaii: A quick, efficient process of igneous differentiation.  Geological Society of America Bulletin - GEOL SOC AMER BULL101578-594.  10.1130/0016-7606(1989)101¡0578:DTOMIK¿2.3.CO;2
  * Jaupart . (2007) Jaupart07Jaupart, C., Molnar, P. Cottrell, E.  2007\.  Instability of a chemically dense layer heated from below and overlain by a deep less viscous fluid Instability of a chemically dense layer heated from below and overlain by a deep less viscous fluid.  Journal of Fluid Mechanics572433-469.
  * Jellinek Kerr (2001) Jellinek01Jellinek, A. Kerr, R.  2001\.  Magma dynamics, crystallization, and chemical differentiation of the 1959 Kilauea Iki lava lake, Hawaii, revisited Magma dynamics, crystallization, and chemical differentiation of the 1959 kilauea iki lava lake, hawaii, revisited.  Journal of Volcanology and Geothermal Research1103235 - 263.  http://www.sciencedirect.com/science/article/pii/S0377027301002128 https://doi.org/10.1016/S0377-0273(01)00212-8
  * Kaminski . (2020) Kaminski20Kaminski, E., Limare, A., Kenda, B. Chaussidon, M.  2020\.  Early accretion of planetesimals unraveled by the thermal evolution of parent bodies of magmatic iron meteorites Early accretion of planetesimals unraveled by the thermal evolution of parent bodies of magmatic iron meteorites.  Earth and Planetary Science Letters548116469.  https://doi.org/10.1016/j.epsl.2020.116469
  * Kerr . (19901) Kerr90aKerr, R., Woods, A., Worster, M. Huppert, H.  19901\.  Solidification of an alloy cooled from above. Part 1. Equilibrium growth. Solidification of an alloy cooled from above. part 1. equilibrium growth.  Journal of Fluid Mechanics216323-342.  10.1017/S0022112090000453
  * Kerr . (19902) Kerr90bKerr, R., Woods, A., Worster, M. Huppert, H.  19902\.  Solidification of an alloy cooled from above. Part 2. Non-equilibrium interfacial kinetics. Solidification of an alloy cooled from above. part 2. non-equilibrium interfacial kinetics.  Journal of Fluid Mechanics217331-348.  10.1017/S002211209000074X
  * Lajeunesse . (2010) Lajeunesse10Lajeunesse, E., Malverti, L. Charru, F.  2010\.  Bed load transport in turbulent flow at the grain scale: Experiments and modeling Bed load transport in turbulent flow at the grain scale: Experiments and modeling.  Journal of Geophysical Research: Earth Surface115F4.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/2009JF001628 10.1029/2009JF001628
  * Lavorel Le Bars (2009) Lavorel09Lavorel, G. Le Bars, M.  200910\.  Sedimentation of particles in a vigorously convecting fluid Sedimentation of particles in a vigorously convecting fluid.  Physical review. E, Statistical, nonlinear, and soft matter physics80046324.  10.1103/PhysRevE.80.046324
  * Limare . (2019) Limare19Limare, A., Jaupart, C., Kaminski, E., Fourel, L. Farnetani, C.  2019\.  Convection in an internally heated stratified heterogeneous reservoir Convection in an internally heated stratified heterogeneous reservoir.  Journal of Fluid Mechanics87067–105.  10.1017/jfm.2019.243
  * Limare . (2021) Limare21Limare, A., Kenda, B., Kaminski, E., Surducan, E., Surducan, V. Neamtu, C.  2021\.  Transient convection experiments in internally-heated systems Transient convection experiments in internally-heated systems.  MethodsX8101224.
  * Limare . (2013) Limare13Limare, A., Surducan, E., Surducan, V., Neamtu, C., di Giuseppe, E., Vilella, K.Jaupart, C.  2013\.  Microwave-based laboratory experiments for internally-heated mantle convection Microwave-based laboratory experiments for internally-heated mantle convection.  AIP Conference Proceedings1565114-18.  https://aip.scitation.org/doi/abs/10.1063/1.4833687 10.1063/1.4833687
  * Limare . (2015) Limare15Limare, A., Vilella, K., Di Giuseppe, E., Farnetani, CG., Kaminski, E., Surducan, E.et al.  2015\.  Microwave-heating laboratory experiments for planetary mantle convection Microwave-heating laboratory experiments for planetary mantle convection.  Journal of Fluid Mechanics77750–67.  10.1017/jfm.2015.347
  * Mandler Elkins-Tanton (2013) Mandler13Mandler, B. Elkins-Tanton, L.  2013\.  The origin of eucrites, diogenites, and olivine diogenites: Magma ocean crystallization and shallow magma chamber processes on Vesta The origin of eucrites, diogenites, and olivine diogenites: Magma ocean crystallization and shallow magma chamber processes on vesta.  Meteoritics & Planetary Science48112333-2349.  https://onlinelibrary.wiley.com/doi/abs/10.1111/maps.12135 10.1111/maps.12135
  * Mark (2007) HandbookPolymerMark, J. ().  2007\.  Physical Properties of Polymers Handbook, 2nd edition Physical properties of polymers handbook, 2nd edition ( 29).  Springer.  10.1007/978-0-387-69002-5
  * Martin Nokes (1988) MartinNokes88Martin, D. Nokes, R.  1988\.  Crystal settling in a vigorously convecting magma chamber Crystal settling in a vigorously convecting magma chamber.  Nature3326164534-536.  https://doi.org/10.1038/332534a0 10.1038/332534a0
  * Martin Nokes (1989) MartinNokes89Martin, D. Nokes, R.  198912\.  A Fluid-Dynamical Study of Crystal Settling in Convecting Magmas A Fluid-Dynamical Study of Crystal Settling in Convecting Magmas.  Journal of Petrology3061471-1500.  https://doi.org/10.1093/petrology/30.6.1471 10.1093/petrology/30.6.1471
  * Maurice . (2020) Maurice20Maurice, M., Tosi, N., Schwinger, S., Breuer, D. Kleine, T.  2020\.  A long-lived magma ocean on a young Moon A long-lived magma ocean on a young moon.  Science Advances628.
  * Namur . (2011) Namur11Namur, O., Charlier, B., Pirard, C., Hermann, J., Liégeois, JP. Auwera, J.  201109\.  Anorthosite formation by plagioclase flotation in ferrobasalt and implications for the lunar crust Anorthosite formation by plagioclase flotation in ferrobasalt and implications for the lunar crust.  Geochimica et Cosmochimica Acta754998.  10.1016/j.gca.2011.06.013
  * Reynolds (1885) Reynolds85Reynolds, O.  1885\.  LVII. On the dilatancy of media composed of rigid particles in contact. With experimental illustrations Lvii. on the dilatancy of media composed of rigid particles in contact. with experimental illustrations.  Philosophical Magazine Series 120469-481.
  * Righter Drake (1997) Righter97Righter, K. Drake, M.  199702\.  A magma ocean on Vesta: Core formation and petrogenesis of eucrites and diogenites A magma ocean on vesta: Core formation and petrogenesis of eucrites and diogenites.  Meteoritics & Planetary Science32929 - 944.  10.1111/j.1945-5100.1997.tb01582.x
  * Roberts (1967) Roberts67Roberts, P.  1967\.  Convection in horizontal layers with internal heat generation. Theory Convection in horizontal layers with internal heat generation. theory.  Journal of Fluid Mechanics30133-49.  10.1017/S0022112067001284
  * Shields (1936) Shields36Shields, A.  1936\.  Anwendung der Aehnlichkeitsmechanik und der Turbulenzforschung auf die Geschiebebewegung, Mitteilung Preussischen Versuchsanstalt Wasserbau Schiffbau [Application of similarity mechanics and turbulence research on shear flow] Anwendung der aehnlichkeitsmechanik und der turbulenzforschung auf die geschiebebewegung, mitteilung preussischen versuchsanstalt wasserbau schiffbau [application of similarity mechanics and turbulence research on shear flow].  Preussischen Versuchsanstalt für Wasserbau und Schiffbau, Berlin26.
  * Solomatov (2000) Solomatov00Solomatov, V.  200001\.  Fluid dynamics of a terrestrial magma ocean Fluid dynamics of a terrestrial magma ocean.  Origin of the Earth and Moon1.
  * Solomatov . (1993) Solomatov93aSolomatov, V., Olson, P. Stevenson, D.  1993\.  Entrainment from a bed of particles by thermal convection Entrainment from a bed of particles by thermal convection.  Earth and Planetary Science Letters1203387 - 393\.  http://www.sciencedirect.com/science/article/pii/0012821X93902525 https://doi.org/10.1016/0012-821X(93)90252-5
  * Solomatov Stevenson (1993) Solomatov93bSolomatov, V. Stevenson, D.  1993\.  Suspension in convective layers and style of differentiation of a terrestrial magma ocean Suspension in convective layers and style of differentiation of a terrestrial magma ocean.  Journal of Geophysical Research: Planets98E35375-5390.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/92JE02948 10.1029/92JE02948
  * Sparks . (2019) Sparks19Sparks, RSJ., Annen, C., Blundy, JD., Cashman, KV., Rust, AC. Jackson, MD.  2019\.  Formation and dynamics of magma reservoirs Formation and dynamics of magma reservoirs.  Philosophical Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences377213920180019.  https://royalsocietypublishing.org/doi/abs/10.1098/rsta.2018.0019 10.1098/rsta.2018.0019
  * Sturtz . (2021) Sturtz21Sturtz, C., Kaminski, E., Limare, A. Tait, S.  2021\.  The fate of particles in a volumetrically heated convective fluid at high Prandtl number. The fate of particles in a volumetrically heated convective fluid at high prandtl number.  arXiv:2101.06706.  submitted
  * Suckale, Elkins-Tanton Sethian (2012) Suckale12bSuckale, J., Elkins-Tanton, L. Sethian, J.  201208\.  Crystals stirred up: 2. Numerical insights into the formation of the earliest crust on the Moon Crystals stirred up: 2. numerical insights into the formation of the earliest crust on the moon.  Journal of Geophysical Research117E08005.  10.1029/2012JE004067
  * Suckale, Sethian . (2012) Suckale12aSuckale, J., Sethian, J., Yu, Jd. Elkins-Tanton, L.  2012\.  Crystals stirred up: 1. Direct numerical simulations of crystal settling in nondilute magmatic suspensions Crystals stirred up: 1. direct numerical simulations of crystal settling in nondilute magmatic suspensions.  Journal of Geophysical Research117E08004.  10.1029/2012JE004066
  * Surducan . (2014) Surducan14Surducan, E., Surducan, V., Limare, A., Neamtu, C. Giuseppe, E.  2014\.  Microwave heating device for internal heating convection experiments, applied to Earth’s mantle dynamics Microwave heating device for internal heating convection experiments, applied to earth’s mantle dynamics.  Review of Scientific Instruments8512124702.  10.1063/1.4902323
  * Vilella Kaminski (2017) Vilella17Vilella, K. Kaminski, E.  2017\.  Fully determined scaling laws for volumetrically heated convective systems, a tool for assessing habitability of exoplanets Fully determined scaling laws for volumetrically heated convective systems, a tool for assessing habitability of exoplanets.  Physics of the Earth and Planetary Interiors26618 \- 28.  http://www.sciencedirect.com/science/article/pii/S0031920117300092 https://doi.org/10.1016/j.pepi.2017.02.001
  * Vilella . (2018) Vilella18Vilella, K., Limare, A., Jaupart, C., Farnetani, C., Fourel, L. Kaminski, E.  2018\.  Fundamentals of laminar free convection in internally heated fluids at values of the Rayleigh-Roberts number up to $10^{9}$ Fundamentals of laminar free convection in internally heated fluids at values of the rayleigh-roberts number up to $10^{9}$.  Journal of Fluid Mechanics846966–998.  10.1017/jfm.2018.316
  * Warren (1985) Warren85Warren, P.  1985\.  The magma ocean concept and lunar evolution The magma ocean concept and lunar evolution.  Ann. Rev. Earth Planet. Sci.13201-240.
  * Wood (1970) Wood70aWood, J.  1970\.  Petrology of the lunar soil and geophysical implications Petrology of the lunar soil and geophysical implications.  Journal of Geophysical Research (1896-1977)75326497-6513.  https://agupubs.onlinelibrary.wiley.com/doi/abs/10.1029/JB075i032p06497 10.1029/JB075i032p06497
  * Wood . (1970) Wood70bWood, J., Dickey, J., Marvin, U. Powell, B.  1970\.  Lunar anorthosites and a geophysical model of the moon Lunar anorthosites and a geophysical model of the moon.  Proceedings of the Apollo 11 Lunar Science Conference Proceedings of the apollo 11 lunar science conference ( 1, 965-988).


* Gómez González  _et al._ (2023) R. Gómez González, E. Abad, S. Bravo Yuste, and V. Garzó, “Diffusion of intruders in granular suspensions: Enskog theory and random walk interpretation,” arXiv:2305.09259 (2023).
  * Bird (1994) G. A. Bird,  _Molecular Gas Dynamics and the Direct Simulation Monte Carlo of Gas Flows_ (Clarendon, Oxford, 1994).
  * Gómez González and Garzó (2022a) R. Gómez González and V. Garzó, “Kinetic theory of binary granular suspensions at low density. Thermal diffusion segregation,” in  _Nonequilibrium Thermodynamics and Fluctuation Kinetics_ , edited by L. Brenig, N. Brilliantov, and M. Tlidi (Fundamental Theories of Physics 208, Springer, 2022).
  * Williams and MacKintosh (1996) D. R. M. Williams and F. C. MacKintosh, “Driven granular media in one dimension: Correlations and equation of state,” Phys. Rev. E 54, R9–R12 (1996).
  * van Noije and Ernst (1998) T. P. C. van Noije and M. H. Ernst, “Velocity distributions in homogeneous granular fluids: the free and heated case,” Granular Matter 1, 57–64 (1998).
  * Van der Hoef, Beetstra, and Kuipers (2005) M. A. Van der Hoef, R. Beetstra, and J. A. M. Kuipers, “Lattice-Boltzmann simulations of low-Reynolds number flow past mon- and bidisperse arrays of spheres: resutls for the permeability drag force,” J. Fluid Mech. 528, 233–254 (2005).
  * Beetstra, Van der Hoef, and Kuipers (2007) R. Beetstra, M. A. Van der Hoef, and J. A. M. Kuipers, “Erratum: Drag force of intermediate Reynolds number flow past mon- and bidisperse arrays of spheres [AIChE J. 53, 489-591 (2007)],” AIChE J. 53, 3020 (2007).
  * Yin and Sundaresan (2009) X. Yin and S. Sundaresan, “Fluid-particle drag in low-Reynolds-number polydisperse gas-solid suspensions,” AIChE 55, 1352–1368 (2009).
  * van Kampen (2007) N. G. van Kampen,  _Stochastic Processes in Physics and Chemistry_ (North Holland, Amsterdam, 2007).
  * Gómez González and Garzó (2021) R. Gómez González and V. Garzó, “Time-dependent homogeneous states of binary granular suspensions,” Phys. Fluids 33, 093315 (2021).
  * Garzó (2019) V. Garzó,  _Granular Gaseous Flows_ (Springer Nature, Cham, 2019).
  * Note (1) Note that there is a mistake in Eq. (35) for the coefficient $D^{*}$. Equation (77) corrects this mistake.
  * Mason (1954) E. A. Mason, “Transport properties of gases obeying a modified Buckingham potential,” J. Chem. Phys. 22, 169–192 (1954).
  * López de Haro and Cohen (1984) M. López de Haro and E. G. D. Cohen, “The Enskog theory for multicomponent mixtures. III. Transport properties of dense binary mixtures with one tracer component,” J. Chem. Phys. 80, 408–415 (1984).
  * Carnahan and Starling (1969) N. F. Carnahan and K. E. Starling, “Equation of state for nonattracting rigid spheres,” J. Chem. Phys. 51, 635–636 (1969).
  * Boublik (1970) T. Boublik, “Hard-sphere equation of state,” J. Chem. Phys. 53, 471 (1970).
  * Reed and Gubbins (1973) T. M. Reed and K. E. Gubbins,  _Applied Statistical Mechanics_ (MacGraw-Hill, New York, 1973).
  * McLennan (1989) J. A. McLennan,  _Introduction to Nonequilibrium Statistical Mechanics_ (Prentic–Hall, New Yersey, 1989).
  * Montanero and Garzó (2002) J. M. Montanero and V. Garzó, “Monte Carlo simulation of the homogeneous cooling state for a granular mixture,” Granular Matter 4, 17–24 (2002).
  * Khalil and Garzó (2014) N. Khalil and V. Garzó, “Homogeneous states in driven granular mixtures: Enskog kinetic theory versus molecular dynamics simulations,” J. Chem. Phys. 140, 164901 (2014).
  * Rosato  _et al._ (1987) A. Rosato, K. J. Strandburg, F. Prinz, and R. H. Swendsen, “Why the Brazil nuts are on top: Size segregation of particulate matter by shaking,” Phys. Rev. Lett. 58, 1038–1040 (1987).
  * Duran, Rajchenbach, and Clément (1993) J. Duran, J. Rajchenbach, and E. Clément, “Arching effect model for particle size segregation,” Phys. Rev. Lett. 70, 2431–2434 (1993).
  * Cooke  _et al._ (1996) W. Cooke, S. Warr, J. M. Huntley, and R. C. Ball, “Particle size segregation in a two-dimensional bed undergoing vertical vibration,” Phys. Rev. E 53, 2812 (1996).
  * Shinbrot and Muzzio (1998) T. Shinbrot and F. J. Muzzio, “Reverse buoyancy in shaken granular beds,” Phys. Rev. Lett. 81, 4365–4368 (1998).
  * Hong, Quinn, and Luding (2001) D. C. Hong, P. V. Quinn, and S. Luding, “Reverse Brazil nut problem: Competition between percolation and condensation,” Phys. Rev. Lett. 86, 3423–3426 (2001).
  * Arnarson and Willits (1998) B. Arnarson and J. T. Willits, “Thermal diffusion in binary mixtures of smooth, nearly elastic spheres with and without gravity,” Phys. Fluids 10, 1324–1328 (1998).
  * Jenkins and Yoon (2002) J. T. Jenkins and D. K. Yoon, “Segregation in binary mixtures under gravity,” Phys. Rev. Lett. 88, 194301 (2002).
  * Trujillo, Alam, and Herrmann (2003) L. Trujillo, M. Alam, and H. J. Herrmann, “Segregation in a fluidized binary granular mixture: Competition between buoyancy and geometric forces,” Europhys. Lett. 64, 190–196 (2003).
  * Alam, Trujillo, and Herrmann (2006) M. Alam, L. Trujillo, and H. J. Herrmann, “Hydrodynamic theory for reverse Brazil nut segregation and the non-monotonic ascension dynamics,” J. Stat. Phys. 124, 587–623 (2006).
  * Yoon and Jenkins (2006) D. K. Yoon and J. T. Jenkins, “The influence of different species’ granular temperatures on segregation in a binary mixture of dissipative grains,” Phys. Fluids 18, 073303 (2006).
  * Garzó (2006) V. Garzó, “Segregation in granular binary mixtures: Thermal diffusion,” Europhys. Lett. 75, 521–527 (2006).
  * Garzó (2008) V. Garzó, “Brazil-nut effect versus reverse Brazil-nut effect in a moderately granular dense gas,” Phys. Rev. E 78, 020301 (R) (2008).
  * Garzó (2009) V. Garzó, “Segregation by thermal diffusion in moderately dense granular mixtures,” Eur. Phys. J. E 29, 261–274 (2009).
  * Galvin, Dahl, and Hrenya (2005) J. E. Galvin, S. R. Dahl, and C. M. Hrenya, “On the role of non-equipartition in the dynamics of rapidly flowing granular mixtures,” J. Fluid Mech. 528, 207–232 (2005).
  * Brey, Ruiz-Montero, and Moreno (2005) J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, “Energy partition and segregation for an intruder in a vibrated granular system under gravity,” Phys. Rev. Lett. 95, 098001 (2005).
  * Brey, Ruiz-Montero, and Moreno (2006) J. J. Brey, M. J. Ruiz-Montero, and F. Moreno, “Hydrodynamic profiles for an impurity in an open vibrated granular gas,” Phys. Rev. E 73, 031301 (2006).
  * Brito  _et al._ (2008) R. Brito, H. Enríquez, S. Godoy, and R. Soto, “Segregation induced by inelasticity in a vibrofluidized granular mixture,” Phys. Rev. E 77, 061301 (2008).
  * Brito and Soto (2009) R. Brito and R. Soto, “Competition of Brazil nut effect, buoyancy, and inelasticity induced segregation in a granular mixture,” Eur. Phys. J. Special Topics 179, 207–219 (2009).
  * Kincaid, Cohen, and López de Haro (1987) J. M. Kincaid, E. G. D. Cohen, and M. López de Haro, “The Enskog theory for multicomponent mixtures. iv. Thermal diffusion,” J. Chem. Phys. 86, 963–975 (1987).
  * Wildman, Huntley, and Parker (2001) R. D. Wildman, J. M. Huntley, and D. J. Parker, “Granular temperature profiles in three-dimensional vibrofluidized granular beds,” Phys. Rev. E 63, 061311 (2001).
  * Breu  _et al._ (2003) A. P. J. Breu, H. M. Ensner, C. A. Kruelle, and I. Rehberg, “Reversing the Brazil-nut effect: Competition between percolation and condensation,” Phys. Rev. Lett. 90, 014302 (2003).
  * Schautz  _et al._ (2005) T. Schautz, R. Brito, C. A. Kruelle, and I. Rehberg, “A horizontal Brazil-nut effect and its reverse,” Phys. Rev. Lett. 95, 028001 (2005).
  * Qiao  _et al._ (2021) J. Qiao, C. Duan, K. Dong, W. Wang, H. Jiang, H. Zhu, and Y. Zhao, “Dem study of segregation degree and velocity of binary granular mixtures subject to vibration,” Powder Technol. 382, 107–117 (2021).
  * Shao  _et al._ (2023) Y. Shao, A. Li, Z. Wang, M. Sun, and D. Huang, “Resonance-induced acceleration of the RBNE-BNE segregation inversion,” AIChE J. e18101 (2023).
  * Garzó (2011) V. Garzó, “Thermal diffusion segregation in granular binary mixtures described by the Enskog equation,” New J. Phys. 13, 055020 (2011).
  * Gómez González and Garzó (2022b) R. Gómez González and V. Garzó, “Kinetic theory of granular particles immersed in a molecular gas,” J. Fluid Mech. 943, A9 (2022b).
  * Khalil and Garzó (2013) N. Khalil and V. Garzó, “Transport coefficients for driven granular mixtures at low-density,” Phys. Rev. E 88, 052201 (2013).