<html><head><meta name="color-scheme" content="light dark"></head><body><pre style="word-wrap: break-word; white-space: pre-wrap;">\makeatletter
\@ifundefined{HCode}
{
\documentclass[CRPHYS, Unicode, screen, Thematic]{cedram}
\newenvironment{noXML}{}{}
\newenvironment{Table}{\begin{table}}{\end{table}}
\def\thead{\noalign{\relax}\hline}
\def\endthead{\noalign{\relax}\hline}
\def\tabnote#1{\vskip4pt\parbox{.85\linewidth}{#1}}
\def\tsup#1{$^{{#1}}$}
\def\tsub#1{$_{{#1}}$}
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\usepackage{amssymb}
\RequirePackage{etoolbox}
\usepackage{amsmath}
\def\botline{\\\hline}
\csdef{Seqnsplit}{\\}
\def\xxboldsymbol#1{\boldsymbol{#1}}
\def\tbody{\noalign{\relax}\hline}
\def\jobid{crphys20240722}
%\graphicspath{{/tmp/\jobid_figs/web/}}
\graphicspath{{./figures/}}
\newcounter{runlevel}
\let\MakeYrStrItalic\relax
\def\refinput#1{}
\def\0{\phantom{0}}
\usepackage{multirow}
\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\sfrac#1#2{{#1}/{#2}}
\def\stfrac#1#2{({#1})/{#2}}
\def\stofrac#1#2{(({#1})/{#2})}
\def\back#1{}
\DOI{10.5802/crphys.227}
\datereceived{2024-09-07}
\daterevised{2024-12-17}
\dateaccepted{2024-12-17}
\ItHasTeXPublished
\def\ndash{\text{--}}
\def\mathbi#1{\text{\textit{\textbf{#1}}}}
}
{\documentclass[crphys]{article}
\def\CDRdoi{10.5802/crphys.227}
\let\newline\break
\def\sfrac#1#2{{#1}/{#2}}
\def\sofrac#1#2{({#1}/{#2})}
\makeatletter
\def\selectlanguage#1{} 
\def\href#1#2{\url[#1]{#2}}
}
\makeatother

\usepackage{mathrsfs}

\def\xxboldsymbol#1{#1}
\usepackage{upgreek}

\dateposted{2025-01-13}
\begin{document}

\begin{noXML}

