Small-Scale Metal/Silicate Equilibration During Core Formation: The Influence of Stretching Enhanced Diffusion on Mixing

TL;DR

This study explores how turbulent stretching of metal during core formation enhances chemical mixing efficiency, revealing that stretching significantly reduces equilibration times and depends weakly on diffusivity at high Péclet numbers.

Contribution

It introduces scaling laws and numerical simulations demonstrating the impact of stretching enhanced diffusion on metal-silicate equilibration during core formation.

Findings

Stretching reduces equilibration times at high Péclet numbers.

Equilibration time depends weakly on diffusivity in vigorous stretching regimes.

Vigorous stretching leads to a stirring regime with efficient mixing.

Abstract

Geochemical data provide key information on the timing of accretion and on the prevailing physical conditions during core/mantle differentiation. However, their interpretation depends critically on the efficiency of metal/silicate chemical equilibration, which is poorly constrained. Fluid dynamics experiments suggest that, before its fragmentation, a volume of liquid metal falling into a magma ocean undergoes a change of topology from a compact volume of metal toward a collection of sheets and ligaments. We investigate here to what extent the vigorous stretching of the metal phase by the turbulent flow can increase the equilibration efficiency through what is known as stretching enhanced diffusion. We obtain scaling laws giving the equilibration times of sheets and ligaments as functions of a P\'eclet number based on the stretching rate. At large P\'eclet, stretching drastically decreases the equilibration time, which in this limit depends only weakly on the diffusivity. We also perform 2D numerical simulations of the evolution of a volume of metal falling into a magma ocean, from which we identify several equilibration regimes depending on the values of the P\'eclet (Pe), Reynolds (Re), and Bond (Bo) numbers. At large Pe, Re and Bo, the metal phase is vigorously stretched and convoluted in what we call a stirring regime. The equilibration time is found to be independent of viscosity and surface tension and depends weakly on diffusivity. Equilibration is controlled by an efficient thermochemical stretching enhanced diffusion mechanism developing from the mean flow and entraining the surrounding silicate phase.

Tables

Table 1. Table 1: Typical parameters of a metallic impactor’s core falling into a magma ocean. Modified after Solomatov ( 2015 ) and Deguen et al. ( 2011 ) . Dimensionless numbers are calculated for a 100 km impactor’s core.

	Magma ocean
Magma ocean depth, $L$	${10}^6$		m
Gravity, $g$	10		m.s-2
	Silicate phase	Metal phase
Density, ${\rho }_i$	$4\times {10}^3$	$7.8\times {10}^3$	kg.m-3
Thermal expansion, ${\alpha }_i$	$5\times {10}^{-5}$	${10}^{-5}$	K-1
Specific heat capacity, $c_{p_i}$	${10}^3$	$5\times {10}^2$	J.kg-1.K-1
Dynamic viscosity, ${\eta }_i$	${10}^{-1}$	${10}^{-2}$	Pa.s
Diffusion coefficient, ${\kappa }_i^c$	${10}^{-9}$	${10}^{-9}$	m2.s-1
Thermal conductivity, ${\lambda }_i$	4	${10}^2$	W.m-1.K-1
Thermal diffusivity, ${\kappa }_i$	${10}^{-6}$	${10}^{-5}$	m2.s-1
Surface energy, $\gamma$	0.5	0.5	J.m-2
Reynolds number, Re	${10}^{12}$	${10}^{14}$	-
Thermal Péclet number, ${\mathrm{Pe}}_T$	${10}^{14}$	${10}^{12}$	-
Compositional Péclet number, ${\mathrm{Pe}}_C$	${10}^{17}$	${10}^{17}$	-
Thermal Ohnesorge number, OhT	${10}^{-12}$	${10}^{-10}$	-
Compositional Ohnesorge number, OhC	${10}^{-15}$	${10}^{-15}$	-
Bond number, Bo	${10}^{15}$	${10}^{15}$	-
Weber number, We	${10}^{15}$	${10}^{15}$	-

Table 2. Table 2: Summary of the correspondence between thermal and compositional cases. δ i ​ j subscript 𝛿 𝑖 𝑗 \delta_{ij} is the Kronecker symbol. The general scalar field χ i subscript 𝜒 𝑖 \chi_{i} is associated with modified diffusivity K i subscript 𝐾 𝑖 K_{i} and conductivity k i subscript 𝑘 𝑖 k_{i} in each phase i 𝑖 i . q i subscript 𝑞 𝑖 q_{i} is also defined as the ratio between k i subscript 𝑘 𝑖 k_{i} and K i subscript 𝐾 𝑖 K_{i} .

Thermal	$T_i$	${\kappa }_i$	${\lambda }_i$	${\rho }_i{c}_{p_i}$
Compositional	$C_i/D_{1/2}^{{\delta }_{i1}}$	${\kappa }_i^c$	${\kappa }_i^c{D}_{1/2}^{{\delta }_{i1}}$	$D_{1/2}^{{\delta }_{i1}}$
General	${\chi }_i$	$K_i$	$k_i$	$q_i$

Equations

$$t_h\sim \frac{1}{2\dot{ϵ}}\ln \left[2\frac{\dot{ϵ}{s}_0^2}{K}+1\right],$$

$$T_1$$

$${\lambda }_1\nabla T_1\cdot \mathbf{n}$$

$$D_{1/2}$$

$${\kappa }_1^c\nabla C_1\cdot \mathbf{n}$$

$${\chi }_i=T_i\text{or}{\chi }_i=\frac{C_i}{D_{1/2}^{{\delta }_{i1}}},$$

$$\frac{\mathrm{D}{\chi }_i}{\mathrm{D}t}=K_i{\nabla }^2{\chi }_i,$$

$${\chi }_1$$

$$k_1\nabla {\chi }_1\cdot \mathbf{n}$$

$${\rho }_i\left(\frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \nabla \mathbf{u}\right)$$

$$\nabla \cdot \mathbf{u}$$

$$\left[{\sigma }_1-{\sigma }_2\right]\cdot \mathbf{n}=\gamma (\nabla \cdot \mathbf{n})\mathbf{n},$$

$${\mathrm{Re}}_1=\frac{s_0^2\dot{ϵ}}{{\nu }_1},{\mathrm{Pe}}_1=\frac{s_0^2\dot{ϵ}}{K_1},{\mathrm{Oh}}_{{\kappa }_1}=\frac{K_1\sqrt{{\rho }_1}}{\sqrt{\gamma s_0}},\frac{K_1}{K_2},\frac{k_1}{k_2},\frac{{\rho }_1}{{\rho }_2},\frac{{\nu }_1}{{\nu }_2},$$

$${\mathrm{We}}_1=\frac{{\rho }_1{s}_0^3{\dot{ϵ}}^2}{\gamma },$$

$${\mathrm{Re}}_1=\sqrt{\frac{\Delta \rho }{{\rho }_1}\frac{g{s}_0^3}{{\nu }_1^2}},{\mathrm{Pe}}_1=\sqrt{\frac{\Delta \rho }{{\rho }_1}\frac{g{s}_0^3}{K_1^2}},\mathrm{Bo}=\frac{\Delta \rho g{s}_0^2}{\gamma },\frac{K_1}{K_2},\frac{k_1}{k_2},\frac{{\rho }_1}{{\rho }_2},\frac{{\nu }_1}{{\nu }_2}.$$

$$\mathrm{R}{\mathrm{e}}_1={\left(\frac{{\gamma }^3{\rho }_1^2}{\Delta \rho g{\eta }_1^4}\right)}^{1/4}{\mathrm{Bo}}^{3/4}.$$

$$\frac{\partial {\chi }_i(x,t)}{\partial t}+\dot{ϵ}x\frac{\partial {\chi }_i(x,t)}{\partial x}=K_i\frac{{\partial }^2{\chi }_i(x,t)}{\partial x^2}.$$

$$\xi$$

$$\tau$$

$$\frac{\partial {\chi }_i(\xi ,\tau )}{\partial \tau }=\frac{K_i}{K_1}\frac{{\partial }^2{\chi }_i(\xi ,\tau )}{\partial {\xi }^2}.$$

$${\chi }_1(|\xi |\le 1,\tau =0)$$

$${\chi }_2(|\xi |>1,\tau =0)$$

$$\chi (\xi =0,\tau )=\Delta \chi (1+p)\sum _{n=0}^{\infty }(-p{)}^n\mathrm{erf}\left(\frac{2n+1}{4\sqrt{\tau }}\right).$$

