# Fluid dynamics of planetary differentiation

Renaud Deguen, Ludovic Huguet, Maylis Landeau, Victor Lherm, Augustin Maller, Jean-Baptiste Wacheul

► To cite this version:

Renaud Deguen, Ludovic Huguet, Maylis Landeau, Victor Lherm, Augustin Maller, et al. Fluid dynamics of planetary differentiation. Comptes Rendus. Physique, 2025, 25 (S3), pp.1-45. 10.5802/cr-phys.227. hal-04990367

HAL Id: hal-04990367

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

Submitted on 13 Mar 2025

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

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

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

Distributed under a Creative Commons Attribution 4.0 International License
Comptes Rendus
Physique

Renaud Deguen, Ludovic Huguet, Maylis Landeau, Victor Lherm, Augustin Maller and Jean-Baptiste Wacheul

**Fluid dynamics of planetary differentiation**

Published online: 13 January 2025

**Part of Special Issue:** Geophysical and astrophysical fluid dynamics in the laboratory

**Guest editors:** Stephan Fauve (Laboratoire de Physique de l'ENS, CNRS, PSL Research University, Sorbonne Université, Université Paris Cité, Paris, France) and Michael Le Bars (CNRS, Aix Marseille Univ, Centrale Marseille, IRPHE, Marseille, France)

https://doi.org/10.5802/crphys.227

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

This article is licensed under the
CREATIVE COMMONS ATTRIBUTION 4.0 INTERNATIONAL LICENSE.

http://creativecommons.org/licenses/by/4.0/
Review article / *Article de synthèse*

Geophysical and astrophysical fluid dynamics in the laboratory /
*Dynamique des fluides géophysiques et astrophysiques au*
*laboratoire*

Fluid dynamics of planetary differentiation

*La différenciation des planètes telluriques vue par la mécanique des fluides*

Renaud Deguen<sup>®,*,a</sup>, Ludovic Huguet<sup>®,a</sup>, Maylis Landeau<sup>®,b</sup>, Victor Lherm<sup>®,c</sup>, Augustin Maller<sup>®,b</sup> and Jean-Baptiste Wacheul<sup>d</sup>

**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.

**Résumé.** La structure des planètes terrestres — un noyau métallique riche en fer entouré d'un manteau silicaté — a été établie pendant leur accrétion, lorsque la fusion généralisée a permis la séparation des phases métalliques et silicatées. Le transfert d'éléments chimiques et de chaleur entre métal et silicates qui s'est

<footnote><sup>a</sup> ISTERRE, Université Grenoble Alpes, Université Savoie Mont Blanc, CNRS, IRD, Université Gustave Eiffel, 38000 Grenoble, France</footnote>

<footnote><sup>b</sup> Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, 75005 Paris, France</footnote>

<footnote><sup>c</sup> Université Paris-Saclay, CNRS, FAST, 91405, Orsay, France</footnote>

<footnote><sup>d</sup> Univ Lyon, ENSL, UCBL, UJM, CNRS, LGL-TPE, F-69007 Lyon, France</footnote>

<footnote>*E-mails*: renaud.deguen@univ-grenoble-alpes.fr (R. Deguen), ludovic.huguet@univ-grenoble-alpes.fr (L. Huguet), landeau@ipgp.fr (M. Landeau), victor.lherm@cnrs.fr (V. Lherm), maller@ipgp.fr (A. Maller), jb.wacheul@orange.fr (J.-B. Wacheul)</footnote>

<footnote>* Corresponding author</footnote>
produit pendant cette période est déterminant pour la composition et la température initiale du noyau et du manteau, et a des implications importantes pour l'évolution et la dynamique à long terme des planètes. Après avoir résumé les principales contraintes observationnelles sur la différenciation noyau-manteau, l'article suit la séquence des processus qui ont conduit à la formation des noyaux planétaires, en se concentrant sur les contributions de l'approche de dynamique des fluides expérimentale à notre compréhension 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étaires, en utilisant des expériences de laboratoire pour illustrer et quantifier les processus d'impact et de cratérisation, et la dispersion de la phase métallique qui en résulte. Nous examinons ensuite l'écoulement diphasique qui suit un impact, lorsque le noyau d'un impacteur chute dans un océan de magma. Nous introduisons le modèle de thermique turbulent, qui constitue selon nous un bon modèle de référence pour l'écoulement post-impact. Nous discutons ensuite des facteurs additionnels — immiscibilité et fragmentation, inertie héritée de l'impact, force de Coriolis, sédimentation — pouvant affecter les prédictions de ce modèle. Enfin, la dernière partie de l'article est consacrée à la migration de la phase métallique à travers une partie solide du manteau.

**Keywords.** Early earth, Core formation, Fluid dynamics experiments.

**Mots-clés.** Terre primitive, formation du noyau, dynamique des fluides expérimentale.

**Funding.** French Academy of Science, Electricité de France, European Research Council (ERC) (grant number 716429), Labex OSUG@2020 (ANR10 LABX56), Programme National de Planétologie (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.

Manuscript received 7 September 2024, revised and accepted 17 December 2024.

# 1. Introduction

## 1.1. 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 [1]. Following the accretion of dust from the solar nebula into planetesimals with diameters ranging from several tenths of km to a few hundred km [2,3], 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 [3]. 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 [4].

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 dis-equilibrium 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 [5]. 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 [6–8], 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 $fO_2$ [5, 9, 10]. This chemical fractionation also has major geodynamical implications. In particular, the nature and abundance in the core of radioactive elements [11–13], magnesium [14, 15], or light elements (e.g. O, S, Si) depend on the pressure *P*, temperature *T*, and redox conditions under which metal and silicates equilibrated [16–18]. 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 [19, 20]. 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 [21], the existence of an early dynamo [22, 23], 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 *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 1a 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 1, the accretion time estimated from Hf/W is about 10 My if perfect equilibration is assumed [6–8], but can be much larger if equilibration is imperfect [8, 24, 25]. 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 ≈40% [8, 24–26]. This implies a significant flux of chemical elements between metal and silicate during differentiation.

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 [28]. 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 *excess siderophile problem* [29]. This can be seen in Figure 1b which compares the apparent core/mantle concentration ratio $D_{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* > 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* < 1, e.g. Mn and V). The most plausible explanation for this discrepancy is chemical equilibration of the metal and silicate at
<figure><img src="image_3.png" /><figcaption>**Figure 1.** 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 [8]). (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 [27] and partitioning parameterization of [10] for a range of oxygen fugacity (IW–1.5 to IW–3.5).</figcaption></figure>

conditions of high temperature and high pressure (>2000 K and 30–60 GPa), under which the partitioning of these elements is more consistent with the silicate Earth composition [30–32]. 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 [5, 8, 33–36].

Geochemical (as summarised above) and physical arguments have led to the scenario of planetary cores formation sketched in Figure 2. Radioactive and impact heating led to widespread melting of the growing terrestrial planets, providing favourable conditions for metal-silicate separation [38–40]. 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 [41]. The metal delivered by each impact is thought to have fragmented into droplets that eventually settled at the base of the magma ocean [42–44]. The estimated size of these droplets (mm to cm) is small enough to have allowed them to equilibrate efficiently with the surrounding magma ocean [42–44]. 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 [38]. 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 [39, 43, 45, 46]. Diking (fracturation) or percolation are also possible mechanisms [38]. 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
<figure><img src="image_4.png" /><figcaption>**Figure 2.** Sketch illustrating the various metal-silicate segregation processes that can occur during the differentiation of terrestrial planets. From [37], adapted from [4].</figcaption></figure>

these conditions [30–32,47]. 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 [48, 49], and it is far from obvious that volumes of molten iron of size ~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 [8, 24–26]. The equilibration efficiency is constrained by Hf–W systematics (Figure 1a), but the range of permissible values allows for a number of different accretion and differentiation scenarios [8]. 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.

## 1.2. 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 fragmentation of the metal phase, as well as stirring and mixing processes. Phase separation, fragmentation, and mixing are of importance in many industrial and natural contexts, and have received significant interest from various communities [50,51]. 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 $ρ_m$, $ρ_s$ and $v_m$, $v_s$, respectively (where the subscripts m and s denote metal and silicates). We denote by $Δρ = ρ_m - ρ_s$ the density contrast. The fluids being immiscible, their boundary has a surface tension $σ$. The time and length scales at stake can be such that the rotation of the planet may affect the flow; the rotation rate $Ω$ 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 $κ$ needs to be included in the set of relevant parameters ($κ$ 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

$$Re = \frac{UL}{v_s}, \quad We = \frac{\rho_s U^2 L}{\sigma}, \quad Fr = \frac{U^2}{gL}, \quad Ro = \frac{U}{\Omega L}, \quad Pe = \frac{UL}{\kappa}, \quad M = \frac{U}{c}, \quad \frac{\Delta\rho}{\rho_s}, \quad \frac{v_m}{v_s}. \quad (1)$$

The Reynolds number Re, Weber number We, Froude number Fr, and Rossby number Ro compare inertia to viscous forces, surface tension, gravity forces, and Coriolis force, respectively. The Péclet number 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 1):

**(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 Re, We, and 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 \text{ km} \cdot \text{s}^{-1}$ has Re ~ We ~ $10^{14}$ and Pe ~ $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 ~1 for large impacts and $\gg 1$ for small-sized impactors (see discussion in Section 2.1). 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 2.3.1). 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 3.4.2).

The huge values of Re and We imply highly turbulent conditions, with viscosity and surface tension playing a role only at the smallest scales [65, 66]. Post-impact flows must have involved a wide range of length scales, from the size of the impactor to less than ~mm scale at which
<table><caption>**Table 1.** Representative values of relevant physical parameters</caption><tbody><tr><td>Impactor radius</td><td></td><td>R<sub>i</sub></td><td>50–3000 km</td><td>[2,3]</td></tr><tr><td>Impact velocity</td><td></td><td>v<sub>i</sub></td><td>5–40 km·s<sup>-1</sup></td><td>[52]</td></tr><tr><td>Earth rotation rate</td><td></td><td>Ω</td><td>2–5 present-day rate</td><td></td></tr><tr><td>Earth's gravity</td><td></td><td>g</td><td>≤ 10 m·s<sup>-2</sup></td><td></td></tr><tr><td rowspan="3">Density</td><td>Metal</td><td>ρ<sub>m</sub></td><td>5000–10 000 kg·m<sup>-3</sup></td><td>[53]</td></tr><tr><td>Silicate</td><td>ρ<sub>s</sub></td><td>3000–5000 kg·m<sup>-3</sup></td><td>[54]</td></tr><tr><td>Molten metal</td><td>η<sub>m</sub></td><td>10<sup>-3</sup>–10<sup>-2</sup> Pa·s</td><td>[55, 56]</td></tr><tr><td rowspan="2">Dynamic viscosity</td><td>Molten silicate</td><td>η<sub>s</sub></td><td>3 × 10<sup>-3</sup>–10<sup>-1</sup> Pa·s</td><td>[57, 58]</td></tr><tr><td>Solid silicate</td><td>η<sub>s</sub></td><td>10<sup>12</sup>–10<sup>21</sup> Pa·s</td><td>[59, 60]</td></tr><tr><td rowspan="2">Thermal conductivity</td><td>Metal</td><td>k<sub>m</sub></td><td>30–100 W·m<sup>-1</sup>·K<sup>-1</sup></td><td>[61]</td></tr><tr><td>Silicate</td><td>k<sub>s</sub></td><td>1–2 W·m<sup>-1</sup>·K<sup>-1</sup></td><td>[62]</td></tr><tr><td rowspan="2">Compositional diffusivity</td><td>Metal</td><td>κ<sub>c,m</sub></td><td>10<sup>-12</sup>–10<sup>-8</sup> m<sup>2</sup>·s<sup>-1</sup></td><td>[44]</td></tr><tr><td>Silicate</td><td>κ<sub>c,s</sub></td><td>10<sup>-12</sup>–10<sup>-8</sup> m<sup>2</sup>·s<sup>-1</sup></td><td>[44, 62]</td></tr><tr><td>Interfacial tension</td><td>Metal-silicate</td><td>σ</td><td>0.1–1 J·m<sup>-2</sup></td><td>[63]</td></tr><tr><td rowspan="2">Heat capacity</td><td>Metal</td><td>c<sub>p,m</sub></td><td>500–1000 J·kg<sup>-1</sup>·K<sup>-1</sup></td><td>[64]</td></tr><tr><td>Silicate</td><td>c<sub>p,s</sub></td><td>500–1000 J·kg<sup>-1</sup>·K<sup>-1</sup></td><td>[64]</td></tr><tr><td>Sound velocity</td><td>Silicates</td><td>c</td><td>3–10 km·s<sup>-1</sup></td><td>[54]</td></tr></tbody></table>

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.