%\makeatletter
%\def\TITREspecial{\relax}
%\def\cdr@specialtitle@english{Geophysical and astrophysical fluid dynamics in the laboratory}
%\def\cdr@specialtitle@french{Dynamique des fluides g\'eophysiques et astrophysiques aulaboratoire}
%\makeatother

\CDRsetmeta{articletype}{review}

\title{Fluid dynamics of planetary differentiation}

\alttitle{La diff\'erenciation des plan\`etes telluriques vue par la
m\'ecanique des fluides}

\author{\firstname{Renaud} \lastname{Deguen}\CDRorcid{0000-0002-6883-0655}\IsCorresp}
\address{ISTerre, Universit\'e Grenoble Alpes, Universit\'e Savoie Mont
Blanc, CNRS, IRD, Universit\'e Gustave Eiffel, 38000 Grenoble, France}
\email[R. Deguen]{renaud.deguen@univ-grenoble-alpes.fr}

\author{\firstname{Ludovic} \lastname{Huguet}\CDRorcid{0000-0002-5532-6302}}
\addressSameAs{1}{ISTerre, Universit\'e Grenoble Alpes, Universit\'e
Savoie Mont Blanc, CNRS, IRD, Universit\'e Gustave Eiffel, 38000
Grenoble, France}
\email[L. Huguet]{ludovic.huguet@univ-grenoble-alpes.fr}

\author{\firstname{Maylis} \lastname{Landeau}\CDRorcid{0000-0002-8836-4903}} 
\address{Universit\'e Paris Cit\'e, Institut de Physique du Globe de
Paris, CNRS, 75005 Paris, France}
\email[M. Landeau]{landeau@ipgp.fr}

\author{\firstname{Victor} \lastname{Lherm}\CDRorcid{0000-0001-5814-0637}}
\address{{Universit\'e Paris-Saclay, CNRS, FAST, 91405, Orsay, France}}
\email[V. Lherm]{victor.lherm@cnrs.fr}

\author{\firstname{Augustin} \lastname{Maller}\CDRorcid{0000-0001-9788-6239}}
\addressSameAs{2}{Universit\'e Paris Cit\'e, Institut de Physique du
Globe de Paris, CNRS, 75005 Paris, France}
\email[A. Maller]{maller@ipgp.fr}

\author{\firstname{Jean-Baptiste} \lastname{Wacheul}}
\address{Univ Lyon, ENSL, UCBL, UJM, CNRS, LGL-TPE, F-69007 Lyon, France}
\email[J.-B. Wacheul]{jb.wacheul@orange.fr}

\shortrunauthors

\thanks{French Academy of Science, Electricit\'e de France, European
Research Council (ERC) (grant number 716429), Labex OSUG@2020 (ANR10
LABX56), Programme National de Plan\'etologie (PNP) of CNRS-INSU
co-funded by CNES, Center for Matter at Atomic Pressure (CMAP), 
National Science Foundation (NSF) Physics Frontier Center under award
PHY-2020249.}

\begin{abstract}
The basic structure of the terrestrial planets---an iron-rich metallic
core surrounded by a silicate mantle---was established during their
accretion, when widespread melting allowed the metal and silicate
phases to separate.  The transfer of chemical elements and heat between
the metal and silicate that occurred during this period is critical for
the composition and initial temperature of the core and mantle, and has
important implications for the long-term evolution and dynamics of the
planets. After having summarised the main observational constraints on
core-mantle differentiation,  the article follows the sequence of
processes that led to the formation of planetary cores, focusing on the
contributions of laboratory fluid dynamics experiments to our
understanding of these processes, and discussing the relevance and
limitations of this approach to this problem.  We first focus on the
dynamics of planetary impacts, using laboratory experiments to
illustrate and quantify the impact and cratering processes and the
resulting metal phase dispersion. We then consider the two-phase flow
that follows an impact, when a molten impactor core falls by buoyancy
in a magma ocean.  The model of miscible turbulent thermal, which we
argue is a good reference model for the post-impact flow, is presented.
We then discuss how additional factors---immiscibility and
fragmentation, inertia inherited from the impact, Coriolis force,
sedimentation---affect the predictions of this model, and discuss the
extent of chemical equilibration. Finally, a last part of the article
is devoted to the migration of the metal phase through a solid part of
the mantle. 
\end{abstract}

\begin{altabstract}
La structure des plan\`etes terrestres --- un noyau m\'etallique riche en
fer entour\'e d'un manteau silicat\'e --- a \'et\'e \'etablie pendant
leur accr\'etion, lorsque la fusion g\'en\'eralis\'ee a permis la
s\'eparation des phases m\'etalliques et silicat\'ees.  Le transfert
d'\'el\'ements chimiques et de chaleur entre m\'etal et silicates qui
s'est produit pendant cette p\'eriode est d\'eterminant pour la
composition et la temp\'erature initiale du noyau et du manteau, et a
des implications importantes pour l'\'evolution et la dynamique \`a long
terme des plan\`etes. Apr\`es avoir r\'esum\'e les principales
contraintes observationnelles sur la diff\'erenciation noyau-manteau, 
l'article suit la s\'equence des processus qui ont conduit \`a la
formation des noyaux plan\'etaires, en se concentrant sur les
contributions de l'approche de dynamique des fluides exp\'erimentale \`a
notre compr\'ehension de ces processus, et en discutant de la
pertinence et des limites de cette approche.  Nous nous concentrons
d'abord sur la dynamique des impacts plan\'etaires, en utilisant des
exp\'eriences de laboratoire pour illustrer et quantifier les processus
d'impact et de crat\'erisation, et la dispersion de la phase
m\'etallique qui en r\'esulte. Nous examinons ensuite l'\'ecoulement
diphasique qui suit un impact, lorsque le noyau d'un impacteur chute
dans un oc\'ean de magma.  Nous introduisons le mod\`ele de thermique
turbulent, qui constitue selon nous un bon mod\`ele de r\'ef\'erence
pour l'\'ecoulement post-impact.  Nous discutons ensuite des facteurs
additionnels --- immiscibilit\'e et fragmentation, inertie h\'erit\'ee
de l'impact, force de Coriolis, s\'edimentation --- pouvant affecter les
pr\'edictions de ce mod\`ele. Enfin, la derni\`ere partie de l'article
est consacr\'ee \`a la migration de la phase m\'etallique \`a travers une
partie solide du manteau. 
\end{altabstract}

\keywords{\kwd{Early earth}
\kwd{Core formation}
\kwd{Fluid dynamics experiments}}

\altkeywords{\kwd{Terre primitive}
\kwd{formation du noyau}
\kwd{dynamique des fluides exp\'erimentale}}

\maketitle

\end{noXML}

\let\rmmu\upmu

\def\eqref#1{(\ref{#1})}

\section{Introduction}\label{sec1}

\subsection{Geological context}\label{sec:geological_context}

The terrestrial, or ``Earth-like'', planetary bodies all share the same
internal structure: an iron-rich metallic core surrounded by a silicate
(rocky) mantle. This is an expected outcome of planetary formation: the
metal and silicate phases being immiscible and of different densities,
they separated under the action of gravity to form the core and mantle.
Efficient gravitational separation requires that at least one of the
two phases is molten, a condition likely met during the accretion of
the terrestrial planets, about 4.56~Gy ago.  Planetary accretion models
indeed strongly suggest that the Earth started
hot~\cite{Stevenson2015}. Following the accretion of dust from the
solar nebula into planetesimals with diameters ranging from several
tenths of km to a few hundred 
km~\cite{morbidelli2009asteroids,Lunine2011}, the next 10 to 100
million years of the Solar System are thought to have been dominated by
increasingly energetic collisions between massive planetary bodies
several thousand~km in diameter~\cite{Lunine2011}. The dissipation into
heat of the kinetic energy of these collisions likely resulted in
widespread melting, leading to the formation of magma oceans in which
the metal and silicate phases separated.  A significant fraction of
these collisions likely involved bodies that were already
differentiated into a core and mantle, as a result of previous impact
heating and radioactive heating from short-lived isotopes such as
$^{26}$Al and $^{60}$Fe~\cite{Rubie2015a}.

Importantly, the metal and silicate gradually added to the Earth were
generally not in thermodynamic equilibrium with the silicate mantle of
the growing Earth. Indeed, the temperature and composition of the
bodies impacting the Earth very likely showed significant variability,
and were therefore in general different from Earth's. In addition, the
metal-silicate partition coefficients (the ratio of the concentrations
of a given chemical element in metal and silicate phases at equilibrium
with each other) depend on pressure, temperature, and composition,
which makes the conditions of thermodynamic equilibrium time-dependent
and spatially variable.  This disequilibrium must have driven fluxes of
heat and chemical elements between the metal and silicate phases, which
set the initial compositions and heat content of the core and mantle. 
Once the core has formed and the mantle has crystallised, the flux of
chemical elements between the core and the mantle becomes small (it is
limited by the slow solid state convective motions of the mantle) and
its effect on the average compositions of the core and mantle is
insignificant. This implies that the present-day composition of the
core and mantle can in principle inform us about the differentiation
process~\cite{Wood2006}.  The chemical fractionation between the core
and the mantle has been used to constrain the timing of
differentiation, with radio-chronometers like the Hafnium--Tungsten
(Hf--W) system~\cite{Yin2002,Kleine2002,Rudge2010}, and the physical
conditions under which it occurred, by comparing the estimated
core/mantle chemical partitioning with the thermodynamic metal-silicate
partition coefficients, which depend on the conditions of pressure $P$,
temperature $T$, and oxygen fugacity
$f{\mathrm{O}_2}$~\cite{Wood2006,Corgne2008,Siebert2011}. This chemical
fractionation also has major geodynamical implications. In particular,
the nature and abundance in the core of radioactive
elements~\cite{Gessmann2002,Corgne2007,Bouhifd2007},
magnesium~\cite{ORourkeStevenson2016,Badro2016}, or light elements
(e.g.\ O, S, Si) depend on the pressure $P$, temperature $T$, and
redox conditions under which metal and silicates
equilibrated~\cite{Siebert2012,Badro2015,Fischer2015}. Like the
chemical elements, the heat released during accretion was partitioned
between the mantle and the core in a way that depended on the physical
mechanisms of phase separation~\cite{king2011heat,clesi2023linking}. 
The distribution of this heat defined the value of the temperature
contrast between the mantle and the core at the end of accretion, a key
parameter for the dynamics of the planet, and in particular for the
existence of a basal magma ocean~\cite{Labrosse2007}, the existence of
an early dynamo~\cite{Williams2004,monteux2011might}, and more
generally for the thermal and magnetic evolution of the planet.

How long did it take to form the terrestrial planets? How were chemical
elements and heat partitioned between the mantle and core?  What was
the state of the Earth at the end of accretion?  Geochemistry provides
critical data to address these questions, but  their interpretation is
not unique and depends on a number of assumptions. In particular, one
key parameter is the \textit{efficiency of metal-silicate
equilibration}---a measure of how efficiently metal and silicates have
been able to exchange chemical elements during core formation---which
depends on the details of the physical mechanisms of metal-silicate
segregation. This is of particular importance for the interpretation of
core formation chronometers. There is an inverse trade-off between the
estimated differentiation timescale and the efficiency of
equilibration. Figure~\ref{fig:geochemistry}a shows the estimated core
formation time according to Hf--W systematics, as a function of the
assumed equilibration efficiency $\mathcal{E}_i$, here defined as the
total mass of element $i$ exchanged between metal and silicates
normalized by their maximum possible value, had all the metal
equilibrated with an infinitely larger silicate reservoir. As shown in
Figure~\ref{fig:geochemistry}, the accretion time estimated from Hf/W
is about $10$~My if perfect equilibration is
assumed~\cite{Yin2002,Kleine2002,Rudge2010}, but can be much larger if
equilibration is imperfect \cite {Rudge2010,Halliday2004,Kleine2004}. 
There is a limit where Hf--W systematics does not constrain the timing
of core-formation anymore, but instead give constraints on the
efficiency of metal-silicate chemical equilibration, which must be
larger than ${\simeq} 40$\%~\cite{Halliday2004,Kleine2004,Nimmo2010,Rudge2010}.  This implies
a significant flux of chemical elements between metal and silicate
during differentiation.

\begin{figure}
\includegraphics{fig01}
\caption{\label{fig:geochemistry}Geochemical constraints on core/mantle differentiation:
{(a)} Accretion time $\tau_a$ as a function of the equilibration
efficiency, assuming an exponentially decreasing accretion rate
(modified from~\cite{Rudge2010}). {(b)} Estimated core/mantle
partitioning {vs.}\ low pressure, high temperature partition
coefficient for a few low volatility elements, according to core and
mantle model composition of~\cite{HallidayWood2015} and partitioning
parameterization of~\cite{Siebert2011} for a range of oxygen fugacity
(IW$-1.5$ to IW$-3.5$).}
\end{figure}

The abundance of siderophile elements (e.g.\ Ni, Co) in Earth's mantle
provides further evidence of metal-silicate chemical equilibration, at
rather high temperature and pressure~\cite{Murthy1991}.  The
concentration of most siderophile elements is much greater in the
silicate Earth than expected from their low-temperature and
low-pressure partitioning behaviour,  a situation known as the
\textit{excess siderophile problem}~\cite{Ringwood1966}.  This can be
seen in Figure~\ref{fig:geochemistry}b which compares the apparent
core/mantle concentration ratio $D_{\mathrm{c/m}}$ of several chemical elements
(as estimated from geochemical models of mantle and core compositions)
to the corresponding equilibrium partition coefficient $D$ at low
pressure (as measured in the laboratory).  The estimated core/mantle
concentration ratio is lower than expected from low-pressure data for
siderophile elements ($D&gt;1$, e.g.\ Co and Ni), which means that these
elements are more abundant than expected in the mantle. The opposite
observation holds for lithophile elements ($D&lt;1$, e.g.\ Mn and V). The
most plausible explanation for this discrepancy is chemical
equilibration of the metal and silicate at \mbox{conditions} of high
temperature and high pressure (${&gt;}2000$ K and 30--60 GPa), under which
the partitioning of these elements is more consistent with the silicate
Earth composition~\cite{Li1996,Chabot2005,Siebert2013}. Though the
basic implications are clear---significant metal-silicate
equilibration at high temperature and pressure---, a more precise
interpretation of these data in terms of physical conditions during
accretion is difficult. It must often rely on many assumptions,
including on the degree of equilibration, starting compositions of
primitive bodies, and evolution of the oxidation conditions in Earth's
mantle~\cite{Wade2005,Wood2006,wood2008accretion,Rudge2010,Rubie2011,Rubie2015}.

\begin{figure}
\includegraphics{fig02}
\caption{\label{fig:intro}Sketch illustrating the various
metal-silicate segregation processes that can occur during the
differentiation of terrestrial planets. From~\cite{Lherm2021}, adapted
from~\cite{Rubie2015a}.}
\end{figure}

Geochemical (as summarised above) and physical arguments have led to
the scenario of planetary cores formation sketched in
Figure~\ref{fig:intro}.  Radioactive and impact heating led to
widespread melting of the growing terrestrial planets, providing
favourable conditions for metal-silicate
separation~\cite{stevenson1981models,tonks1992core,Tonks1993}. The
extent of mantle melting is likely to have varied during accretion,
possibly in a non-monotonic way, but the tendency was probably towards
increasingly deep melting (due to increasingly energetic impacts),
possibly down to the core-mantle boundary~\cite{gu2023comparisons}. 
The metal delivered by each impact is  thought to have fragmented into
droplets that eventually settled at the base of the magma
ocean~\cite{stevenson1990fluid,karato1997core,Rubie2003}. The estimated
size of these droplets (mm to cm) is small enough to have allowed them
to equilibrate efficiently with the surrounding magma
ocean~\cite{stevenson1990fluid,karato1997core,Rubie2003}. These
droplets either coalesced with the core if the molten mantle region
extended down to the core-mantle boundary, or otherwise accumulated at
the bottom of the magma ocean to form an iron
pool~\cite{stevenson1981models}. In the latter case, a Rayleigh--Taylor
instability at the boundary between the iron pool and the underlying
mantle would have produced metal diapirs that could have migrated down
to the core~\cite{tonks1992core,karato1997core,samuel2008dynamics,monteux2009model}. 
Diking (fracturation) or percolation are also possible
mechanisms~\cite{stevenson1981models}. In this scenario, metal-silicate
chemical exchange was likely negligible during the migration through
the solid portions of the mantle.  Geochemical models of core-formation
therefore usually assume that the metal last equilibrated with the
silicate at the base of the magma ocean, with its composition set by
the equilibrium partitioning coefficient at these
conditions~\cite{Li1996,thibault1995influence,Chabot2005,Siebert2013}. 
Each metal addition to the core would therefore have a composition
reflecting the conditions of pressure, temperature, and oxygen fugacity
at the base of the magma ocean.

Although this scenario is in qualitative agreement with the geochemical
data, some of the assumptions are questionable.  For example, a
significant proportion of the impactors were probably
differentiated~\cite{taylor1990accretion,greenwood2005widespread}, and
it is far from obvious that volumes of molten iron of size ${\sim}100$~km 
or more delivered by differentiated impactors could indeed fully
fragment into cm-sized droplets. Furthermore, and as discussed above,
the interpretation of geochemical data in terms of chronological
constraints and physical conditions during Earth formation depends
critically on assumptions on the efficiency of chemical transfer
between metal and
silicate~\cite{Halliday2004,Kleine2004,Nimmo2010,Rudge2010}. The
equilibration efficiency is constrained by Hf--W systematics
(Figure~\ref{fig:geochemistry}a), but the range of permissible values
allows for a number of different accretion and differentiation
scenarios~\cite{Rudge2010}.  It is therefore desirable to develop
physical models of metal-silicate chemical exchange that could
supplement geochemical data and possibly improve the geochemically
based constraints, or at least help interpret the geochemical
constraints in terms of physical processes. Heat transfer between metal
and silicates has received comparatively less attention, but, as
explained above, is important for the thermal state and evolution of
the planets.


\subsection{A fluid dynamicist's point of view}

From a fluid dynamics point of view, studying the differentiation of
terrestrial planets comes down to investigating gravity-induced phase
separation, and heat and mass transfer between the two immiscible
phases. In practice, this also requires considering the question of the
possible \textit{fragmentation} of the metal phase, as well as
\textit{stirring} and \textit{mixing} processes.  Phase separation,
fragmentation, and mixing are of importance in many industrial and
natural contexts, and have received significant interest from various
communities~\cite{villermaux2019mixing,Villermaux2020}.  However, phase
separation in the context of planetary differentiation has some
important specificities, as will be shown below, that warrant dedicated
studies.


Consider a two-phase flow with characteristic length and velocity
scales $L$ and $U$.  The separation of the two phases is induced by the
gravity $g$.  The two fluids have different densities and kinematic
viscosities, which we denote by $\rho_{m},\ \rho_{s}$ and $\nu_{m},\
\nu_{s}$, respectively (where the subscripts $m$ and $s$ denote metal
and silicates). We denote by $\Delta\rho=\rho_{m}-\rho_{s}$ the density
contrast.  The fluids being immiscible, their boundary has a surface
tension $\sigma$. The time and length scales at stake can be such that
the rotation of the planet may affect the flow; the rotation rate
$\Omega$ of the planet is therefore a possibly relevant parameter. In
addition, most of the planetary collisions that deliver the two phases
are supersonic, and the sound velocity $c$ is yet another parameter to
consider. Finally, we will be interested in heat or mass exchange
between the two phases, and therefore at least one diffusivity $\kappa$
needs to be included in the set of relevant parameters ($\kappa$ is
either a thermal diffusivity or compositional diffusivity).  Overall,
this forms a set of 11 parameters with 3 independent units (kg, m, s). 
Applying the Vaschy--Buckingham theorem, one finds that 8 independent
dimensionless numbers can be formed from these parameters, one possible
set being
{\begin{equation}
 \mathit{Re}=\frac{U L}{\nu_{s}},\quad \mathit{We}=\frac{\rho_{s} U^{2} L}{\sigma},\quad 
 \mathit{Fr}=\frac{U^{2}}{g L},\quad \mathit{Ro}=\frac{U}{\Omega L},\quad 
 \mathit{Pe}=\frac{U L}{\kappa},\quad M=\frac{U}{c},\quad 
 \frac{\Delta \rho}{\rho_{s}}, \quad \frac{ \nu_{m}}{\nu_{s}}.
 \label{Eq:NonDimensionalNumbers_1}
\end{equation}}\unskip
The Reynolds number \textit{Re}, Weber number \textit{We}, Froude
number \textit{Fr}, and Rossby number \textit{Ro} compare inertia to
viscous forces, surface tension, gravity forces, and Coriolis force,
respectively. The P\'eclet number \textit{Pe} compares the advection of
a scalar field (here either heat or concentration) to diffusion.  The
Mach number $M$ is the ratio of the velocity scale to the sound
velocity, and informs us about the possible importance of compression
effects.  The two last parameters are density and viscosity ratios. 

We end up with a eight-dimension parameter space to explore. Only parts
of this parameter space are relevant to planetary differentiation,
depending on the state (solid or molten) and degree of dispersion of
the phases (representative ranges of values for relevant physical
properties and parameters are given in  Table~\ref{Table:Properties}):

\begin{table}
\caption{Representative values of relevant physical parameters\label{Table:Properties}}
\begin{tabular}{lllll}
\tbody
Impactor radius  &amp; &amp; $R_{i}$ &amp; 50--3000~km &amp; \cite{morbidelli2009asteroids,Lunine2011} \\
Impact velocity &amp; &amp; $v_{i}$ &amp; 5--40~km{${\cdot}$}s$^{-1}$ &amp;~\cite{raymond2009building}\vspace*{8pt} \\
Earth rotation rate &amp; &amp; $\Omega$ &amp; 2--5 present-day rate &amp; \\
Earth's gravity &amp; &amp; $g$ &amp; $\leq 10$ m${\cdot}$s$^{-2}$ &amp;\vspace*{8pt} \\
\morerows{1}{Density} &amp; Metal  &amp; $\rho_{m}$ &amp; 5000--10\,000~kg${\cdot}$m$^{-3}$ &amp;~\cite{morard2013earth} \\
&amp; Silicate &amp; $\rho_{s}$ &amp; 3000--5000~kg${\cdot}$m$^{-3}$ &amp; ~\cite{Spera2010} \\
\morerows{2}{Dynamic viscosity} &amp; Molten metal  &amp; $\eta_{m}$  &amp; $10^{-3}$--$10^{-2}$  Pa${\cdot}$s &amp;\cite{poirier1988transport,de1998viscosity} \\
&amp; Molten silicate  &amp; $\eta_{s}$  &amp; $3\times 10^{-3}$--$10^{-1}$  Pa${\cdot}$s &amp;\cite{karki2010viscosity,xie2021direct}\\
&amp; Solid silicate  &amp; $\eta_{s}$   &amp; $10^{12}$--$10^{21}$ Pa${\cdot}$s &amp;\cite{scott2006effect,hansen20152} \\
\morerows{1}{Thermal conductivity} &amp; Metal &amp; $k_{m}$ &amp; 30--100 W${\cdot}$m$^{-1}\cdot$K$^{-1}$ &amp;\cite{williams2018thermal} \\
&amp; Silicate &amp; $k_{s}$ &amp; 1--2 W${\cdot}$m$^{-1}\cdot$K$^{-1}$ &amp;\cite{ni2015transport} \\
\morerows{1}{Compositional diffusivity} &amp; Metal &amp; $\kappa_{c,m}$ &amp; $10^{-12}$--$10^{-8}$ m$^{2}{\cdot}$s$^{-1}$ &amp;\cite{Rubie2003} \\
&amp;Silicate &amp; $\kappa_{c,s}$ &amp; $10^{-12}$--$10^{-8}$ m$^{2}{\cdot}$s$^{-1}$ &amp;\cite{Rubie2003,ni2015transport} \\
Interfacial tension &amp; Metal-silicate &amp; $\sigma$ &amp; $0.1$--$1$~J${\cdot}$m$^{-2}$ &amp;\cite{terasaki2012interfacial} \\
\morerows{1}{Heat capacity} &amp; Metal &amp; $c_{p,m}$ &amp; 500--1000~J${\cdot}$kg$^{-1}{\cdot}$K$^{-1}$ &amp;  \\
&amp; Silicate &amp; $c_{p,s}$ &amp; 500--1000~J${\cdot}$kg$^{-1}{\cdot}$K$^{-1}$ &amp;\cite{solomatov2015magma} \\
Sound velocity  &amp; Silicates &amp; $c$ &amp; 3--10~km{${\cdot}$}s$^{-1}$ &amp;\cite{Spera2010}
\botline
\end{tabular}
\tabnote{The ranges given here take into account the pressure and
temperature ranges that are representative of the Earth's mantle and
magma oceans, and also take into account the effect of the composition
of the metal and/or silicate phases.  The range given for the viscosity
of solid silicates accounts for the effect of temperature, melt fraction,
and grain size.}
\end{table}


\paragraph*{(i) Metal-silicate separation in a magma ocean.}
During the
migration of the metal phase in a magma ocean ({i.e.}\ when the metal
and silicates are both molten), the possibly very large length and
velocity scales involved ensure that \textit{Re}, \textit{We}, and
\textit{Pe} are all very large. For example, an iron volume 100~km in
diameter falling in a magma ocean with an initial velocity of
1~km${\cdot}$s$^{-1}$ has $\mathit{Re}\sim \mathit{We} \sim 10^{14}$
and $\mathit{Pe} \sim 10^{14}$ for heat exchange or ${\gtrsim} 10^{16}$
for a chemical species.  In the impact stage, the Froude number based
on impact velocity is ${\sim} 1$ for large impacts and $\gg 1$ for 
small-sized impactors (see discussion in
Section~\ref{subsec:impacts_rationale}). Later, when the metal falls
under its own weight, the typical velocity is set by a balance between
inertia and buoyancy, and the Froude number is set by the dynamics. 
The Mach number is in general larger than 1 during the impact (and thus
compressibility effects can be important) but then decreases quickly
below 1 (see discussion in Section~\ref{sec:compressibilty}).  The
Rossby number is large during the impact (and thus rotational effects
can be neglected) but can later decrease below 1, depending on the
sizes of the impactor and magma ocean (Section~\ref{subsec:rotation}).


The huge values of \textit{Re} and \textit{We} imply highly turbulent conditions,
with viscosity and surface tension playing a role only at the smallest
scales~\cite{Dahl2010,Deguen2011c}. Post-impact flows must have
involved a wide range of length scales, from the size of the impactor
to less than ${\sim}$mm scale at which surface tension and viscosity
eventually come into play. Numerical simulations of planetary-scale
impacts~\cite{Benz1988,Canup2000,cuk2012making,nakajima2015melting}
can include relevant effects such as compressibility, phase change,
self-gravity, and spherical geometry, that cannot be investigated with
fluid dynamics experiments. However, they have spatial resolutions
(${\sim} 14$~km in~\cite{kegerreis2022immediate}) which are usually too
rough to model the small turbulent scales at which chemical
equilibration can occur (though useful information can be extracted,
for example by calculating the Lagrangian stretching of advected
structures~\cite{kendall2016}). Direct numerical simulations of  the
fragmentation of the metal phase in a magma
ocean~\cite{Ichikawa2010,Samuel2012} have been limited so far to
laminar conditions. In contrast, fluid dynamics laboratory experiments
can more easily reach turbulent conditions, and provide insight on the
effect of turbulence on fragmentation and mass and heat transfer
between the two phases. 

\paragraph*{(ii) Metal-silicate separation in a solid mantle}  When the
silicate phase is solid, the migration of the metal phase can occur
either through diapirism, percolation, or diking. Diking, whose
importance is not well understood in the context of core formation,
will be discussed in mostly qualitative terms in
Section~\ref{subsec:diapirs_complications}.  We will focus our
discussion on the diapirism and percolation scenarios. In these cases
the Reynolds number \textit{Re} is small, either because the flow intensity is
controlled by the deformation of the very viscous silicate phase
(diapirism), or because of the small scale of the flow (percolation).
Since the typical velocity is set by a balance between buoyancy and
viscous forces, buoyancy forces also dominates over inertia ($\mathit{Fr}\ll
1$). Because of the negligible importance of inertia, the importance of
surface tension and Coriolis force is better measured by comparing them
to either buoyancy or viscous forces rather than inertia. Instead of
\textit{We} and \textit{Ro}, it is therefore more relevant to use the Bond number
\textit{Bo} and the Ekman number \textit{Ek} defined as 
{\begin{equation} 
\mathit{Bo} = \frac{\Delta \rho g L^{2}}{\sigma} =
\frac{\Delta\rho}{\rho_{s}}\frac{\mathit{We}}{\mathit{Fr}}, 
\quad \mathit{Ek} =
\frac{\nu_{s}}{\Omega L^{2}} = \frac{\mathit{Ro}}{\mathit{Re}}. 
\end{equation}}\unskip
The Bond
number is large in the case of metal diapirs, but possibly as small as
${\sim} 1$ in the case of percolation.  The Coriolis force is negligible
compared to the viscous forces ($\mathit{Ek}\gg 1$) in all of these cases.
Typical velocities are much smaller than the speed of sound ($M\ll 1$)
and thus dynamic compressibility effects (shock waves) are unimportant.
However compression heating can be of importance for the thermal
evolution of the metal phase, due to the large pressure variations
experienced by the metal as it travels through the mantle
(Section~\ref{subsec:diapir_ref_model}). The P\'eclet number \textit{Pe} is
large in the case of diapirism because of the large length scale
involved (Section~\ref{subsec:diapir_ref_model}), indicating that
transport of heat toward the core is faster than diffusion toward the
surrounding mantle. In contrast diffusion is more important in
percolation flow, due to the smallest scale of the flow ({i.e.}\ grain
size of the silicate matrix).

\subsection{A note on metal-silicate chemical equilibration efficiency}

A particular focus of  the experiments that will be discussed is the
question of chemical equilibration between the metal and silicate
phases. As discussed in Section~\ref{sec:geological_context}, this is
of primary importance for the interpretation of geochemical data in
terms of core formation timing and conditions. The equilibration
efficiency $\mathscr{E}_{i}$ can be defined as the total mass of a
given element $i$ exchanged between metal and silicates, normalized by
its maximum possible value, had all the metal re-equilibrated with an
infinitely larger silicate reservoir~\cite{Deguen2014}. Assuming that
the mass of metal fully equilibrates with a mass of silicate $\Delta$
times larger, one finds that 
{\begin{equation}
\mathscr{E}_{i} = \frac{1}{1+D_{i}/\Delta},
\label{eq:efficiency}
\end{equation}}\unskip
where $D_{i}$ is the metal/silicate partitioning coefficient of element
$i$.  Equation \eqref{eq:efficiency} highlights the importance of the
parameter $\Delta$, which has been termed \textit{metal dilution}. 
This is a convenient parameter for geochemical models because it
appears quite naturally in the conservation equations on which
geochemical models of core formation are based, but it is not
straightforward to estimate in an experiment.  Rather, what has
typically been measured is the volume of ``silicate'' liquid  in which
the ``metal'' is
dispersed~\cite{Deguen2014,landeau2021,lherm2022rayleigh,Kriaa2024}. In
the following sections we will refer to this type of measurement as
\textit{metal dispersion}. Equating the volume of equilibrated
silicates to the volume in which the metal is dispersed gives an upper
bound on the mass transfer between metal and silicates, and therefore
an upper bound on the equilibration efficiency $\mathscr{E}_{i}$.  A
better estimate of the mass transfer would require quantifying the
variations of composition in both phases, possibly through their
Probability Density Functions (Sections~\ref{sec:stat_c} 
and~\ref{subsection:ChemicalEquilibration}).

\subsection{Content}

The structure of the article follows the sequence of processes that led
to the formation of planetary cores.  Section~\ref{sec:impacts} is
devoted to the dynamics of planetary impacts. We present the governing
dimensionless numbers and discuss the relevance and limitations of
fluid dynamics laboratory experiments to the understanding of planetary
impacts. Laboratory experiments are used to illustrate and quantify the
impact and cratering processes and the resulting metal phase
dispersion. Section~\ref{sec:magmaocean} focuses on the two-phase flow
that follows an impact, when a molten impactor core falls by buoyancy
in a magma ocean.  We start by presenting the model of miscible
turbulent thermal, which we argue to be a good reference model for the
post-impact flow. We then discuss how additional
factors---immiscibility and fragmentation, inertia inherited from the
impact, Coriolis force, sedimentation---affect the predictions of this
model, and discuss the extent of chemical equilibration. Finally,
Section~\ref{sec:solidmantle} is devoted to the migration of the metal
phase through a solid part of the mantle, focusing on diapirism and
percolation.  Given the topic of this special issue on
\textit{Geophysical and Astrophysical Fluid Dynamics in the lab}, the
paper is biased toward the experimental fluid dynamics contributions to
the understanding of planetary differentiation.  Results from numerical
simulations and experimental petrology will be mentioned, but the
review is by no means intended to be comprehensive on these approaches.

\section{Planetary impacts}\label{sec:impacts}
\subsection{Rationale}\label{subsec:impacts_rationale}

Planetary impacts have been classified as being either in the so-called
\textit{strength regime} or \textit{gravity regime}, depending on the
importance of the resistance to deformation of the target
material~\cite{Melosh1989}. This can be evaluated by comparing the
level of impact-induced stresses to either the yield stress for a solid
target, or viscous stresses in the case of an impact onto a magma
ocean. On Earth, impact-induced stresses exceed the yield stress for
impactors that are larger than a few tens of
meters~\cite{holsapple1993scaling,Melosh1989}. Given the expected range
of sizes of planetesimals (${\gtrsim} 100$ km in diameter) and embryos
(${\gtrsim} 1000$ km in diameter), most of the impacts that formed the
Earth were thus likely in the gravity regime in which the level of
induced stresses is significantly larger than the yield
stress~\cite{holsapple1993scaling,Melosh1989}. In the case of an
impactor striking a pre-existing magma ocean, the impact-induced stress
is much larger than viscous stress {i.e.}\ $\mathit{Re}\gg 1$. The resistance to
deformation is thus usually of negligible importance in this context. 
During an impact, \textit{We}, \textit{Pe}, and \textit{Ro} are all much larger than
unity.  The Mach number $M$ based on the impact velocity is usually
larger than 1, and possibly up to ${\sim} 10$.

In this limit, the dominant forces are inertia, pressure forces, and
gravity, and the non-dimensional numbers of most relevance are \textit{Fr},
$M$, and the density ratio, with \textit{Fr} and $M$ defined in Equation
\eqref{Eq:NonDimensionalNumbers_1} using the impactor radius and impact
velocity.  An additional parameter of importance is the impact angle,
{i.e.}\ the angle between the impact velocity vector and the horizontal.

It will be useful to also introduce a modified Froude number defined as 
{\begin{equation}
\mathit{Fr}^{*} = \frac{\rho_i}{\rho_t} \frac{U_i^2}{g\, R_i},
\end{equation}}\unskip
where $R_i$ is the diameter of the impactor, $U_i$ its velocity at
impact, and $\rho_{i}$ and $\rho_{t}$ are the densities of the impactor
and target material,  which may be mixtures of metal and silicate
(which is why we do not use the notations $\rho_{m}$ and $\rho_{s}$ in
this section). With this definition, $\mathit{Fr}^{*}$ is the ratio of the
kinetic energy of the impacting body (${\sim} \rho_{i} U_{i}^{2}
R_i^{3}$) to the gravitational potential energy (${\sim} \rho_{t} g
R_i^{4}$) released by displacing a volume ${\sim} R_i^3$ of target
material over a height ${\sim} R_i$.  Assuming that the impact velocity
scales as the escape velocity, $U_{i}\simeq \sqrt{g D_p}$, where $D_p$
is the diameter of the growing planet and $g$ the acceleration of
gravity at its surface, we find that 
{\begin{equation}
\mathit{Fr}^{*} \simeq \frac{\rho_{i}}{\rho_{t}} \frac{D_p}{R_{i}}.
\label{eq:Fr_escape_v}
\end{equation}}\unskip
Since $\rho_{i}/\rho_{t}$ varies within a limited range, the Froude
number $\mathit{Fr}^{*}$ depends primarily on the impactor-to-target size ratio.
$\mathit{Fr}^{*}$ ranges from ${\sim} 1$ in the case of a giant impact to  several
hundreds in the case of ${\sim} 100$~km-diameter planetesimals impacting
the Earth. Recent meteorite impacts have much larger Froude numbers. As
an example, Meteor Crater (Arizona) is thought to have been formed by
the impact of an iron mass about 50 m in diameter with a velocity of
more than 10~km${\cdot}$s$^{-1}$, which gives a Froude number ${\sim} 2\times
10^{5}$.

The full range of $\mathit{Fr}^{*}$ values relevant to the formation of
terrestrial planets can be explored with impact experiments involving
millimetre- to centimetre-sized drops released from different heights
above a bath of a less dense liquid.  This approach differs from
experiments involving supersonic impacts of solid projectiles onto a
solid target, more typically carried out in a planetological
context~\cite{pierazzo2008validation,nagaoka2014degree}.  Both
approaches have their advantages and limitations: the possibility of
reaching supersonic impact speeds comes at the expense of the value of
$\mathit{Fr}^{*}$, which  is typically on the order of $10^9$ in this type of
experiments, {i.e.}\ seven orders of magnitude larger than during the
impact of a 100~km-diameter body; in contrast, fluid dynamics
laboratory experiments as discussed below allow to explore the relevant
range of Froude number, but ignore any possible effects of
compressibility.

To be relevant to the context of planetary impacts, experiments need to
be done at high \textit{Re} (to be in the ``gravity'' regime) and high \textit{We}
(if involving two immiscible liquids). In an experiment in which the
impactor is produced by letting a volume of liquid fall under its own
weight toward the target, the impact velocity $U_i$ is controlled by
the drop height. It can be varied between near zero and the terminal
velocity $U_{i}^t\sim
\sqrt{({\rho_i-\rho_{\mathrm{air}}})/{\rho_{\mathrm{air}}}g R_i}\sim
\mathcal{B}^{1/2}/R_i$ of a drop falling into air, where
$\rho_{\mathrm{air}}$ is the density of air and $\mathcal{B} =
({\rho_i-\rho_{\mathrm{air}}})/{\rho_{\mathrm{air}}}g
\sofrac{4\uppi}{3}R_i^3$ is the buoyancy of the drop. The Reynolds, Weber,
and Froude numbers based on the impactor radius $R_i$ and velocity
$U_i$ can then be written and bounded as follows
{\begin{eqnarray}
\mathit{Re} &amp;=&amp; \frac{\mathcal{B}^{1/2}}{\nu} \frac{U_i}{U_i^t}\lesssim
\frac{\mathcal{B}^{1/2}}{\nu} , \Seqnsplit
\mathit{We} &amp;=&amp; \frac{\rho_i-\rho_{\mathrm{air}}}{\rho_{\mathrm{air}}}
\mathit{Bo} \left( \frac{U_i}{U_i^t} \right)^2 \lesssim
\frac{\rho_i-\rho_{\mathrm{air}}}{\rho_{\mathrm{air}}} \mathit{Bo},
\Seqnsplit
\mathit{Fr} &amp;=&amp; \frac{\rho_i-\rho_{\mathrm{air}}}{\rho_{\mathrm{air}}}
\left( \frac{U_i}{U_i^t} \right)^2  \lesssim
\frac{\rho_i-\rho_{\mathrm{air}}}{\rho_{\mathrm{air}}}, 
\end{eqnarray}}\unskip
where the Bond number is defined as
{\begin{equation}
 \mathit{Bo} = \frac{\rho_i g R_i^2}{\sigma}
\end{equation}}\unskip
(where here $\sigma$ is the interfacial tension between the two
immiscible liquids, not between the drop and air).  From this, it is
found that Reynolds number values ${\gtrsim} 10^4$ and Weber number
values ${\gtrsim} 10^3$ can easily be reached with drops of a few mm in
diameters, at relatively large values of
\textit{Fr}~\cite{engel1966crater,lherm2022rayleigh,lherm2023velocity}.
However, reaching high values of \textit{Re} and \textit{We} at
moderate values of \textit{Fr} requires larger impactors. We can indeed
notice that $\mathit{We} = \mathit{Bo} \mathit{Fr}$, which implies that
reaching a high value of \textit{We} while maintaining a moderate value
of \textit{Fr} requires a high value of \textit{Bo}, and therefore a
larger impactor, typically of a few cm in diameter at least. Drops with
a diameter up to ${\sim} 5$ mm can be easily formed by dripping liquid
from the tip of a syringe or
capillar~\cite{Lhuissier2013,lherm2022rayleigh,lherm2023velocity}.  At
this size, the drop is kept coherent by surface tension as it falls
(though its shape starts to deviate from spherical, and can oscillate).
Laboratory experiments with drops produced this way reach Froude
numbers in the range $[10^2\ndash
10^3]$~\cite{lherm2022rayleigh,lherm2023velocity}. Producing liquid
impactors of a few centimeters with moderate deformation is more
complicated, but necessary to explore smaller values of \textit{Fr}
while still keeping $\mathit{Re},\mathit{We}\gg
1$~\cite{landeau2021,MallerInRev}. One method consists in enclosing the
impactor liquid in a balloon which is pierced just before
impact~\cite{Wacheul2014,Wacheul2018,landeau2021,MallerInRev}.


\subsection{Phenomenology and scaling laws}
\label{sec:ImpactPhenomenology}

The impact of a drop on a deep pool has been studied
extensively~\cite{worthington1883impact,worthington1897v,worthington1900iv,engel1966crater,
engel1967initial,prosperetti1993impact,liow2001splash,fedorchenko2004some,bisighini2010crater,ray2015regimes,lherm2023velocity} 
since the pioneering work of Worthington at the turn of the 19th and
20th
centuries~\cite{worthington1883impact,worthington1897v,worthington1900iv}.
Although the analogy with planetary impacts and its limitations were
not fully formalised, the idea of using drop impacts as models for
meteorite impacts was already explicit in Engel's 1966
paper~\cite{engel1966crater}.  Most studies focused on the impact of a
drop on a pool of the same liquid. The effect of a difference of
density between the drop and
pool~\cite{lherm2021rayleigh,lherm2022rayleigh,landeau2021} and of
immiscibility~\cite{Lhuissier2013,Jain2019,Hadj2021,MallerInRev} have
been the subject of fewer studies, and not always in regimes relevant
to planetary differentiation. Experiments on oblique impacts are scarce
as well~\cite{gielen2017oblique}.

The main results from these studies will be discussed in the following
subsections, and illustrated with the experiments shown in
Figure~\ref{fig:impact_images}. Figure~\ref{fig:impact_images}a shows
successive snapshots from an experiment in which a drop of Caesium
Chloride ($\rho_i=1789$~kg${\cdot}$m$^{-3}$) strikes a pool of
water~\cite{lherm2022rayleigh}. The drop has a radius $R_i=2.36$ mm and
its velocity at impact is $U_i=3.54$ m${\cdot}$s$^{-1}$, which results in
$\mathit{Fr}=540$ and $\rho_i/\rho_t=1.789$. Figure~\ref{fig:impact_images}b
shows snapshots from an experiment in which a volume of
Perfluorodecalin ($\rho_i=1940$~kg${\cdot}$m$^{-3}$) strikes a pool of
water~\cite{MallerInRev}, with $R_i=2$ cm and $U_i=4.3$ m${\cdot}$s$^{-1}$,
which results in $\mathit{Fr}=94$, $\rho_i/\rho_t=1.94$.  Perfluorodecalin and
water are immiscible and have an interfacial tension $\sigma =50$
mJ${\cdot}$m$^{-2}$, giving $\mathit{We}=1.4 \times 10^4$.  
Figures~\ref{fig:impact_images}c and~\ref{fig:impact_images}d are composite
images from the experiment shown in Figure~\ref{fig:impact_images}a
(left part of the images) and similar experiments in which velocity
field (PIV) and concentration (PLIF) measurements have been
done~\cite{lherm2023velocity}. Figure~\ref{fig:impact_images}e shows a
close-up of image {iv} of Figure~\ref{fig:impact_images}b.

We will discuss the three main stages of the impact and cratering
process that can be identified from these experiments at large \textit{Re} and
\textit{We}:  (i) the opening and growth of the crater
(Section~\ref{sec:CraterOpening}), with a particular focus on
instabilities and fragmentation that are observed during this stage;
(ii) the collapse of the crater and formation of a central peak
(Section~\ref{sec:CraterClosing}); and\break (iii)~the collapse of the
central peak (Section~\ref{sec:collapse}).

\begin{figure}[p!]
{\vspace*{10pt}}
\includegraphics{fig03}
{\vspace*{4pt}}
\caption{\label{fig:impact_images}{(a)} Time sequence from a liquid impact experiment at
$\mathit{Fr}=540$ and $\mathit{Re}=8.3\times 10^3$ with miscible drop and target liquids
(using a CsCl aqueous solution with $\rho_i=1789$~kg${\cdot}$m$^{-3}$
for the drop, water as the target liquid, $R_i=2.36$ mm, and $U_i=3.54$
m${\cdot}$s$^{-1}$).  The CsCl solution appears blue on the pictures. {(b)}
Time sequence from a liquid impact experiment at $\mathit{Fr}=94$,
$\mathit{We}=1.4\times 10^4$, and $\mathit{Re}= 8.6\times 10^4$ with immiscible drop
and target liquids (using Perfluorodecalin with
$\rho_i=1940$~kg${\cdot}$m$^{-3}$ for the drop, water as the target
liquid, $R_i=2$ cm, and $U_i=4.3$ m${\cdot}$s$^{-1}$)~\cite{MallerInRev}. {(c)}
and {(d)} Composite images with pictures from the experiments shown on
subfigure (a) on the left, concentration fields obtained using Planar
Laser Induced Fluorescence on the right, and the velocity field
obtained from Particle Image Velocinetry shown with black arrows.  The
concentration measurements have been performed in an experiment with
parameter values similar but slightly different from (a) ($\mathit{Fr}=572$,
$\rho_i/\rho_t=1.754$); the velocity measurements have been performed
in an experiment involving a water drop (to avoid strong variations in
optical index) and has $\rho_i/\rho_t=1$ and $\mathit{Fr}=446$. {(e)} Close-up
on image (b)iv.}
\end{figure}

\subsubsection{The opening phase: Rayleigh--Taylor instability and
fragmentation} \label{sec:CraterOpening}

When it strikes the surface of the pool, the drop is abruptly
decelerated, and flattens onto the surface of the pool. Part of the
momentum of the drop is imparted to the pool, which induces an outward
flow that opens a cavity---the crater 
(Figures~\ref{fig:impact_images}ai and~\ref{fig:impact_images}bi). The
crater then expands and quickly becomes near spherical, until reaching
a maximum in size. The liquid ejected from the crater forms a crown
(images i and ii of Figures~\ref{fig:impact_images}a and
\ref{fig:impact_images}b) which expands concomitantly with the
crater. 

\paragraph*{Evolution of crater size} 
The evolution of the crater mean radius $R_{c}(t)$ has been modelled by
considering the energy balance of the
system~\cite{engel1966crater,engel1967initial,liow2001splash,lherm2022rayleigh}.
Neglecting viscous dissipation (the flow around the crater is close to
potential~\cite{lherm2023velocity}), one can assume that the sum of
kinetic energy of the flow, potential energy, and surface tension
energy is constant. This can be written as
{\begin{equation}
 \xi \uppi \rho_t R_{c}^3 \dot R_{c}^2 + \uppi \sigma R_{c}^2 + 
 \frac{\uppi}{4} \rho_t g R_{c}^4 = \phi \frac{2 \uppi}{3} \rho_i R_i^3 U_i^2.
 \label{eq:EnergyBalance}
\end{equation}}\unskip
The first term on the LHS is the kinetic energy of the flow around the
crater, with $\xi$ a correction parameter accounting for the unknown
spatial variations of the flow. $\xi=1$ if the flow is radially outward
and incompressible~\cite{liow2001splash,lherm2022rayleigh}.  PIV
velocity measurements~\cite{lherm2023velocity} have shown that this is
actually a poor approximation (as shown in
Figure~\ref{fig:impact_images}c): the effective value is smaller than
1 and time-dependent. The second and third terms account for changes of
surface energy and potential energy associated with the opening of the
crater (here neglecting for simplicity the contributions of the crown).
The term on the RHS is the kinetic energy of the impact, corrected by a
factor $\phi$ which accounts for possible energy loss. Fitting the
solution of \eqref{eq:EnergyBalance} to experimental data requires
$\xi\simeq 0.34$ and $\phi \simeq 0.4$~\cite{lherm2022rayleigh}.
Surface tension does not affect much the shape of the crater during the
opening phase if \textit{We} is large enough ($\mathit{We}\gtrsim 54 \mathit{Fr}^{0.25}$
according to~\cite{liow2001splash}), and its only effect is to slow
down the opening.

Despite its simplicity, the energy approach gives the correct scaling
relationships for the maximum size of the crater and the time at which
this maximum is reached.  The maximum crater size can easily be
obtained by setting $\dot R_{c}=0$ in Equation
\eqref{eq:EnergyBalance}. Neglecting the surface tension term for
simplicity, this yields,
{\begin{equation}
R_{c,\mathrm{max}} = R_{i} \left( \frac{8}{3}\phi\right)^{1/4} 
{\mathit{Fr}^{*}}^{1/4},
\label{eq:scaling_Rmax}
\end{equation}}\unskip
which is in good agreement with drop impact
experiments~\cite{lherm2022rayleigh}. Importantly, this scaling law is
also in good agreement with numerical simulations of planetary impacts
until relatively large values of the Mach number (see
Section~\ref{sec:compressibilty}).

The time $\tau_{\mathrm{max}}$ at which the crater reaches its maximum
size can be predicted by integrating in time equation
\eqref{eq:EnergyBalance}.  This predicts
{\begin{equation}
\tau_{\mathrm{max}} \sim \frac{R_{i}}{U_{i}} 
\left( \frac{\rho_i}{\rho_t} \right)^{-1/2}{\mathit{Fr}^{*}}^{5/8},
\label{eq:tau_max}
\end{equation}}\unskip
which again is consistent with experimental observations
(see~\cite{lherm2022rayleigh} for an analytical expression and
extension to include the effect of surface tension). For a planetary
impact with $R_{i}=100$~km and $U_{i}=10$~km{${\cdot}$}s$^{-1}$, Equation
\eqref{eq:tau_max} predicts a crater opening time of ${\sim} 3$ min.

\paragraph*{Instabilities}

Inspection of Figures~\ref{fig:impact_images}a,b,c shows that the
liquid carried by the impactor forms a layer around the opening crater,
on which small-scale corrugations and plumes quickly develop (see
snapshots {i} and {ii} of
Figures~\ref{fig:impact_images}a and b, and~\ref{fig:impact_images}a).
In addition, measurements of the thickness of the layer show that it
increases with time, while conservation of mass predicts a decrease of
thickness in the absence of incorporation of ambient liquid, as a
result of the spreading of the liquid around the growing
crater~\cite{lherm2022rayleigh}.

This is interpreted as being the result of a Rayleigh--Taylor
instability triggered by the deceleration of the crater
boundary~\cite{lherm2021rayleigh,lherm2022rayleigh,landeau2021}. In the
opening phase, the crater radius $R_{c}$ is an increasing function of
time and $\dot R_{c} = \mathrm{d}R_{c}/\mathrm{d}t$ 
is a positive but decreasing function
of time (it is equal to 0 when the crater reaches its maximum size).
The acceleration of the crater boundary $\ddot R_{c} =
\mathrm{d}^{2}R_{c}/\mathrm{d}t^{2}$ is therefore negative: the boundary decelerates. In
a reference frame moving with the boundary,  this is equivalent to an
outward apparent gravity of magnitude $ {-} \ddot R_{c}$. If the impactor
is denser than the ambient, the situation is analogous to having a
dense layer above a less dense liquid, a situation which is known to be
unstable against the Rayleigh--Taylor instability.  In practice, $\ddot
R_{c}$ can be significantly larger than the acceleration of gravity
$g$. For a significant fraction of the crater opening phase, the
kinetic energy of the flow surrounding the crater dominates over the
contributions of the potential and surface energies.  Neglecting these
two terms in Equation \eqref{eq:EnergyBalance}, taking the derivative
with respect to time and re-arranging gives $\ddot R_{c}/g \sim \mathit{Fr}^*
(R_i/R_{c})^4$. When $\mathit{Fr}^* \gg 1$, the magnitude of the deceleration is
therefore initially much larger than $g$. $\ddot R_{c}$ becomes ${\sim}
g$ when the crater approaches its maximum size.


The total increase of volume of the mixing layer depends on the
magnitude of $\sfrac{\Delta \rho}{\rho_t} |\ddot R|$, but depends also
on the time interval in which the mixing layer can grow, which is close
to the opening time $\tau_{\mathrm{max}}$. Lherm 
et al.~\cite{lherm2022rayleigh} found that the volume of the mixing
layer, normalized by the impactor's volume, is ${\simeq} 0.07
(\rho_i/\rho_t)^{2.8}\mathit{Fr}^{0.67}$. 


\paragraph*{Fragmentation}
Maller et~al.~\cite{MallerInRev} recently investigated the impact of a
large, centimetre-sized, spherical liquid volume onto a pool of a
lighter immiscible liquid.  When the impact speed is large enough, they
observed that the impacting liquid fragments into millimetric drops
(last two images in Figure~\ref{fig:impact_images}b, and
Figure~\ref{fig:impact_images}e).  Varying systematically the impactor
size and speed, the authors quantified the degree of fragmentation of
the impacting liquid by measuring the number of fragments prior to the
collapse of the central peak. They found that the entire impacting
volume fragments into drops when the Froude number \textit{Fr} exceeds 40,
regardless of the value of the Weber number \textit{We}, and hence the
strength of surface tension.  This result is unexpected as surface
tension, through the Weber number, often controls the fragmentation of
a liquid. 

To explain this observation, the authors assumed that, during the
opening of the crater, the Rayleigh--Taylor instabilities described
above stretch the impacting liquid into ligaments, which then break
into drops by capillary instabilities.  The timescale of capillary
instabilities is much smaller than the cratering time in these
experiments. This means that, as soon as a ligament forms, surface
tension breaks it into drops almost instantaneously. The breakup of
ligaments is therefore not the limiting process here. 


In contrast, the timescale of nonlinear Rayleigh--Taylor instabilities
is on the same order of magnitude as the cratering time. This suggests
that fragmentation is possible only when Rayleigh--Taylor instabilities
have time to grow and stretch the impacting liquid into ligaments. 
Because $\mathit{We}&gt;300$ in these experiments, Rayleigh--Taylor instabilities
and impact cratering do not depend on surface tension, which explains
why the degree of fragmentation is independent of \textit{We}. 

During the formation of planets, most impactors are differentiated into
a metallic core and a silicate mantle.  Maller et al.~\cite{MallerInRev}
also conducted experiments with two-phase impactors using a balloon
enclosed into another. They found that the impact Froude number still
controls the fragmentation, which again occurs when $\mathit{Fr}&gt;40$. 

Applied to planetary impacts, the above results predict that the larger
the impact speed, the larger the Froude number \textit{Fr}, and hence the more
fragmented the impactor core.  As discussed in
Section~\ref{subsec:impacts_rationale}, the Froude number of a
planetary impact scales as the ratio of the planet to impactor
diameters (Equation \eqref{eq:Fr_escape_v}). This leads to the
following result: the larger the impactor size, the smaller the Froude
number, and hence the degree of fragmentation. For impacts onto Earth
at the escape speed, Maller et al.~\cite{MallerInRev} predict that the
cores of impactors with a radius $R_i&lt;300$~km fully fragment into drops
during the impact process, while the cores of giant impactors with
$R_{i}&gt;1600$~km remain unbroken. For impact velocities higher than 1.5
times the escape speed, a significant fraction of the impactor core
fragments for any impactor size.


\subsubsection{Crater closing and formation of a central peak}
\label{sec:CraterClosing}

The crater reaches its maximum radius when all the kinetic energy of
the impact has been consumed by the change in gravitational potential
energy associated with the opening of the crater.  This potential
energy is then converted back into kinetic energy when the crater
closes.  Qualitative~\cite{engel1966crater} and
quantitative~\cite{lherm2023velocity} velocity measurements show that
the flow is convergent toward the axis of symmetry of the crater
(Figure~\ref{fig:impact_images}d). This convergence leads to the
formation of an upward-growing central peak, or jet
(Figures~\ref{fig:impact_images}a{iv}
and~\ref{fig:impact_images}b{iv} and {v}), which is made
of a mixture of drop and pool liquids~\cite{lherm2022rayleigh}.  In
millimetre drop impact experiments
(Figures~\ref{fig:impact_images}a{iv} and {v}), the
central peak is often observed to fragment due to the development of a
capillary instability~\cite{worthington1883impact}, a mechanism that is
irrelevant in the case of a planetary impact. The central peak always
remains coherent (no fragmentation) in experiments with
multi-centimetre drops~\cite{landeau2021,MallerInRev}. The flow
associated with the closure of the crater and the formation of the
central peak strongly stretches the mixing zone formed by the
Rayleigh--Taylor instability~\cite{landeau2021}.

\subsubsection{Collapse of the central peak and further mixing}
\label{sec:collapse}

The central peak reaches a maximum height before collapsing, thus
re-injecting into the pool the drop/pool liquid mixture it carries. At
moderate values of $\mathit{Fr}^{*}$, this mixture forms a well-defined
compact cloud with little residual velocity~\cite{landeau2021}, as can
be seen for example on image {v} of Figure~\ref{fig:impact_images}a. 
The dispersion of the drop liquid at this stage has been quantified
experimentally by estimating the volume over which the drop liquid is
spread~\cite{landeau2021}. This volume has been found to scale as
${\sim} R_i^{3} (({\rho_i-\rho_t})/{\rho_t})^{-0.57}
\mathit{Fr}^{1.2}$~\cite{landeau2021}: it is an increasing function of
\textit{Fr} but, perhaps surprisingly, a decreasing function of the
buoyancy ratio $(\rho_i-\rho_t)/\rho_t$.  These experiments have been
performed with density ratios $\rho_i/\rho_t$ larger than but close to
1;  it remains to be checked that the above scaling law still applies
at density ratios approaching the metal-silicate value of ${\simeq} 2$.

\subsection{Limitations and future work}

\subsubsection{On the effect of compressibility}
\label{sec:compressibilty}

An obvious limitation of fluid mechanics experiments applied to
planetary differentiation is their low Mach numbers $M$, and therefore
absence of compressibility effects. This seems unavoidable; the effect
of compressibility needs to be investigated with numerical models.

One well-identified effect of compressibility is the formation, when
$M&gt;1$, of a shock wave which propagates away from the impact, followed
and eventually caught up by a rarefaction wave (or pressure relief
wave).  The compression/decompression cycle associated with the shock
and rarefaction waves is an irreversible process that results in
melting, and possibly vaporisation, of the shocked
material~\cite{Melosh1989}. The impact crater opens in the
decompressed, molten region behind the rarefaction wave (see for
example Figure~4 from~\cite{allibert2023planetary} showing snapshots
from a numerical model of an impact at $M=5$, or Figure~2
from~\cite{kendall2016}). This, in fact, justifies to some extent the
use of liquids in laboratory experiments to model planetary impacts on
solid targets~\cite{allibert2023planetary}. However, compressibility
does affect the maximum size of the crater, as shown by the systematic
investigation of the effect of \textit{Fr} and $M$ on the maximum crater
depth~\cite{allibert2023planetary}.  The effect of the Mach number on
the maximum crater depth is moderate up to $M \sim 4$ (less than ${\sim}
20$\%), but becomes more significant for $M\gtrsim 4$, a limit that
is also relevant to planetary formation~\cite{allibert2023planetary}.
At large $M$, the maximum crater depth becomes independent of the
impact velocity and depends on a Froude number defined from the speed
of sound rather than the impact velocity~\cite{allibert2023planetary},
with the maximum crater depth scaling as ${\sim} R_{i} (\mathit{Fr}
M^{-2})^{1/4} = R_{i} (\sfrac{c^{2}}{g R_{i}})^{1/4}$.

Compressibility can also affect the smaller scales of the flow around
the crater, and therefore the dispersion and fragmentation of the
metal. It is indeed know that compressibility can affect the coherence
of large scale turbulent eddies in turbulent flows, which  can be
understood as follows.  In a fluid, the pressure field adjusts to time
variations of the flow through pressure waves ({i.e.}\ sound
waves). At low Mach number, pressure waves travel much faster than the
flow and the pressure field adjusts almost instantaneously, thus
allowing long-range interactions between different parts of the flow.
When the sound velocity is on the same order of magnitude or larger
than the flow velocity, the finite speed of sound introduces a time
delay in the transmission of pressure signals from one point to
another, which can disrupt the coherence of large turbulent
eddies~\cite{Breidenthal1992,Freund2000,Pantano2002}.  This has been
observed experimentally to have a significant effect on dispersion and
mixing in jets and mixing layers when $M\gtrsim
1$~\cite{Brown1974,Freund2000}.

However, there are reasons to believe that the effect of
compressibility may not be as important as the value of the Mach number
at impact would suggest, because the velocity of the flow after the
impact is smaller than the impact velocity. The maximum velocity is
abruptly decreased when the impactor comes into contact with the target
because part of the momentum is transmitted to the target material.
This results in a decrease of the maximum velocity, and therefore of
the instantaneous Mach number, by roughly a factor of
2~\cite{Melosh1989}. In addition, the velocity is further decreased
during crater opening because: (i) the flow around the crater occupies
a volume significantly larger than the volume of the impactor, which
means that the kinetic energy is distributed over an increasingly large
volume; (ii) the total kinetic energy decreases in proportion to the
increase of the potential energy associated with the excavation of the
crater, as suggested by the energy balance (Equation
\eqref{eq:EnergyBalance}). This reasoning suggests that the effect of
compressibility may be limited to a short interval of time after
contact between the impactor and the target.  This is supported by the
analysis of snapshots from numerical simulations of planetary impacts,
which shows that the flow around the crater is subsonic for a
significant part of the crater opening
phase~\cite{kendall2016,allibert2023planetary}.


\subsubsection{Oblique impacts}

\begin{figure}
\includegraphics{fig04}
\caption{\label{fig:ObliqueImpacts}Time sequences from three oblique impact experiments with
impact angles $\theta$ (with respect to the horizontal) equal to
59\textdegree~(left column), 43\textdegree~(middle column), and
25\textdegree~(right column).  In all three experiments, the impactor
is made of 290 mL of dyed water ($R_{i}=4.1$~cm) and is released from a
height 70 cm above the surface of the tank, which is filled with fresh
water.  The impacting liquid is initially enclosed in a balloon which
is allowed to slide with respect to a tilted guiding rail, giving it an
oblique trajectory.  The balloon falls under its own weight and is
ruptured by a scalpel blade just before impact (first images of the
time sequences). On these images, the two dark horizontal lines
correspond to the contact lines of the free surface with the container.
Two lines rather than one appear because the free surface is seen
slightly from above.}
\end{figure}

The experiments discussed in Section~\ref{sec:ImpactPhenomenology}
focused on vertical impacts. This is an important limitation: planetary
impacts are more likely to be oblique, with the most probable impact
angle being 45\textdegree~\cite{pierazzo2000understanding}.  This is
known from numerical simulations and hyper-velocity or drop impact
experiments to strongly affect the shape of the crater and the
dispersion of the
impactor~\cite{gault1978experimental,pierazzo2000understanding,pierazzo2000melt,leneweit2005regimes,okawa2008effect,kendall2016,gielen2017oblique}. 


Figure~\ref{fig:ObliqueImpacts} shows snapshots from yet unpublished
experiments~\cite{Wacheul_in_prep} in which a volume of dyed water with
radius $R_{i}=4.1$ cm strikes a water bath with an impact angle of
either 59\textdegree, 43\textdegree, or 25\textdegree. In contrast with
vertical impacts, the crater, crown, and central peak are asymmetric,
which is consistent with numerical
models~\cite{pierazzo2000understanding,kendall2016} and fluid dynamics
experiments using mm-sized
drops~\cite{leneweit2005regimes,okawa2008effect,gielen2017oblique}. The
snapshots shown on the second line of Figure~\ref{fig:ObliqueImpacts}
show that a large-scale instability develops at the boundary between
the impactor and pool liquid, in particular at the lowest impact angles
of 43\textdegree\ and 25\textdegree.  The shape of the observed
structures is reminiscent of Kelvin--Helmoltz billows, which suggests
that they originate from the shear associated with the anisotropic
growth of the crater.  Qualitatively, this instability appears to
increase the amount of dispersion of the impactor material, but this
remains to be quantified in detail. The effect of the impact angle on
the dispersion and mixing of the impactor material remains to be
quantified and modelled. 

\section{Metal-silicate separation in a magma ocean}\label{sec:magmaocean}

\subsection{Rationale}\label{sec:magma_oceean_rationale}

The experiments discussed in Section~\ref{sec:impacts} show that the
impact leaves the impactor material as a coherent mass without much
momentum, distributed over a volume $V_{0}$ larger than the impactor's.
This volume $V_{0}$ is a mixture of the impactor and target mantle
material, which may or may not be fragmented.  It is negatively buoyant
and therefore falls under its own weight. Given its expected size (say
a radius $R_{0}$ equal to 10 to 1000~km) and the low viscosity of the
magma ocean (${\sim} 10^{-2}$~Pa${\cdot}$s), its velocity initially scales as
$\sqrt{\sofrac{\Delta\rho_0}{\rho_s} g R_{0}}$, where $\Delta \rho_{0}$
is the effective density difference between the impactor/magma mixture
and surrounding magma. Using this scaling to estimate \textit{Re} and \textit{We}
gives
{\begin{eqnarray}
 \mathit{Re} &amp;\sim&amp; \frac{\mathcal{B}_{0}^{1/2}}{\nu_s}, \Seqnsplit
 \mathit{We} &amp;\sim&amp; \frac{\Delta \rho_{0} g R_{0}^2}{\sigma} = \mathit{Bo},
\end{eqnarray}}\unskip
where $\mathcal{B}_{0} = \sofrac{\Delta\rho_{0}}{\rho_s} g V_{0}$ is the
total buoyancy of the metal volume. For falling impactors, very large
values are expected for these two dimensionless numbers.  For example,
taking $R_{0}=10$~km (a relatively small impact) and $\Delta
\rho_{0}=1000$ kg.m$^3$ (which corresponds to a significant dispersion
upon impact) gives $\mathit{Re} \sim \mathit{We} \sim 10^{12}$ (taking $\sigma \sim 1$
J${\cdot}$m$^2$ for the metal-silicate surface tension, $g=10$ m${\cdot}$s$^{-2}$,
$\rho_s=4000$~kg${\cdot}$m$^{-3}$, $\nu_s=10^{-6}$ m$^2{\cdot}$s$^{-1}$). The
large value of \textit{We} suggests that surface tension is unimportant for
large-scale dynamics. As a first step, the effect of surface tension
can be ignored ({i.e.}\ taking the limit $\mathit{We} \rightarrow \infty$),
which can be done experimentally by using \textit{miscible} liquids
instead of immiscible liquids. At large Reynolds numbers, a buoyant
mass of fluid forms what is known as a \textit{turbulent thermal}. This
is discussed in Section~\ref{sec:miscible_thermals}, and provides a
useful reference on which additional physical effects can be added
incrementally. The effects of immiscibility, rotation, and
impact-induced dispersion, which are all relevant to planetary
formation, are discussed in the following subsections. We will also
discuss experimental measurement of the evolution of the concentration
field within a miscible thermal.

\subsection{The infinite \textit{We} limit: miscible turbulent thermals}
\label{sec:miscible_thermals}

Figure~\ref{fig:thermals}a shows snapshots from an experiment in which
a negatively buoyant volume of an aqueous solution of sodium iodide
(NaI) is released into fresh water. The volume of dyed liquid increases
as it falls, which indicates that it entrains and incorporates ambient
liquid, resulting in its gradual dilution. Measurements on many similar
experiments and numerical simulations have shown that the mean radius
of the negatively buoyant mixture increases linearly with the distance
from the point of
release~\cite{Batchelor1954,Morton1956,Scorer1957,Woodward1959,Turner1969,Lecoanet2019}. 
PIV measurements on similar experiments and numerical simulations show
that the velocity field has a vortex ring structure, with most of the
entrainment of ambient liquid occurring from the rear of the
cloud~\cite{Lecoanet2019,Lherm2021}. This is what is known as a
\textit{turbulent thermal} in the fluid mechanics and atmospheric
science
communities~\cite{Batchelor1954,Scorer1957,Woodward1959,Turner1969}, a
name inherited from the usage of glider pilots, for whom a thermal is
an isolated mass of warm air rising through the lower atmosphere. 
Though the buoyancy in atmospheric thermals is due to temperature
differences, the nature of the source of buoyancy (thermal or
compositional) happens to be of secondary importance and introduces no
qualitative difference. The term thermal has since been used to denote
an isolated buoyant mass of fluid rising or falling, regardless of the
nature of the source of buoyancy.  In the following, we will
interchangeably use the terms \textit{turbulent thermal},
\textit{buoyant cloud}, \textit{metal-silicate cloud}, or simply
\textit{cloud}. 

\begin{figure}
\includegraphics{fig05}
\caption{\label{fig:thermals}{(a)} A 169 mL volume of an aqueous solution of NaI
($\rho_m=1502$~kg${\cdot}$m$^{-3}$) falling into fresh water, at
$\mathit{Re}=2.9\times 10^{4}$ (as defined by Equation
\eqref{eq:Re_thermal}). The time interval between each image is 0.3~s. 
{(b)} A 169 mL volume of an aqueous solution of NaI
($\rho_m=1601$~kg${\cdot}$m$^{-3}$) falling into a low viscosity
silicon oil ($\rho_s =821$~kg${\cdot}$m$^{-3}$, $\eta = 8.2\times
10^{-4}$~Pa${\cdot}$s),  at $Bo=9.3\times10^{3}$ and $\mathit{Re}=4.0\times
10^{4}$. The time interval between each image is 0.2 s.  {(c)} The
\textit{turbulent thermal} model: a volume of fluid with initial radius
$R_{0}$ and excess density $\Delta\rho$ is released at $z=0$ in an
initially motionless fluid. The thermal has a mean vertical velocity
$u_{z}$.  Its mean radius $R$ increases with the distance $z$ from the
source due to the entrainment of ambient fluid at a rate $u_{e}=\alpha
u_{z}$. {(d)} and {(e)} Close-ups on the second and forth images of
subfigure (b).  Figures (a), (b), and (c) are modified
from~\cite{deguen2020fluid}.}
\end{figure}

\subsubsection{Large-scale asymptotic behaviour}\label{sec:thermal_asymptotic}


A \textit{miscible} turbulent thermal falling in an ambient liquid of
density $\rho_s$ under the the action of a gravity field $g$ is fully
defined by its initial volume $V_{0}$ (or equivalent radius $R_0$) and
shape, solute concentration $c_{0}$ (here measured in kilogram of
solute per unit of volume), and its kinematic viscosity $\nu$ and
solute diffusivity $\kappa$ (the two latter properties being here
assumed independent of composition). Experiments such as shown in
Figure~\ref{fig:thermals} show that  ambient fluid is continuously
entrained in the thermal as it moves, such that the volume $V$ of the
thermal increases with time, while its average concentration $\langle c
\rangle$ decreases, together with the mean density difference $\Delta
\rho$ between the thermal and
surrounding~\cite{Batchelor1954,Morton1956}. Assuming that the loss of
material from the thermal (by ``detrainment'') is small~\cite{lj2019}
implies that the mass of solute within the thermal, which we denote by
$\mathcal{M} = \langle c \rangle V$, is conserved and equal to its
initial value
{\begin{equation}
 \mathcal{M}_0 = c_0 V_0.
\end{equation}}\unskip
This implies that the total buoyancy $\mathcal{B}=\sofrac{\Delta
\rho}{\rho_s} g V$ of the thermal is constant as well, and equal to
{\begin{equation}
 \mathcal{B}_0 = \frac{\Delta \rho_0}{\rho_s} g V_0,
 \label{eq:def_B0}
\end{equation}}\unskip
where $\Delta\rho_0$ is the density difference at release between the
dense and ambient liquids.

The asymptotic time evolution of the bulk properties of the thermal in
the Boussinesq limit can be predicted using dimensional
analysis\cite{Batchelor1954,Scorer1957,Turner1969}.  Far from the
source, the state of the thermal is expected to be independent of the
details of its release: it should depend on the conserved quantities
$\mathcal{M}_0$ and $\mathcal{B}_0$, but not on its initial radius or
shape. If in addition we assume that the large-scale properties are
independent of the momentum and solute diffusivity, on account of the
turbulent nature of the flow~\cite{Batchelor1954}, then dimensional
analysis predicts that the size $R$ of the thermal, vertical velocity
$u_{z}$, and average concentration $\langle c\rangle$ evolve with time
or distance from the source according to 
{\begin{eqnarray}
 R &amp;\sim&amp; \mathcal{B}_{0}^{1/4} t^{1/2} \sim z, \label{eq:r} \Seqnsplit
 u_{z} &amp;\sim&amp; \mathcal{B}_{0}^{1/4} t^{-1/2} \sim \mathcal{B}_{0}^{1/2} z^{-1}, \label{eq:w} \Seqnsplit
 \langle c\rangle &amp;\sim&amp; \mathcal{M}_0 \mathcal{B}_{0}^{-\frac{3}{4}} t^{-\frac{3}{2}} \sim \mathcal{M}_0 z^{-3}, \label{eq:c}
\end{eqnarray}}\unskip
when far from the release point.
The proportionality between $r$ and $z$ predicted by Equation
\eqref{eq:r} is often written as
{\begin{equation}
R = \alpha z, 
\label{eq:radius_depth}
\end{equation}}\unskip
where the factor $\alpha$ is the so-called \textit{entrainment
coefficient}. $\alpha$ is found to be typically in the range 0.1--0.3,
and seems to be independent of the viscosity and diffusivity as long as
the flow is turbulent~\cite{Turner1969} . One of the main factors
affecting the value of $\alpha$ appears to be the initial aspect ratio
of the thermal, with $\alpha$ being higher in oblate thermals than in
prolate thermals~\cite{lai2015numerical,morrison2023controls}. Note
also that the value of $\alpha$ depends on how the boundary of the
thermal is defined, so care should be taken when comparing estimates
from different studies. Equations \eqref{eq:r} to \eqref{eq:c} have
been extensively verified in
experiments~\cite{Scorer1957,Turner1969,zhao2013internal} and numerical
simulations~\cite{lj2019,morrison2023controls}. Close to the source,
the time evolution of the thermal radius is well described by a
relation of the form $R = R_0 + \alpha z$ as soon as $z \gtrsim
R_0$~\cite{Landeau2014}. The asymptotic
expression~\ref{eq:radius_depth} is therefore valid for $z \gtrsim
R_0/\alpha$. 

Interestingly, one can note that the Reynolds and P\'eclet numbers
estimated using the above predictions for $R(t)$ and $u_z(t)$ are
independent of time, since $R\propto t^{1/2}$ and $u_z\propto
t^{-1/2}$. Ignoring any unnecessary dimensionless coefficients, this
offers a convenient way to define the Reynolds and P\'eclet number of a
turbulent thermal as
{\begin{eqnarray}
Re &amp;=&amp; \frac{\mathcal{B}_{0}^{1/2}}{\nu}, \label{eq:Re_thermal}\Seqnsplit 
Pe_{c} &amp;=&amp; \frac{\mathcal{B}_{0}^{1/2}}{\kappa_{c}} \label{eq:Pe_thermal}. 
\end{eqnarray}}\unskip
This will be further discussed in Section~\ref{sec:self-similarity}.

\subsubsection{Turbulent entrainment model}
\label{sec:TurbulentEntrainmentModel}


Alternatively, models for the time evolution of thermals have been
built from the conservation of mass and momentum by assuming that the
ambient liquid is entrained in the thermal at a rate that is
proportional to the mean velocity of the thermal~\cite{Morton1956}.
This assumption implies that the time derivative of the volume of the
thermal is given by
{\begin{equation}
 \frac{\mathrm{d}}{\mathrm{d}t}\left( \frac{4\uppi}{3} R^3\right) = 4\uppi R^2 \alpha u_{z},
 \label{eq:thermal_mass_balance}
\end{equation}}\unskip
where $\alpha$ is the \textit{entrainment coefficient}. Since
$\sfrac{\mathrm{d}}{\mathrm{d}t}=u_{z} 
\sfrac{\mathrm{d}}{\mathrm{\mathrm{d}}z}$, integration of Equation
\eqref{eq:thermal_mass_balance} gives $R=R_0 + \alpha z$, which shows
that the coefficient introduced in Equation \eqref{eq:radius_depth} to
relate the size of the thermal to the distance it has travelled is in
fact the entrainment coefficient. Conservation of momentum gives a
prognostic equation for the vertical velocity $u_{z}$ of the centre of
mass of the cloud. Ignoring fluid drag on the thermal and a possible
loss of buoyancy in the wake, conservation of momentum simply states
that
{\begin{equation}
\frac{\mathrm{d}}{\mathrm{d}t} \left( \frac{4\uppi}{3} R^{3} \bar\rho u_{z} \right) = \frac{4\uppi}{3} R_{0}^{3} \Delta\rho_0 g,
\label{Eq:ThermalMomentumConservation1}
\end{equation}}\unskip
where $\bar\rho(t)$ is the mean density of the thermal. 
Conservation of mass implies that
{\begin{equation}
\frac{\mathrm{d}}{\mathrm{d}t}\left( \frac{4\uppi}{3} R^{3} \bar\rho
\right) = 4\uppi \rho_s R^{2} \alpha u_{z}.
\end{equation}}\unskip
This set of equations, in general, requires to be solved numerically
from initial conditions.  It can be shown that the solution has an
asymptotic limit consistent with the prediction of dimensional analysis
(Equations \eqref{eq:r} to \eqref{eq:c}). The asymptotic laws are valid
at times large compared to the timescale
{\begin{equation}
t_{g} = \frac{R_0^{2}}{\mathcal{B}_{0}^{1/2}} \sim \sqrt{\frac{\rho_s}{\Delta \rho_0}\frac{R_0}{g}}
\label{eq:t_g}
\end{equation}}\unskip
characteristic of the acceleration phase which follows the release from
rest of the thermal ($t_{g}$ is the time needed to reach from rest the
free fall velocity $\sqrt{\sofrac{\Delta \rho_{0}}{\rho_s}g R_0}$ with an
acceleration $\sofrac{\Delta \rho_0}{\rho_s} g$).

\begin{figure}
{\vspace*{4pt}}
\includegraphics{fig06}
{\vspace*{4pt}}
\caption{\label{fig:Thermal_PLIF} {(a)} Negatively buoyant mass (sodium
chloride aqueous solution with rhodamine B) falling in an unstratified
environment (pure water) and forming a turbulent thermal with a
Reynolds number of $\mathit{Re} = 4175$ and a P\'eclet number of
$\mathit{Pe}_{c}=1.0\times10^7$. The red line denotes the contour of the thermal
with a threshold equal to $c/\langle c\rangle = 0.5$, with $\langle
c\rangle$ the mean of the concentration inside the contour.  {(b)}
Probability density function (PDF) of the concentration $c/c_0$ as a
function of the concentration $c/c_0$ for the snapshot (a) at
$t/t_g=51$.  {(c)} PDF of concentration normalized by the mean of the
concentration $c/\langle c\rangle$ as a function of $c/\langle
c\rangle$ for different times between $t/t_g=20$ and $t/t_g=75$. Here,
red dashed lines show the temporal mean of the $c/\langle c\rangle$.
Inserts in (b) and (c) are the semi-log version of the PDF. Panel (a)
adapted with permission from Lherm~(2021)~\cite{Lherm2021}.}
{\vspace*{4pt}}
\end{figure}

\subsubsection{Statistics of the concentration field}
\label{sec:stat_c}

Measurements of the concentration field within turbulent miscible
thermals have been carried out by~\cite{Lherm2021} using Planar Laser
Induced Fluorescence, and further analysed by~\cite{huguet_in_prep}. 
Figure~\ref{fig:Thermal_PLIF}(a) shows a snapshot of the concentration
field in a vertical plane within a thermal having $\mathit{Re}=4175$ and
$\mathit{Pe}_{c}=1.0 \times 10^{7}$. The red line is an isocontour of the
concentration field which is used to define the boundary of the
thermal. Figure~\ref{fig:Thermal_PLIF}(a) shows that the concentration
within the thermal is far from being homogeneous, the concentration
field exhibiting small-scale, highly corrugated elongated structures. 
The standard deviation of the concentration field (${\simeq} 2.5 \times
10^{-3} c_{0}$ for this particular snapshot) is about half of the
concentration average (${\simeq} 4.2 \times 10^{-3} c_{0}$), which shows
that the concentration has a relatively wide distribution.
Concentration variations can be further quantified by calculating the
probability density function (PDF) of the concentration within the
thermal, as shown in Figure~\ref{fig:Thermal_PLIF}(b). The probability
density function is a strongly decreasing function of $c$, with a
distribution dominated by the small concentrations.

Dimensional analysis predicts that far from the source the moment of
order $k$ of the concentration is of the form
{\begin{equation}
\langle c^k \rangle = \left(\mathcal{M}_{0} \mathcal{B}_{0}^{-\frac{3}{4}} 
t^{-\frac{3}{2} }\right)^k f_{c,k}(\mathit{Re},\mathit{Pe}_{c}),
\label{eq:moments_c}
\end{equation}}\unskip
where \textit{Re} and $\mathit{Pe}_{c}$ are defined according to Equations
\eqref{eq:Re_thermal} and \eqref{eq:Pe_thermal}, and $f_{c,k}$
non-dimensional factors which may be functions of \textit{Re} and
$\mathit{Pe}_{c}$~\cite{huguet_in_prep}. It follows that the moments of
$c/\langle c \rangle$, 
{\begin{equation}
 \langle (c/\langle c \rangle)^k \rangle = 
 \frac{\langle c^k \rangle}{\langle c \rangle^{k}} = 
 \frac{f_{c,k}(\mathit{Re},\mathit{Pe}_{c})}{f_{c,1}^{k}(\mathit{Re},\mathit{Pe}_{c})},
\end{equation}}\unskip
are predicted to be independent of time, which implies that the
probability density function of $c/\langle c \rangle$ is independent of
time as well. Figure~\ref{fig:Thermal_PLIF}(c) shows that the PDF of
$c/\langle c \rangle$ calculated at different times in the asymptotic
regime  collapse on a well-defined master curve, thereby validating the
above prediction. The PDFs of $c/\langle c \rangle$ are well modelled
by a decreasing exponential given by
{\begin{equation}
p(c/\langle c \rangle) = \exp\left( - \frac{c}{\langle c \rangle}\right)
\label{eq:pdf_c_norm}
\end{equation}}\unskip
with $\langle c \rangle$ given by Equation \eqref{eq:c} (with a
dimensionless prefactor adjusted to measurements of $\langle c
\rangle$). The time average of the PDFs (dashed red line) indeed lies
very close to the prediction of Equation~\eqref{eq:pdf_c_norm} (grey
line). Qualitatively, the stationary nature of the distribution of
$c/\langle c \rangle$ can be understood as resulting from the competing
effects of homogenization within the thermal and continuous addition of
ambient liquid through entrainment. This is well illustrated by the
fact that the standard deviation is always commensurate with the
average concentration, $\sqrt{\langle c^2 \rangle - \langle c
\rangle^2} \sim \langle c \rangle$. The details of the mechanisms by
which the stationary distribution is maintained, and the reasons why
the resulting distribution is exponential, have yet to be understood
and modelled.

\subsubsection{A note on self-similarity}
\label{sec:self-similarity}

The state at time $t$ of a miscible thermal can be characterised by its
Reynolds and P\'eclet numbers.  As noted at the end of
Section~\ref{sec:thermal_asymptotic}, it so happens that the
instantaneous values of \textit{Re} and $\mathit{Pe}_{c}$ estimated with the asymptotic
predictions for $R(t)$ and $u_z(t)$ are independent of time. It can be
argued that the self-similar behaviour of turbulent thermals is in fact
linked with the time invariance of the instantaneous \textit{Re} and $\mathit{Pe}_{c}$.
This allows the velocity and composition fields to be self-similar in
time, in the sense that, when normalised by their spatial average $u_z$
and $\langle c \rangle$, the statistical distribution of the velocity
and concentration can be constant in time.

Adding ingredients like surface tension or rotation might impede
self-similarity. One can note that the Weber number based on the
instantaneous size and velocity in the asymptotic regime (Equations
\eqref{eq:r} and \eqref{eq:w}) is a decreasing function of time and
distance given by
{\begin{equation}
\mathit{We} \sim \frac{\rho_s \mathcal{B}_{0}}{\sigma R_{0}} \left(\frac{t}{t_{g}} \right)^{-1/2}
\sim \mathit{Bo} \left(\frac{t}{t_{g}} \right)^{-1/2} \sim \alpha \mathit{Bo} \frac{R_0}{z} ,
\label{eq:We_Bo}
\end{equation}}\unskip
where here the Bond number is defined as 
{\begin{equation}
 \mathit{Bo} = \frac{\Delta\rho_0 g R_0^2}{\sigma } \sim \frac{\rho_{0} \mathcal{B}_{0}}{\sigma R_{0}} .
 \label{eq:Bo_thermal}
\end{equation}}\unskip
The proportionality between \textit{We} and \textit{Bo} reflects the
force balance (inertia ${\sim}$ buoyancy) which controls the velocity of
the thermal.

Similarly, the Rossby number, which measures the relative importance of
inertia and Coriolis force, writes 
{\begin{equation}
\mathit{Ro} = \frac{u_z(t)}{2 \Omega R(t)} \sim \frac{1/\Omega}{t} \sim 
\left( \frac{\mathcal{B}_{0}^{1/4} \Omega^{-1/2}}{ z} \right)^{2}
\end{equation}}\unskip
when estimated with the asymptotic expressions for $R(t)$ and
$u_{z}(t)$. This predicts that the Coriolis force becomes comparable to
inertia after a time interval ${\sim} 1$ day, or, equivalently, when the
distance travelled by the thermal is larger than ${\sim}
\mathcal{B}_{0}^{1/4}
\Omega^{-1/2}$~\cite{ayotte1994motion,helfrich1994thermals,kriaa2022effects}. 

Finally, if the thermal is made of buoyant droplets or particles, then
an additional relevant number is the Rouse number $\mathcal{R}$ defined
as the ratio of the sedimentation velocity of the droplets or
particles, $u_{\mathrm{sed}}$, to a characteristic velocity scale of the
surrounding fluid.  Taking for this velocity scale the asymptotic
scaling of $u_z$, the Rouse number writes 
{\begin{equation}
\mathcal{R} = \frac{u_{\mathrm{sed}}}{u_{z}} \sim \frac{u_{\mathrm{sed}}}
{\mathcal{B}_{0}^{1/4} } t^{1/2} \sim \frac{u_{\mathrm{sed}}}{\mathcal{B}_{0}^{1/2} } z,
\label{eq:Rouse_number}
\end{equation}}\unskip
which is again time-dependent.

The fact that some of the relevant dimensionless parameters are
time-dependent must prevent true self-similarity. The implications are
further discussed in the following subsections. 

\subsection{Immiscible turbulent thermals}

Section~\ref{sec:impacts} suggests two possible scenarios for the fate
of the impactor core. For projectiles smaller than $300$~km in radius,
or falling faster than twice the escape speed, the core fragments at
the top of the magma ocean during the cratering
process~\cite{MallerInRev}. In this case, the metal drops behave
collectively and generate a turbulent cloud of metal drops and
entrained silicates that fall in the magma ocean~\cite{Deguen2011}.
Such particle-laden clouds behave as the miscible thermals described in
the previous section~\cite{Bush2003,kriaa2022effects,Kriaa2024}. In
contrast, the cores of larger impactors coming onto Earth at the escape
speed may not fragment during the impact~\cite{MallerInRev}. In this
scenario, the impactor core, still intact after the cratering stage,
falls by buoyancy in an immiscible pool of liquid silicates under
extremely large Reynolds and Weber numbers, as discussed in
Section~\ref{sec:magma_oceean_rationale}.

\subsubsection{Experimental test of the turbulent entrainment model}

Figure~\ref{fig:thermals}b shows snapshots from an experiment in which
a negatively buoyant volume of sodium iodide is released into a
low-viscosity silicon oil, under conditions otherwise similar to the
experiment shown on Figure~\ref{fig:thermals}a.  The NaI solution and
the silicon oil are \textit{immiscible}, so this experiment is one
important step closer to the core-mantle differentiation configuration.
The experimental fluids and configuration have been chosen so as to
maximize the values of the Bond and Reynolds numbers, which are
$\mathit{Bo}=9.3\times10^{3}$ and $\mathit{Re}=4.0\times 10^{4}$ (with \textit{Bo} and \textit{Re} as
defined by Equations \eqref{eq:Bo_thermal} and \eqref{eq:Re_thermal}). 
We use here a Bond number rather than a Weber number, which is
time-dependent in these experiments; as shown by Equation
\eqref{eq:We_Bo}, \textit{Bo} offers a convenient measure of the order of
magnitude of $We(t)$.  The density ratio is $\rho_{m}/\rho_{s}=1.95$
(close to metal-silicate), and the viscosity ratio is
$\eta_{m}/\eta_{s}=1.2$.  The large-scale evolution of the negatively
buoyant volume is strikingly similar to what has been observed in the
\textit{miscible} experiment: the volume of the negatively buoyant
fluid increases with distance, which indicates that it entrains and
incorporates silicon oil. The analogy with miscible thermals also holds
quantitatively: the evolution of the vertical position, size, and
velocity of immiscible thermals are well predicted by the turbulent
entrainment model.  This has been verified in experiments similar to
those shown in Figure~\ref{fig:thermals}b~\cite{Deguen2014,Landeau2014}
as well as in experiments in which a volume of galinstan
($\rho=6440$~kg${\cdot}$m$^{-3}$) is released in a pool of UCON oil
($\rho=1050$~kg${\cdot}$m$^{-3}$) of different
viscosities~\cite{Wacheul2014,wacheul2017fall,Wacheul2018}: when 
$\mathit{Bo} \gtrsim 100$, the entrainment coefficient $\alpha=0.25\pm0.1$ is
indistinguishable from that of miscible thermals.

\subsubsection{Fragmentation}

Immiscible thermals strongly differ from miscible thermals when looking
at the small scales: in the former, the released liquid eventually
breaks into drops. In the limit of large Weber numbers, the interface
between the released and ambient liquids is stretched by the turbulence
until it becomes fractal~\cite{Deguen2014}. As this fractal interface
builds up, ligaments emerge everywhere inside the thermal. All the
ligaments break up almost simultaneously, leading to the sudden
formation of droplets everywhere in the thermal
(Figure~\ref{fig:thermals}b)~\cite{Deguen2014,Landeau2014}. This sudden
and global fragmentation strikingly contrasts with the fragmentation of
a liquid in air, which is often gradual and localised at a well-defined
interface~\cite{Marmottant2004,Villermaux2020}. 

At the very large values of \textit{We} or \textit{Bo} estimated in
Section~\ref{sec:magma_oceean_rationale}, the immiscible thermal is
expected to be in  this asymptotic regime, in which the impactor core
forms a turbulent immiscible thermal and suddenly breaks up into drops.
In experiments, the breakup into drops occurs at a constant depth of
$4$ to $8$ initial radii when $\mathit{Bo}\gtrsim 100$ in the explored range
of \textit{Bo} ($\lesssim 10^4$)~\cite{Landeau2014}. Coupling this result with
the disruption caused by the impact stage,~\cite{MallerInRev} predicts
the depth of breakup as a function of the impactor size and velocity.
They found that, on Earth, only impactors larger than about $2000$~km
in radius, and coming slowly at the escape speed, reach the bottom of
the magma ocean unbroken. Any smaller or faster impactor is fragmented
at the bottom of the ocean. 

Once fragmented, the degree of chemical equilibration of the impactor
core may depends in part on the drop size. In an isotropic and
homogeneous turbulent flow with typical Weber number \textit{We}, the drop
size $d$ is comparable to the Hinze scale~\cite{Hinze1955}, which is
obtained assuming a Kolmogorov cascade such that the speed at scale $l$
is $u_l\sim (l/L)^{1/3}\,U$, where $L$ and $U$ are the typical length
and velocity at the largest scales. Then, at scale $l$, the local Weber
number is 
{\begin{equation}
 \mathit{We}_l=\frac{\rho_s\,u_l^2\,l}{\sigma} \sim We\,\left(\frac{l}{L}\right)^{5/3}.\label{eq:local_Weber}
\end{equation}}\unskip
The Hinze scale is such that inertia balances surface tension, {i.e.}\ 
$\mathit{We}_l \sim 1$. From Equation \eqref{eq:local_Weber}, one obtains 
{\begin{equation}
 l_H\sim L\, \mathit{We}^{-3/5}.\label{eq:Hinze}
\end{equation}}\unskip
Assuming that the thermal speed arises from the buoyancy force only, we
can take $L=R(t)$ and use Equation \eqref{eq:We_Bo} for \textit{We} to find
{\begin{equation}
 l_H\sim R_0 (\alpha^2 Bo)^{-3/5} \left( \frac{R}{R_0} \right)^{8/5} \sim 
 \left(\frac{\alpha^2 \Delta\rho_0 g}{\sigma}\right)^{-3/5} 
 \left( \frac{R}{R_0} \right)^{8/5} R_0^{-1/5},
\end{equation}}\unskip
which is only weakly dependent on the initial size $R_0$ of the metal
cloud.  Estimated at $R/R_0$ between 4 and 8 (the fragmentation
distance experimentally observed) with $R_0$ between 10~km and 1000~km,
the Hinze scale is between a few millimeters and a few centimeters. 
Drops of this size are small enough to fully equilibrate with the
surrounding silicates~\cite{Rubie2003,Ichikawa2010,Qaddah2019}. 

However, scaling \eqref{eq:Hinze} remains to be tested experimentally
for immiscible thermals.  The relevance of the Kolmogorov cascade to
turbulent thermals is questionable (turbulence is likely affected by
buoyancy forces and by the shear associated with the vortex ring
structure of thermals). In addition, when fragmentation occurs during
the cratering stage (Section~\ref{sec:impacts}), multiple choices are
possible for the length scale $L$ and the typical velocity $U$ in the
Weber number and the Hinze scale (\ref{eq:Hinze}).

Experiments using galinstan and UCON oil have been used to investigate
the effect of the viscosity contrast on the fragmentation
behaviour~\cite{Wacheul2014,wacheul2017fall,Wacheul2018,wacheul2019etude}, though at smaller values of the Bond number ($Bo&lt;100$).
Fragmentation has been observed to happen in a single global event, as
described above for a viscosity ratio ${\sim} 1$, in experiments with a
viscosity ratio $\eta_s/\eta_m$ equal to $0.383$, $8.37$, and $72.5$
(where $\eta$ denotes the dynamics viscosity).  The distribution of
drop sizes does not seem to be affected by the viscosity ratio in this
range, and the average drop size seems to be in agreement with Hinze
scaling (Equation~(\ref{eq:Hinze}))~\cite{Wacheul2018}. At the highest
viscosity ratio investigated ($\eta_s/\eta_m=691$), fragmentation has
been observed to result from a different regime, the
\textit{Rayleigh--Taylor regime}~\cite{wacheul2019etude}, in which the
initial deformation of the dense liquid results from large scale
Rayleigh--Taylor instabilities. However, the Reynolds number of these
experiments at $\eta_s/\eta_m=691$ is moderate ($\mathit{Re}&lt;80$)\footnote{The
Reynolds number given here is calculated using
Equation~(\ref{eq:Re_thermal}), which gives a value twice larger than
quoted in~\cite{Wacheul2018,wacheul2019etude}.}, and it is not clear
whether this regime persists at higher values of \textit{Re}. In fact \textit{Re} is
inversely correlated with the viscosity ratio in this suite of
experiments (because the viscosity ratio can only be increased by
increasing the ambient liquid viscosity), and it is therefore not
straightforward to discriminate between the effects of the viscosity
ratio and Reynolds number.

\subsubsection{Post-fragmentation dynamics}
If the metal phase does fragment before reaching the base of the magma
ocean, then the dynamics of the resulting cloud of droplets depends on
whether the droplets have a collective behaviour or fall without
significant interaction between them.  This is controlled by the Rouse
number $\mathcal{R}$ defined as the ratio of the settling velocity of
individual drops to the magnitude of the flow and, to a lesser extent,
by the drop particle number $Re_{d}$ (based on drop size and settling
velocity) which affects the flow around each drop and the interactions
between\break drops. 

In experiments on immiscible turbulent thermals such as shown in
Figure~\ref{fig:thermals} or described
by~\cite{Deguen2014,Landeau2014}, the size of the droplets resulting
from fragmentation was such that, given the height of the tank, the
cloud velocity remained larger than the settling velocity for most of
the duration of the experiments, except when the cloud was approaching
the base of the tank. There it was observed that the largest drops
separated from the bulk of the thermal while smaller drops were still
carried by the vortex ring flow, but the limited observations available
did not allow any quantitative and systematic characterisation of the
transition.  This observation is consistent with the fact that the
Rouse number estimated using the asymptotic scaling law for the
turbulent thermal velocity (Equation \eqref{eq:w}) is a decreasing
function of time, as shown by\break Equation \eqref{eq:Rouse_number}.

Using clouds made of solid particles rather than drops allows to
characterize the effect of the settling velocity of the suspended phase
in a much simpler and systematic way. Experiments have shown that
particle clouds behave as turbulent thermals when $\mathcal{R}\lesssim
1$ but form bowl-shaped \textit{particle swarms} when
$\mathcal{R}\gtrsim 1$, falling at a velocity near the settling
velocity of individual
particles~\cite{Bush2003,Deguen2011,kriaa2022effects}. When estimated
using the asymptotic scaling law for the turbulent thermal velocity
(Equation \eqref{eq:w}), the Rouse number is given by Equation
\eqref{eq:Rouse_number}. It is a decreasing function of time, due to
the decrease of the cloud velocity resulting from the gradual
entrainment of ambient liquid with no momentum. It is found from
Equation \eqref{eq:Rouse_number} that $\mathcal{R}$ reaches the value
of 1 when the cloud has travelled a distance of the order of
${\mathcal{B}_0^{1/2}}/{u_{\mathrm{sed}}}$. \cite{Bush2003} found
experimentally that the fallout distance  $z_{\mathrm{fall}}$ at which
particle-ladden turbulent thermals transition to particle swarms indeed
scale as ${\mathcal{B}_0^{1/2}}/{u_{\mathrm{sed}}}$, with a prefactor
${\simeq} 9\pm 2$, when the particle Reynolds number is larger than
${\simeq} 4.4$.

When the drop Reynolds number is large, then the settling velocity is
$u_{s} \sim \sqrt{\sfrac{\Delta\rho}{\rho_{s}}g a}$ where $a$ is the
particule or drop radius. Using this scaling,~\cite{Bush2003} found
experimentally that the ratio of the fallout distance to particle
radius is only a function of the number of particles $N_{p}$,
$z_{\mathrm{fall}}/a\simeq (9\pm2) N_{p}^{1/2}$. Re-arranging, and
denoting by $\phi_0$ the initial metal fraction in the cloud, we can
write the fallout distance as a function of the initial radius $R_{0}$
of the metal and of the drop\break size as
{\begin{equation}
 \frac{z_{\mathrm{fall}}}{R_{0}} \simeq (9\pm 2) \left(\phi_0\frac{R_{0}}{a}\right)^{1/2},
\end{equation}}\unskip
which is more convenient for application to metal-silicate separation
in a magma ocean. Taking $a\sim 1$ cm, $z_{\mathrm{fall}}$ would only be
smaller than ${\sim} 3000$~km (the thickness of the mantle) if $R_{0}$ is
smaller than ${\sim} 1$~km. The critical value of $R_{0}$ decreases to
${\sim}500$ m if we rather assume $a \sim 1$ mm. This suggests that
fallout of drops from a metal-silicate thermal was unimportant, except
for very small metal volumes.

\subsection{Additional effects}
\subsubsection{Impact-generated thermals}

The experiments discussed in Sections~\ref{sec:miscible_thermals}
to~\ref{subsec:rotation} have been performed by releasing volumes of
dense liquid \textit{within} the host liquid, without accounting for
the impact and cratering phase described in Section~\ref{sec:impacts}.
As discussed above, the assumption is that the impact leaves the
impactor material as a coherent and localized mass which then evolves
as a turbulent thermal.

Experiments using both miscible liquid pairs~\cite{landeau2021} and
immiscible liquid pairs (Figure~\ref{fig:themals_impacts}a) show that
after the collapse of the crater the liquid released by the impact does
indeed evolve as a turbulent thermal, at least as long as the Froude
number is not too large.  This was the case for all experiments studied
by~\cite{landeau2021}, with $\mathit{Fr}\in [4,100]$. Landeau 
{et~al.}~\cite{landeau2021} quantified the amount of dispersion following
the cratering phase, and found that the evolution of the equivalent
radius of the dense liquid cloud evolves as
{\begin{equation}
 R = R_0 + \alpha z,
\end{equation}}\unskip
as predicted by the turbulent entrainment model, if the radius $R_0$ of
the cloud at $z=0$ is taken to be
{\begin{equation}
 R_0 = 0.29 R_i \left( \frac{\rho_m-\rho_s}{\rho_s} \right)^{-0.19} \mathit{Fr}^{0.4}
\end{equation}}\unskip
instead of the impactor radius $R_i$. This accounts for the dispersion
induced by the impact which is, as expected, an increasing function of
\textit{Fr}.

\begin{figure}
\includegraphics{fig07}
\caption{\label{fig:themals_impacts}(a) Time sequence from an
experiment in which a volume of water ($R_i=2.3$ cm, $\nu=10^{-6}$
m$^2{\cdot}$s$^{-1}$, $\rho_i=1000$~kg${\cdot}$m$^{-3}$) strikes a bath
of low-viscosity silicon oil ($\nu=1.2\times 10^{-6}$
m$^2{\cdot}$s$^{-1}$, $\rho_t=821$~kg${\cdot}$m$^{-3}$) at a velocity
$U_i=1.1$ m${\cdot}$s$^{-1}$. The time step is $0.1$ s between the
three first images, and then $0.3$ s.  (b) Time sequence from an
experiment in which a volume of NaCl solution ($R_i=3.6$~cm,
$\nu=10^{-6}$ m$^2{\cdot}$s$^{-1}$, $\rho_i=1088$~kg${\cdot}$m$^{-3}$)
strikes a bath of water at a velocity $U_i=1.6$ m${\cdot}$s$^{-1}$ with
an angle $\theta=52.5$\textdegree~with  respect to the horizontal.}
\end{figure}

Figure~\ref{fig:themals_impacts}(b) shows snapshots from an experiment
in which a volume of NaCl solution ($R_i=3.6$~cm,
$\rho_i=1088$~kg${\cdot}$m$^{-3}$, appearing blue) strikes the surface
of a water pool with an impact angle $\theta=52.5$\textdegree, at
$\mathit{Fr}=7.1$ and $\mathit{Re}=5.8\times 10^4$. In contrast with the snapshots shown
in Figure~\ref{fig:ObliqueImpacts}, which focused on the cratering
dynamics, we focus here on the longer term evolution of the liquid
delivered by the impact. Superimposed on the images are predictions
from a turbulent entrainment model allowing for a non-vertical
trajectory of the dense liquid (an extension of the model described in
Section~\ref{sec:TurbulentEntrainmentModel}),  with the initial
conditions (size of the cloud and velocity of its center of mass)
fitted to match the observations. The solid red line shows the model
prediction for the trajectory of the center of mass of the cloud, the
position of which is shown with a red cross at the times of the
snapshots;  the dashed red lines show the positions of the left and
right edges of the cloud according to the turbulent entrainment model.
The evolution of the dense liquid cloud is well described by the
turbulent entrainment model, though one can notice that a significant
amount of impactor liquid is dispersed at impact and left behind the
leading cloud. Preliminary investigations shows that the initial
conditions to be given to the turbulent entrainment model are functions
of $\theta$ and \textit{Fr}. A parameterization of these initial conditions,
coupled to the turbulent entrainment model, will allow to build a
predictive model for the dispersion following an oblique impact. 


\subsubsection{Effect of planetary rotation}\label{subsec:rotation}

As discussed in Section~\ref{sec:self-similarity}, the Rossby number
based on the instantaneous size and velocity predicted in the
asymptotic regime is a decreasing function of time or, equivalently,
distance travelled by the cloud. Rotation starts to have a significant
effect when the thermal reaches a depth $z_\Omega =
\mathcal{B}_{0}^{1/4} \Omega^{-1/2} \propto R_0^{3/4}
\Omega^{-1/2}$~\cite{ayotte1994motion,helfrich1994thermals,kriaa2022effects}.
This critical depth can also be written as
{\begin{equation}
 z_\Omega \sim \left( \frac{1\ \text{day}}{t_g}\right)^{1/2} R_0 ,
 \label{eq:z_Omega}
\end{equation}}\unskip
where $t_g$ is the acceleration time scale defined in Equation
\eqref{eq:t_g}. The prefactor in Equation \eqref{eq:z_Omega} can be
estimated from the asymptotic solutions of the turbulent entrainment
model (e.g.~\cite{Escudier1973,Deguen2014}), and is found to be close
to 1, depending on the parameters, including $\alpha$, of the
entrainment model.

During core formation, the rotation of the Earth was likely faster than
it is today. Taking $R_0$ between 10~km and 1000~km and a rotation rate
twice as fast as today gives $z_\Omega$ between ${\sim} 400$~km and ${\sim}
12\,000$~km, which needs to be compared to the thickness of the magma
ocean, which is at most ${\sim} 3000$~km. $z_\Omega$ is equal to the
thickness $H$ of the magma ocean for  $R_0 \sim H^{4/3} \Omega^{2/3}
(\rho_s /( g \Delta \rho))^{1/3}$. Taking $H= 2800$~km and $\Omega$
twice the present day value gives $R_0\sim 500$~km, in agreement with a
more complete calculation made by~\cite{Kriaa2024}: impactors with a
core radius larger than 500~km will not be affected much by rotation;
smaller impactors can be affected, depending on the thickness of the
magma ocean.  Taking a rotation rate 5 times larger than today
increases $z_\Omega$ by a factor $\sqrt{5/2}\simeq 1.6$, while
increasing by a factor $(5/2)^{2/3}\simeq 1.85$ the maximum radius at
which rotation can affect the flow. In contrast, assuming a magma ocean
thickness of 1400~km instead of 2800~km reduces the critical radius to
${\simeq} 200$~km.

Experiments on miscible turbulent
thermals~\cite{ayotte1994motion,helfrich1994thermals} and cloud
particles~\cite{kriaa2022effects,Kriaa2024} with a rotation axis
parallel to the gravity have shown that the thermal stops expanding
laterally after reaching the critical depth $z_{\Omega}$, and then
grows as a column which keeps a constant radius ${\sim} \alpha
z_{\Omega}$ but expands in the vertical direction. The amount of
dispersion in this regime has been quantified and modelled
by~\cite{Kriaa2024}. For small impactors falling in deep magma oceans,
the amount of dispersion is significantly lowered compared to the
prediction of the turbulent thermal model. In details, the results
depends on the initial conditions. Assuming that the impactor is left
with no vertical momentum at the end of the impact phase, and taking
into account the dispersion produced by the impact~\cite{landeau2021},
Kriaa {et~al.}~\cite{Kriaa2024} found that in a 2800~km thick
magma ocean with a planetary rotation rate twice that of the Earth
today, the evolution of all impactors with a radius less than 450~km is
affected by rotation.  For these impactors, rotation reduces the mass
of silicate that is entrained and equilibrated with the impactor core. 
However, because dilution of the impactor core is very efficient for
such small impactors, the degree of equilibration is already close to
one and the effect of rotation on equilibration is quite moderate. It
only becomes significant for highly siderophile elements with a
metal/silicate partition coefficient $D \gtrsim 10^3$, for which the
degree of equilibration is reduced by a factor of up to 3 when rotation
is included. The effect of rotation is stronger at faster rotation
rates: the evolution of all impactors is affected by rotation if the
rotation rate was 5 times that of Earth today, again with a significant
effect on the degree of equilibration for elements with a
metal/silicate partition coefficient $D \gtrsim
10^3$~\cite{Kriaa2024}. In contrast, the effect of rotation is less
important in a shallower magma ocean: the critical impactor radius is
reduced to 250~km in a 1400~km thick magma ocean~\cite{Kriaa2024}.

The experiments discussed
above~\cite{ayotte1994motion,helfrich1994thermals,kriaa2022effects,Kriaa2024} 
have been done with a rotation axis parallel to the gravity. When
transposed to planetary settings, this corresponds to impacts in the
polar regions, which probably account for only a moderate fraction of
the total mass accreted.  Experiments with a rotation axis tilted with
respect to the gravity field will be more challenging.  In principle,
this could be done in experiments with a rotation rate high enough to
give a centrifuge acceleration close to or greater than $g$. 

\subsection{Chemical equilibration within the cloud}
\label{subsection:ChemicalEquilibration}

Published models of metal-silicate equilibration assume that the metal
phase fully equilibrates with the entrained
silicates~\cite{Deguen2014,landeau2021}, in proportions predicted by
the turbulent entrainment and impact mixing models. This allows to
build simple models of mass transfer but ignores the concentration
heterogeneities that must be continuously produced by the entrainment
of fresh, unequilibrated material.  For a chemical element that does
not partition between metal and silicates (metal/silicate partition
coefficient $D=1$), the assumption of full equilibration amounts to
assuming that the concentration within the thermal is homogeneous. This
is inconsistent with the probability density function of the
concentration field discussed in Section~\ref{sec:stat_c}. This is also
inconsistent with measurements made by~\cite{Wacheul2018} on the
temperature of the metal phase in experiments in which a volume of
galinstan falls as a turbulent immiscible thermal in a pool of UCON
oil. It is indeed found in these experiments that the temperature of
the galinstan is higher than predicted by full equilibration, which
corresponds to less efficient heat transfer~\cite{Wacheul2018}.

The implications for the degree of equilibration of the concentration
density function given in Equation \eqref{eq:pdf_c_norm} need to be
further quantified and generalized to two-phase flows and chemical
elements with a partition coefficient $D \neq 1$. We can already note
that the median of the distribution of $c$ corresponding to Equation
\eqref{eq:pdf_c_norm} is equal to $(\ln 2) \langle c \rangle \simeq
0.69 \langle c \rangle$, which means that half of the thermal volume
has concentration levels smaller than 69\% of the value $\langle c
\rangle$ assumed in perfect equilibration models. Taking into account
the distribution of $c$ within the thermal will therefore result in a
decrease in the efficiency of equilibration.


\subsection{Reaching the bottom of the magma ocean}
\label{sec:sedimentation}

The metal cloud may still have a velocity up to several hundreds of
m${\cdot}$s$^{-1}$ when it reaches the bottom of the magma ocean, depending on
the size of the impactor and depth of the magma ocean. This could cause
a ``secondary impact'' which, while subsonic, can be energetic enough to
result in significant additional dispersion~\cite{Deguen2011c}. 
Depending on the extent of mantle melting, several scenarios can be
envisaged~\cite{Deguen2014}. If the magma ocean extends down to the
core, the metal cloud can entrain mantle silicates into the core,
leading to further chemical exchange at high pressure and
temperature~\cite{landeau2016core}. If instead a solid portion of
mantle subsists, the metal-silicate mixture may spread laterally at the
base of the magma ocean as a turbulent gravity current if the bottom of
the magma ocean is flat, or may possibly be driven back upward by
inertia and re-suspended if the magma ocean has an hemispherical
shape~\cite{Deguen2011c}. If this is the case, then one may wonder how
long the metal drops can be kept in suspension, given that typical
convective velocity in a magma ocean are estimated to be 
${\sim} 10$ m${\cdot}$s$^{-1}$~\cite{solomatov2015magma}, significantly larger than the
settling velocity $v_{s}$ of a mm-to-cm sized drops, which is around
0.3--0.5~m${\cdot}$s$^{-1}$~\cite{Rubie2003,Samuel2012}.

Laboratory experiments have shown that sedimentation from a convective
layer happens on a timescale ${\sim} H/v_{s}$, where $H$ is the thickness
of the convecting layer, even if the settling velocity is much smaller
than the convective
velocity~\cite{martin1988crystal,martin1989fluid,lavorel2009sedimentation}.
This results from the presence of a boundary layer at the bottom of the
convective layer in which the vertical velocity of the convecting
liquid goes to zero, thus allowing the particles to sediment in spite
of vigorous convection in the bulk~\cite{martin1988crystal}.  With
$v_{s} \sim 0.3\ndash 0.5$~m${\cdot}$s$^{-1}$~\cite{Rubie2003,Samuel2012}, the
sedimentation time from a 2000~km deep magma ocean would be $H/v_{s}
\sim 2$ months. This is small compared to the lifetime of a magma ocean
(${\gtrsim} 1000$ years~\cite{solomatov2015magma}). 

Laboratory experiments using solid particles have also shown that the
volumetric fraction of suspended particles reach a non-zero equilibrium
value, set by a balance between particle sedimentation and
re-entrainment of particles from the layer of deposited
particles~\cite{solomatov1993entrainment,lavorel2009sedimentation}. The
processes behind entrainment from a dense liquid layer likely differs
from those behind particles entrainment, because of the necessity to
first produce new drops that could then be entrained by the convection.
It is likely to be less efficient as it requires the convective
stresses to overcome the surface tension in addition to the
gravitational forces. As an illustration, we can however use the model
developed for solid
particles~\cite{solomatov1993entrainment,lavorel2009sedimentation} to
estimate an upper bound on the equilibrium metal volume fraction in a
convecting magma ocean.  Using Equation~(12)
from~\cite{lavorel2009sedimentation} with a magma ocean thickness of
2000~km and a Rayleigh number in the range $10^{28}$--$10^{32}$, we find
the predicted equilibrium metal volume fraction to be in the range
$10^{-10}$--$10^{-5}$, suggesting that only a very small fraction of metal
could have remained in the mantle.

Planetary rotation may affect the settling and dispersion of the metal
phase, in particular through its effect on the regime and strength of
magma ocean convection.  This has not been investigated experimentally,
but numerical simulations of the gravitational settling of particles in
a vigorously convecting layer suggest that the effect of rotation can
be strong if convection is strongly influenced by the Coriolis force,
which depends on the state of the magma ocean~\cite{maas2019dynamics}. 
In this case, Maas {et~al.}~\cite{maas2021fate} found that the
degree of dispersion of the metal depends on the latitude of the
impact: the degree of dispersion is rather limited in the case of
impacts in the polar regions, and much higher in the case of impacts in
the equatorial band. In these simulations, this is due to the presence
of a strong zonal flow that transports and disperses the metal
throughout the equatorial band~\cite{maas2021fate}.

\section{Metal-silicate separation in a solid mantle}\label{sec:solidmantle}

\subsection{Rayleigh--Taylor instability and diapirism}

\subsubsection{A simple reference model}
\label{subsec:diapir_ref_model}


Let us consider a situation in which a global magma ocean or a more
localized magma pond is separated from the core by an at least
partially crystallised silicate layer. As discussed in
Section~\ref{sec:sedimentation}, most of the metal is expected to
sediment from the magma ocean in a timescale of a few months.  This is
small compared to the time needed to crystallise the magma ocean
(${\gtrsim} 1000$ years~\cite{solomatov2015magma}). A layer of metal
would therefore accumulate at the base of the magma
ocean~\cite{stevenson1981models}. Molten iron being denser than the
underlying silicates, this situation is gravitationally unstable and
gives rise to Rayleigh--Taylor
instabilities~\cite{elsasser1963early,stevenson1981models,sasaki1986metal,ida1987earth,stevenson1990fluid}.

\begin{figure}
{\vspace*{-2pt}}
\includegraphics{fig08}
\caption{\label{fig:RT_diapirs}(a) Snapshots extracted and modified from 
Figure~2
of~\cite{whitehead1975}, with permission from Wiley (CCC license number
5781330644178).  The dark liquid is glycerin ($\eta=175$~Pa${\cdot}$s,
$\rho=1250$~kg${\cdot}$m$^{-3}$) and the ambient liquid (transparent)
is silicon oil ($\eta=5520$~Pa${\cdot}$s, $\rho=920$~kg${\cdot}$m$^{-3}$). 
Note that in the original paper the images were flipped upside-down to
be consistent with the problem studied (plume formation from a thin
layer at the base of the mantle).  Here the images have been put back
to their original orientation. {(b)} A diapir formed in a three-layer
experiment with initial conditions consisting of a layer of gallium
(appearing dark) sandwiched between two layers of glucose solutions
with different viscosities and densities, the lower layer (s2) being
denser and more viscous than the upper layer (s1). {(c)} Similar to (b)
except that the intermediate layer is a mixture of gallium and liquid
from the upper layer.  The density and viscosity of the working liquids
are $\rho=5910$~kg${\cdot}$m$^{-3}$, $\eta_m=10^{-3}$~Pa${\cdot}$s for gallium,
and $\rho=1256$~kg${\cdot}$m$^{-3}$, $\eta_s=0.1$~Pa${\cdot}$s and
$\rho=1417$~kg${\cdot}$m$^{-3}$, $\eta_s=11.75$~Pa${\cdot}$s for the s1 and s2
glucose solutions.  (b) and (c) are parts of Figures~1 and~2
of~\cite{fleck2018iron}, under the CC BY 4.0 license
(\url{http://creativecommons.org/licenses/by/4.0/}).}
\end{figure}


In the relevant limit in which the lower layer is much more viscous
than the thin upper layer, linear stability analysis shows that the
most unstable mode has a wavelength $\lambda$ that is commensurate with
the thickness $H$ of the lower viscous layer if it can be approximated
as isoviscous~\cite{chandrasekhar2013hydrodynamic}. The timescale of
instability growth is
{\begin{equation}
\tau_{RT} \sim \frac{\eta_s}{\Delta\rho g \lambda} \sim \left(\frac{\eta_{s}}
{10^{18}~\mathrm{Pa}{\cdot}\mathrm{s}} \right) 
\left(\frac{10^{3}~\mathrm{km}}{\lambda} \right) \times 1  \text{ year},
 \label{eq:RT_timescale}
\end{equation}}\unskip
where $\eta_s$ is here the effective viscosity of the solid mantle
layer, and $\Delta \rho$ the density difference between the metal and
silicate layers. With a viscosity of $10^{18}$~Pa${\cdot}$s, which corresponds
to the estimate by~\cite{karato1997core} just below the solidus, and a
wavelength of ${\sim} 10^{3}$~km, $\tau_{RT}$ is on the order of 1 year. 
A viscosity of $10^{18}$~Pa${\cdot}$s corresponds to the estimate
by~\cite{karato1997core} just below the solidus.  The viscosity may be
orders of magnitude smaller if the mantle is partially molten.
$\tau_{RT}$ is small enough compared to the lifetime of a magma ocean
(${\gtrsim} 1000$ years~\cite{solomatov2015magma}) to allow the
development of the instability before magma ocean crystallisation. 

The non-linear development of the Rayleigh--Taylor instability at the
interface between a thin low viscosity layer and much thicker and more
viscous layer has been investigated experimentally as a model for
mantle plume~\cite{whitehead1975,Bercovici1997}.
Figure~\ref{fig:RT_diapirs}a shows snapshots from one such experiment,
in which a thin layer of glycerin with a density
$\rho_m=1250$~kg${\cdot}$m$^{-3}$ (appearing black on the image) is
placed at the top of a deep layer of a 32 times more viscous silicon
oil of density $\rho_s=920$~kg${\cdot}$m$^{-3}$~\cite{whitehead1975}.
Figure~\ref{fig:RT_diapirs}a shows that the interface between the two
layers develops protrusions of low viscosity fluid that subsequently
separate into several nearly spherical pockets of low-viscosity fluid,
or diapirs~\cite{whitehead1975}. These diapirs can be modelled as
spherical inviscid spheres falling into a viscous liquid.  A reasonable
prediction for the velocity of a metal diapir of radius $R$ migrating
through solid silicate is therefore provided by the Stokes velocity
(for an inviscid sphere)
{\begin{equation}
u_{z} = \frac{\Delta \rho g R^{2}}{3 \eta_{s}},
\label{eq:Stokes_law}
\end{equation}}\unskip
which gives a transit time
{\begin{equation}
\tau_{tr} = \frac{3 \eta_{s} H}{\Delta \rho g R^{2}} \simeq 
\left( \frac{\eta_{s}}{10^{18}~\mathrm{Pa}{\cdot}\mathrm{s}}
\right)\left( \frac{R}{100~\mathrm{km}}\right)^{2}
\left( \frac{H}{1000~\mathrm{km}}\right)\times 250 \ \text{year}.
\end{equation}}\unskip

Of particular interest is the thermal evolution of such a diapir and
the amount of heat it exchanges with the surrounding mantle before
reaching the core. The calculation detailed in Appendix~\ref{AppendixA}
suggests that the amplitude of temperature variations within the diapir
is small compared to the temperature contrast with the surrounding
mantle.  We therefore assume that the temperature within the diapir is
uniform. Accounting for the compression of the diapir, the heat flow
$Q_{m\rightarrow s}$ from the metal diapir to the surrounding, and the
dissipation $\psi_{m}$ in the metal, the energy balance of a diapir
writes
{\begin{equation}
M c_{p,m} \frac{\mathrm{d}T}{\mathrm{d}t} = 
M \frac{\alpha_m T}{\rho_m} \frac{\mathrm{d}P}{\mathrm{d}t} - 
Q_{m\rightarrow s} + \psi_{m}
\end{equation}}\unskip
where $M$ is the diapir mass, $T$ and $P$ its average temperature and
pressure, and $c_{p,m}$ and $\alpha_m$ the specific heat capacity and
thermal expansion coefficient of the metal. This can be rewritten,
using the transformation $\sfrac{\mathrm{d}(\cdots)}{\mathrm{d}t} = u_{z}
\sfrac{\mathrm{d}(\cdots)}{\mathrm{d}z}$, as
{\begin{equation}
 \frac{\mathrm{d}T}{\mathrm{d}z} = \frac{\alpha_m g T}{c_{p,m}} - 
 \frac{Q_{m\rightarrow s} }{M c_{p,m} u_{z}} + \frac{\psi_{m}}{M c_{p_,m} u_{z}}.
 \label{eq:diapir_dTdz}
\end{equation}}\unskip
The first term on the RHS is the isentropic gradient, which gives the
rate of change of the temperature of the diapir if it evolves in an
adiabatic ($Q_{m\rightarrow s}=0$) and reversible ($\psi_{m}=0$) way.

Assuming the diapir keeps a nearly spherical shape, the heat flux from
the diapir to its surroundings is
{\begin{equation}
 Q_{m\rightarrow s} = a 4\uppi R^2 k_s \frac{T-T_s}{R} \mathit{Pe}^{1/2},
 \label{eq:Q_diapir}
\end{equation}}\unskip
with $a \simeq 0.46$ \cite{clift2005bubbles,ulvrova2011compositional},
$T_s$ the temperature of the surrounding mantle, and $k_s$ the thermal
conductivity of the silicate mantle. The P\'eclet number \textit{Pe} is here
defined with the diapir velocity and mantle thermal diffusivity.
Assuming the diapir falls at the Stokes velocity,
{\begin{equation}
 \mathit{Pe} = \frac{R u_z}{\kappa_s} = \frac{\Delta \rho g R^3}{3 \kappa_s \eta_s} 
 \sim \frac{\mathcal{B}}{\kappa_s \nu_s},
 \label{eq:Pe_diapir}
\end{equation}}\unskip
where $\mathcal{B} = \sofrac{4\uppi}{3}R^3\sfrac{\Delta \rho}{\rho_s} g$ is
the total buoyancy of the diapir, and $\nu_s=\eta_s/\rho_s$ and
$\kappa_s$ the kinematic viscosity and thermal diffusivity of the
silicate. As defined here, \textit{Pe} is akin to a Rayleigh number.  Using
Equation \eqref{eq:Q_diapir}, Equation \eqref{eq:diapir_dTdz} can be
recast as
{\begin{equation}
 \frac{\mathrm{d}T}{\mathrm{d}z} + \left(\frac{1}{\ell_\kappa} - 
 \frac{1}{\ell_c} \right) T = \frac{T_s}{\ell_\kappa} + 
 \frac{\psi_{m}}{M c_{p,m}}.
\end{equation}}\unskip
where 
{\begin{eqnarray}
 \ell_c &amp;=&amp; \frac{c_{p,m}}{\alpha_m g},\Seqnsplit
 \ell_\kappa &amp;=&amp; \frac{\mathit{Pe}^{1/2}}{3a} \frac{\rho_m c_{p,m}}
 {\rho_s c_{p,s}} R,
\end{eqnarray}}\unskip
are the characteristic length scales for compression heating and heat
exchange with the surrounding, respectively.

The length scale $\ell_\kappa$ can be interpreted as a \textit{thermal
equilibration distance} with the surrounding mantle. Replacing \textit{Pe}
with its expression \eqref{eq:Pe_diapir}, we find that $\ell_\kappa$
is larger than the solid mantle thickness $H$ when the diapir radius is
larger than
{\begin{equation} 
 R_{\mathrm{eq},1} = H\, 27^{1/5} \left( a \frac{\rho_s c_{p,s}}{\rho_m c_{p,m}} \right)^{2/5} 
 \left(\frac{\Delta \rho g H^3}{\kappa_s \eta_s } \right)^{-1/5} \simeq 
 \left( \frac{\eta_s}{10^{18}~\mathrm{Pa}{\cdot}\mathrm{s}} \right)^{1/5} 
 \left( \frac{H}{10^{3}~\mathrm{km}} \right)^{2/5} \times 10~\mathrm{km}.
\end{equation}}\unskip
Note that $R_{\mathrm{eq},1}$ is weakly dependent on the mantle
viscosity and solid mantle thickness. Only diapirs smaller than about
10~km could have been affected by heat transfer to the mantle, a
prediction which is in agreement with numerical
models~\cite{monteux2009model,samuel2010heat}.

The ratio of the compression and heat exchange terms (1st and 2nd terms
on the RHS of \eqref{eq:diapir_dTdz}) is given by the ratio of
$\ell_\kappa$ and $\ell_c$. It is larger than 1 if the diapir radius is
larger than a critical radius $R_{\mathrm{eq},2}$ given by
{\begin{equation}
 R_{\mathrm{eq},2} = \left(3 a \frac{c_{p,m}}{\alpha_m g} \frac{\rho_s c_{p,s}}{\rho_m c_{p,m}}\right)^{2/5}
 \left(\frac{3 \kappa_s \eta_s}{\Delta \rho g}\right)^{1/5} \sim 
 \left( \frac{\eta_s}{10^{18}~\mathrm{Pa}{\cdot}\mathrm{s}} \right)^{1/5}\times 20 ~\mathrm{km}.
\end{equation}}\unskip
The ratio $\ell_\kappa/\ell_c = (R/R_{\mathrm{eq},2})^{5/2}$ strongly
increases with $R$, which implies that heat exchange quickly becomes
negligible compared to compression heating when $R$ exceeds
$R_{\mathrm{eq},2}$.

The total viscous dissipation $\psi$ associated with the motion of the
diapir is equal to the rate of change of gravitational energy. Of this,
only a small part $\psi_m$ is expected to be dissipated in the diapir,
the remaining part $\psi_s$ being dissipated in the surrounding mantle.
The ratio $\psi_m/\psi_s$ is indeed ${\sim} \eta_{m}/\eta_{s}$ if the
flow within the diapir is laminar (consistent
with~\cite{monteux2009model}), and ${\sim} \sofrac{\rho_m}{\rho_s}
\sofrac{\Delta \rho g R^3}{\eta_s^2}$ if the flow within the diapir is
turbulent~\cite{clesi2023linking}.  Because of the very large viscosity
of the solid silicate, the contribution to the thermal evolution of the
diapir is negligible in both cases, and dissipative heating can be
safely neglected in the diapir energy balance.

Overall, the above simple calculations suggest that: (i) migration by
diapirism of the metal through a solid mantle is fast (it happens on a
timescale which is short compared to the duration of accretion); (ii)
heat exchange between the diapir and surrounding mantle is negligible
if the diapir radius is ${\gtrsim} 10$~km; (iii) dissipation of the
potential energy released during the motion of the diapir is localized
in the mantle, and has a negligible effect on the thermal evolution of
the metal diapir. Together, points (ii) and (iii) suggest that the
diapir temperature evolves along an isentrope.

\subsubsection{Additional complexities}
\label{subsec:diapirs_complications}

Several factors can alter the simple picture we have just described:

\paragraph*{Three-layer systems.} 
The experiments shown in
Figure~\ref{fig:RT_diapirs} have been done in a two-layer system: a
thin dense liquid layer above a thicker, less dense but much more
viscous layer. A more realistic set-up would include three layers with
a dense layer representing the metal pond sandwiched between two layers
representing the (partially) solidified mantle and magma ocean. This
configuration has been investigated using gallium for the metal layer,
and glucose solutions with different viscosities and densities to model
the solid mantle and magma
ocean~\cite{olson2008experiments,fleck2018iron}. Like in the two-layer
system, the interface between the dense metal layer and the underlying
viscous liquid destabilizes through a Rayleigh--Taylor instability. 
Figure~\ref{fig:RT_diapirs}b shows a diapir that developed from this
instability~\cite{fleck2018iron}. The structure differs markedly from
the diapirs obtained in two-layer experiments
(Figure~\ref{fig:RT_diapirs}a). The diapir has a shape close to
hemispherical, with a flat trailing edge, and is followed by a trailing
conduit filled with liquid from the upper layer (liquid s1 on the
figure) corresponding to the magma
ocean~\cite{olson2008experiments,fleck2018iron}. The persistence of the
trailing conduit is allowed by the slow viscous relaxation of the
surrounding viscous liquid (s2).  Compared to a spherical diapir, these
``half diapirs'' are faster. The drag is indeed reduced by the presence
of the low viscosity trailing conduit on the trailing side of the half
diapir~\cite{fleck2018iron}. 

Previous studies~\cite{olson2008experiments,fleck2018iron} have also
investigated the gravitational instability of metal emulsions, under
the hypothesis that the metal droplets resulting from the fragmentation
of impactor cores did not fully coalesce before migrating through the
solid part of the mantle (either because coalescence is slow, or
because the residence time of droplets in the magma ocean is larger
than $\tau_{RT}$).  Figure~\ref{fig:RT_diapirs}c shows snapshots from
an experiment of this type.  The metal migrates as a two-phase plume
consisting of a coalesced half-diapir followed by a trailing conduit
made of metal droplets and low-viscosity liquid. Previous
investigations~\cite{olson2008experiments,king2011heat,fleck2018iron}
argue that dissipation heating in such metal-silicate plume could heat
up both phases and increase the temperature of the metal by ${\sim} 1000$
K.  This contrasts with the reference model of
Section~\ref{subsec:diapir_ref_model} which predicts negligible viscous
heating of the metal. During a single differentiation event, the diapir
transit time is much larger than the sedimentation time of metal
droplets, at least with the nominal parameter values that we used.  In
this limit, the amount of metal droplets in suspension in the trailing
conduit and therefore dissipation heating are likely small.
However,~\cite{olson2008experiments,fleck2018iron} argue that the
conduit can be kept open for a significant time if additional metal is
added before the conduit closure by viscous relaxation. 

\paragraph*{Temperature-dependent viscosity.}
The effective viscosity of the (partially) solidified mantle depends
strongly on temperature and melt fraction. This can have at least two
effects:

(i)  A vertical stratification in viscosity can affect the wavelength
of the most unstable mode of the Rayleigh--Taylor
instability~\cite{stevenson1981models}. The wavelength $\lambda$ can be
smaller if the effective viscosity increases with depth, though likely
not to the degree ($\lambda \sim$ a few~km) proposed in 1981 by
Stevenson~\cite{stevenson1981models}, who considered the
Rayleigh--Taylor instability of a molten metal layer above a cold and
highly viscous non-differentiated layer.

(ii)  If the diapir is hotter than the mantle, then the silicate would
have a reduced viscosity in the thermal boundary layer surrounding the
diapir, which can decrease the drag on the diapir and allow it to fall
more rapidly than predicted by Stokes law \eqref{eq:Stokes_law}. This
has been investigated
theoretically~\cite{marsh1982mechanics,Morris1982,ansari1985effects}
and experimentally~\cite{Ribe1983,ansari1985effects} in the context of
magma intrusions (plutons) in Earth's crust.  It is found that the drag
is a function of \textit{Pe} and of the ratio between the viscosity
$\eta_{s,0}$ of the silicate in contact with the diapir and the
viscosity $\eta_{s,\infty}$ far from it. Qualitatively, two end-members
have been identified: (i) the \textit{Stokes} limit as described in
Section~\ref{subsec:diapir_ref_model}, in which the drag is controlled
by the viscosity far from the diapir; and (ii) the \textit{lubrication}
limit, in which deformation is localized in a thin low viscosity layer
coating the forward side of the diapir~\cite{ansari1985effects}. Using
the prediction of~\cite{ansari1985effects} for the boundary between the
lubrication and Stokes regimes (and equating the drag force with the
buoyancy), we find that a diapir is in the Stokes regime if
{\begin{equation}
 \frac{\Delta\rho g R^3}{\eta_{s,\infty} \kappa_s} \left( \frac{\eta_{s,0}}{\eta_{s,\infty}}\right)^{1/3} 
 \ln \left( \frac{\eta_{s,\infty}}{\eta_{s,0}}\right) \gtrsim 1.
\end{equation}}\unskip
Taking for example $R=10$~km and $\eta_{s,\infty}=10^{18}$~Pa${\cdot}$s, this
would be verified if $\eta_{s,0} \gtrsim 1$~Pa${\cdot}$s. Any realistic
viscosity value for partially molten silicates above the rheological
transition brings the diapir well within the Stokes regime. The model
described in Section~\ref{subsec:diapir_ref_model} then gives a good
approximation of the diapir dynamics. The low sensitivity of the diapir
velocity to the temperature effect on viscosity results from the large
\textit{Pe} value of metal diapirs:  since the boundary layer thickness is
${\sim} R \mathit{Pe}^{-1/2}$, only a small volume around the diapir would have a
reduced viscosity when $\mathit{Pe}\gg 1$.  This could only be compensated if
$\eta_{s,\infty}/\eta_{s,0}\gg 1$. 

Dissipation heating is negligible in laboratory scale experiment, but
could possibly affect the diapir velocity by softening a volume of
silicates extending outside the thermal boundary layer. As discussed in
Section~\ref{subsec:diapir_ref_model}, dissipative heating is expected
to be localized in the adjacent mantle, in agreement with numerical
simulations having $\eta_s/\eta_m \gg
1$~\cite{monteux2009model,samuel2010heat}. In addition to be of
importance for the thermal evolution of the mantle, this can in
principle also affect the velocity of diapirs by softening the
surrounding silicates, thus decreasing the transit time.  Numerical
modelling~\cite{samuel2008dynamics} suggests that the effect on the
diapir velocity is relatively weak in the parameters range investigated
(a 1.7-fold increase for $\eta_{s,\infty}/\eta_{s,0}=10^4$). This can
be qualitatively understood along the same lines as in the case without
dissipation (previous paragraph).  In the Stokes regime, the volumetric
dissipation rate $\tau :\nabla\mathbf{u}$ (where
$\tau$ is the deviatoric stress tensor) evolves with distance $r$ from
the diapir center as ${\sim} r^{-4}$ (since the velocity around the
diapir decreases as $1/r$). Dissipative heating is therefore localized
in the close vicinity of the diapir, and the volume of softened
silicate is moderate, and at least not sufficiently large to allow the
diapir to be in the lubrication regime.

\paragraph*{Diapirism {vs.}\ diking}
The possibility that the metal migrated through the solid mantle as
dikes has been proposed by~\cite{stevenson1981models}. In the 1981's
paper in which this mechanism was proposed~\cite{stevenson1981models},
the inner part of the planet was hypothesized to be cold and
undifferentiated, with an effective viscosity possibly in excess of
$10^{26}$~Pa${\cdot}$s~\cite{karato1997core}. In this limit the migration of
metal diapirs would have been very slow~\cite{stevenson1981models},
while the brittle nature of cold mantle rocks would have allowed the
propagation of iron-filled cracks (analog to magma dikes in the
present-day crust), in a possibly very efficient
way~\cite{stevenson2003mission}. However, whether this mechanism was
operative in a hot and soft mantle near or above the solidus is not
clear.

The competition between diapirism and diking has been a central
question in the context of magma transport in the Earth's crust.
Rubin~\cite{rubin1993dikes} has developed a model of dike propagation
in a viscoelastic rock,  which predicts that the shape of dikes is not
significantly affected by viscous deformation of the host rock if the
viscosity of the melt is smaller than some fraction of the viscosity of
the host rock. This fraction depends on the ratio of excess pressure in
the dike to elastic stiffness of the host rock. Applied to an
iron-filled crack, the model predicts negligible effect of viscous
deformation on the shape of the dike if the mantle viscosity exceeds
${\sim} 10^{11}$--$10^{14}$~Pa${\cdot}$s~\cite{Rubie2015a}. However, the model is
``dike-oriented'': by construction, it is a model of dike propagation in
which the host rock has a viscoelastic rheology, focusing \pagebreak on how
viscosity of the host rock affects the shape of a volume of magma
propagating along a pre-existing crack.\break It does not consider in a
quantitative way the initiation of the crack.  Application of the model
therefore does not provide a definitive answer to the diapir
{vs.}\ dike question. Experiments on buoyancy driven crack
propagation in a viscoelastic material have been performed by Sumita
and Ota~\cite{sumita2011experiments} around the brittle--ductile
transition, using a low concentration agar as the host material and
aqueous solutions for the dike.  They found that near the transition,
the liquid migrates as a hybrid of diapir and dike, with a diapir-like
head and a dike-like tail. Whether these diapir-dike hybrids can be of
relevance for core/mantle differentiation remains to be explored.

\subsection{Percolation}

When the silicate phase is solid, phase separation may also result from
the percolation of molten metal through solid or partially molten
silicate. This is possible if the metal forms an interconnected network
within the crystalline silicate matrix, which depends on the metal
volume fraction and on the capacity of the metal to wet silicate
grains.  Phase separation through percolation may also be relevant to
the separation of liquid metal and silicate at the base of magma oceans
during the sedimentation of metal drops, as suggested
by~\cite{sato2007experiments,olson2008experiments,fleck2018iron}. 

Most of the laboratory work relevant to metal-silicate separation by
percolation has been focused on describing the wetting properties and
mobility of iron-rich melts in mixtures of silicate crystals and iron
alloys at various temperatures and pressures, using an experimental
petrology
approach~\cite{shannon1998percolation,yoshino2003core,terasaki2007interconnectivity,terasaki2008percolative,neri2019metal}. 
As this falls somewhat outside of the scope of this paper, we will not
give much details about this approach (a relatively recent review can
be found in Section 2.3.1 of~\cite{Rubie2015a}), and only list a few
salient points. Iron-rich melts with geologically relevant compositions
do not wet silicates, except at low pressure (${\lesssim} 3$~GPa) at high
O and S content and perhaps in the lower mantle, though this is
debated~\cite{Rubie2015a}. As a result, the percolation of iron-rich
melt is usually only possible if its volume fraction exceeds a
threshold volume fraction (${\sim} 5$~vol\%) above which it is
interconnected. The threshold melt fraction may be smaller 
(${\sim} 1$~vol\%) if the material undergoes shear
deformation~\cite{hustoft2006metal}, though this is also
debated~\cite{Rubie2015a}. Partial melting of the silicates further
complicates the picture, with experiments suggesting the formation of
aggregates of silicate crystals and molten
iron~\cite{neri2021textural}.

In most situations familiar to fluid dynamicists, percolation occurs
through a solid, undeformable porous medium, and is described by
Darcy's law. In the context of core formation, percolation requires in
addition that the porous medium---here the solid silicate matrix---can
deform.  Indeed, conservation of mass implies that any net radial flux
of metal within a planetary body must be accompanied by changes of
silicate volume fraction corresponding to either compaction of the
silicate matrix (where metal is removed) or decompaction (where metal
is added). In other words, the migration of the metal toward the center
is only possible if the silicate matrix deforms to make room for the
metal at the center of the planet. Models of core formation by
percolation must therefore account for the deformation of the matrix.

Consider a two-phase medium made of a silicate matrix with effective
dynamic viscosity $\eta_{s}$ through which a liquid metal of viscosity
$\eta_{m}$ is percolating. We denote by $k_0$ a representative value of
the matrix permeability, and by $\Delta \rho$ the density difference
between the metal and silicate. We assume that melt fraction, and
therefore the melt and matrix velocities, vary spatially over a length
scale $\ell$.  Focusing on the relevant limit of {gravity-driven} phase
separation, one can introduce two velocity scales characteristic of
Darcy flow and matrix deformation (or compaction) (see
also~\cite{monnereau2023differentiation}): the Darcy velocity 
{\begin{equation}
V_D = \frac{\Delta \rho g k_{0}}{\eta_m },
\end{equation}}\unskip
which is independent of $\ell$, and the compaction velocity 
{\begin{equation}
 V_C = \frac{\Delta \rho g \ell^2}{\eta_s}.
\end{equation}}\unskip
The ratio of these two velocity scales is 
{\begin{equation}
\frac{V_D}{V_C} = \frac{\eta_s}{\eta_m}\frac{k_{0}}{\ell^2} = \left(\frac{\delta}{\ell}\right)^2,
 \label{eq:matric_Darcy_timescales_ratio}
\end{equation}}\unskip
where 
{\begin{equation}
\delta = \sqrt{\frac{\eta_s}{\eta_m} k_{0}}
\end{equation}}\unskip
is the so-called \textit{compaction length}. This name can be somewhat
misleading and care must be taken when interpreting the value of
$\delta$.  As shown by Equation
\eqref{eq:matric_Darcy_timescales_ratio}, $\delta$ gives the length
scale  at which Darcy friction and viscous stress associated with
matrix deformation affect phase separation to a similar extent. Phase
separation is rate-limited by Darcy friction if $\ell\gg \delta$, and
by matrix deformation if $\ell\ll \delta$. Both limits are relevant to
the core formation problem, depending on the size and state of the
planetary body~\cite{monnereau2023differentiation}.

\begin{figure}
\includegraphics{fig09}
\caption{\label{fig:SatoSumita2003}Snapshot from an experiment in which a mixture of silicon oil
with viscosity $12.1$~Pa${\cdot}$s (appearing yellow on the pictures) and
hydroxyethylcellulose (HEC) solution with viscosity $0.023$~Pa${\cdot}$s
(appearing blue) separate under the action of gravity. The density
difference between the HEC solution and the silicon oil is
$31$~kg${\cdot}$m$^{-3}$.  On the left is shown a vertical section
(with a height of $12.8$ cm) of the fluid mixture. Subfigures (a), (b),
and (c) are close-ups corresponding to the boxes shown on the left
figure. The scale bar in (a) is 1 cm. Details on the experiment can be
found in~\cite{sato2007experiments} Reproduced from Figure~9
of~\cite{sato2007experiments}, 2007 \copyright~Cambridge University Press,
reproduced with permission (Licence number 5797030122437).}
\end{figure}

The rich dynamics coming from the interplay between matrix deformation
and Darcy flow is difficult to study experimentally, and has so far
been almost exclusively studied from theoretical and numerical point of
views~\cite{vsramek2010multiphase,vsramek2012thermal,lichtenberg2019magma,monnereau2023differentiation}.
To our knowledge the only fluid dynamics experiments devoted to
gravity-driven percolation in the context of core formation are the
experiments of Sato and Sumita~\cite{sato2007experiments}, in which two
immiscible liquids separate under the action of gravity after having
been dispersed by mechanical stirring.  Their purpose was to model the
separation of molten silicate and metal following the sedimentation of
the metal at the base of a magma ocean. They therefore focused on
moderate viscosity ratios $\eta_s/\eta_m$ (smaller than ${\simeq} 600$),
using silicon oil and hydroxyethylcellulose (HEC) solutions as working
fluids. Figure~\ref{fig:SatoSumita2003} shows a snapshot of an
experiment using a silicon oil with viscosity $12.1$~Pa${\cdot}$s (appearing
yellow on the pictures) and a HEC solution with viscosity $0.023$~Pa${\cdot}$s
(appearing blue). Mixtures of silicon oil and HEC appear green. 
Interestingly, it was found in these experiments that the most viscous
phase forms drops, while the less viscous forms a continuous phase
between the drops (Figure~\ref{fig:SatoSumita2003}).  Applied to
metal-silicate separation at the base of a magma ocean, this would
predict that the silicate phase would form drops surrounded by the
metal, which is opposite to the standard view of an iron rain. It
contrasts with the gallium-glucose mixtures of the experiments
discussed in Section~\ref{subsec:diapirs_complications} and shown in
Figure~\ref{fig:RT_diapirs} in which the less viscous phase (gallium)
was observed to form drops~\cite{olson2008experiments,fleck2018iron}.

These experiments can be used to illustrate the difficulty of designing
fluid dynamics experiments able to reach regimes in which matrix
deformation is of importance for the separation dynamics. The
compaction length in these experiments can be estimated from the
viscosity values of the two liquids as given
by~\cite{sato2007experiments} and the permeability law they used to
model the experiment, namely the Blake--Kozeny--Carman permeability
law 
{\begin{equation}
 k = \frac{a^2}{K_0} \frac{\phi^3}{(1-\phi)^2},
 \label{eq:permeability}
\end{equation}}\unskip
where $a$ is the droplet size and $K_0\simeq 15$ is an empirical
constant. With $a\sim 0.1$ mm as estimated in the paper and a maximum
viscosity ratio of ${\simeq} 600$, one finds that $\delta$ is at most
${\simeq} 0.5$ mm, which is not much larger than the droplets size. 
Given the thickness of the tank ($H=14$ cm $\gg \delta$), phase
separation is therefore expected to be controlled by Darcy flow, which
is consistent with the measured separation times which are of the order
of ${\sim} H/V_D$ in this set of experiments~\cite{sato2007experiments}.
Increasing the importance of matrix deformation necessitates to
increase the ratio of the compaction length to tank thickness, which
can be done by either decreasing the tank thickness or 
increasing~$\delta$. To make sure the mixture can be modelled as a continuum, the
tank size has to be much larger than the drop size, which suggests a
tank size of at least a few centimeters as in the experiments of Sato
\&amp; Sumita. To increase the value of $\delta$, one could use a pair of
liquids with a higher viscosity ratio. One may for example use a
silicon oil with viscosity $10^3$~Pa${\cdot}$s and water ($10^{-3}$~Pa${\cdot}$s) to
reach a viscosity ratio of ${\sim} 10 ^6$. With a drop size $a \sim 0.1$
mm and assuming the permeability law \eqref{eq:permeability}, this
would give $\delta \sim 2$ cm. It is therefore in principle possible to
design a fluid dynamics experimental set-up that could allow to study
gravity-driven separation in a regime in which the matrix deformation
velocity and Darcy velocity are of similar magnitudes. However, a
number of practical difficulties are to be expected. Among others, one
could note that producing emulsions with small drop size may be an
issue with highly viscous liquids, or that longer separation times may
increase the importance of drops coalescence, thus complicating the
interpretation of the experiments.

\section{Conclusion} \label{sec:conclusion}


The experimental fluid mechanics approach has proved relevant to the
study of the differentiation of terrestrial planets, complementing the
more established numerical modelling approach to the subject. Its main
interest lies in its ability to incorporate physical effects and
processes that are difficult to account for in numerical simulations,
such as turbulence, surface tension and fragmentation, heat and mass
transfer at large P\'eclet numbers, or large viscosity contrasts.
Perhaps unexpectedly, experiments using common working fluids at
ambient conditions can achieve dynamical conditions that are closer to
early Earth conditions than current numerical simulations in terms of
several (but not all) of the key governing non-dimensional parameters.
In particular, laboratory experiments have allowed to model
metal-silicate separation and fragmentation in turbulent conditions, at
\textit{Re} and \textit{We} values around $10^{4}$ that are still difficult to
reach with numerical simulations.


However, the experimental approach also has important limitations. A
number of important effects cannot be investigated easily in fluid
dynamics experiments.  The effects of compressibility (including shock
waves propagation), phase change (melting and vaporisation), or
self-gravitation on planetary impacts and differentiation are best
studied using numerical models.  Even when ignoring these effects, the
results of the experiments still must be used with caution, because of
the wild extrapolations required to apply them to the early Earth.  For
example, the metal-silicate equilibration model built from experiments
such as described in Section~\ref{sec:magmaocean} has been incorporated
in several geochemical models of core-mantle formations with some
success~\cite{Rubie2015,Badro2016,Jacobson2017}. Yet, the model is
incomplete and needs to be refined and further tested.  The mechanisms
leading to fragmentation and homogenisation in turbulent thermals or
during crater opening are not fully elucidated yet, and additional
experiments and analysis are needed. A better understanding of
fragmentation and mixing in turbulent thermals is essential to give
more weight to the extrapolation of experimental results.

\thanks{The authors thank M. Le Bars and S. Fauve for offering them the
opportunity to write this review article, and two anonymous reviewers
for their very detailed and constructive comments. RD thanks the French
Academy of Science and Electricit\'e de France for granting their
``Amp\`ere Prize'' to the Geodynamo team at ISTerre.  This work was
supported by the European Research Council (ERC) under the European
Unions Horizon 2020 research and innovation programme (grant number
716429). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). This work
was supported by the Programme National de Plan\'etologie (PNP) of
CNRS-INSU co-funded by CNES. Funding for this research was supported by
the Center for Matter at Atomic Pressure (CMAP), a National Science
Foundation (NSF) Physics Frontier Center, under award PHY-2020249. Any
opinions, findings, conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect
those of the National Science Foundation.}

\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}
The authors thank M. Le Bars and S. Fauve for offering them the
opportunity to write this review article, and two anonymous reviewers
for their very detailed and constructive comments. RD thanks the French
Academy of Science and Electricit\'e de France for granting their
``Amp\`ere Prize'' to the Geodynamo team at ISTerre.  This work was
supported by the European Research Council (ERC) under the European
Unions Horizon 2020 research and innovation programme (grant number
716429). ISTerre is part of Labex OSUG@2020 (ANR10 LABX56). This work
was supported by the Programme National de Plan\'etologie (PNP) of
CNRS-INSU co-funded by CNES. Funding for this research was supported by
the Center for Matter at Atomic Pressure (CMAP), a National Science
Foundation (NSF) Physics Frontier Center, under award PHY-2020249. Any
opinions, findings, conclusions or recommendations expressed in this
material are those of the author(s) and do not necessarily reflect
those of the National Science Foundation.

\back{}

\appendix{}

\def\appendixlabel{Appendix \Alph{section}}

\section{Temperature inside a metal diapir}\label{AppendixA}

The model of diapir thermal evolution presented in Section
\ref{subsec:diapir_ref_model} assumes that the temperature in the
diapir is uniform.  This assumes that convection within the diapir is
strong enough to ensure that temperature variations within the diapir
are small compared to the temperature contrast between the diapir and
the surrounding mantle. The magnitude of temperature variations within
a diapir can be estimated by equating the heat flux associated with
convection within the diapir with the mantle-side heat flux controlled
by the slow viscous flow around the diapir.  We denote by $\Delta T$
the temperature difference between the centre of the diapir and the
surrounding mantle, and by $\Delta T_d$ the difference between the
temperature at the centre of the diapir and the mean boundary
temperature.  The heat flux from the diapir to the mantle, when
estimated from the diapir side, can be related to $\Delta T_d$ using a
scaling law derived from experimental and numerical studies of natural
convection in a sphere cooled from its boundary in a uniform
(non-radial) gravity field~\cite{schmidt1956,chow1975,hutchins1989pseudosteady}.
According to these studies, the heat flux
associated with internal convection is given by 
{\begin{equation}
q =  a_d k_m \frac{\Delta T_d}{R} Ra_d^\beta =  
a_d k_m \frac{\Delta T}{R} \left(\frac{\Delta T_d}{\Delta T}
\right)^{1+\beta} Ra_{\Delta T}^\beta 
\end{equation}}\unskip
where $ a_d=0.553$ and $\beta=0.2521$~\cite{hutchins1989pseudosteady}, and
{\begin{equation}
Ra_d = \frac{\alpha_m g R^3 \Delta T_d}{\nu_m \kappa_m}, 
\quad Ra_{\Delta T} = \frac{\alpha_m g R^3 \Delta T}{\nu_m \kappa_m}.
\end{equation}}\unskip
When estimated on the mantle side, the heat flux is given by
{\begin{equation}
    q = a k_s \frac{\Delta T - \Delta T_d}{R} Pe^{1/2}.
\end{equation}}\unskip
Equating the two expressions for $q$ and re-arranging, we find that 
{\begin{equation}
\frac{(\Delta T_d/\Delta T)^{1+\beta}}{1-\Delta T_d/\Delta T} = 
\frac{a_o}{ a_d} \frac{k_s}{k_m} Pe^{1/2} Ra_{\Delta T}^{-\beta} 
\propto R^{0.74} \Delta T^{-0.25}
\end{equation}}\unskip
Taking $R=100$~km and $\Delta T=100$ K gives $\Delta T_d/\Delta T\simeq
1.3$\%. Taking $R=500$~km instead gives $\Delta T_d/\Delta T\simeq
3.5$\%. The ratio is smaller when $\Delta T$ is greater.  This shows
that the amplitude of temperature variations within the diapir is small
compared to the temperature difference with the surrounding mantle,
justifying the assumption of a uniform temperature within the diapir.

\bibliographystyle{crunsrt}
\bibliography{crphys20240722}
\refinput{crphys20240722-reference.tex}

\end{document}
</pre></body></html>