$$p=\frac{1-\sqrt{\frac{k_1}{k_2}\frac{q_1}{q_2}}}{1+\sqrt{\frac{k_1}{k_2}\frac{q_1}{q_2}}}.$$

$${\tau }_{1/2}(p)\sim \frac{1}{(p+1{)}^2},$$

$$s(t)=s_0\exp \left(-\dot{ϵ}t\right).$$

$$t_{1/2}=\frac{1}{2\dot{ϵ}}\ln \left[2\mathrm{P}{\mathrm{e}}_1{\tau }_{1/2}(p)+1\right]=\frac{1}{2\dot{ϵ}}\ln \left[\frac{2s_0^2\dot{ϵ}}{K_1}{\tau }_{1/2}(p)+1\right],$$

$$t_{1/2}\sim \frac{1}{2\dot{ϵ}}\ln \left[\mathrm{P}{\mathrm{e}}_1{\tau }_{1/2}(p)\right].$$

$$t_{1/2}\sim \frac{1}{2\dot{ϵ}}\ln \left[\mathrm{P}{\mathrm{e}}_1\frac{1}{(p+1{)}^2}\right]\sim \frac{1}{2\dot{ϵ}}\ln \left[\mathrm{P}{\mathrm{e}}_2{\left(\frac{q_1}{q_2}\right)}^2\right],$$

$$s_{1/2}\sim \frac{1}{\sqrt{{\tau }_{1/2}(p)}}{\left(\frac{K_1}{\dot{ϵ}}\right)}^{1/2},$$

$$s_{1/2}\sim \frac{2}{D_{1/2}}{\left(\frac{{\kappa }_2^c}{\dot{ϵ}}\right)}^{1/2}.$$

$$t_{1/2}=\frac{1}{2\dot{ϵ}}\ln \left[2\mathrm{P}{\mathrm{e}}_1{\tau }_{1/2}(K,k)+1\right].$$

$$s_{1/2}\sim \frac{1}{\sqrt{{\tau }_{1/2}(K,k)}}{\left(\frac{K_1}{\dot{ϵ}}\right)}^{1/2}$$

Full text

\journalname

JGR-Solid Earth

Université de Lyon, UCBL, ENSL, CNRS, Laboratoire de Géologie de Lyon : Terre, Planètes, Environnement, 69622 Villeurbanne, France

\correspondingauthor

Victor Lhermvictor.lherm@ens-lyon.fr

{keypoints}

Stretching enhanced diffusion accelerates equilibration and mixing in stretched sheets and ligaments

Equilibration happens when the metal is stretched down to a scale (a modified Batchelor scale) smaller than $\sim 5$ cm

Impactors’ cores are vigorously stretched and convoluted in a stirring regime

1 Introduction

The initial conditions of the Earth and other terrestrial planets of the Solar System are inherited from their concomitant accretion and differentiation. Terrestrial planets differentiate into an iron core and a silicate mantle, inducing chemical fractionation and heat partitioning between the metal phase migrating toward the core and the surrounding silicate phase (Rubie et al., 2015). The partitioning of heat provided by the impacts and the reduction of the potential gravitational energy sets the initial conditions of the temperature contrast between the core and the mantle (Rubie et al., 2015). This is crucial for the early thermal and magnetic evolution of planetary bodies (Monteux et al., 2011; Williams and Nimmo, 2004), and the formation and evolution of primitives magma oceans (Sun et al., 2018), with in particular the possibility of forming a basal magma ocean (Labrosse et al., 2007). The partitioning of chemical elements has also profound geodynamical implications. For example, the identity and abundance of light and radioactive elements in the core (Badro et al., 2015; Corgne et al., 2007), which are key parameters for the dynamics and evolution of the core, depend on ($P$,$T$,$fO_{2}$) conditions at which metal and silicate have equilibrated.

Geochemical data such as geochronometers, estimates of the mantle and core composition, meteorite geochemistry, and partition coefficients (Li and Agee, 1996; Righter, 2011; Siebert et al., 2011) constitute an inverse problem for the timing of accretion, the ($P$,$T$,$fO_{2}$) segregation conditions and the degree of chemical equilibration of the impactors’ cores with the mantle of the Earth. This inverse problem is under-determined since there are more degree of freedom than constraints. The interpretation of core formation chronometers well illustrates this issue since there is a strong trade-off between the estimated accretion time and the equilibration efficiency: the Hf-W chronometer predicts an accretion time around 30 Myr if perfect equilibration is assumed (Kleine et al., 2002), but that time is significantly increased in case of imperfect equilibration (Rudge et al., 2010).

The interpretation of geochemical data crucially depends on the equilibration efficiency of the segregation mechanisms. The scenario typically used in geochemical core formation models involves the impactor’s core equilibration in a deep magma ocean (cf. table 1 for typical physical parameters) (Rubie et al., 2015). The metal phase is assumed to fragment into centimetric drops, before separating quickly from the silicates as an iron rain and sedimenting to the base of the magma ocean (Ichikawa et al., 2010; Karato and Murthy, 1997; Rubie et al., 2003; Stevenson, 1990). The centimeter scale corresponds to the size at which surface tension prevents further fragmentation of the drops (Rubie et al., 2003). It is small enough to ensure efficient chemical equilibration between the droplets and the surrounding magma ocean (Rubie et al., 2003). Iron finally migrates through the solid part of the mantle toward the forming core by diapirism, initiated by Rayleigh-Taylor instabilities within the liquid pond (Karato and Murthy, 1997; Monteux et al., 2009; Samuel et al., 2010; Stevenson, 1990), diking (Stevenson, 2003) or percolation (Stevenson, 1990). No further equilibration is expected at this stage owing to the solid state of the underlying mantle and the large size of the diapirs.

The iron rain scenario is based on the assumption of rapid fragmentation of the impactor cores into centimetric drops. Whether this is realistic remains an open question. Differentiation of terrestrial planets started early and most of the Earth was accreted from large differentiated planetesimals and embryos (100 to 1000 km) (Rubie et al., 2015). Whether the huge volumes of metal delivered by these impacts were indeed able to fragment into cm-scale drops is unclear (Dahl and Stevenson, 2010; Deguen et al., 2014; Kendall and Melosh, 2016; Landeau et al., 2014; Wacheul et al., 2014). Liquid fragmentation is a well studied problem of fluid mechanics and the key mechanisms are reasonably well understood: though the details depends on the configuration of the flow, the route to fragmentation necessarily involves the formation of elongated liquid ligaments (e.g. Marmottant and Villermaux, 2004a; Villermaux, 2004, 2007), which is the only geometry unstable against the Rayleigh-Plateau capillary instability leading to fragmentation into drops. Fragmentation must therefore be preceded by a change of topology of the metal phase, from a compact volume (the core of the impactor) toward metal ligaments. Ligaments can be produced either directly by hydrodynamic instabilities or turbulence, or indirectly by the bursting of metal sheets, the liquid from the sheets then collecting into ligaments (e.g. Villermaux, 2007). In the context of metal fragmentation in a magma ocean, laboratory experiments and numerical simulations suggest that the production of both metal ligaments and sheets can result from a combination of Rayleigh-Taylor and shear instabilities, turbulent fluctuations, and interactions with the shear related to the crater opening and the following flow (Deguen et al., 2011, 2014; Kendall and Melosh, 2016; Landeau et al., 2014; Wacheul et al., 2014; Wacheul and Le Bars, 2018). Once formed, ligaments can fragment into drops as a result of the Rayleigh-Plateau capillary instability (Eggers and Villermaux, 2008; Marmottant and Villermaux, 2004a, b). Figure 1 illustrates the fragmentation sequence just described, which is to be considered as generic: if fragmentation of planetesimal’s core indeed happens, it corresponds to the typical and necessary sequence of events required for the fragmentation of an initially spherical metallic impactor into stable droplets.