surface tension and viscosity eventually come into play. Numerical simulations of planetary-scale impacts [67–70] 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 (~14 km in [71]) 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 [72]). Direct numerical simulations of the fragmentation of the metal phase in a magma ocean [73, 74] 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.

**(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 4.1.2. We will focus our discussion on the diapirism and percolation scenarios. In these cases the Reynolds number *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 (*Fr* ≪ 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 *We* and *Ro*, it is therefore more relevant to use the Bond number *Bo* and the Ekman number *Ek*
defined as

$$Bo = \frac{\Delta \rho g L^2}{\sigma} = \frac{\Delta \rho We}{\rho_s Fr}, \quad Ek = \frac{\nu_s}{\Omega L^2} = \frac{Ro}{Re} \quad (2)$$

The Bond number is large in the case of metal diapirs, but possibly as small as ~1 in the case of percolation. The Coriolis force is negligible compared to the viscous forces (Ek ≫ 1) in all of these cases. Typical velocities are much smaller than the speed of sound (M ≪ 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 4.1.1). The Péclet number Pe is large in the case of diapirism because of the large length scale involved (Section 4.1.1), 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).

## 1.3. 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 1.1, this is of primary importance for the interpretation of geochemical data in terms of core formation timing and conditions. The equilibration efficiency $\mathcal{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 [75]. Assuming that the mass of metal fully equilibrates with a mass of silicate $\Delta$ times larger, one finds that

$$\mathcal{E}_i = \frac{1}{1 + D_i / \Delta}, \quad (3)$$

where $D_i$ is the metal/silicate partitioning coefficient of element i. Equation (3) highlights the importance of the parameter $\Delta$, which has been termed *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 [75–78]. In the following sections we will refer to this type of measurement as *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 $\mathcal{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 3.2.3 and 3.5).

## 1.4. Content

The structure of the article follows the sequence of processes that led to the formation of planetary cores. Section 2 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 3 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 4 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 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.

## 2. Planetary impacts

### 2.1. Rationale

Planetary impacts have been classified as being either in the so-called strength regime or gravity regime, depending on the importance of the resistance to deformation of the target material [79]. 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 [79, 80]. Given the expected range of sizes of planetesimals (≥100 km in diameter) and embryos (≥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 [79, 80]. In the case of an impactor striking a pre-existing magma ocean, the impact-induced stress is much larger than viscous stress i.e. Re ≫ 1. The resistance to deformation is thus usually of negligible importance in this context. During an impact, We, Pe, and 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 ~10.

In this limit, the dominant forces are inertia, pressure forces, and gravity, and the non-dimensional numbers of most relevance are Fr, M, and the density ratio, with Fr and M defined in Equation (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

$$Fr^* = \frac{\rho_i U_i^2}{\rho_t g R_i}, \qquad (4)$$

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, $Fr^*$ is the ratio of the kinetic energy of the impacting body (~$\rho_i U_i^2 R_i^3$) to the gravitational potential energy (~$\rho_t g R_i^4$) released by displacing a volume ~$R_i^3$ of target material over a height ~$R_i$. Assuming that the impact velocity scales as the escape velocity, $U_i \approx \sqrt{gD_p}$, where $D_p$ is the diameter of the growing planet and g the acceleration of gravity at its surface, we find that

$$Fr^* \approx \frac{\rho_i D_p}{\rho_t R_i}. \qquad (5)$$

Since $\rho_i/\rho_t$ varies within a limited range, the Froude number $Fr^*$ depends primarily on the impactor-to-target size ratio. $Fr^*$ ranges from ~1 in the case of a giant impact to several hundreds in the case of ~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·s⁻¹, which gives a Froude number ~2 × 10⁵.
The full range of $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 planetary context [81, 82]. Both approaches have their advantages and limitations: the possibility of reaching supersonic impact speeds comes at the expense of the value of $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 $Re$ (to be in the “gravity” regime) and high $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_{air})/\rho_{air}gR_i \sim \mathcal{B}^{1/2}/R_i}$ of a drop falling into air, where $\rho_{air}$ is the density of air and $\mathcal{B} = (\rho_i - \rho_{air})/\rho_{air}g(4\pi/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

$$Re = \frac{\mathcal{B}^{1/2}}{\nu} \frac{U_i}{U_i^t} \lesssim \frac{\mathcal{B}^{1/2}}{\nu}, \qquad (6)$$

$$We = \frac{\rho_i - \rho_{air}}{\rho_{air}} Bo \left( \frac{U_i}{U_i^t} \right)^2 \lesssim \frac{\rho_i - \rho_{air}}{\rho_{air}} Bo, \qquad (7)$$

$$Fr = \frac{\rho_i - \rho_{air}}{\rho_{air}} \left( \frac{U_i}{U_i^t} \right)^2 \lesssim \frac{\rho_i - \rho_{air}}{\rho_{air}}, \qquad (8)$$

where the Bond number is defined as

$$Bo = \frac{\rho_i g R_i^2}{\sigma} \qquad (9)$$

(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 $Fr$ [77, 83, 84]. However, reaching high values of $Re$ and $We$ at moderate values of $Fr$ requires larger impactors. We can indeed notice that $We = BoFr$, which implies that reaching a high value of $We$ while maintaining a moderate value of $Fr$ requires a high value of $Bo$, and therefore a larger impactor, typically of a few cm in diameter at least. Drops with a diameter up to ~5 mm can be easily formed by dripping liquid from the tip of a syringe or capillar [77, 84, 85]. 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-10^3$] [77, 84]. Producing liquid impactors of a few centimeters with moderate deformation is more complicated, but necessary to explore smaller values of $Fr$ while still keeping $Re, We \gg 1$ [76, 86]. One method consists in enclosing the impactor liquid in a balloon which is pierced just before impact [76, 86–88].

## 2.2. Phenomenology and scaling laws

The impact of a drop on a deep pool has been studied extensively [83, 84, 89–97] since the pioneering work of Worthington at the turn of the 19th and 20th centuries [89–91]. 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 [83]. 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 [76, 77, 98] and of immiscibility [85, 86, 99, 100] have been the subject of fewer studies, and not always in regimes relevant to planetary differentiation. Experiments on oblique impacts are scarce as well [101].

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

We will discuss the three main stages of the impact and cratering process that can be identified from these experiments at large $Re$ and $We$: (i) the opening and growth of the crater (Section 2.2.1), 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 2.2.2); and (iii) the collapse of the central peak (Section 2.2.3).

### 2.2.1. *The opening phase: Rayleigh-Taylor instability and fragmentation*

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 3a and 3b). 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 3a and 3b) which expands concomitantly with the crater.

**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 [77, 83, 92, 94]. Neglecting viscous dissipation (the flow around the crater is close to potential [84]), 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

$$ \xi\pi\rho_t R_c^3 \dot{R}_c^2 + \pi\sigma R_c^2 + \frac{\pi}{4}\rho_t g R_c^4 = \phi \frac{2\pi}{3}\rho_i R_i^3 U_i^2. \quad (10) $$

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 [77, 94]. PIV velocity measurements [84] have shown that this is actually a poor approximation (as shown in Figure 3c): 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 (10) to experimental data requires $\xi \approx 0.34$ and $\phi \approx 0.4$ [77]. Surface tension does not affect much the shape of the crater during the opening phase if $We$ is large enough ($We \gtrsim 54Fr^{0.25}$ according to [94]), 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
<figure><img src="image_5.png" /><figcaption>(b) Impact and crater opening</figcaption></figure>

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

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

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

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

<figure><img src="image_10.png" /><figcaption>Crater collapse and jet rise</figcaption></figure>

<caption>Jet collapse</caption>

<figure><img src="image_11.png" /><figcaption>(d)</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>(e)</figcaption></figure>

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