In our conceptual fragmentation model (figure 1) and in laboratory experiments on which it is based (Deguen et al., 2014; Landeau et al., 2014), the liquid metal phase is vigorously stirred and stretched by the turbulent flow following the impact, before fragmentation actually happens. This suggests that mass and heat transfer between metal and silicates may be aided, or even controlled, by stretching enhanced diffusion, a mechanism identified as crucial in the context of the mixing of a stirred diffusive heterogeneity (e.g. Coltice and Schmalzl, 2006; Duplat and Villermaux, 2008; Kellogg and Turcotte, 1987; Olson et al., 1984; Ottino, 1989; Ranz, 1979; Ricard, 2015; Venaille and Sommeria, 2008; Villermaux, 2004). Stretching enhanced diffusion works as follows: take an initially compact patch of a diffusing tracer (temperature or concentration of solute) advected by a flow. Unless the flow is uniform, the patch will be deformed according to the local strain tensor. In 2D, it will be stretched in one of the principal strain directions, and compressed in the other one (thus leading to the formation of 2D sheets). In 3D, it will be either stretched into one of the principal strain direction and compressed in the two others (thus leading to the formation of a filament), or compressed into one direction and elongated in the two others (leading to the formation of sheets). This stretching can drastically accelerate the homogenization of the tracer because the flow has the effect of increasing the exchange surface and maintaining strong concentration gradients in the compression direction(s), thus enabling efficient diffusive transport of the tracer. The effect is important if the stretching time (the inverse of the stretching rate) is smaller than the diffusion timescale based on the initial blob size. For example, if an heterogeneity is stretched into a sheet of thickness $s(t)$ at a constant rate $\dot{\epsilon}=\mathrm{d}(\ln s)/\mathrm{d}t$, the homogenization timescale $t_{h}$ is given by

[TABLE]

where $s_{0}$ the initial thickness of the heterogeneity, and $K$ the tracer diffusivity (e.g. Ranz (1979); see also Kellogg and Turcotte (1987)). If the stretching time $1/\dot{\epsilon}$ is large compared to the diffusion time $s_{0}^{2}/K$ (i.e. if ${\dot{\epsilon}s_{0}^{2}}/{K}\ll 1$), a Taylor expansion of equation 1 shows that $t_{h}\to s_{0}^{2}/K$, which shows that in this limit the homogenization is controlled solely by diffusion. However, if ${\dot{\epsilon}s_{0}^{2}}/{K}\gg 1$, equation 1 shows that $t_{h}$ depends predominantly on the stretching rate $\dot{\epsilon}$, and only weakly (logarithmically) on the tracer diffusivity $K$. Stretching enhanced diffusion is a kinematic theory, and makes no assumption on the nature of the flow responsible for the stretching. Though in deep Earth geodynamics this formalism has been used in the context of low-Reynolds, laminar flows (mixing of compositional heterogeneities in the mantle) (e.g. Coltice and Schmalzl, 2006; Kellogg and Turcotte, 1987; Olson et al., 1984), it has been developed with turbulent flows in mind (Ranz, 1979), and has been used for example to characterize mixing in turbulent jets (e.g. Duplat and Villermaux, 2008; Villermaux, 2004).

Our goal is to test whether stretching enhanced diffusion can be an efficient mean of equilibrating (thermally and chemically) the metal phase of the impactors’ cores with the immiscible molten silicates of the magma ocean, and to generalize its formalism (and equation 1 in particular) to heat and mass transfer between two stirred immiscible liquid phases. After presenting the governing equations and introducing a change of variable allowing to consider heat and mass transfer with the same set of equations (section 2), we build a regime diagram for the deformation of an initially round drop falling under its own weight (section 3). Our numerical simulations and laboratory experiments (Deguen et al., 2014; Landeau et al., 2014; Wacheul et al., 2014; Wacheul and Le Bars, 2018) show that a mass of metal falling into a magma ocean is most probably in a regime which we call stirring regime, in which the metal is stretched into convoluted sheets and ligaments, before fragmenting into drops. We then study in section 4 the equilibration of what we consider to be the “building blocks” of the fragmentation and equilibration sequence, i.e. sheets, ligaments and droplets. In section 5, we finally come back to the stirring regime which we characterize in light of our results on the equilibration of stretched sheets.

2 Equations and Numerical Model

2.1 Geometries

In sections 3, 4, and 5, we will carry out analytical and numerical calculations (cf. table S1 for numerical parameters) on the temperature and concentration evolution of sheets, ligaments and drops. In each geometry (figure 2), we consider a two-phase flow involving a liquid silicate-like outer phase (phase 2) and a liquid metallic-like inner phase (phase 1). Sheets and ligaments are stretched by a constant pure shear flow at a strain rate $\dot{\epsilon}$. We don’t consider the effect of gravity. The two phases have different densities, viscosities, and diffusivities, both for analytical (sections 4.1 and 4.2.1 for sheets and ligaments, respectively) and numerical calculations (supplementary information and section 4.2.2 for sheet and ligaments, respectively). In numerical simulations, ligaments are free to fragment. We also consider the case of a free-falling, initially round metal mass. Numerical calculations on deformed drops assume identical viscosities and diffusivities between the two phases (sections 3, 4.3 and 5). Analytical calculations on non-deformed drops do take into account the conductivity and diffusivity contrasts (section 4.3).

2.2 Governing equations, and equivalence between the heat and mass transfer problems

We study here the evolution of either temperature $T_{i}$ or concentration $C_{i}$ in a two-phase flow, where $i\in\{1,2\}$ refers to the metallic ($i=1$) or the silicate phase ($i=2$). In the Boussinesq approximation and in the absence of volumetric heat sources (radioactive heating or latent heat), the evolution of temperature and composition within both phases is governed by two transport equations of identical mathematical form. Mathematically speaking, the only difference between heat and mass transfer problems lies in the boundary conditions at the metal/silicates interface, which can be written

[TABLE]

for temperature, and

[TABLE]

for composition, where $\lambda_{i}$ is the thermal conductivity, $\kappa_{i}^{c}$ the mass diffusivity, $D_{1/2}$ the partition coefficient and ${\bf n}$ the unit vector normal to the interface pointing toward phase 2. Equations 3 and 5 express the continuity of heat and mass flux at the interface, respectively.

Temperature is continuous across the interface, but concentration is not. Yet the two problems can be made mathematically equivalent by introducing a new variable $\chi_{i}$ defined as either

[TABLE]

where $\delta_{ij}$ is the Kronecker delta. Rather than investigating separately heat and mass transfer, we will thus consider the scalar field $\chi_{i}$, which evolution is governed by the transport equation

[TABLE]

where $\mathrm{D\bullet}/\mathrm{D}t$ is the Lagrangian derivative, with boundary conditions

[TABLE]

at the metal/silicate interface. $\chi_{i}$ has a diffusivity $K_{i}$ (i.e. thermal diffusivity $\kappa_{i}$ or mass diffusivity $\kappa_{i}^{c}$) and a conductivity $k_{i}$ (equal to either $\lambda_{i}$ or $D_{1/2}^{\delta_{i1}}\kappa_{i}^{c}$). We will also make use of the ratio $q_{i}=k_{i}/K_{i}$, which is equal to either $\rho_{i}c_{p_{i}}$ or $D_{1/2}^{\delta_{i1}}$, where $\rho_{i}$ and $c_{p_{i}}$ are respectively density and specific heat capacity of the phase $i$. The correspondence between the material properties of $\chi_{i}$ and their thermal and compositional counterparts are summarized in table 2.

In the following calculations, the velocity field will either be imposed, or obtained by solving the Navier-Stokes equation under the assumption of incompressibility,

[TABLE]

where ${\bf u}$ is the flow velocity, $P$ is pressure, $\eta_{i}$ is viscosity, and $g$ is the acceleration of gravity. The non-linear term on the left-hand side of equation 10 is the source of turbulence. It is included in our numerical simulations, but since our simulations are 2D, we will not be able to reach fully developed turbulent regimes. At the metal/silicates interface, the velocity field and tangential stress are continuous, while the normal stress is discontinuous due to interfacial tension,

[TABLE]

where ${\bf\sigma_{i}}$ is the stress tensor, ${\bf\nabla}\cdot{\bf n}$ is the local curvature, and $\gamma$ is interfacial tension.

2.3 Dimensionless numbers

From these equations, we define a diffusion time $t_{\kappa}=s_{0}^{2}/K_{1}$, an advection time related to stretching $t_{\dot{\epsilon}}=1/\dot{\epsilon}$, a free-fall time related to gravity $t_{g}=\sqrt{(\rho_{1}/\Delta\rho)s_{0}/g}$, a viscous time $t_{\nu}=s_{0}^{2}/\nu_{1}$ and a capillary time $t_{\gamma}=\sqrt{\rho_{1}s_{0}^{3}/\gamma}$, with $s_{0}$ the typical length (initial thickness of the sheet, or diameter of the ligament or drop), $K_{1}$ the thermal or mass diffusivity of phase 1, $\dot{\epsilon}$ the stretching rate, $\rho_{1}$ the density of phase 1, $\Delta\rho$ the density difference, and $\nu_{1}$ the kinematic viscosity of phase 1.

Based on these timescales and the remaining parameters of the set of equations, we build two sets of dimensionless numbers depending on the geometry of the flow we consider. In the case of a stretched sheet or ligament in which we ignore the effect of gravity, using an advection timescale equal to the stretching timescale allows to build the following set of independent dimensionless parameters:

[TABLE]

where Re1 is the Reynolds number, Pe1 the Péclet number, and Oh${}_{\kappa_{1}}$ the thermal/compositional Ohnesorge number. We will also make use of the Weber number defined as

[TABLE]

which is equal to $\left(\mathrm{Oh}_{\kappa_{1}}\mathrm{Pe}_{1}\right)^{2}$. In the case of a free-falling mass of metal, we use an advection timescale equal to the free-fall time and the problem is now described by the following set of independent dimensionless parameters:

[TABLE]

The Ohnesorge number can be obtained from Bo and Pe1 as $\mathrm{Oh}_{\kappa_{1}}=\mathrm{Bo}^{1/2}/\mathrm{Pe}_{1}$.

The Reynolds number defined either as $\mathrm{Re}_{1}=t_{\nu}/t_{\dot{\epsilon}}$ or $\mathrm{Re}_{1}=t_{\nu}/t_{g}$ compares inertia to viscous forces. The Péclet number $\mathrm{Pe}_{1}=t_{\kappa}/t_{\dot{\epsilon}}$ or $\mathrm{Pe}_{1}=t_{\kappa}/t_{g}$ accounts for the relative importance of advective transport compared to diffusive transport. The thermal/compositional Ohnesorge number $\mathrm{Oh}_{\kappa}=t_{\gamma}/t_{\kappa}$ compares the influence of surface tension to thermal or mass diffusion. The number $\mathrm{Oh}_{\kappa_{1}}$ thus defined has the same form as the classical Ohnesorge number ($\nu_{1}\sqrt{\rho_{1}}/\sqrt{\gamma s_{0}}$), except here the $t_{\gamma}$ is compared to the scalar diffusion time rather than to the momentum diffusion time $t_{\nu}$. The Weber number $\mathrm{We}=(t_{\gamma}/t_{\dot{\epsilon}})^{2}$ accounts for the relative importance of fluid’s inertia compared to its interfacial tension, and the Bond number $\mathrm{Bo}=(t_{\gamma}/t_{g})^{2}$ for the relative importance of gravitational forces compared to tension forces.

2.4 Numerical simulations with Basilisk

The equations introduced in section 2.2 are solved in 2D with direct numerical simulations using the free software Basilisk (see basilisk.fr and Popinet (2009)). This partial differential equations solver uses an adaptive Cartesian mesh, making simulations possible up to Reynolds number of $10^{4}$. Mesh refinement is performed when the discretization error of a selected field (e.g. temperature, concentration, velocity) is larger than an arbitrary value. It allows a good spatial resolution when and where it is needed, particularly in the vicinity of the stretched metal phase. On the contrary, spatial resolution is reduced far from the metal-silicate interface, where velocity variations are spatially smoother.

The incompressible Navier-Stokes equations (eq. 10 and 11) are solved with a Courant-Friedrichs-Lewy (CFL) condition limited time step, using a Bell-Collela-Glaz (BCG) advection scheme and a multigrid Poisson-Helmholtz solver for the viscous term. The BCG solver uses a second order upwind scheme. The diffusion term of the advection-diffusion equation (eq. 7) is solved using a time implicit backward Euler discretization. The non-linear advection term of the equation is solved with the BCG scheme.

The two-phase flow is implemented using a level-set method. Each phase is characterized by its physical properties (e.g. diffusivity, viscosity, etc) and the initial value of the scalar and vector fields (e.g. temperature, concentration, velocity). The interface between the two phases is described by a level set function $\phi$ determining its shape. The time evolution of the interface between the phases is computed with a Volume-Of-Fluid (VOF) advection scheme. First, the interface is reconstructed defining lines (2D) in each cell corresponding to the interface. Then, these lines and planes are moved with a geometrical flux computation. Surface tension (eq. 12) is finally implemented by computing the curvature of the interface.

3 The Path Toward Equilibration: Regime Diagram and Timescales

3.1 Dynamical shape regimes

We start here by building a regime diagram (figure 3) for the deformation of an initially round metal volume falling from rest into molten silicates, from 2D numerical simulations. We choose to construct a regime diagram as a function of the Reynolds number in the inner phase Re1 (eq. 15, build with the inertial terminal velocity scaling) and of the Bond number Bo (eq. 15), which measures the relative importance of buoyancy and surface tension. Since we use $\nu_{1}=\nu_{2}$ in these simulations, $\mathrm{Re_{2}}=\sqrt{\rho_{1}/\rho_{2}}~{}\mathrm{Re_{1}}$ and the regime diagram is thus a function of $\mathrm{Re_{2}}$ divided by 1.4. Besides, viscosities in the metal and the magma ocean also vary, with a potential influence on the regime diagram (Wacheul and Le Bars, 2018). Since in the case of a free-falling drop buoyancy is balanced by the largest of inertia and viscous forces, a small value of Bo implies that surface tension always dominates over the other forces, irrespectively of Re1. Conversely, a large value of Bo implies that surface tension is small compared to either viscous forces (low Re1) or inertia (high Re1).

Though based on 2D simulations, our results are broadly consistent with published results from experiments and numerical simulations (in axisymmetric configurations) (e.g. Clift et al., 1978; Han and Tryggvason, 1999; Landeau et al., 2014).

At low Re1, the drop remains undeformed, irrespectively of the value of Bo. The flow around the drop is laminar and stationary, and the drop develops an internal circulation. Keeping Bo small and increasing Re1, the velocity field eventually becomes time-dependent, with the magnitude of the velocity fluctuations increasing with Re1. The velocity fluctuations do deform the drop surface and modify its trajectory, leading to a wobbling regime (Clift et al., 1978), but surface tension still remains large compared to the time-dependent inertial stresses, which prevents any significant stretching of the metal and limits the amount of deformation. At moderate values of Bo ($\sim 10$) and Re1 ($\sim 10-10^{2}$), we identify an ellipsoidal-cap regime (Clift et al., 1978) similar to the spherical regime except for the flattened shape of the droplet.

Deformation of the drop is significant only at large values of Bo and moderate-to-large values of Re1. At $\mathrm{Bo}\gtrsim 10^{2}$, we find that the initially round drop deforms into a backward-facing bag shape for Re1 on the order of 10 (Baumann et al., 1992; Clift et al., 1978; Han and Tryggvason, 1999; Samuel, 2012; Thomson and Newall, 1886), and into a forward-facing bag shape (similar to the skirted mode of Clift et al. (1978) and to the jellyfish mode of Landeau et al. (2014)) for Re1 on the order of $10^{2}$. At even larger values of Re1, we identify a regime which we will refer to as stirring regime where the velocity field consists in two counter-rotating vortices, with smaller scale velocity fluctuations super-imposed. The metal phase is vigorously stretched and deformed into sheets. This regime is qualitatively similar to the immiscible turbulent thermal regime observed experimentally by Deguen et al. (2014) and Landeau et al. (2014), even though the 2D geometry of our simulations does not allow as much turbulence to develop.

From the dimensionless numbers definition of section 2.3, we write Re1 as a function of Bo as

[TABLE]

The factor in between parenthesis is the inverse of a Morton number (e.g. Clift et al., 1978). It only involves material properties in addition to $g$, and typically varies by a factor of 10 mainly due to uncertainties on the viscosity $\eta_{1}$, which is relatively small compared to Re1 and Bo varying by orders of magnitude depending on the size of the metal. Equation 16 (with uncertainties) thus define the part of the regime diagram which is the most relevant for metal segregation in a magma ocean. This is shown in figure 3 (dashed blue lines) for parameters values given in table 1 and a viscosity of the metal phase in the range $10^{-3}-10^{-1}$ Pa.s (Rubie et al., 2003, 2015). If starting from a volume of metal on the large Re and Bo side of the regime diagram (i.e. Re${}_{1}\gtrsim 10^{3}$ and Bo$\gtrsim 10^{2}$, which corresponds to a diameter of $\gtrsim 5$ cm), the falling metal would be expected to be in the stirring regime (turbulent thermal). If fragmentation produces droplets at moderate Re1 and Bo, these droplets will most likely be in either the wobbling or spherical regimes, or possibly in the ellipsoidal-cap regime.

3.2 Equilibration times