**Figure 3.** (a) Time sequence from a liquid impact experiment at *Fr* = 540 and *Re* = 8.3 × 10³ with miscible drop and target liquids (using a CsCl aqueous solution with *ρ*ᵢ = 1789 kg·m⁻³ for the drop, water as the target liquid, *R*ᵢ = 2.36 mm, and *U*ᵢ = 3.54 m·s⁻¹). The CsCl solution appears blue on the pictures. (b) Time sequence from a liquid impact experiment at *Fr* = 94, *We* = 1.4 × 10⁴, and *Re* = 8.6 × 10⁴ with immiscible drop and target liquids (using Perfluorodecalin with *ρ*ᵢ = 1940 kg·m⁻³ for the drop, water as the target liquid, *R*ᵢ = 2 cm, and *U*ᵢ = 4.3 m·s⁻¹) [86]. (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) (*Fr* = 572, *ρ*ᵢ/*ρ*ᵣ = 1.754); the velocity measurements have been performed in an experiment involving a water drop (to avoid strong variations in optical index) and has *ρ*ᵢ/*ρ*ᵣ = 1 and *Fr* = 446. (e) Close-up on image (b)iv.
size can easily be obtained by setting $\dot{R}_c = 0$ in Equation (10). Neglecting the surface tension term for simplicity, this yields,

$$R_{c,\max} = R_i \left(\frac{8}{3}\phi\right)^{1/4} Fr^{*1/4}, \qquad (11)$$

which is in good agreement with drop impact experiments [77]. 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 2.3.1).

The time $\tau_{\max}$ at which the crater reaches its maximum size can be predicted by integrating in time equation (10). This predicts

$$\tau_{\max} \sim \frac{R_i}{U_i} \left(\frac{\rho_i}{\rho_t}\right)^{-1/2} Fr^{*5/8}, \qquad (12)$$

which again is consistent with experimental observations (see [77] 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 \text{ km} \cdot \text{s}^{-1}$, Equation (12) predicts a crater opening time of ~3 min.

**Instabilities.** Inspection of Figures 3a,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 3a and b, and 3a). 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 [77].

This is interpreted as being the result of a Rayleigh–Taylor instability triggered by the deceleration of the crater boundary [76, 77, 98]. In the opening phase, the crater radius $R_c$ is an increasing function of time and $\dot{R}_c = dR_c/dt$ 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 = d^2R_c/dt^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, $\dot{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 (10), taking the derivative with respect to time and re-arranging gives $\ddot{R}_c/g \sim Fr^*(R_i/R_c)^4$. When $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 $\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_{\max}$. Lherm et al. [77] found that the volume of the mixing layer, normalized by the impactor's volume, is $\approx 0.07(\rho_i/\rho_t)^{2.8} Fr^{0.67}$.

**Fragmentation.** Maller et al. [86] 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 3b, and Figure 3e). 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 $Fr$ exceeds 40, regardless of the value of the Weber number $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 $We > 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 $We$.

During the formation of planets, most impactors are differentiated into a metallic core and a silicate mantle. Maller et al. [86] 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 $Fr > 40$.

Applied to planetary impacts, the above results predict that the larger the impact speed, the larger the Froude number $Fr$, and hence the more fragmented the impactor core. As discussed in Section 2.1, the Froude number of a planetary impact scales as the ratio of the planet to impactor diameters (Equation (5)). 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. [86] predict that the cores of impactors with a radius $R_i < 300$ km fully fragment into drops during the impact process, while the cores of giant impactors with $R_i > 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.

### 2.2.2. Crater closing and formation of a central peak

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 [83] and quantitative [84] velocity measurements show that the flow is convergent toward the axis of symmetry of the crater (Figure 3d). This convergence leads to the formation of an upward-growing central peak, or jet (Figures 3aiv and 3biv and v), which is made of a mixture of drop and pool liquids [77]. In millimetre drop impact experiments (Figures 3aiv and v), the central peak is often observed to fragment due to the development of a capillary instability [89], 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 [76,86]. 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 [76].

### 2.2.3. Collapse of the central peak and further mixing

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 $Fr^*$, this mixture forms a well-defined compact cloud with little residual velocity [76], as can be seen for example on image v of Figure 3a. The dispersion of the drop liquid at this stage has been quantified experimentally by estimating the volume over which the drop liquid is spread [76]. This volume has been found to scale as $\sim R_i^3 ((\rho_i - \rho_t)/\rho_t)^{-0.57} Fr^{1.2}$ [76]: it is an increasing function of $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 $ρ_i/ρ_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 ≈2.

## 2.3. Limitations and future work

### 2.3.1. On the effect of compressibility

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 > 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 [79]. The impact crater opens in the decompressed, molten region behind the rarefaction wave (see for example Figure 4 from [102] showing snapshots from a numerical model of an impact at $M = 5$, or Figure 2 from [72]). This, in fact, justifies to some extent the use of liquids in laboratory experiments to model planetary impacts on solid targets [102]. However, compressibility does affect the maximum size of the crater, as shown by the systematic investigation of the effect of $Fr$ and $M$ on the maximum crater depth [102]. The effect of the Mach number on the maximum crater depth is moderate up to $M \sim 4$ (less than ~20%), but becomes more significant for $M \gtrsim 4$, a limit that is also relevant to planetary formation [102]. 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 [102], with the maximum crater depth scaling as $\sim R_i(FrM^{-2})^{1/4} = R_i(c^2/gR_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 known 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 [103–105]. This has been observed experimentally to have a significant effect on dispersion and mixing in jets and mixing layers when $M \gtrsim 1$ [104, 106].

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 [79]. 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 (10)). 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,
<figure><img src="image_14.png" /><figcaption>**Figure 4.** Time sequences from three oblique impact experiments with impact angles $\theta$ (with respect to the horizontal) equal to 59° (left column), 43° (middle column), and 25° (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.</figcaption></figure>

which shows that the flow around the crater is subsonic for a significant part of the crater opening phase [72, 102].

### 2.3.2. Oblique impacts

The experiments discussed in Section 2.2 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° [107]. 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 [72, 101, 107–111].
Figure 4 shows snapshots from yet unpublished experiments [112] 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°, 43°, or 25°. In contrast with vertical impacts, the crater, crown, and central peak are asymmetric, which is consistent with numerical models [72, 107] and fluid dynamics experiments using mm-sized drops [101, 110, 111]. The snapshots shown on the second line of Figure 4 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° and 25°. 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.

# 3. Metal-silicate separation in a magma ocean

## 3.1. Rationale

The experiments discussed in Section 2 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 (~$10^{-2}$ Pa·s), its velocity initially scales as $\sqrt{(\Delta\rho_0/\rho_s)gR_0}$, where $\Delta\rho_0$ is the effective density difference between the impactor/magma mixture and surrounding magma. Using this scaling to estimate $Re$ and $We$ gives

$$Re \sim \frac{\mathcal{B}_0^{1/2}}{v_s}, \qquad (13)$$

$$We \sim \frac{\Delta\rho_0 g R_0^2}{\sigma} = Bo, \qquad (14)$$

where $\mathcal{B}_0 = (\Delta\rho_0/\rho_s)gV_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 $Re \sim We \sim 10^{12}$ (taking $\sigma \sim 1$ J·m$^2$ for the metal-silicate surface tension, $g = 10$ m·s$^{-2}$, $\rho_s = 4000$ kg·m$^{-3}$, $v_s = 10^{-6}$ m$^2$·s$^{-1}$). The large value of $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 $We \to \infty$), which can be done experimentally by using miscible liquids instead of immiscible liquids. At large Reynolds numbers, a buoyant mass of fluid forms what is known as a *turbulent thermal*. This is discussed in Section 3.2, 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.

## 3.2. The infinite $We$ limit: miscible turbulent thermals

Figure 5a 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 [113–118]. 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 [37, 118]. This is what is known as a *turbulent thermal* in the fluid mechanics and atmospheric science communities [113, 115–117], 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 *turbulent thermal*, *buoyant cloud*, *metal-silicate cloud*, or simply *cloud*.

### 3.2.1. Large-scale asymptotic behaviour

A *miscible* turbulent thermal falling in an ambient liquid of density $ρ_s$ under 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 $ν$ and solute diffusivity $κ$ (the two latter properties being here assumed independent of composition). Experiments such as shown in Figure 5 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 ρ$ between the thermal and surrounding [113, 114]. Assuming that the loss of material from the thermal (by “detrainment”) is small [120] implies that the mass of solute within the thermal, which we denote by $M = \langle c \rangle V$, is conserved and equal to its initial value

$$M_0 = c_0 V_0. \qquad (15)$$

This implies that the total buoyancy $B = (\Delta ρ / ρ_s) gV$ of the thermal is constant as well, and equal to

$$B_0 = \frac{\Delta \rho_0}{\rho_s} g V_0, \qquad (16)$$

where $\Delta ρ_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 [113, 115, 117]. 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 $M_0$ and $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 [113], 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

$$R \sim B_0^{1/4} t^{1/2} \sim z, \qquad (17)$$

$$u_z \sim B_0^{1/4} t^{-1/2} \sim B_0^{1/2} z^{-1}, \qquad (18)$$

$$\langle c \rangle \sim M_0 B_0^{-3/4} t^{-3/2} \sim M_0 z^{-3}, \qquad (19)$$

when far from the release point. The proportionality between $r$ and $z$ predicted by Equation (17) is often written as

$$R = \alpha z, \qquad (20)$$

where the factor $\alpha$ is the so-called *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
<figure><img src="image_15.png" /><figcaption>**Figure 5.** (a) A 169 mL volume of an aqueous solution of NaI ($\rho_m = 1502 \text{ kg} \cdot \text{m}^{-3}$) falling into fresh water, at $Re = 2.9 \times 10^4$ (as defined by Equation (21)). The time interval between each image is 0.3 s. (b) A 169 mL volume of an aqueous solution of NaI ($\rho_m = 1601 \text{ kg} \cdot \text{m}^{-3}$) falling into a low viscosity silicon oil ($\rho_s = 821 \text{ kg} \cdot \text{m}^{-3}$, $\eta = 8.2 \times 10^{-4} \text{ Pa} \cdot \text{s}$), at $Bo = 9.3 \times 10^3$ and $Re = 4.0 \times 10^4$. The time interval between each image is 0.2 s. (c) The *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 fourth images of subfigure (b). Figures (a), (b), and (c) are modified from [119].</figcaption></figure>

is turbulent [117]. 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 [121, 122]. 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 (17) to (19) have
been extensively verified in experiments [115, 117, 123] and numerical simulations [120, 122]. 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$ [124]. The asymptotic expression 20 is therefore valid for $z \gtrsim R_0/\alpha$.

Interestingly, one can note that the Reynolds and Péclet 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éclet number of a turbulent thermal as

$$Re = \frac{\mathcal{B}_0^{1/2}}{\nu}, \qquad (21)$$

$$Pe_c = \frac{\mathcal{B}_0^{1/2}}{\kappa_c}. \qquad (22)$$

This will be further discussed in Section 3.2.4.

### 3.2.2. Turbulent entrainment model

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 [114]. This assumption implies that the time derivative of the volume of the thermal is given by

$$\frac{d}{dt}\left(\frac{4\pi}{3}R^3\right) = 4\pi R^2 \alpha u_z, \qquad (23)$$

where $\alpha$ is the entrainment coefficient. Since $d/dt = u_z d/dz$, integration of Equation (23) gives $R = R_0 + \alpha z$, which shows that the coefficient introduced in Equation (20) 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

$$\frac{d}{dt}\left(\frac{4\pi}{3}R^3 \bar{\rho} u_z\right) = \frac{4\pi}{3} R_0^3 \Delta \rho_0 g, \qquad (24)$$

where $\bar{\rho}(t)$ is the mean density of the thermal. Conservation of mass implies that

$$\frac{d}{dt}\left(\frac{4\pi}{3}R^3 \bar{\rho}\right) = 4\pi \rho_s R^2 \alpha u_z. \qquad (25)$$

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 (17) to (19)). The asymptotic laws are valid at times large compared to the timescale

$$t_g = \frac{R_0^2}{\mathcal{B}_0^{1/2}} \sim \sqrt{\frac{\rho_s R_0}{\Delta \rho_0 g}} \qquad (26)$$

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{(\Delta\rho_0/\rho_s)gR_0}$ with an acceleration $(\Delta\rho_0/\rho_s)g$).

### 3.2.3. Statistics of the concentration field

Measurements of the concentration field within turbulent miscible thermals have been carried out by [37] using Planar Laser Induced Fluorescence, and further analysed by [125]. Figure 6(a) shows a snapshot of the concentration field in a vertical plane within a thermal having $Re = 4175$ and $Pe_c = 1.0 \times 10^7$. The red line is an isocontour of the concentration field which is used
<figure><img src="image_16.png" /><figcaption>**Figure 6.** (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 Re = 4175 and a Péclet number of Pe<sub>c</sub> = 1.0 × 10<sup>7</sup>. The red line denotes the contour of the thermal with a threshold equal to c/⟨c⟩ = 0.5, with ⟨c⟩ the mean of the concentration inside the contour. (b) Probability density function (PDF) of the concentration c/c<sub>0</sub> as a function of the concentration c/c<sub>0</sub> for the snapshot (a) at t/t<sub>g</sub> = 51. (c) PDF of concentration normalized by the mean of the concentration c/⟨c⟩ as a function of c/⟨c⟩ for different times between t/t<sub>g</sub> = 20 and t/t<sub>g</sub> = 75. Here, red dashed lines show the temporal mean of the c/⟨c⟩. Inserts in (b) and (c) are the semi-log version of the PDF. Panel (a) adapted with permission from Lherm (2021) [37].</figcaption></figure>

to define the boundary of the thermal. Figure 6(a) shows that the concentration within the ther-
mal is far from being homogeneous, the concentration field exhibiting small-scale, highly corru-
gated elongated structures. The standard deviation of the concentration field (≈2.5 × 10⁻³c₀ for
this particular snapshot) is about half of the concentration average (≈4.2 × 10⁻³c₀), which shows
that the concentration has a relatively wide distribution. Concentration variations can be fur-
ther quantified by calculating the probability density function (PDF) of the concentration within
the thermal, as shown in Figure 6(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 concen-
tration is of the form

$$
\langle c^k \rangle = \left( M_0 \mathcal{B}_0^{-\frac{3}{4}} t^{-\frac{3}{2}} \right)^k f_{c,k}(Re, Pe_c), \quad (27)
$$
where Re and Pe<sub>c</sub> are defined according to Equations (21) and (22), and f<sub>c,k</sub> non-dimensional factors which may be functions of Re and Pe<sub>c</sub> [125]. It follows that the moments of c/⟨c⟩,

$$ \langle (c/\langle c \rangle)^k \rangle = \frac{\langle c^k \rangle}{\langle c \rangle^k} = \frac{f_{c,k}(Re, Pe_c)}{f_{c,1}^k(Re, Pe_c)}, \quad (28) $$

are predicted to be independent of time, which implies that the probability density function of c/⟨c⟩ is independent of time as well. Figure 6(c) shows that the PDF of c/⟨c⟩ calculated at different times in the asymptotic regime collapse on a well-defined master curve, thereby validating the above prediction. The PDFs of c/⟨c⟩ are well modelled by a decreasing exponential given by

$$ p(c/\langle c \rangle) = \exp\left(-\frac{c}{\langle c \rangle}\right) \quad (29) $$

with ⟨c⟩ given by Equation (19) (with a dimensionless prefactor adjusted to measurements of ⟨c⟩). The time average of the PDFs (dashed red line) indeed lies very close to the prediction of Equation (29) (grey line). Qualitatively, the stationary nature of the distribution of c/⟨c⟩ 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.

### 3.2.4. A note on self-similarity

The state at time t of a miscible thermal can be characterised by its Reynolds and Péclet numbers. As noted at the end of Section 3.2.1, it so happens that the instantaneous values of Re and Pe<sub>c</sub> estimated with the asymptotic predictions for R(t) and u<sub>z</sub>(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 Re and Pe<sub>c</sub>. This allows the velocity and composition fields to be self-similar in time, in the sense that, when normalised by their spatial average u<sub>z</sub> and ⟨c⟩, 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 (17) and (18)) is a decreasing function of time and distance given by

$$ We \sim \frac{\rho_s \mathcal{B}_0}{\sigma R_0} \left(\frac{t}{t_g}\right)^{-1/2} \sim Bo \left(\frac{t}{t_g}\right)^{-1/2} \sim \alpha Bo \frac{R_0}{z}, \quad (30) $$

where here the Bond number is defined as

$$ Bo = \frac{\Delta \rho_0 g R_0^2}{\sigma} \sim \frac{\rho_0 \mathcal{B}_0}{\sigma R_0}. \quad (31) $$

The proportionality between We and Bo reflects the force balance (inertia ~ buoyancy) which controls the velocity of the thermal.

Similarly, the Rossby number, which measures the relative importance of inertia and Coriolis force, writes

$$ Ro = \frac{u_z(t)}{\frac{1}{2\Omega R(t)}} \sim \frac{1}{t} \sim \left( \frac{\mathcal{B}_0^{1/4} \Omega^{-1/2}}{z} \right)^2 \quad (32) $$

when estimated with the asymptotic expressions for R(t) and u<sub>z</sub>(t). This predicts that the Coriolis force becomes comparable to inertia after a time interval ~1 day, or, equivalently, when the distance travelled by the thermal is larger than ~B<sub>0</sub><sup>1/4</sup>Ω<sup>-1/2</sup> [126–128].

Finally, if the thermal is made of buoyant droplets or particles, then an additional relevant number is the Rouse number R defined as the ratio of the sedimentation velocity of the droplets
or particles, $u_{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

$$ \mathcal{R} = \frac{u_{sed}}{u_z} \sim \frac{u_{sed}}{\mathcal{B}_0^{1/4}} t^{1/2} \sim \frac{u_{sed}}{\mathcal{B}_0^{1/2}} z, \quad (33) $$

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.

### 3.3. Immiscible turbulent thermals

Section 2 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 [86]. 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 [129]. Such particle-laden clouds behave as the miscible thermals described in the previous section [78, 128, 130]. In contrast, the cores of larger impactors coming onto Earth at the escape speed may not fragment during the impact [86]. 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 3.1.

#### 3.3.1. Experimental test of the turbulent entrainment model

Figure 5b 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 5a. The NaI solution and the silicon oil are *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 $Bo = 9.3 \times 10^3$ and $Re = 4.0 \times 10^4$ (with Bo and Re as defined by Equations (31) and (21)). We use here a Bond number rather than a Weber number, which is time-dependent in these experiments; as shown by Equation (30), 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 *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 5b [75, 124] as well as in experiments in which a volume of galinstan ($\rho = 6440 \text{ kg} \cdot \text{m}^{-3}$) is released in a pool of UCON oil ($\rho = 1050 \text{ kg} \cdot \text{m}^{-3}$) of different viscosities [87, 88, 131]: when $Bo \gtrsim 100$, the entrainment coefficient $\alpha = 0.25 \pm 0.1$ is indistinguishable from that of miscible thermals.

#### 3.3.2. 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 [75]. 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 5b) [75, 124]. 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 [51, 132].

At the very large values of We or Bo estimated in Section 3.1, 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 Bo ≳ 100 in the explored range of Bo ( ≲ 10⁴) [124]. Coupling this result with the disruption caused by the impact stage, [86] 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 We, the drop size d is comparable to the Hinze scale [133], 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

$$We_l = \frac{\rho_s u_l^2 l}{\sigma} \sim We \left(\frac{l}{L}\right)^{5/3}. \qquad (34)$$

The Hinze scale is such that inertia balances surface tension, i.e. $We_l \sim 1$. From Equation (34), one obtains

$$l_H \sim L We^{-3/5}. \qquad (35)$$

Assuming that the thermal speed arises from the buoyancy force only, we can take $L = R(t)$ and use Equation (30) for We to find

$$l_H \sim R_0 (\alpha^2 B o)^{-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}, \qquad (36)$$

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 [44, 73, 134].

However, scaling (35) 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 2), multiple choices are possible for the length scale L and the typical velocity U in the Weber number and the Hinze scale (35).

Experiments using galinstan and UCON oil have been used to investigate the effect of the viscosity contrast on the fragmentation behaviour [87, 88, 131, 135], though at smaller values of the Bond number ($Bo < 100$). Fragmentation has been observed to happen in a single global event, as described above for a viscosity ratio ~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 (35)) [88]. At the highest viscosity ratio investigated ($\eta_s/\eta_m = 691$), fragmentation has been observed to result from a different regime, the *Rayleigh–Taylor regime* [135], 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 ($Re < 80$)<fnref n="1" />, and it is not clear whether this regime persists at higher

<footnote n="1">¹The Reynolds number given here is calculated using Equation (21), which gives a value twice larger than quoted in [88, 135].</footnote>
values of Re. In fact 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.

### 3.3.3. 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 drops.

In experiments on immiscible turbulent thermals such as shown in Figure 5 or described by [75, 124], 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 (18)) is a decreasing function of time, as shown by Equation (33).

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 particle swarms when $\mathcal{R} \gtrsim 1$, falling at a velocity near the settling velocity of individual particles [128–130]. When estimated using the asymptotic scaling law for the turbulent thermal velocity (Equation (18)), the Rouse number is given by Equation (33). 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 (33) 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_{sed}$. [130] found experimentally that the fallout distance $z_{fall}$ at which particle-ladden turbulent thermals transition to particle swarms indeed scale as $\mathcal{B}_0^{1/2}/u_{sed}$, with a prefactor $\approx 9 \pm 2$, when the particle Reynolds number is larger than $\approx 4.4$.

When the drop Reynolds number is large, then the settling velocity is $u_s \sim \sqrt{\Delta\rho/\rho_s ga}$ where a is the particle or drop radius. Using this scaling, [130] found experimentally that the ratio of the fallout distance to particle radius is only a function of the number of particles $N_p$, $z_{fall}/a \approx (9 \pm 2)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 size as

$$ \frac{z_{\text{fall}}}{R_0} \approx (9 \pm 2) \left( \phi_0 \frac{R_0}{a} \right)^{1/2}, \qquad (37) $$

which is more convenient for application to metal-silicate separation in a magma ocean. Taking $a \sim 1$ cm, $z_{fall}$ would only be smaller than ~3000 km (the thickness of the mantle) if $R_0$ is smaller than ~1 km. The critical value of $R_0$ decreases to ~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.
<figure><img src="image_17.png" /><figcaption>(b) Oblique impact involving miscible liquids | $Fr = 7.1$, $Re = 5.8 \times 10^4$, $\rho_i / \rho_t = 1.088$, $\theta = 52.5^\circ$</figcaption></figure>

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

**Figure 7.** (a) Time sequence from an experiment in which a volume of water ($R_i = 2.3$ cm, $v = 10^{-6}$ m$^2 \cdot s^{-1}$, $\rho_i = 1000$ kg$\cdot$m$^{-3}$) strikes a bath of low-viscosity silicon oil ($v = 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, $v = 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^\circ$ with respect to the horizontal.

## 3.4. Additional effects

### 3.4.1. Impact-generated thermals

The experiments discussed in Sections 3.2 to 3.4.2 have been performed by releasing volumes of dense liquid within the host liquid, without accounting for the impact and cratering phase described in Section 2. 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 [76] and immiscible liquid pairs (Figure 7a) 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 [76], with $Fr \in [4, 100]$. Landeau et al. [76] 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

$$R = R_0 + \alpha z, \qquad (38)$$

as predicted by the turbulent entrainment model, if the radius $R_0$ of the cloud at $z = 0$ is taken to be

$$R_0 = 0.29 R_i \left( \frac{\rho_m - \rho_s}{\rho_s} \right)^{-0.19} Fr^{0.4} \qquad (39)$$

instead of the impactor radius $R_i$. This accounts for the dispersion induced by the impact which is, as expected, an increasing function of $Fr$.
Figure 7(b) shows snapshots from an experiment in which a volume of NaCl solution ($R_i = 3.6$ cm, $\rho_i = 1088$ kg·m⁻³, appearing blue) strikes the surface of a water pool with an impact angle $\theta = 52.5^\circ$, at $Fr = 7.1$ and $Re = 5.8 \times 10^4$. In contrast with the snapshots shown in Figure 4, 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 3.2.2), 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 $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.

### 3.4.2. Effect of planetary rotation

As discussed in Section 3.2.4, 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}$ [126–128]. This critical depth can also be written as

$$z_{\Omega} \sim \left( \frac{1 \text{ day}}{t_{g}} \right)^{1/2} R_{0}, \qquad (40)$$

where $t_g$ is the acceleration time scale defined in Equation (26). The prefactor in Equation (40) can be estimated from the asymptotic solutions of the turbulent entrainment model (e.g. [75, 136]), 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 ~400 km and ~12000 km, which needs to be compared to the thickness of the magma ocean, which is at most ~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 [78]: 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} \approx 1.6$, while increasing by a factor $(5/2)^{2/3} \approx 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 $\approx 200$ km.

Experiments on miscible turbulent thermals [126, 127] and cloud particles [78, 128] 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 ~$\alpha z_\Omega$ but expands in the vertical direction. The amount of dispersion in this regime has been quantified and modelled by [78]. 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 [76], Kriaa et al. [78] 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 of the 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$ [78]. 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 [78].

The experiments discussed above [78, 126-128] 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$.

### 3.5. *Chemical equilibration within the cloud*

Published models of metal-silicate equilibration assume that the metal phase fully equilibrates with the entrained silicates [75, 76], 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 3.2.3. This is also inconsistent with measurements made by [88] 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 [88].

The implications for the degree of equilibration of the concentration density function given in Equation (29) 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 (29) is equal to $(\ln 2)\langle c \rangle \approx 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.

### 3.6. *Reaching the bottom of the magma ocean*

The metal cloud may still have a velocity up to several hundreds of m·s⁻¹ 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 [66]. Depending on the extent of mantle melting, several scenarios can be envisaged [75]. 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 [137]. 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 [66]. 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 ~10 m·s⁻¹ [64], significantly larger than the settling velocity *v*<sub>s</sub> of a mm-to-cm sized drops, which is around 0.3–0.5 m·s⁻¹ [44, 74].

Laboratory experiments have shown that sedimentation from a convective layer happens on a timescale ~*H*/*v*<sub>s</sub>, where *H* is the thickness of the convecting layer, even if the settling velocity is much smaller than the convective velocity [138–140]. 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 [138]. With *v*<sub>s</sub> ~ 0.3–0.5 m·s<sup>-1</sup> [44, 74], the sedimentation time from a 2000 km deep magma ocean would be *H*/*v*<sub>s</sub> ~ 2 months. This is small compared to the lifetime of a magma ocean (≥1000 years [64]).

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 [140, 141]. 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 [140, 141] to estimate an upper bound on the equilibrium metal volume fraction in a convecting magma ocean. Using Equation (12) from [140] with a magma ocean thickness of 2000 km and a Rayleigh number in the range 10<sup>28</sup>–10<sup>32</sup>, we find the predicted equilibrium metal volume fraction to be in the range 10<sup>-10</sup>–10<sup>-5</sup>, 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 [142]. In this case, Maas et al. [143] 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 [143].

# 4. Metal-silicate separation in a solid mantle

## 4.1. Rayleigh–Taylor instability and diapirism

### 4.1.1. *A simple reference 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 3.6, 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 (≥1000 years [64]). A layer of metal would therefore accumulate at the base of the magma ocean [38].
Molten iron being denser than the underlying silicates, this situation is gravitationally unstable and gives rise to Rayleigh–Taylor instabilities [38, 42, 144–146].

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 [149]. The timescale of instability growth is

$$ \tau_{RT} \sim \frac{\eta_s}{\Delta\rho g \lambda} \sim \left( \frac{\eta_s}{10^{18} \text{ Pa}\cdot\text{s}} \right) \left( \frac{10^3 \text{ km}}{\lambda} \right) \times 1 \text{ year}, \quad (41) $$

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} \text{ Pa}\cdot\text{s}$, which corresponds to the estimate by [43] just below the solidus, and a wavelength of $\sim 10^3 \text{ km}$, $\tau_{RT}$ is on the order of 1 year. A viscosity of $10^{18} \text{ Pa}\cdot\text{s}$ corresponds to the estimate by [43] 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 \text{ years}$ [64]) 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 [147, 150]. Figure 8a shows snapshots from one such experiment, in which a thin layer of glycerin with a density $\rho_m = 1250 \text{ kg}\cdot\text{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 \text{ kg}\cdot\text{m}^{-3}$ [147]. Figure 8a 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 [147]. 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)

$$ u_z = \frac{\Delta\rho g R^2}{3\eta_s}, \quad (42) $$

which gives a transit time

$$ \tau_{tr} = \frac{3\eta_s H}{\Delta\rho g R^2} \approx \left( \frac{\eta_s}{10^{18} \text{ Pa}\cdot\text{s}} \right) \left( \frac{R}{100 \text{ km}} \right)^2 \left( \frac{H}{1000 \text{ km}} \right) \times 250 \text{ year}. \quad (43) $$

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 A 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\to s}$ from the metal diapir to the surrounding, and the dissipation $\psi_m$ in the metal, the energy balance of a diapir writes

$$ Mc_{p,m} \frac{dT}{dt} = M \frac{\alpha_m T}{\rho_m} \frac{dP}{dt} - Q_{m\to s} + \psi_m \quad (44) $$

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 $d(\cdots)/dt = u_z d(\cdots)/dz$, as

$$ \frac{dT}{dz} = \frac{\alpha_m g T}{c_{p,m}} - \frac{Q_{m\to s}}{Mc_{p,m} u_z} + \frac{\psi_m}{Mc_{p,m} u_z}. \quad (45) $$

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\to s} = 0$) and reversible ($\psi_m = 0$) way.
<figure><img src="image_19.png" /><figcaption>**Figure 8.** (a) Snapshots extracted and modified from Figure 2 of [147], with permission from Wiley (CCC license number 5781330644178). The dark liquid is glycerin ($\eta = 175$ Pa·s, $\rho = 1250$ kg·m⁻³) and the ambient liquid (transparent) is silicon oil ($\eta = 5520$ Pa·s, $\rho = 920$ kg·m⁻³). 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·m⁻³, $\eta_m = 10^{-3}$ Pa·s for gallium, and $\rho = 1256$ kg·m⁻³, $\eta_s = 0.1$ Pa·s and $\rho = 1417$ kg·m⁻³, $\eta_s = 11.75$ Pa·s for the s1 and s2 glucose solutions. (b) and (c) are parts of Figures 1 and 2 of [148], under the CC BY 4.0 license (http://creativecommons.org/licenses/by/4.0/).</figcaption></figure>
Assuming the diapir keeps a nearly spherical shape, the heat flux from the diapir to its surroundings is

$$Q_{m \to s} = a4\pi R^2 k_s \frac{T - T_s}{R} Pe^{1/2}, \quad (46)$$

with $a \approx 0.46$ [151, 152], $T_s$ the temperature of the surrounding mantle, and $k_s$ the thermal conductivity of the silicate mantle. The Péclet number $Pe$ is here defined with the diapir velocity and mantle thermal diffusivity. Assuming the diapir falls at the Stokes velocity,

$$Pe = \frac{Ru_z}{\kappa_s} = \frac{\Delta\rho g R^3}{3\kappa_s \eta_s} \sim \frac{\mathcal{B}}{\kappa_s v_s}, \quad (47)$$

where $\mathcal{B} = (4\pi/3)R^3\Delta\rho/\rho_s g$ is the total buoyancy of the diapir, and $v_s = \eta_s/\rho_s$ and $\kappa_s$ the kinematic viscosity and thermal diffusivity of the silicate. As defined here, $Pe$ is akin to a Rayleigh number. Using Equation (46), Equation (45) can be recast as

$$\frac{dT}{dz} + \left(\frac{1}{\ell_{\kappa}} - \frac{1}{\ell_{c}}\right) T = \frac{T_s}{\ell_{\kappa}} + \frac{\psi_m}{Mc_{p,m}}. \quad (48)$$

where

$$\ell_c = \frac{c_{p,m}}{\alpha_m g}, \quad (49)$$

$$\ell_{\kappa} = \frac{Pe^{1/2}}{3a} \frac{\rho_m c_{p,m}}{\rho_s c_{p,s}} R, \quad (50)$$

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 thermal equilibration distance with the surrounding mantle. Replacing $Pe$ with its expression (47), we find that $\ell_{\kappa}$ is larger than the solid mantle thickness $H$ when the diapir radius is larger than

$$R_{\text{eq},1} = H27^{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} \approx \left( \frac{\eta_s}{10^{18} \text{ Pa}\cdot\text{s}} \right)^{1/5} \left( \frac{H}{10^3 \text{ km}} \right)^{2/5} \times 10 \text{ km}. \quad (51)$$

Note that $R_{\text{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 [46, 153].

The ratio of the compression and heat exchange terms (1st and 2nd terms on the RHS of (45)) 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_{\text{eq},2}$ given by

$$R_{\text{eq},2} = \left( 3a \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} \text{ Pa}\cdot\text{s}} \right)^{1/5} \times 20 \text{ km}. \quad (52)$$

The ratio $\ell_{\kappa}/\ell_c = (R/R_{\text{eq},2})^{5/2}$ strongly increases with $R$, which implies that heat exchange quickly becomes negligible compared to compression heating when $R$ exceeds $R_{\text{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 [46]), and $\sim (\rho_m/\rho_s)(\Delta\rho g R^3/\eta_s^2)$ if the flow within the diapir is turbulent [20]. 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 ≥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.

### 4.1.2. Additional complexities

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

**Three-layer systems.** The experiments shown in Figure 8 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 [148, 154]. 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 8b shows a diapir that developed from this instability [148]. The structure differs markedly from the diapirs obtained in two-layer experiments (Figure 8a). 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 [148, 154]. 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 [148].

Previous studies [148, 154] 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 8c 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 [19, 148, 154] argue that dissipation heating in such metal-silicate plume could heat up both phases and increase the temperature of the metal by ~1000 K. This contrasts with the reference model of Section 4.1.1 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, [148, 154] argue that the conduit can be kept open for a significant time if additional metal is added before the conduit closure by viscous relaxation.

**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 [38]. The wavelength $\lambda$ can be smaller if the effective viscosity increases with depth, though likely not to the degree ($\lambda$ ~ a few km) proposed in 1981 by Stevenson [38], 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 (42). This has been investigated theoretically [155-157] and experimentally [157, 158] in the context of magma intrusions (plutons)
in Earth's crust. It is found that the drag is a function of 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 Stokes limit as described in Section 4.1.1, in which the drag is controlled by the viscosity far from the diapir; and (ii) the lubrication limit, in which deformation is localized in a thin low viscosity layer coating the forward side of the diapir [157]. Using the prediction of [157] 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

$$ \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. \qquad (53) $$

Taking for example $R = 10$ km and $\eta_{s,\infty} = 10^{18}$ Pa·s, this would be verified if $\eta_{s,0} \gtrsim 1$ Pa·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 4.1.1 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 Pe value of metal diapirs: since the boundary layer thickness is $\sim RPe^{-1/2}$, only a small volume around the diapir would have a reduced viscosity when 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 4.1.1, dissipative heating is expected to be localized in the adjacent mantle, in agreement with numerical simulations having $\eta_s/\eta_m \gg 1$ [46, 153]. 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 [45] 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 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.

**Diapirism vs. diking.** The possibility that the metal migrated through the solid mantle as dikes has been proposed by [38]. In the 1981's paper in which this mechanism was proposed [38], the inner part of the planet was hypothesized to be cold and undifferentiated, with an effective viscosity possibly in excess of $10^{26}$ Pa·s [43]. In this limit the migration of metal diapirs would have been very slow [38], 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 [159]. 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 [160] 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·s [4]. 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 on how viscosity of
the host rock affects the shape of a volume of magma propagating along a pre-existing crack. 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 [161] 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.

## 4.2. 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 [148, 154, 162].

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 [163–167]. 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 [4]), and only list a few salient points. Iron-rich melts with geologically relevant compositions do not wet silicates, except at low pressure (≈3 GPa) at high O and S content and perhaps in the lower mantle, though this is debated [4]. As a result, the percolation of iron-rich melt is usually only possible if its volume fraction exceeds a threshold volume fraction (~5 vol%) above which it is interconnected. The threshold melt fraction may be smaller (~1 vol%) if the material undergoes shear deformation [168], though this is also debated [4]. Partial melting of the silicates further complicates the picture, with experiments suggesting the formation of aggregates of silicate crystals and molten iron [169].

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 [170]): the Darcy velocity

$$V_D = \frac{\Delta\rho g k_0}{\eta_m}, \qquad (54)$$
which is independent of $\ell$, and the compaction velocity

$$V_C = \frac{\Delta \rho g \ell^2}{\eta_s} \qquad (55)$$

The ratio of these two velocity scales is

$$\frac{V_D}{V_C} = \frac{\eta_s k_0}{\eta_m \ell^2} = \left(\frac{\delta}{\ell}\right)^2, \qquad (56)$$

where

$$\delta = \sqrt{\frac{\eta_s}{\eta_m} k_0} \qquad (57)$$

is the so-called compaction length. This name can be somewhat misleading and care must be taken when interpreting the value of $\delta$. As shown by Equation (56), $\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 [170].

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 [170–173]. 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 [162], 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 ≈600), using silicon oil and hydroxyethylcellulose (HEC) solutions as working fluids. Figure 9 shows a snapshot of an experiment using a silicon oil with viscosity 12.1 Pa·s (appearing yellow on the pictures) and a HEC solution with viscosity 0.023 Pa·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 9). 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 4.1.2 and shown in Figure 8 in which the less viscous phase (gallium) was observed to form drops [148,154].

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 [162] and the permeability law they used to model the experiment, namely the Blake–Kozeny–Carman permeability law

$$k = \frac{a^2 \phi^3}{K_0 (1 - \phi)^2}, \qquad (58)$$

where $a$ is the droplet size and $K_0 \approx 15$ is an empirical constant. With $a \sim 0.1$ mm as estimated in the paper and a maximum viscosity ratio of ≈600, one finds that $\delta$ is at most ≈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 [162]. 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
<figure><img src="image_20.png" /><figcaption>**Figure 9.** Snapshot from an experiment in which a mixture of silicon oil with viscosity 12.1 Pa·s (appearing yellow on the pictures) and hydroxyethylcellulose (HEC) solution with viscosity 0.023 Pa·s (appearing blue) separate under the action of gravity. The density difference between the HEC solution and the silicon oil is 31 kg·m⁻³. 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 [162] Reproduced from Figure 9 of [162], 2007 © Cambridge University Press, reproduced with permission (Licence number 5797030122437).</figcaption></figure>

larger than the drop size, which suggests a tank size of at least a few centimeters as in the experiments of Sato & 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·s and water ($10^{-3}$ Pa·s) to reach a viscosity ratio of $\sim 10^6$. With a drop size $a \sim 0.1$ mm and assuming the permeability law (58), 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.

**5. 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éclet 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 *Re* and *We* values around 10⁴ 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 compress-
ibility (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 cau-
tion, because of the wild extrapolations required to apply them to the early Earth. For exam-
ple, the metal-silicate equilibration model built from experiments such as described in Section 3
has been incorporated in several geochemical models of core-mantle formations with some
success [15, 36, 174]. 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.

**Declaration of interests**

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

**Acknowledgments**

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é de France for granting their “Ampère 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étologie (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.
## Appendix A. Temperature inside a metal diapir

The model of diapir thermal evolution presented in Section 4.1.1 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 [175–177]. According to these studies, the heat flux associated with internal convection is given by

$$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} \quad (59)$$

where $a_d = 0.553$ and $\beta = 0.2521$ [177], and

$$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}. \quad (60)$$

When estimated on the mantle side, the heat flux is given by

$$q = a k_s \frac{\Delta T - \Delta T_d}{R} Pe^{1/2}. \quad (61)$$

Equating the two expressions for $q$ and re-arranging, we find that

$$\frac{(\Delta T_d / \Delta T)^{1+\beta}}{1 - \Delta T_d / \Delta T} = \frac{a_o k_s}{a_d k_m} Pe^{1/2} Ra_{\Delta T}^{-\beta} \propto R^{0.74} \Delta T^{-0.25} \quad (62)$$

Taking $R = 100$ km and $\Delta T = 100$ K gives $\Delta T_d / \Delta T \approx 1.3\%$. Taking $R = 500$ km instead gives $\Delta T_d / \Delta T \approx 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.

## References

[1] D. J. Stevenson, "Evolution of the Earth: An introduction and overview", in *Treatise on Geophysics* (G. Schubert, ed.), vol. 9, Elsevier, Oxford, second ed., 2015, ISBN 978-0-444-53803-1, p. 1-9.

[2] A. Morbidelli, W. F. Bottke, D. Nesvorný, H. F. Levison, "Asteroids were born big", *Icarus* **204** (2009), no. 2, p. 558-573.

[3] J. I. Lunine, D. P. O'Brien, S. N. Raymond, A. Morbidelli, T. Quinn, A. L. Graps, "Dynamical models of terrestrial planet formation", *Adv. Sci. Lett.* **4** (2011), no. 2, p. 325-338.

[4] D. C. Rubie, F. Nimmo, H. J. Melosh, "Formation of the Earth's core", in *Treatise on Geophysics* (G. Schubert, ed.), vol. 9, Elsevier, Oxford, second ed., 2015, ISBN 978-0-444-53803-1, p. 43-79.

[5] B. Wood, M. Walter, J. Wade, "Accretion of the Earth and segregation of its core", *Nature* **441** (2006), no. 7095, p. 825-833.

[6] Q. Yin, S. B. Jacobsen, K. Yamashita, J. Blichert-Toft, P. Telouk, F. Albarede, "A short timescale for terrestrial planet formation from Hf/W chronometry of meteorites", *Nature* **418** (2002), no. 6901, p. 949-952.

[7] T. Kleine, C. Münker, K. Mezger, H. Palme, "Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry", *Nature* **418** (2002), p. 952-955.

[8] J. F. Rudge, T. Kleine, B. Bourdon, "Broad bounds on Earth's accretion and core formation constrained by geochemical models", *Nat. Geosci.* **3** (2010), p. 439-443.

[9] A. Corgne, S. Keshav, B. J. Wood, W. F. McDonough, Y. Fei, "Metal silicate partitioning and constraints on core composition and oxygen fugacity during Earth accretion", *Geochim. Cosmochim. Acta* **72** (2008), p. 574-589.
[10] J. Siebert, A. Corgne, F. J. Ryerson, "Systematics of metal-silicate partitioning for many siderophile elements applied to Earth's core formation", *Geochim. Cosmochim. Acta* **75** (2011), no. 6, p. 1451-1489.

[11] C. K. Gessmann, B. J. Wood, "Potassium in the Earth's core?", *Earth Planet. Sci. Lett.* **200** (2002), p. 63-78.

[12] A. Corgne, S. Keshav, Y. Fei, W. F. McDonough, "How much potassium is in the Earth's core? New insights from partitioning experiments", *Earth Planet. Sci. Lett.* **256** (2007), p. 567-576.

[13] M. A. Bouhifd, L. Gautron, N. Bolfan-Casanova, V. Malavergne, T. Hammouda, D. Andrault, A. P. Jephcoat, "Potassium partitioning into molten iron alloys at high-pressure: Implications for Earth's core", *Phys. Earth Planet. Inter.* **160** (2007), p. 22-33.

[14] J. G. O'Rourke, D. J. Stevenson, "Powering Earth's dynamo with magnesium precipitation from the core", *Nature* **529** (2016), no. 7586, p. 387-389.

[15] J. Badro, J. Siebert, F. Nimmo, "An early geodynamo driven by exsolution of mantle components from Earth's core", *Nature* **536** (2016), no. 7616, p. 326-328.

[16] J. Siebert, J. Badro, D. Antonangeli, F. J. Ryerson, "Metal-silicate partitioning of Ni and Co in a deep magma ocean", *Earth Planet. Sci. Lett.* **321** (2012), p. 189-197.

[17] J. Badro, J. P. Brodholt, H. Piet, J. Siebert, F. J. Ryerson, "Core formation and core composition from coupled geochemical and geophysical constraints", *Proc. Natl. Acad. Sci.* **112** (2015), no. 40, p. 12310-12314.

[18] R. A. Fischer, Y. Nakajima, A. J. Campbell *et al.*, "High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O", *Geochim. Cosmochim. Acta* **167** (2015), p. 177-194.

[19] C. King, P. Olson, "Heat partitioning in metal-silicate plumes during Earth differentiation", *Earth Planet. Sci. Lett.* **304** (2011), no. 3-4, p. 577-586.

[20] V. Clesi, R. Deguen, "Linking the core heat content to Earth's accretion history", *Geochem. Geophys. Geosyst.* **24** (2023), no. 5, article no. e2022GC010661.

[21] S. Labrosse, C. Jaupart, "Thermal evolution of the Earth: Secular changes and fluctuations of plate characteristics", *Earth Planet. Sci. Lett.* **260** (2007), p. 465-481.

[22] J.-P. Williams, F. Nimmo, "Thermal evolution of the Martian core: Implications for an early dynamo", *Geology* **32** (2004), p. 97-100.

[23] J. Monteux, A. M. Jellinek, C. L. Johnson, "Why might planets and moons have early dynamos?", *Earth Planet. Sci. Lett.* **310** (2011), no. 3, p. 349-359.

[24] A. Halliday, "Mixing, volatile loss and compositional change during impact-driven accretion of the Earth", *Nature* **427** (2004), no. 6974, p. 505-509.

[25] T. Kleine, K. Mezger, H. Palme, C. Münker, "The W isotope evolution of the bulk silicate Earth: Constraints on the timing and mechanisms of core formation and accretion", *Earth Planet. Sci. Lett.* **228** (2004), p. 109-123.

[26] F. Nimmo, D. P. O'Brien, T. Kleine, "Tungsten isotopic evolution during late-stage accretion: Constraints on Earth-Moon equilibration", *Earth Planet. Sci. Lett.* **292** (2010), no. 3-4, p. 363-370.

[27] A. N. Halliday, B. J. Wood, "The composition and major reservoirs of the earth around the time of the moon-forming giant impact", in *Treatise on Geophysics* (G. Schubert, ed.), vol. 9, Elsevier, Oxford, second ed., 2015, ISBN 978-0-444-53803-1, p. 11-42.

[28] V. R. Murthy, "Early differentiation of the earth and the problem of mantle siderophile elements: A new approach", *Science* **253** (1991), no. 5017, p. 303-306.

[29] A. Ringwood, "Chemical evolution of the terrestrial planets", *Geochim. Cosmochim. Acta* **30** (1966), no. 1, p. 41-104.

[30] J. Li, C. Agee, "Geochemistry of mantle-core differentiation at high pressure", *Nature* **381** (1996), no. 6584, p. 686-689.

[31] N. Chabot, D. Draper, C. Agee, "Conditions of core formation in the Earth: Constraints from Nickel and Cobalt partitioning", *Geochim. Cosmochim. Acta* **69** (2005), no. 8, p. 2141-2151.

[32] J. Siebert, J. Badro, D. Antonangeli, F. J. Ryerson, "Terrestrial accretion under oxidizing conditions", *Science* **339** (2013), no. 6124, p. 1194-1197.

[33] J. Wade, B. J. Wood, "Core formation and the oxidation state of the Earth", *Earth Planet. Sci. Lett.* **236** (2005), p. 78-95.

[34] B. J. Wood, "Accretion and core formation: Constraints from metal-silicate partitioning", *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* **366** (2008), no. 1883, p. 4339-4355.

[35] D. C. Rubie, D. J. Frost, U. Mann, "Heterogeneous accretion, composition and core-mantle differentiation of the Earth", *Earth Planet. Sci. Lett.* **301** (2011), no. 1, p. 31-42.

[36] D. C. Rubie, S. A. Jacobson, D. P. O'Brien, E. D. Young, J. de Vries, F. Nimmo, H. Palme, D. J. Frost, "Accretion and differentiation of the terrestrial planets with implications for the compositions of early-formed Solar System bodies and accretion of water", *Icarus* **248** (2015), p. 89-108.

[37] V. Lherm, *Dynamique du fractionnement thermique et chimique lors de la différenciation des planètes telluriques*, Phd thesis, Université de Lyon, 2021, https://theses.hal.science/tel-03377150.

[38] D. J. Stevenson, "Models of the Earth's core", *Science* **214** (1981), no. 4521, p. 611-619.

[39] W. B. Tonks, H. J. Melosh, "Core formation by giant impacts", *Icarus* **100** (1992), no. 2, p. 326-346.
[40] W. B. Tonks, H. J. Melosh, "Magma ocean formation due to giant impacts", *J. Geophys. Res.* **98** (1993), p. 5319-5333.

[41] J. T. Gu, R. A. Fischer, M. C. Brennan, M. S. Clement, S. A. Jacobson, N. A. Kaib, D. P. O'Brien, S. N. Raymond, "Comparisons of the core and mantle compositions of earth analogs from different terrestrial planet formation scenarios", *Icarus* **394** (2023), article no. 115425.

[42] D. J. Stevenson, "Fluid dynamics of core formation", in *Origin of the Earth* (H. E. Newsom, J. H. Jones, eds.), Oxford University Press, New York, 1990, p. 231-249.

[43] S.-i. Karato, V. R. Murthy, "Core formation and chemical equilibrium in the Earth—I. Physical considerations", *Phys. Earth Planet. Inter.* **100** (1997), no. 1-4, p. 61-79.

[44] D. C. Rubie, H. J. Melsoh, J. E. Reid, C. Liebske, K. Righter, "Mechanisms of metal-silicate equilibration in the terrestrial magma ocean", *Earth Planet. Sci. Lett.* **205** (2003), no. 3-4, p. 239-255.

[45] H. Samuel, P. J. Tackley, "Dynamics of core formation and equilibration by negative diapirism", *Geochem. Geophys. Geosyst.* **9** (2008), no. 6, article no. Q06011.

[46] J. Monteux, Y. Ricard, N. Coltice, F. Dubuffet, M. Ulvrova, "A model of metal-silicate separation on growing planets", *Earth Planet. Sci. Lett.* **287** (2009), no. 3-4, p. 353-362.

[47] Y. Thibault, M. J. Walter, "The influence of pressure and temperature on the metal-silicate partition coefficients of nickel and cobalt in a model C1 chondrite and implications for metal segregation in a deep magma ocean", *Geochim. Cosmochim. Acta* **59** (1995), no. 5, p. 991-1002.

[48] S. R. Taylor, M. D. Norman, "Accretion of differentiated planetesimals to the Earth", in *Origin of the Earth* (H. E. Newsom, J. H. Jones, eds.), 1990, p. 29-43.

[49] R. C. Greenwood, I. A. Franchi, A. Jambon, P. C. Buchanan, "Widespread magma oceans on asteroidal bodies in the early solar system", *Nature* **435** (2005), no. 7044, p. 916-918.

[50] E. Villermaux, "Mixing versus stirring", *Annu. Rev. Fluid Mech.* **51** (2019), p. 245-273.

[51] E. Villermaux, "Fragmentation versus Cohesion", *J. Fluid Mech.* **898** (2020), article no. P1.

[52] S. N. Raymond, D. P. OBrien, A. Morbidelli, N. A. Kaib, "Building the terrestrial planets: Constrained accretion in the inner solar system", *Icarus* **203** (2009), no. 2, p. 644-662.

[53] G. Morard, J. Siebert, D. Andrault, N. Guignot, G. Garbarino, F. Guyot, D. Antonangeli, "The Earth's core composition from high pressure density measurements of liquid iron alloys", *Earth Planet. Sci. Lett.* **373** (2013), p. 169-178.

[54] F. J. Spera, M. S. Ghiorso, D. Nevins, "Structure, thermodynamic and transport properties of liquid MgSiO3: Comparison of molecular models and laboratory results", *Geochim. Cosmochim. Acta* **75** (2010), p. 1272-1296.

[55] J. P. Poirier, "Transport properties of liquid metals and viscosity of the Earth's core", *Geophys. J. Int.* **92** (1988), no. 1, p. 99-105.

[56] G. A. De Wijs, G. Kresse, L. Vovcadlo, D. Dobson, D. Alfe, M. J. Gillan, G. D. Price, "The viscosity of liquid iron at the physical conditions of the Earth's core", *Nature* **392** (1998), no. 6678, p. 805-807.

[57] B. B. Karki, L. P. Stixrude, "Viscosity of MgSiO3 liquid at Earths mantle conditions: Implications for an early magma ocean", *Science* **328** (2010), no. 5979, p. 740-742.

[58] L. Xie, A. Yoneda, T. Katsura, D. Andrault, Y. Tange, Y. Higo, "Direct viscosity measurement of peridotite melt to lower-mantle conditions: A further support for a fractional magma-ocean solidification at the top of the lower mantle", *Geophys. Res. Lett.* **48** (2021), no. 19, article no. e2021GL094507.

[59] T. Scott, D. L. Kohlstedt, "The effect of large melt fraction on the deformation behavior of peridotite", *Earth Planet. Sci. Lett.* **246** (2006), no. 3-4, p. 177-187.

[60] L. Hansen, D. L. Kohlstedt, "2.18 constitutive equations, rheological behavior, and viscosity of rocks", *Treatise Geophys.* **2** (2015), p. 441-472.

[61] Q. Williams, "The thermal conductivity of Earth's core: A key geophysical parameter's constraints and uncertainties", *Annu. Rev. Earth Planet. Sci.* **46** (2018), no. 1, p. 47-66.

[62] H. Ni, H. Hui, G. Steinle-Neumann, "Transport properties of silicate melts", *Rev. Geophys.* **53** (2015), no. 3, p. 715-744.

[63] H. Terasaki, S. Urakawa, D. C. Rubie, K.-i. Funakoshi, T. Sakamaki, Y. Shibazaki, S. Ozawa, E. Ohtani, "Interfacial tension of Fe-Si liquid at high pressure: Implications for liquid Fe-alloy droplet size in magma oceans", *Phys. Earth Planet. Inter.* **202** (2012), p. 1-6.

[64] V. Solomatov, "Magma oceans and primordial mantle differentiation", in *Treatise on Geophysics* (G. Schubert, ed.), vol. 9, Elsevier, Oxford, second ed., 2015, ISBN 978-0-444-53803-1, p. 81-104.

[65] T. W. Dahl, D. J. Stevenson, "Turbulent mixing of metal and silicate during planet accretion — And interpretation of the Hf-W chronometer", *Earth Planet. Sci. Lett.* **295** (2010), p. 177-186.

[66] R. Deguen, P. Olson, P. Cardin, "Experiments on turbulent metal-silicate mixing in a magma ocean", *Earth Planet. Sci. Lett.* **310** (2011), p. 303-313.

[67] W. Benz, W. L. Slattery, A. G. W. Cameron, "Collisional stripping of Mercury's mantle", *Icarus* **74** (1988), p. 516-528.

[68] R. M. Canup, C. B. Agnor, "Accretion of the terrestrial planets and the Earth-Moon system", in *Origin of the Earth and Moon* (R. M. Canup, K. Righter, eds.), University of Arizona Press, Tucson, AZ, 2000, p. 113-129.
[69] M. Ćuk, S. T. Stewart, "Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning", *Science* **338** (2012), no. 6110, p. 1047-1052.

[70] M. Nakajima, D. J. Stevenson, "Melting and mixing states of the Earth's mantle after the Moon-forming impact", *Earth Planet. Sci. Lett.* **427** (2015), p. 286-295.

[71] J. A. Kegerreis, S. Ruiz-Bonilla, V. R. Eke, R. J. Massey, T. D. Sandnes, L. F. Teodoro, "Immediate origin of the moon as a post-impact satellite", *Astrophys. J. Lett.* **937** (2022), no. 2, article no. L40.

[72] J. D. Kendall, H. J. Melosh, "Differentiated planetesimal impacts into a terrestrial magma ocean: fate of the iron core", *Earth Planet. Sci. Lett.* **448** (2016), p. 24-33.

[73] H. Ichikawa, S. Labrosse, K. Kurita, "Direct numerical simulation of an iron rain in the magma ocean", *J. Geophys. Res.* **115** (2010), article no. B01404.

[74] H. Samuel, "A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans", *Earth Planet. Sci. Lett.* **313** (2012), p. 105-114.

[75] R. Deguen, M. Landeau, P. Olson, "Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans", *Earth Planet. Sci. Lett.* **391** (2014), p. 274-287.

[76] M. Landeau, R. Deguen, D. Phillips, J. A. Neufeld, V. Lherm, S. B. Dalziel, "Metal-silicate mixing by large Earth-forming impacts", *Earth Planet. Sci. Lett.* **564** (2021), article no. 116888.

[77] V. Lherm, R. Deguen, T. Alboussière, M. Landeau, "Rayleigh-Taylor instability in impact cratering experiments", *J. Fluid Mech.* **937** (2022), article no. A20.

[78] Q. Kriaa, M. Landeau, M. Le Bars, "Influence of planetary rotation on metal-silicate mixing and equilibration in a magma ocean", *Phys. Earth Planet. Inter.* **352** (2024), article no. 107168.

[79] H. J. Melosh, in *Impact Cratering: A Geologic Process* (H. J. Melosh, ed.), Oxford Monographs on Geology and Geophysics, No. 11, Oxford University Press, New York, 1989, p. 253.

[80] K. A. Holsapple, "The scaling of impact processes in planetary sciences", *Annu. Rev. Earth Planet. Sci.* **21** (1993), p. 333-373.

[81] E. Pierazzo, N. Artemieva, E. Asphaug *et al.*, "Validation of numerical codes for impact and explosion cratering: Impacts on strengthless and metal targets", *Meteorit. Planet. Sci.* **43** (2008), no. 12, p. 1917-1938.

[82] H. Nagaoka, S. Takasawa, A. M. Nakamura, K. Sangen, "Degree of impactor fragmentation under collision with a regolith surface-Laboratory impact experiments of rock projectiles", *Meteorit. Planet. Sci.* **49** (2014), no. 1, p. 69-79.

[83] O. G. Engel, "Crater depth in fluid impacts", *J. Appl. Phys.* **37** (1966), no. 4, p. 1798-1808.

[84] V. Lherm, R. Deguen, "Velocity field and cavity dynamics in drop impact experiments", *J. Fluid Mech.* **962** (2023), article no. A21.

[85] H. Lhuissier, C. Sun, D. Lhose, "Drop fragmentation at impact onto a bath of an immiscible liquid", *Phys. Rev. Lett.* **110** (2013), no. 26, article no. 264503.

[86] A. Maller, M. Landeau, L. Allibert, S. Charnoz, "Condition for metal fragmentation during Earth-forming collisions", *Phys. Earth Planet. Inter.* **352** (2024), article no. 107199.

[87] J. B. Wacheul, M. Le Bars, J. Monteux, J. M. Aurnou, "Laboratory experiments on the breakup of liquid metal diapirs", *Earth Planet. Sci. Lett.* **403** (2014), p. 236-245.

[88] J.-B. Wacheul, M. Le Bars, "Experiments on fragmentation and thermo-chemical exchanges during planetary core formation", *Phys. Earth Planet. Inter.* **276** (2018), p. 134-144.

[89] A. M. Worthington, "On impact with a liquid surface", *Proc. R. Soc. Lond.* **34** (1883), no. 220-223, p. 217-230.

[90] A. M. Worthington, R. S. Cole, "V. Impact with a liquid surface, studied by the aid of instantaneous photography", *Philos. Trans. R. Soc. Lond. Ser. A-Contain. Pap. Math. Phys. Character.* **189** (1897), p. 137-148.

[91] A. M. Worthington, R. S. Cole, "IV. Impact with a liquid surface studied by the aid of instantaneous photography. Paper II", *Philos. Trans. R. Soc. Lond. Ser. A-Contain. Pap. Math. Phys. Character.* **194** (1900), no. 252-261, p. 175-199.

[92] O. G. Engel, "Initial pressure, initial flow velocity, and the time dependence of crater depth in fluid impacts", *J. Appl. Phys.* **38** (1967), no. 10, p. 3935-3940.

[93] A. Prosperetti, H. N. Oguz, "The impact of drops on liquid surfaces and the underwater noise of rain", *Annu. Rev. Fluid Mech.* **25** (1993), no. 1, p. 577-602.

[94] J.-L. Liow, "Splash formation by spherical drops", *J. Fluid Mech.* **427** (2001), p. 73-105.

[95] A. I. Fedorchenko, A.-B. Wang, "On some common features of drop impact on liquid surfaces", *Phys. Fluids* **16** (2004), no. 5, p. 1349-1365.

[96] A. Bisighini, G. E. Cossali, C. Tropea, I. V. Roisman, "Crater evolution after the impact of a drop onto a semi-infinite liquid target", *Phys. Rev. E* **82** (2010), no. 3, article no. 036319.

[97] B. Ray, G. Biswas, A. Sharma, "Regimes during liquid drop impact on a liquid pool", *J. Fluid Mech.* **768** (2015), p. 492-523.

[98] V. Lherm, R. Deguen, T. Alboussière, M. Landeau, "Rayleigh-Taylor instability in drop impact experiments", *Phys. Rev. Fluids* **6** (2021), no. 11, article no. 110501.

[99] U. Jain, M. Jalaal, D. Lohse, D. van der Meer, "Deep pool water-impacts of viscous oil droplets", *Soft Matter* **15** (2019), no. 23, p. 4629-4638.
[100] M. Hadj-Achour, N. Rimbert, M. Gradeck, R. Meignen, "Fragmentation of a liquid metal droplet falling in a water pool", *Phys. Fluids* **33** (2021), no. 10, article no. 103315.

[101] M. V. Gielen, P. Sleutel, J. Benschop *et al.*, "Oblique drop impact onto a deep liquid pool", *Phys. Rev. Fluids* **2** (2017), no. 8, article no. 083602.

[102] L. Allibert, M. Landeau, R. Röhlen, A. Maller, M. Nakajima, K. Wünnemann, "Planetary impacts: Scaling of crater depth from subsonic to supersonic conditions", *J. Geophys. Res.: Planets* **128** (2023), no. 8, article no. e2023JE007823.

[103] R. E. Breidenthal, "Sonic eddy-a model for compressible turbulence", *AIAA J.* **30** (1992), no. 1, p. 101-104.

[104] J. B. Freund, S. K. Lele, P. Moin, "Compressibility effects in a turbulent annular mixing layer. Part 1. Turbulence and growth rate", *J. Fluid Mech.* **421** (2000), p. 229-267.

[105] C. Pantano, S. Sarkar, "A study of compressibility effects in the high-speed turbulent shear layer using direct simulation", *J. Fluid Mech.* **451** (2002), no. 1, p. 329-371.

[106] G. L. Brown, A. Roshko, "On density effects and large structure in turbulent mixing layers", *J. Fluid Mech.* **64** (1974), p. 775-816.

[107] E. Pierazzo, H. J. Melosh, "Understanding oblique impacts from experiments, observations, and modeling", *Annu. Rev. Earth Planet. Sci.* **28** (2000), no. 1, p. 141-167.

[108] D. E. Gault, J. A. Wedekind, "Experimental studies of oblique impact", in *Lunar and Planetary Science Conference, 9th, Houston, Tex., March 13-17, 1978, Proceedings. Volume 3.* (A79-39253 16-91) *New York*, vol. 9, Pergamon Press, Inc., 1978, p. 3843-3875.

[109] E. Pierazzo, H. J. Melosh, "Melt production in oblique impacts", *Icarus* **145** (2000), no. 1, p. 252-261.

[110] G. Leneweit, R. Koehler, K. G. Roesner, G. Schäfer, "Regimes of drop morphology in oblique impact on deep fluids", *J. Fluid Mech.* **543** (2005), p. 303-331.

[111] T. Okawa, T. Shiraishi, T. Mori, "Effect of impingement angle on the outcome of single water drop impact onto a plane water surface", *Exp. Fluids* **44** (2008), p. 331-339.

[112] J.-B. Wacheul, R. Deguen, "A laboratory model for oblique planetary impacts" (in preparation).

[113] G. K. Batchelor, "Heat convection and buoyancy effects in fluids", *Q. J. R. Met. Soc.* **80** (1954), p. 339-358.

[114] B. R. Morton, G. Taylor, J. S. Turner, "Turbulent gravitational convection from maintained and instantaneous sources", *Proc. R. Soc. Lond. Ser. A, Math. Phys. Sci.* **234** (1956), p. 1-23.

[115] R. S. Scorer, "Experiments on convection of isolated masses of buoyant fluid", *J. Fluid Mech. Dig. Arch.* **2** (1957), no. 06, p. 583-594.

[116] B. Woodward, "The motion in and around isolated thermals", *Q. J. R. Meteorol. Soc.* **85** (1959), no. 364, p. 144-151.

[117] J. S. Turner, "Buoyant plumes and thermals", *Annu. Rev. Fluid Mech.* **1** (1969), no. 1, p. 29-44.

[118] D. Lecoanet, N. Jeevanjee, "Entrainment in resolved, dry thermals", *J. Atmosph. Sci.* **76** (2019), no. 12, p. 3785-3801.

[119] R. Deguen, M. Lasbleis, "Fluid dynamics of Earth's core: geodynamo, inner core dynamics, core formation", in *Fluid Mechanics of Planets and Stars* (M. Le Bars, D. Lecoanet, eds.), Springer, 2020, p. 129-212.

[120] D. Lecoanet, N. Jeevanjee, "Entrainment in resolved, dry thermals", *J. Atmosph. Sci.* **76** (2019), no. 12, p. 3785-3801.

[121] A. C. Lai, B. Zhao, A. W.-K. Law, E. E. Adams, "A numerical and analytical study of the effect of aspect ratio on the behavior of a round thermal", *Environ. Fluid Mech.* **15** (2015), p. 85-108.

[122] H. Morrison, N. Jeevanjee, D. Lecoanet, J. M. Peters, "What controls the entrainment rate of dry buoyant thermals with varying initial aspect ratio?", *J. Atmos. Sci.* **80** (2023), no. 11, p. 2711-2728.

[123] B. Zhao, A. W. K. Law, A. C. H. Lai, E. E. Adams, "On the internal vorticity and density structures of miscible thermals", *J. Fluid Mech.* **722** (2013), article no. R5.

[124] M. Landeau, R. Deguen, P. Olson, "Experiments on the fragmentation of a buoyant liquid volume in another liquid", *J. Fluid Mech.* **749** (2014), p. 478-518.

[125] L. Huguet, V. Lherm, R. Deguen, T. Le Borgne, J. Heyman, "On the self-similarity of turbulent thermals", *J. Fluid Mech.* (in preparation).

[126] B. A. Ayotte, H. J. Fernando, "The motion of a turbulent thermal in the presence of background rotation", *J. Atmos. Sci.* **51** (1994), no. 13, p. 1989-1994.

[127] K. R. Helfrich, "Thermals with background rotation and stratification", *J. Fluid Mech.* **259** (1994), p. 265-280.

[128] Q. Kriaa, E. Subra, B. Favier, M. Le Bars, "Effects of particle size and background rotation on the settling of particle clouds", *Phys. Rev. Fluids* **7** (2022), no. 12, article no. 124302.

[129] R. Deguen, P. Olson, P. Cardin, "Experiments on turbulent metal-slicate mixing in a magma ocean", *Earth Planet. Sci. Lett.* **310** (2011), no. 3-4, p. 303-313.

[130] J. W. M. Bush, B. A. Thurber, F. Blanchette, "Particle clouds in homogeneous and stratified environments", *J. Fluid Mech.* **489** (2003), p. 29-54.

[131] J.-B. Wacheul, M. Le Bars, "Fall and fragmentation of liquid metal in a viscous fluid", *Phys. Rev. Fluids* **2** (2017), no. 9, article no. 090507.

[132] P. Marmottant, E. Villermaux, "On spray formation", *J. Fluid Mech.* **498** (2004), p. 73-111.
[133] J. O. Hinze, "Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes", *AIChE. J.* **1** (1955), no. 3, p. 289-295.

[134] B. Qaddah, J. Monteux, V. Clesi, M. A. Bouhifd, M. Le Bars, "Dynamics and stability of an iron drop falling in a magma ocean", *Phys. Earth Planet. Inter.* **289** (2019), p. 75-89.

[135] J.-B. Wacheul, *Étude de la fragmentation de métal liquide en chute libre dans un environnement visqueux: application à la formation des planètes*— *Thèses*, fr, Phd thesis, Aix-Marseille, 2019.

[136] M. P. Escudier, T. Maxworthy, "On the motion of turbulent thermals", *J. Fluid Mech.* **61** (1973), no. 03, p. 541-552.

[137] M. Landeau, P. Olson, R. Deguen, B. H. Hirsh, "Core merging and stratification following giant impact", *Nat. Geosci.* **9** (2016), no. 10, p. 786-789.

[138] D. Martin, R. Nokes, "Crystal settling in a vigorously convecting magma chamber", *Nature* **332** (1988), no. 6164, p. 534-536.

[139] D. Martin, R. Nokes, "A fluid-dynamical study of crystal settling in convecting magmas", *J. Petrol.* **30** (1989), no. 6, p. 1471-1500.

[140] G. Lavorel, M. Le Bars, "Sedimentation of particles in a vigorously convecting fluid", *Phys. Rev. E* **80** (2009), no. 4, article no. 046324.

[141] V. S. Solomatov, P. Olson, D. J. Stevenson, "Entrainment from a bed of particles by thermal convection", *Earth Planet. Sci. Lett.* **120** (1993), no. 3-4, p. 387-393.

[142] C. Maas, U. Hansen, "Dynamics of a terrestrial magma ocean under planetary rotation: A study in spherical geometry", *Earth Planet. Sci. Lett.* **513** (2019), p. 81-94.

[143] C. Maas, L. Manske, K. Wünnemann, U. Hansen, "On the fate of impact-delivered metal in a terrestrial magma ocean", *Earth Planet. Sci. Lett.* **554** (2021), article no. 116680.

[144] W. M. Elsasser, "Early history of the Earth", in *Earth Science and Meteorites* (J. Geiss, E. D. Goldberg, eds.), North-Holland Pub. Co., Amsterdam, 1963, p. 113-129.

[145] S. Sasaki, K. Nakazawa, "Metal-silicate fractionation in the growing Earth: Energy source for the terrestrial magma ocean", *J. Geophys. Res.: Solid Earth* **91** (1986), no. B9, p. 9231-9238.

[146] S. Ida, Y. Nakagawa, K. Nakazawa, "The Earth's core formation due to the Rayleigh-Taylor instability", *Icarus* **69** (1987), no. 2, p. 239-248.

[147] J. A. Whitehead Jr, D. S. Luther, "Dynamics of laboratory diapir and plume models", *J. Geophys. Res.* **80** (1975), no. 5, p. 705-717.

[148] J. R. Fleck, C. L. Rains, D. S. Weeraratne et al., "Iron diapirs entrain silicates to the core and initiate thermochemical plumes", *Nat. Commun.* **9** (2018), no. 1, article no. 71.

[149] S. Chandrasekhar, *Hydrodynamic and Hydromagnetic Stability*, Dover Publications, Inc., Mineola, New York, 2013.

[150] D. Bercovici, A. Kelly, "The non-linear initiation of diapirs and plume heads", *Phys. Earth Planet. Inter.* **101** (1997), p. 119-130.

[151] R. Clift, J. R. Grace, M. E. Weber, *Bubbles, Drops, and Particles*, Dover Publications, Inc., Mineola, New York, 2005.

[152] M. Ulvrová, N.oltice, Y. Ricard, S. Labrosse, F. Dubuffet, J. Velímský, O. Šrámek, "Compositional and thermal equilibration of particles, drops, and diapirs in geophysical flows", *Geochem. Geophys. Geosyst.* **12** (2011), no. 10, article no. Q10014.

[153] H. Samuel, P. J. Tackley, M. Evonuk, "Heat partitioning in terrestrial planets during core formation by negative diapirism", *Earth Planet. Sci. Lett.* **290** (2010), no. 1-2, p. 13-19.

[154] P. Olson, D. Weeraratne, "Experiments on metal-silicate plumes and core formation", *Philos. Trans. R. Soc. A: Math. Phys. Eng. Sci.* **366** (2008), no. 1883, p. 4253-4271.

[155] B. D. Marsh, "On the mechanics of igneous diapirism, stoping, and zone melting", *Am. J. Sci.* **282** (1982), no. 6, p. 808-855.

[156] S. Morris, "The effects of a strongly temperature-dependent viscosity on slow flow past a hot sphere", *J. Fluid Mech.* **124** (1982), p. 1-26.

[157] A. Ansari, S. Morris, "The effects of a strongly temperature-dependent viscosity on Stokes's drag law: Experiments and theory", *J. Fluid Mech.* **159** (1985), p. 459-476.

[158] N. M. Ribe, "Diapirism in the Earth's mantle: experiments on the motion of a hot sphere in a fluid with temperature dependent viscosity", *J. Volcanol. Geotherm. Res.* **16** (1983), p. 221-245.

[159] D. J. Stevenson, "Mission to Earth's core—a modest proposal", *Nature* **423** (2003), no. 6937, p. 239-240.

[160] A. M. Rubin, "Dikes vs. diapirs in viscoelastic rock", *Earth Planet. Sci. Lett.* **117** (1993), no. 3-4, p. 653-670.

[161] I. Sumita, Y. Ota, "Experiments on buoyancy-driven crack around the brittle-ductile transition", *Earth Planet. Sci. Lett.* **304** (2011), no. 3-4, p. 337-346.

[162] M. Sato, I. Sumita, "Experiments on gravitational phase separation of binary immiscible fluids", *J. Fluid Mech.* **591** (2007), p. 289-319.

[163] M. C. Shannon, C. B. Agee, "Percolation of core melts at lower mantle conditions", *Science* **280** (1998), no. 5366, p. 1059-1061.
[164] T. Yoshino, M. J. Walter, T. Katsura, "Core formation in planetesimals triggered by permeable flow", *Nature* **422** (2003), no. 6928, p. 154-157.

[165] H. Terasaki, D. J. Frost, D. C. Rubie, F. Langenhorst, "Interconnectivity of Fe–O–S liquid in polycrystalline silicate perovskite at lower mantle conditions", *Phys. Earth Planet. Inter.* **161** (2007), no. 3-4, p. 170-176.

[166] H. Terasaki, D. J. Frost, D. C. Rubie, F. Langenhorst, "Percolative core formation in planetesimals", *Earth Planet. Sci. Lett.* **273** (2008), no. 1-2, p. 132-137.

[167] A. Néri, J. Guignard, M. Monnereau, M. J. Toplis, G. Quitté, "Metal segregation in planetesimals: Constraints from experimentally determined interfacial energies", *Earth Planet. Sci. Lett.* **518** (2019), p. 40-52.

[168] J. W. Hustoft, D. L. Kohlstedt, "Metal-silicate segregation in deforming dunitic rocks", *Geochem. Geophys. Geosyst.* **7** (2006), no. 2, article no. Q02001.

[169] A. Néri, M. Monnereau, J. Guignard, M. Bystricky, C. Tenailleau, B. Duployer, M. J. Toplis, G. Quitté, "Textural evolution of metallic phases in a convecting magma ocean: A 3D microtomography study", *Phys. Earth Planet. Inter.* **319** (2021), article no. 106771.

[170] M. Monnereau, J. Guignard, A. Néri, M. J. Toplis, G. Quitté, "Differentiation time scales of small rocky bodies", *Icarus* **390** (2023), article no. 115294.

[171] O. Šrámek, Y. Ricard, F. Dubuffet, "A multiphase model of core formation", *Geophys. J. Int.* **181** (2010), no. 1, p. 198-220.

[172] O. Šrámek, L. Milelli, Y. Ricard, S. Labrosse, "Thermal evolution and differentiation of planetesimals and planetary embryos", *Icarus* **217** (2012), no. 1, p. 339-354.

[173] T. Lichtenberg, T. Keller, R. F. Katz, G. J. Golabek, T. V. Gerya, "Magma ascent in planetesimals: Control by grain size", *Earth Planet. Sci. Lett.* **507** (2019), p. 154-165.

[174] S. A. Jacobson, D. C. Rubie, J. Hernlund, A. Morbidelli, M. Nakajima, "Formation, stratification, and mixing of the cores of Earth and Venus", *Earth Planet. Sci. Lett.* **174** (2017), p. 375-386.

[175] E. Schmidt, "Versuche zum Wärmeübergang bei natürlicher Konvektion", *Chemie Ingenieur Technik* **28** (1956), no. 3, p. 175-180.

[176] M. Y. Chow, R. G. Akins, "Pseudosteady-state natural convection inside spheres", *J. Heat Transfer* **97** (1975), no. 1, p. 54-59.

[177] J. Hutchins, E. Marschall, "Pseudosteady-state natural convection heat transfer inside spheres", *Int. J. Heat Mass Transfer* **32** (1989), no. 11, p. 2047-2053.