We now focus on the thermochemical equilibration efficiency, which we characterize with an equilibration time $t_{1/2}$ corresponding to the time at which the mean scalar field in the metal phase $\langle\chi_{1}\rangle$ is half its initial value. Figure 4a shows the equilibration time $t_{1/2}$ normalized by the free-fall time $t_{g}$ obtained from numerical simulations as a function of the Péclet number in the inner phase $\mathrm{Pe_{1}}$, for various values of Re1 and Bo.

At $\mathrm{Pe_{1}}\lesssim 10^{2}$, we find that $t_{1/2}/t_{g}\sim\mathrm{Pe_{1}}$, irrespectively of Re1 and Bo. This implies $t_{1/2}\sim t_{\kappa}$, which means that equilibration is controlled by diffusion. Increasing Pe1 above $\sim 10^{2}$ results in faster equilibration than in the diffusion regime. At Pe${}_{1}\gtrsim 10^{2}$ the equilibration time depends on the deformation regime of the drop, and hence on Re1 and Bo, with the stirring regime being the most efficient at improving equilibration: we find $t_{1/2}/t_{g}\sim\mathrm{Pe_{1}}^{0.77}$ in the high Pe spherical regime (circles), while in the stirring regime, $t_{1/2}/t_{g}$ shows a weak dependency on Pe1 (diamonds) with $t_{1/2}/t_{g}=0.53\,\mathrm{ln(Pe_{1})}$ (the logarithmic Pe1 dependency will be explained in section 5). The effect of the deformation regime can also be seen from figure 4b which shows the effect on $t_{1/2}/t_{g}$ of increasing Re1 while keeping Bo and Pe1 constants at $10^{4}$. At $\mathrm{Bo}=10^{4}$ the deformation regime changes from spherical to backward-facing bag shape, forward-facing bag shape, and finally stirring regime as Re1 is increased. The equilibration time decreases significantly with Re1 (as $t_{1/2}/t_{g}\sim\mathrm{Re_{1}}^{-0.59}$) while $\mathrm{Re_{1}\lesssim 10^{3}}$, but $t_{1/2}/t_{g}$ seems to reach a plateau at Re${}_{1}\gtrsim 10^{3}$ when entering the stirring regime. While higher Re1 simulations would be needed to confirm this point, this suggests that the equilibration time may become independent of Re1 (and hence on viscosity) in the stirring regime at high Re1. Figure 4c shows $t_{1/2}/t_{g}$ as a function of the Bond number, for Pe1 fixed at $10^{4}$ and Re1 fixed at either 10 or $10^{4}$. We find that $t_{1/2}/t_{g}$ decreases as Bo increases, irrespectively of Re1, before reaching a plateau. This is consistent with the fact that increasing Bo allows the metal phase to be more deformed, which helps equilibration. This is particularly drastic in the stirring regime, in which the vigorous stretching and folding of the metal phase leads to very fast equilibration.

4 Equilibration Processes in the “Building Blocks”: Sheets, Ligaments and Droplets

The regime diagram obtained in section 3, and experiments (Deguen et al., 2014; Landeau et al., 2014; Wacheul et al., 2014; Wacheul and Le Bars, 2018), suggest that the most relevant dynamical regime at the beginning of the fragmentation sequence is the stirring regime, in which the metal phase is stretched into thin sheets (in our 2D simulations), and sheets and ligaments in 3D. The metal phase will then eventually fragment into drops. We therefore study here heat and mass transfer in what we consider to be the three “building blocks” of the fragmentation and equilibration sequence: stretched sheets, ligaments, and drops.

4.1 A stretched isolated sheet

We consider here the homogenization of a scalar field $\chi(x,y,t)$ in and around a stretched isolated 2D sheet. Following Ranz (1979), the evolution of $\chi$ is described in a local Lagrangian frame $(x,y)$ which moves and rotates with the fluid such that the $y-$direction is always parallel to the stretching direction, and the $x-$direction perpendicular to it. In this frame of reference, the flow is a pure shear flow (or stagnation flow), and the velocity field is of the form $(u_{x}=-\dot{\epsilon}x,u_{y}=\dot{\epsilon}y)$, where $\dot{\epsilon}$ is the stretching rate. The thickness of the sheet $s(t)$ is given by $s(t)=s_{0}\mathrm{exp}\left(-\int_{0}^{t}\dot{\epsilon}(t^{\prime})\mathrm{d}t^{\prime}\right)$, and is assumed to be much smaller than the radius of curvature of the sheet. Stretching will tend to align the direction of the gradient of $\chi$ with the compression direction, and we can therefore consider the sheet and the scalar field to be locally invariant in the stretching direction. With this assumption and the velocity field given above, the scalar transport equation 7 becomes

[TABLE]

We transform these equations into diffusion equations with the following change of variable (Ranz, 1979):

[TABLE]

This amounts to normalize lengths by $s(t)$, which sets the thickness of the sheet to $\xi=1$, and time by a diffusion time scale based on the diffusivity inside the sheet (phase 1). With this change of variable, equation 17 reduces to

[TABLE]

The scalar field is initialized with a difference $\Delta\chi$ between the phases such as

[TABLE]

The boundary conditions are the continuity of the scalar field and its flux at the phase interface (equations 8 and 9), and $\chi\rightarrow 0$ at infinity.

The full solution of this set of equations is given in appendix A, based on Lovering (1936). In particular, the scalar in the middle of the sheet $\chi(\xi=0)$ evolves in time according to

[TABLE]

The parameter $p\in[-1,1]$ accounts for the diffusivity contrast in the two-phase flow and depends on the diffusivity ratio $k_{1}/k_{2}$ and on the ratio $q_{1}/q_{2}$ (see table 2) as

[TABLE]

Its value dramatically affects the scalar profiles (figure 5b). If $p<0$ (figure 5a), the diffusivity in the inner phase is larger than in the outer phase, leading to a flatter profile in the inner phase. On the contrary, if $p>0$ (figure 5c), the diffusivity in the outer phase is larger than in the inner phase, leading to a flatter profile in the outer phase.

We compute from equation 23 the normalized equilibration time $\tau_{1/2}$ at which $\chi(\xi=0,\tau_{1/2})=\Delta\chi/2$. Equation 23 shows that $\tau_{1/2}$ is a function of $p$ only, and figure 5d shows that $\tau_{1/2}$ decreases when $p$ increases. An approximate expression for $\tau_{1/2}$ can be found in the limit of $p$ close to -1, corresponding to $k_{1}q_{1}\gg k_{2}q_{2}$, as follows. In this limit, we expect equilibration to be limited by diffusion in the outer phase 2. The equilibration time $t_{1/2}$, corresponding to the normalized equilibration time $\tau_{1/2}$, should therefore be independent of $k_{1}$. According to equation 19, $t_{1/2}$ is a function of $\tau_{1/2}/K_{1}=\tau_{1/2}\,q_{1}/k_{1}$, the form of which depends on the evolution of the thickness of the sheet. $\tau_{1/2}(p)$ must therefore be proportional to $k_{1}$. If $p\rightarrow-1$, $p$ simplifies as $p\sim-1+2\sqrt{\frac{k_{2}}{k_{1}}\frac{q_{2}}{q_{1}}}$ from equation 24. It thus satisfies the proportionality with $k_{1}$ only if

[TABLE]

which in the limit $p\to-1$ tends toward $\sim\frac{1}{4}\frac{k_{1}}{k_{2}}\frac{q_{1}}{q_{2}}$. Figure 5d shows that the prediction of equation 25 with a proportionality factor equal to 0.6 is indeed very close to the full solution of equation 23 when $p$ approaches -1. It is still reasonably accurate at higher values of $p$: it overestimates $\tau_{1/2}$ by a factor 2 at most over the all range of $p$.

In the heat transfer case, the value of $p$ in a magma ocean (table 1) is $p\sim-0.6$. In the composition case, $p=(1-\sqrt{\kappa_{1}^{c}/\kappa_{2}^{c}}D_{1/2})/(1+\sqrt{\kappa_{1}^{c}/\kappa_{2}^{c}}D_{1/2})$ depends on the partition coefficient and the mass diffusivity contrast between metal and silicates (figure 5e). In particular, it can be seen that $p$ approaches $-1$ for siderophile elements (e.g. W, Co, Cr, V, Ni), since $D_{1/2}\gg 1$ and $\kappa_{1}^{c}/\kappa_{2}^{c}$ is larger than unity (O’Neill and Palme, 1998).

We now go back to the stretching of the metal sheet and calculate the equilibration time $t_{1/2}$ from equation 19 taken at $\tau=\tau_{1/2}$. For that purpose, we assume that $\dot{\epsilon}$ is constant, as predicted in the case of homogeneous turbulence (Batchelor, 1952), which implies an exponential decrease of the layer thickness

[TABLE]

Using equation 19, we obtain

[TABLE]

which is a generalization of equation 1 to the homogenization of a scalar field between two liquid phases with different transport properties. In the small Péclet limit, a Taylor expansion of equation 27 shows that $t_{1/2}\simeq\tau_{1/2}(p)\,s_{0}^{2}/K_{1}$, as expected. In the large Péclet limit when advection dominates over diffusion, i.e. in the stretching enhanced diffusion regime, equation 27 simplifies to

[TABLE]

In the limit $p\rightarrow-1$ relevant to siderophile elements (figure 5e), using equation 25 gives

[TABLE]

with $\mathrm{Pe}_{2}$ the Péclet number of the phase 2 defined in the same way as $\mathrm{Pe}_{1}$ (eq. 15).

Figure 6a shows the evolution of $\chi$ in the middle of the sheet as a function of the time normalized by $t_{\kappa}=s_{0}^{2}/K_{1}$, as given by equation 23. If $\mathrm{Pe_{1}}\lesssim 1$, the equilibration time is nearly independent of $\mathrm{Pe_{1}}$, whereas if $\mathrm{Pe_{1}}\gtrsim 1$, $\chi(\xi=0)$ decreases over a timescale significantly reduced as $\mathrm{Pe_{1}}$ increases. Figure 6b shows the equilibration time $t_{1/2}$ as a function of $\mathrm{Pe_{1}}$ according to equation 27. If $\mathrm{Pe_{1}}\ll 1$, $t_{1/2}/t_{\kappa}$ is independent of $\mathrm{Pe_{1}}$ and $t_{1/2}$ scales as the diffusion time $t_{\kappa}$, consistently with the low Péclet number approximation of equation 27. We thus obtain the diffusive regime expected in that range of Péclet number (figure S1). On the contrary, if $\mathrm{Pe_{1}}\gg 1$, $t_{1/2}/t_{\kappa}$ is well approximated by equation 28 with a weaker dependence on $p$ than in the diffusive regime (figure 6b, dash-dotted line). If in addition $p\rightarrow-1$, $t_{1/2}/t_{\kappa}$ scales as equation 29 at large $\mathrm{Pe_{1}}$ (figure 6b, dashed line). $t_{1/2}$ depends predominantly on the stretching rate with a comparatively weak (logarithmic) dependency on the diffusivity. If the stretching rate $\dot{\epsilon}$ increases, the equilibration time $t_{1/2}$ decreases dramatically (eq. 27 and figure 6b). In other words, if stretching intensifies, the equilibration efficiency between metal and silicates increases. In that case, we obtain a stretching enhanced diffusion regime (figure S2). The critical $\mathrm{Pe_{1}}$ separating the two regimes depends on $p$. In particular, when $p\to-1$, i.e. for large partition coefficient and diffusivity ratio $\kappa_{1}^{c}/\kappa_{2}^{c}$, the critical $\mathrm{Pe_{1}}$ is significantly decreased below 1.

Using equations 26 and 28, we can also calculate the width of the sheet $s_{1/2}$ at $t=t_{1/2}$. In the $\mathrm{Pe_{1}}\gg 1$ limit, this gives

[TABLE]

where $(K_{1}/\dot{\epsilon})^{1/2}$ is a characteristic diffusion length, known as the Batchelor length scale (Kellogg and Turcotte, 1987). $s_{1/2}$ corresponds to the length at which diffusion is no longer negligible, the sheet being sufficiently thin for diffusion to complete equilibration. A 100 km metallic core falling at 100 m.s*-1* (Deguen et al., 2011) produces a large scale stretching rate of $\mathrm{100~{}m.s^{-1}/100~{}km=10^{-3}~{}s^{-1}}$. Using a thermal diffusivity of $10^{-5}$ m2.s*-1* (table 1), we obtain a thermal equilibration scale around 6 cm. Concerning chemical equilibration, the $p\rightarrow-1$ limit of $\tau_{1/2}(p)$ (eq. 25) gives

[TABLE]

Using a mass diffusivity of $10^{-9}$ m2.s*-1* and $D_{1/2}$ equal to 1, 10 and 100, the width of chemical equilibration is around 2, 0.2 and 0.02 mm, respectively. Thermal equilibration may happen prior to fragmentation because the typical fragmentation scale corresponding to a capillary length is around 3 mm. However, the chemical equilibration length is of the same order as the fragmentation scale, or even smaller regarding siderophile elements. It means that chemical equilibration may happen at the same time, or after, fragmentation of the metallic core occurs.

4.2 A stretched isolated ligament

4.2.1 Analytical solution

Using the same change of variable as for the sheet (equations 18 and 19) allows to obtain an analytical solution of the axisymmetric advection-diffusion equation of the stretched ligament (figure 2b). The scalar field and its flux are continuous across the phase interface. From Carslaw and Jaeger (1959), p. 346, we obtain the solutions given in appendix B.

Using the same method as for the sheet (eq. 19 and 26), we obtain a law for the equilibration time similar to equation 27, except for $\tau_{1/2}$ which now depends independently on $K$ and $k$ (figure S8), leading to

[TABLE]

We obtain a good agreement between the analytical solution from equation 32 (figure 7a, red circles) and numerical simulations (figure 7a, black circles).

As in section 4.1, we obtain the equilibration length, which can be written

[TABLE]

in the large Péclet limit. For a 100 km metallic core falling at 100 m.s*-1* (Deguen et al., 2011) and using typical parameters of table 1, we obtain $\tau_{1/2}(K,k)\sim 5.5$ and a thermal equilibration length around 4 cm. The chemical equilibration length is around 1.6, 0.5 and 0.1 mm for $D_{1/2}$ equal to 1, 10 and 100, respectively. As for the sheet in section 4.1, thermal equilibration may happen prior to fragmentation whereas chemical equilibration may happen during or after fragmentation.

4.2.2 Effect of surface tension

Ligaments are unstable against the Rayleigh-Plateau capillary instability, and will eventually fragment into drops, possibly affecting equilibration. Fragmentation (figure S5) results from the development of the capillary instability, which grows on a timescale on the order of $t_{\gamma}$ (capillary time based on the ligament initial diameter) in the absence of significant stretching, if $t_{\gamma}\ll t_{\dot{\epsilon}}$. It is known that the capillary instability can be damped if the ligament is stretched (e.g. Eggers and Villermaux, 2008; Mikami et al., 1975; Taylor, 1934; Tomotika, 1936); stretching can delay fragmentation if $t_{\dot{\epsilon}}\ll t_{\gamma}$, i.e. if $\mathrm{We_{1}}\gg 1$.

Figure 7 shows the equilibration time $t_{1/2}$, defined here as the time at which the mean value of $\chi$ of the inner phase $\langle\chi_{1}\rangle$ is half its initial value, obtained from numerical calculations as a function of Pe1, We1, and Oh${}_{\kappa_{1}}$. Fragmentation happens after equilibration if $t_{\gamma}$ is large compared to the (no surface tension) equilibration time given by equation 32. This would be the case if $t_{\gamma}$ is large compared to either $t_{\kappa}\tau_{1/2}$ (at small Péclet) or $t_{\dot{\epsilon}}\ln(\mathrm{Pe_{1}}\tau_{1/2})$ (at large Péclet), i.e. if $\mathrm{Oh}_{\kappa_{1}}/\tau_{1/2}\gg 1$ or $\mathrm{We_{1}}\gg\left[\ln(\mathrm{Pe_{1}}\tau_{1/2})\right]^{2}$. Fragmentation has no effect on equilibration in this limit, and the equilibration time obtained from the numerical calculations is consistent with the prediction of equation 32. The flat slopes at large $\mathrm{We_{1}}$ (figure 7b) and $\mathrm{Oh_{\kappa_{1}}}$ (figure 7c) show that equilibration is indeed independent of surface tension. The diffusive (figure S3) and the stretching enhanced diffusion (figure S4) regimes are thus retrieved at small and large Péclet, respectively (figure 7a, circles).

If now fragmentation happens before the equilibration time predicted by equation 32, we find that equilibration is controlled by the time required to develop the Rayleigh-Plateau capillary instability and break the ligament. At low $\mathrm{We_{1}}$ and $\mathrm{Oh_{\kappa_{1}}}$, we indeed find that $t_{1/2}/t_{\dot{\epsilon}}\sim\mathrm{We_{1}}^{1/2}$ (figure 7b) and $t_{1/2}/t_{\kappa}\sim\mathrm{Oh_{\kappa_{1}}}$ (figure 7c), which correspond to $t_{1/2}\sim t_{\gamma}$. This suggest that the flow associated with the fragmentation of the ligament into drops is strong enough to allow for fast equilibration of the drops.

4.3 A free falling droplet

At the end of the fragmentation sequence of the impactors’ cores (figure 1), the fate of the remaining heat and chemical elements in the droplets depends on the Péclet number, i.e. if they equilibrate in the high Pe spherical or in the diffusive regime. In the diffusive regime (figure S6), the equilibration time scales as the typical diffusion time of the drop $t_{\kappa}=s_{0}^{2}/K_{1}$, irrespectively of viscosity and surface tension (figure 4).

In the high Pe spherical regime (figure S7), equilibration is accelerated by the formation of thin thermal or compositinal boundary layers on the leading side of the drop, on both sides of the interface, where advection maintains a strong radial gradient of $\chi$. At a viscosity ratio $\eta_{1}/\eta_{2}$ on the order of 1 or smaller, equilibration is aided by the circulation forced within the drop by the viscous stress at the interface (Ulvrová et al., 2011). The outer boundary layer is connected to a wake tail on the rear side of the drop. A prediction for the equilibration time is obtained from mass conservation on the spherical drop and estimates of the convective flux (see appendix C for details). We find that

[TABLE]

Equation 36 simplify as $t_{1/2}\sim(s_{0}^{2}D_{1/2})/(6\kappa_{2}^{c})\mathrm{Pe_{2}}^{-1/2}$ and $t_{1/2}\sim(s_{0}^{2}/6)(\kappa_{1}^{c}\kappa_{2}^{c})^{-1/2}\mathrm{Pe_{2}}^{-1/2}$ in the limits of siderophile ($D_{1/2}\gg 1$) and lithophile ($D_{1/2}\ll 1$) elements, respectively. The equilibration time we predict for siderophile elements depends linearly on the partition coefficient, which is consistent with the numerical results of Ulvrová et al. (2011). On the contrary, the equilibration time obtained for lithophile elements is independent of the partition coefficient, consistently with the analysis of Samuel (2012).

5 Stirring Regime

5.1 Deformation and stretching dynamic

We now focus on the stirring regime identified in section 3, which we think is the most relevant for metal segregation in a magma ocean (Dahl and Stevenson, 2010; Deguen et al., 2011, 2014; Landeau et al., 2014).

Figure 8 shows snapshots of a simulation in this regime, at $\mathrm{\mathrm{Re_{1}}}=10^{4}$, $\mathrm{Bo}=10^{3}$, and $\mathrm{\mathrm{Pe_{1}}}=10^{4}$. The initially round drop quickly deforms, both because of an instability which we interpret as a combination of Kelvin-Helmholtz (shear-driven) and Rayleigh-Taylor (buoyancy-driven) instabilities, and because of the action of the mean flow, which takes the form of a pair of expanding counter-rotating vortices (the 2D analogue of a vortex ring). The interaction between the vortices and the interface instability leads to strong deformation and stretching of the metal phase, which topology evolves toward a collection of convoluted sheets (see for example figure 8 at $t/t_{g}=4.2$ and figure 10b). In our 2D numerical calculations, the metal phase eventually becomes discontinuous, but this “fragmentation” is an artifact due to the resolution limits of the computation grid. 2D sheets are stable against Rayleigh-Plateau capillary instabilities and are not expected to break. In the simulations, the metal phase becomes discontinuous when the stretched structures become locally too thin to be resolved, but this is not a physical effect. A finer mesh indeed delays this apparent fragmentation. We therefore use the maximum resolution possible, but 2D simulations are fundamentally not designed to address fragmentation. Nonetheless, these calculations are relevant for understanding the interplay between stretching and thermochemical equilibration, considering scales larger than the resolution limit. The Batchelor scale (e.g. equation 30), estimated using a large scale stretching rate related to the pair of counter-rotating vortices, is indeed two order of magnitude larger than the resolution of the grid.

5.2 Equilibration time

The development of sheets from the initial volume of metal can be considered as the superposition of isolated stretched sheets with uniform stretching similar to that of section 4.1. We thus expect the equilibration time to be of the form $t_{1/2}\sim t_{\dot{\epsilon}}\mathrm{ln(Pe_{1}\tau_{1/2})}$ (eq. 28), where $\dot{\epsilon}$ is the stretching rate of the sheets, which needs to be estimated as a function of the flow geometry and strength. $\tau_{1/2}$ is a function of $p$ in 2D (section 4.1 and figure 5d), or a function of $K_{1}/K_{2}$ and $k_{1}/k_{2}$ in 3D (section 4.2.1), owing to the presence of both sheets and ligaments. The flow has a vortex ring structure (figure 8), with a falling velocity $\langle u_{y}\rangle\sim\sqrt{(\Delta\rho/\rho_{1})gs_{0}}$ varying on a scale $s_{0}$. If stretching is governed by the vortices, then $\dot{\epsilon}$ should be on the order of $\langle u_{y}\rangle/s_{0}$, which is equivalent to a stretching timescale equal to the free-fall timescale of the drop $t_{g}$. This leads to

[TABLE]

There is a good agreement between this prediction and the results from the numerical calculations of section 3, which shows that in the stirring regime $t_{1/2}/t_{g}$ is independent of Re1 (figure 4b, diamonds) and Bo (figure 4c, diamonds), and has a Pe1 dependency (figure 4a, dash-dotted line) well fitted by $t_{1/2}=0.53t_{g}\mathrm{ln(Pe_{1})}$. This suggests that equilibration is controlled by stretching enhanced diffusion, with a stretching rate controlled by the large scale component of the flow.

5.3 The role of stretching on equilibration

To address the role of stretching enhanced diffusion, we focus on the rate of gradient production (Ricard, 2015), obtained by applying the operator $2\nabla\chi\cdot\nabla$ to the advection-diffusion equation (eq. 7) and integrating over the 2D numerical domain $\Omega$

[TABLE]

where $\underline{\dot{\epsilon}}$ is the strain rate tensor. The rate of gradient production (term on the left-hand side) is equal to the sum of an advection term (first term on the right-hand side) and a diffusion term (second term on the right-hand side). The advection term is related to the properties of the flow through the strain rate tensor $\underline{\dot{\epsilon}}$ and corresponds to the source of gradient production. It can be either positive or negative depending on the geometry of the flow and scalar field. In contrast, the diffusion term is always negative and damps scalar gradients.

These integrals are calculated in three examples corresponding to small Pe diffusive, high Pe spherical and stirring regime:

(i) In the small Pe diffusive regime (figure 9a), the advection term is equal to zero and the rate of gradient production is equal to the diffusion term because scalar gradients are only dissipated and never generated. The diffusion term converges toward zero over time owing to the progressive equilibration of the drop, on a timescale $\sim t_{\kappa}$.

(ii) In the high Pe spherical regime (figure 9b), the advection term is positive because scalar gradients are produced by advection in the boundary layer in front of the drop during its fall. Starting from 0, the advection term first increases during a short transient during which the thermal boundary layer grows by diffusion up to a thickness at which advection balances diffusion. The advection term then slowly decrease due to the decrease of difference of $\chi$ between the drop and the surrounding, while the thickness of the thermal boundary layer remains constant. The magnitude of the diffusion term follows closely the evolution of the advection term, but is slightly smaller, which shows that equilibration is controlled by the advection and gradients production. The advection term and the diffusion term converge toward zero as the drop equilibrates with the surrounding fluid, with an equilibration time about 100 times smaller than the diffusion time. In the simulation, the shape of the drop slightly varies in time, explaining the fluctuations around the mean exponential decrease obtained on the integrals.

(iii) In the stirring regime (figure 9c), the advection term becomes positive owing to the production of scalar gradients by the stretching. The advection term is larger than the diffusion term as soon as deformation and stretching of the drop develop, leading to a positive rate of gradient production. Equilibration is then controlled by the production of scalar gradients (whose stirring is the source). At first, the rate of gradient production increases, in connection with the development of stretched structures through hydrodynamics instabilities and stirring. Then, the rate of gradient production decreases owing to the diffusive equilibration of the metal phase with the surrounding fluid when the thickness of the sheets become small enough.

Figure 10 shows the negative eigenvalue of the strain rate tensor $\underline{\dot{\epsilon}}$, which we interpret as a stretching magnitude (Ricard, 2015). In the high Pe spherical regime (a), stretching develops on a scale comparable with the size of the drop, which is significantly larger than the thin boundary layer where thermal and compositional exchanges occur. On the contrary, in the stirring regime (b), stretching occurs in the vicinity of the interface, over a very localized deformation area. The source of gradient production responsible for the large advection term (figure 9c, upper) is then located at the same place as the interface deformation. Small-scale stretching processes related to the stirring are responsible for gradient production, arguing in turn for an equilibration controlled by stretching enhanced diffusion.

5.4 Probability Distribution Functions of $\chi$

The qualitative difference between the three regimes is also highlighted by the probability distribution functions (PDF) of the scalar field. In the small Pe diffusive regime (figure 11a), the distribution of $\chi$ tends toward a probability density $\propto 1/\chi$. This is consistent with the radial dependency of $\chi$, which, after the initially discontinuous radial profile of $\chi$ has smoothed out, tends toward a self-similar profile of the form $\chi(r,t)\propto\exp(-r^{2}/(K_{2}t))$. It can be verified that the PDF of a field with this spatial distribution is indeed $\propto 1/\chi$. In the high Pe spherical regime (figure 11b), the PDF slowly shifts toward smaller scalar values over time, conserving its shape consistently with the stable evolution of the scalar field inside the drop during its fall. In the stirring regime (figure 11c), $\chi$ has a peaked distribution, which over time becomes narrower and shifts toward lower scalar levels as a result of the progressive homogenization.

It has been argued (Duplat et al., 2010; Meunier and Villermaux, 2003; Villermaux, 2004) that the interaction between stretched sheets and ligaments lead to the addition of their scalar fields and to the stable self-convolution of their scalar distributions. It is predicted that in this regime the scalar field distribution follows a one parameter gamma distribution,

[TABLE]

where $n=1/\sigma^{2}$ is related to the standard deviation $\sigma$ of the distribution. After a time corresponding to the development of stretched structures, the PDFs from our calculations are indeed well approximated by these one parameters gamma distributions (figure 11c), except in the tail where the probability density is higher. This again is consistent with an equilibration controlled by stretching enhanced diffusion. Self-convolution as a result of scalar fields addition also supports the assumption on the superposition of isolated sheets we made to predict the equilibration time of the drop (eq. 37).

6 Conclusion

Deformation and fragmentation of differentiated impactors’ cores are key mechanisms for heat and chemical elements transfer during the differentiation of planetary bodies. We argue that fragmentation of the metallic core of an impactor must involves a change of topology leading to a collection of sheets and ligaments, which can then fragment into drops (figure 1). This picture is consistent with previous experimental results (Deguen et al., 2011, 2014; Kendall and Melosh, 2016; Landeau et al., 2014; Wacheul et al., 2014; Wacheul and Le Bars, 2018) and our numerical simulations (section 5.1), which suggest that a large volume of metal falling into a magma ocean would be vigorously stirred and stretched. We focus here on the physical concepts associated with the fragmentation and equilibration sequence. The “building blocks” of the fragmentation sequence involve stretched structures (sheets and ligaments) in which stretching enhanced diffusion, generalized here to immiscible fluids, improves significantly chemical and thermal equilibration efficiency between metal and silicates. The equilibration time depends mainly on the stretching rate with, if the stretching rate is constant, a logarithmic dependency on diffusivity and partition coefficient. We find that thermal equilibration may occur before fragmentation whereas chemical equilibration may still be incomplete at this moment. Fragmentation eventually produces droplets whose chemical equilibration depends critically on the partition coefficient (equation 36).

In the stirring regime identified for large volume of metals (in the sense of having large Reynolds and Bond numbers), the metal phase is vigorously stretched and deformed. We find from numerical simulations in that regime that the equilibration time depends mainly on the stretching rate of the metal phase and weakly on diffusivity, irrespectively of viscosity and surface tension. Our results suggest that thermochemical equilibration is controlled by the stretching enhanced diffusion process, which in our numerical simulations develops from the mean flow, as the fall of the drop entrains surrounding silicates. We find that in our numerical calculations the equilibration time is well approximated by $t_{1/2}\sim\sqrt{\frac{\rho_{1}}{\Delta\rho}\frac{s_{0}}{g}}\mathrm{ln\,Pe_{1}}$ (equation 37), where the Péclet number is built on a stretching rate equal to $\langle u_{y}\rangle/s_{0}=\sqrt{(\Delta\rho/\rho_{1})gs_{0}}/s_{0}$. Using this scaling would predict relatively large equilibration distances (50% of the initial temperature/composition difference) in a magma ocean: as an example, the composition and temperature equilibration distances of a 1 km impactor’s core sinking in a magma ocean would be 12 km and 8 km, respectively. Equation 37 could also be used to predict that only impactors with a core smaller than 56 km and 41 km will significantly equilibrate thermally or chemically in a 3000 km deep magma ocean.

However, these predictions are based on the assumption that the stretching rate is governed by the mean flow (eq. 37). This seems to be correct in our 2D calculations, but is unlikely to hold in 3D at higher Re and Bo due to the development of turbulent fluctuations. Turbulent velocity fluctuations would increase the stretching of the metal phase, and thus decrease the equilibration time. The prediction of equation 37 therefore most probably overestimates significantly the equilibration time. We may anticipate that in the limit of high Re1, Bo and Pe1, equilibration of the metal with the surrounding silicates would even be fast enough for the evolution of the mean value of $\chi$ in the metal phase to be limited by the rate of entrainement of fresh silicates rather than by equilibration at the local scale. Our analysis suggests that the key for testing this hypothesis and estimating the equilibration times is to obtain models predicting the stretching rate of the metal phase during the cratering phase and post-cratering flow.

Appendix A Solution of the Diffusion Equation for the Stretched Sheet

The solution of the diffusion equation 20 for a stretched isolated sheet is given in each phase by

[TABLE]

Appendix B Solution of the Diffusion Equation for the Stretched Ligament

The solution of the diffusion equation for a stretched isolated ligament is given in each phase by

[TABLE]

where $K=\sqrt{K_{1}/K_{2}}$, $k=\sqrt{k_{1}/k_{2}}$ and

[TABLE]

$J_{\nu}$ and $Y_{\nu}$ are respectively the Bessel function of the first and the second kind, with $\nu$ their order.

Appendix C Equilibration Time for the Free-Falling Droplet

We consider here the equilibration of a drop of diameter $s_{0}$ and falling velocity $u$, with internal circulation, at high Pe1 and Pe2. In this limit, we expect the formation of two thin thermal/compositional boundary layers, on both sides of the droplet interface, with respective thicknesses $\delta_{1}$ and $\delta_{2}$. Depending on the properties of the two phases, the flux can be limited by one or the other of the boundaries. Denoting by $\left<\chi_{1}\right>$ the mean of $\chi_{1}$ in the drop, by $\chi_{i}$ the value of $\chi$ at the interface, and by $\chi_{2}^{\infty}$ the value of $\chi_{2}$ far from the drop, the heat/mass flux $\Phi$ across the interface is

[TABLE]

The boundary layers thickness on both sides of the interface can be shown to be $\delta_{i}\sim\sqrt{K_{i}s_{0}/u}\sim s_{0}\,\mathrm{Pe}_{i}^{-1/2}$ by balancing advection and diffusion in the vicinity of the interface (see for example Ribe (2015) for the outer boundary layer scaling; a similar reasoning leads to the scaling of the inner boundary layer, which has the same form). Writing $\chi_{i}$ as a function of $\left<\chi_{1}\right>$ and $\chi_{2}^{\infty}$ from equation 46, and noting that $\delta_{1}/\delta_{2}\sim\sqrt{K_{1}/K_{2}}$, the flux becomes

[TABLE]

which is consistent with the heat flux given by Ribe (2015) in the limit $q_{1}k_{1}\gg q_{2}k_{2}$. Using conservation of $\chi$ gives

[TABLE]

from which we obtain the equilibration time

[TABLE]

Acknowledgements.

This project has received funding from the European Research Council (ERC) under the European Unions Horizon 2020 research and innovation programme (grant agreement No 716429). T. Alboussière provided comments that improved the manuscript. We thank Jean-Baptiste Wacheul and an anonymous reviewer for their useful comments. The codes are available at https://figshare.com/articles/lherm_basilisk_stretching/6850271 and basilisk.fr.