<figure><img src="image_76.png" /><figcaption>**Fig. 3.** Evolution of the Earth's mantle Nb/Ta ratio during the planet growth as a function of oxygen fugacity (after Cartier et al., 2014). The model considers the accumulation of impactors, starting with reduced and small-sized bodies, followed by less reduced and larger objects (Rubie et al., 2011). This accretion model also considers the degree of mixing (core–mantle equilibration) after each impact. The first stages are characterized by very low Nb/Ta ratios, due to massive core uptake. The last bodies accreted are more oxidized and their metallic parts do not fully equilibrate with the mantle, yielding to an increase of the Nb/Ta ratio, up to the current estimated value of 14. The arrow indicates the corresponding mass fraction accumulated as accretion progresses.</figcaption></figure>

therefore a good $fO_2$ proxy during core formation. Mars and Vesta-4, for example, are considered to have differentiated under only moderately reducing conditions (IW and IW-1.5, respectively). This explains why SNC and HED meteorites display an almost chondritic Nb/Ta.

In summary, realistic Earth's core formation models are now in our reach that include a combination of a multistage core-mantle differentiation with N-body accretion simulations, or a continuous Earth's core formation where pressure, temperature and oxygen fugacity vary as the Earth's accretion proceed (e.g. Rubie et al., 2015a). In these models, we consider that the accretion of the Earth occurred during a series of large impact events (e.g. Canup, 2008; Wetherill, 1985; O'Brien et al., 2006; Montieux et al., 2009; Nakajima and Stevenson, 2015). The pressure and temperature of metal-silicate equilibrium likely increased as accretion progressed and core formation was therefore a continuous process. Other geochemical arguments based on $^{107}$Pd-$^{107}$Ag and $^{182}$Hf-$^{182}$W (e.g. Halliday and Wood, 2009) provide support for continuous core formation. In addition, Wade and Wood (2005) completed the deep magma ocean model by presenting a model in which the magma ocean deepened in proportion to the size of the proto-Earth, and in which the Earth became more oxidized at the last accretion stages during the first 100 Ma. Results of both models suggest that accretion of Earth was heterogeneous with the initial 60–70% of Earth's mass forming from reducing material and the final 30–40% from relatively oxidized material. These models also give results for metal-silicate partitioning of a significant number of siderophile elements that are consistent with observed mantle abundances provided the oxygen fugacity increases from about IW-5 to IW-2 (5 to 2 log units lower than the Iron-Wüstite buffer) during core separation. This conclusion is consistent with N-body simulations of terrestrial planet formation, which show that the Earth accreted material from farther from the sun at later times (e.g. Raymond et al., 2006; Morbidelli et al., 2012; Fischer and Ciela, 2014; Rubie et al., 2015a). Consequently, Earth began to accrete with a low FeO content in its mantle (in the range of 0.2–0.5 wt% consistent with the reduced chondritic bodies as enstatite chondrites), and at the final stages of accretion we can assume that the Earth's core is 85 wt% Fe and the mantle about 8 wt% FeO (McDonough and Sun, 1995). Finally, for Earth the reducing material can be similar to enstatite chondrites as it is consistent with the $^{17}$O, $^{48}$Ca, $^{50}$Ti, $^{62}$Ni and $^{90}$Mo isotopic study by Dauphas et al. (2014). Enstatite chondrites being the main building blocks of the Earth are also consistent with a recent collisional erosion model (Boujibar et al., 2015) or the separation of the core from a hydrous magma ocean (Clesi et al., 2016).

Other model suggests that the Earth's mantle can start in a highly oxidized state and becomes more reduced as FeO enters the core (Rubie et al., 2004; Siebert et al., 2013). However this type of model is only valid when no Si partitions into the Earth's core. Recent experiments show in fact the simultaneous incorporation of Si and O at the Earth's core formation conditions (Bouhifd and Jephcoat, 2011; Siebert et al., 2012; Tsuno et al., 2013; Fischer et al., 2015). Every mole of silicon that enters into the core releases two moles of FeO that partition into the mantle following the reaction (1), as proposed by Javoy (1995) and Javoy et al. (2010) to argue for an enstatite chondrite model for the Earth:

$$2\text{Fe}_{\text{(metal)}} + \text{SiO}_{2} \text{(silicate)} = \text{Si}_{\text{(metal)}} + 2\text{FeO}_{\text{(silicate)}} \quad (1)$$

Consequently the initially oxidized mantle cannot be reduced to its present state when mass balance is taken into account (e.g. Rubie et al., 2015a).

### 3. Incorporation of highly volatiles in an early magma ocean

Highly volatile elements are defined as elements with low 50% condensation temperature (lower than about 800 K at $10^{-4}$ atm – McDonough and Sun, 1995) and can include hydrogen, carbon, nitrogen, sulfur and noble gases. In the following we will only focus on helium, and how the solubility and partitioning of this element during Earth's core formation places strong constraints on the formation and evolution of our planet. For the other elements their behaviour during the earliest time of our planet can be found in recent papers – even if the origins and abundances of these elements are still much debated – (e.g. Albarède, 2009; Armytage et al., 2013; Bergin et al., 2015; Boujibar et al., 2014; Dasgupta et al., 2013; Gaillard and Scaillet, 2014; Gaillard et al., 2015; Halliday, 2013; Hirschmann, 2016; Hirschmann and Dasgupta, 2009; Hirschmann et al., 2012; Jephcoat et al., 2008; Kadik et al., 2015; Li et al., 2016; Marty, 2012; Marty et al., 2013; Moreira, 2013; Mukhopadhyay, 2012; Roskosz et al., 2013; Wood et al., 2013; and references therein).

#### 3.1. Helium in early magma ocean

To simulate Earth's core formation under conditions of segregation from a deep magma ocean, helium partitioning was determined in a mixture of a molten CI-chondrite model composition and Fe-Ni metal alloy liquid (Bouhifd et al., 2013). The results presented in this study show that the concentrations of He in a CI-chondrite melt increase with increasing pressure up to ~3 GPa, remain roughly constant in the range of 3 to 8 GPa and drop from 2000 ppm to about 500 ppm and again remain constant to the maximum pressure reached of that work (of ~16 GPa). One can note that the shape of He solubility in CI-chondrite melt versus pressure is similar to that observed for Ar with similar silicate composition (Bouhifd and Jephcoat, 2006). We reported in Fig. 4a–b the Ar and He solubilities in CI-chondrite melt at high pressures. Both studies indicate that below 8 GPa (pressure corresponding to a drop of Ar and He solubilities in CI-chondrite melt and corresponding to depths <300 km in the Earth's mantle), Ar and He may dissolve in silicate melts formed by partial melting in the Earth's interior. In a study up to 8 GPa, the partition coefficient of Ar between clinopyroxene and its melt was found to be between $10^{-3}$ and $10^{-4}$, implying that He and Ar would be partitioned in melt and also ultimately degassed (Brooker et al., 2003). This behaviour indicates that Ar and He maintain their geochemically assumed incompatible behaviour with crystalline phases. Similar conclusion can be drawn from experiments at 0.1 GPa where all noble gases show an incompatible behaviour for olivine and clinopyroxene (Heber et al., 2007). However, the high-pressure step decrease of Ar and He in CI-chondrite melt observed at about 8 GPa implies that there may be an increase in compatibility with the crystalline phase (assuming that the crystalline phase does not undergo similar drop in solubility at high pressures). At depths for pressures >8 GPa (~300 km
<figure><img src="image_77.png" /><figcaption>**Fig. 4.** (a) Helium concentrations (ppm) in CI-chondrite melt as a function of pressure (GPa). Figure taken from Bouhifd et al. (2013). (b) Argon concentrations (ppm) in CI-chondrite melt as a function of pressure (GPa). Data are taken from Bouhifd and Jephcoat (2006), Bouhifd et al. (2008) and White et al. (1989). For both figures error bars correspond to 1σ (ranging from 7 to 25%) of the average of the chemical analysis and to 10% of the nominal pressures reported. The dashed lines are used as a guide line only.</figcaption></figure>

in the Earth's mantle) the Ar and He solubility in CI-chondrite melt indicates that it would be possible for primordial and radiogenic Ar and He to be preferentially retained deep in the Earth. One can note here that the solubility of Ar in silicate melts at high pressure is an ongoing debate as the solubility drop of Ar in SiO₂ melt at about 5 GPa (Chamorro-Perez et al., 1996; Bouhifd et al., 2008) was not reproduced in a recent study by Niwa et al. (2013), who instead show a slight decrease of Ar solubility in the range of 5–12 GPa. For basaltic melt, molecular dynamic simulations show also a slight decrease of Ar solubility in the range of 5 to 20 GPa (Guillot and Sator, 2012), but such calculations present large uncertainties. In any case, for silicate melts, it is widely accepted that noble gases dissolve into the interstitial sites, meaning that their solubility is related to the concentration of holes able to accommodate atoms of particular size (Doremus, 1966; Studt et al., 1970; Shackelford et al., 1972; Shelby, 1976). Increasing pressure causes a continuous densification of melts leading to changes in the silicate units/species distribution within the silicate melt. In their molecular dynamics simulations of molten anorthite, Nevins and Spera (1998) show that near 5 GPa the activation volumes of all species present in the melt are higher in absolute value than the partial molar volume of Ar, meaning that Ar can still be accommodated in molten anorthite at this pressure. However, above 20 GPa the activation volumes of all melt species drop below the partial molar volume of argon. This pressure is very consistent with the onset pressure of Ar solubility drop for anorthite liquid of about 17 GPa (Bouhifd and Jephcoat, 2006). Thus, a first-order interpretation of the drop in Ar solubility data is that they represent a closing of a characteristic void space in melts at high pressures and temperatures.

## 3.2. Earth's core as a reservoir for primordial helium

The question now is whether the Earth's core could plausibly be a reservoir for helium, knowing that a large impact that formed the Moon may have depleted the mantle as well as the atmosphere from the noble gases by about 2 orders of magnitude. In fact, Xe isotopes indicate that extensive loss (>97%) occurred from the mantle and atmosphere at two episodes of major mantle degassing, presumably driven by giant impacts, respectively at 20–50 Ma and about 100 Ma after the formation of the solar system which is consistent with the timing of Moon formation derived from Hf-W systematics and simulations of the giant impact (Pepin and Porcelli, 2006; Porcelli et al., 2001). Helium concentration in the early Earth's core was controlled by the initial He abundance in the early mantle at the time of core formation. High mantle concentrations of noble gases may have existed in the early mantle even if the mechanism by which early Earth acquired its noble gases is still debated. Noble gases in the early Earth's mantle can have two main origins: (1) the first origin/mechanism is the dissolution of gases from a dense, gravitationally captured atmosphere (assuming the Earth reached much of its present mass before dissipation of the solar nebula) and convected into the mantle. In this scenario, noble gases of solar composition, in particular solar-type helium and neon are dissolved from the atmosphere (e.g. Harper and Jacobsen, 1996; Honda et al., 1991; Mizuno et al., 1980). (2) The second scenario is that mantle noble gases are already present in the parent bodies of the Earth, due to irradiation by the solar wind of refractory grains (e.g. Moreira, 2013; Moreira and Charnoz, 2016; Porcelli and Halliday, 2001; and references therein). Such mechanisms can provide at a minimum $7 \times 10^{12}$ atoms $^3$He g$^{-1}$ in the mantle (Mizuno et al., 1980; Porcelli and Ballentine, 2002; Porcelli and Halliday, 2001; Trieloff and Kunz, 2005). Considering this value as the initial $^3$He abundance during the Earth's core formation, and the measured molten metal – silicate liquid partition coefficients of He ($D_{He}$) that range from about $5 \times 10^{-3}$ to $2 \times 10^{-2}$ at high pressure and high temperature (Bouhifd et al., 2013; Zhang and Yin, 2012), the early Earth's core would contain between $3 \times 10^{10}$ and $1.2 \times 10^{11}$ $^3$He atoms g$^{-1}$. 1. On the other hand, the $^3$He flux from Oceanic Island Basalts (OIB) ranges from about 38 to 1000 mol yr$^{-1}$ equivalent to values from about $2.3 \times 10^{25}$ to ~6 $\times 10^{26}$ atoms $^3$He yr$^{-1}$ (e.g. Lupton and Craig, 1975; Porcelli and Ballentine, 2002). The ratio of $^3$He flux from OIB to the concentrations of $^3$He in the Earth's core requires that 0.8 to 3% of the mass of the core should have released its primordial helium (Bouhifd et al., 2013). The range of $^3$He/$^4$He ratios observed in the present day OIB samples would then be the result of mixing from the core with radiogenic $^4$He and residual $^3$He from the mantle. In addition, the early core would probably have formed with a high He/U ratio because U is a highly lithophile element (Wheeler et al., 2006; Malavergne et al., 2007; Bouhifd et al., 2013), and the core probably maintained a high, throughout Earth's history, $^3$He/$^4$He ratio. The mantle reservoir debate (at least for helium – more similar work is needed for other noble gases) must now include exchange processes that operate at the core-mantle boundary, especially when molten metal – silicate melt reactions are involved (e.g. Herzberg et al., 2013), and knowing, as reported by Macpherson et al. (1998) that high $^3$He/$^4$He coincides with a domain of ULZV usually associated with partial melting of silicates directly above a molten core. Finally, even if the core helium budget is actually sufficient to provide all mantle $^3$He, it is reasonable to expect that only OIBs originating from the core-mantle boundary can carry the high $^3$He-$^4$He core signature, making the Earth's core a plausible reservoir for these high $^3$He-$^4$He ratios.

## 4. Magma ocean crystallization after core formation

After the core-mantle separation, a second major chemical differentiation event occurred during the magma ocean crystallization. Indeed, during the cooling and the subsequent solidification of the magma
ocean, compatible elements (e.g. Mg, Cr) were preferentially collected in the solid phase while the incompatible elements (e.g. Al, Na, Fe) selectively partitioned into melts. Depending on the density contrast between the solid and the ambient liquid, a major chemical separation occurred within the early mantle. This differentiation event was strongly influenced by the cooling processes that governed the rheology, the melt fraction and hence the efficiency of chemical exchanges within the solidifying early mantle. As the melt fraction depends on the distance between the solidus and liquidus of the mantle, the melting curves of the primitive mantle have major consequences on the cooling dynamics of the early mantle, on the existence of chemical heterogeneities and on their lifetime within primitive mantle reservoirs. Until recently, the determination of melting properties of samples was rather limited in pressure. Now systematic experimental determinations of the solidus and liquidus of geological materials are determined up to the core-mantle boundary (CMB) using the laser-heated diamond anvil cell technique (LHDAC) (e.g. Fiquet et al., 2010; Andrault et al., 2011, 2014; Nomura et al., 2014; Tateno et al., 2014; Pradhan et al., 2015). In Fig. 5 we report the liquidus and solidus of chondritic mantle (the details about the determination of the liquidus and solidus of chondritic material – that are beyond the scope of the present work – are reported in Andrault et al., 2011).

In vigorously convecting systems such as magma oceans, the temperature distribution is nearly adiabatic and isentropic and the crystallization is supposed to start as soon as the adiabatic profile is crossing the liquidus (e.g. Solomatov, 2007). In one-phase systems, such as a completely molten or a completely solid layer, the equation for an adiabat is

$$ \left(\frac{\partial T}{\partial r}\right)_S = -\frac{\alpha g}{C_p} T \qquad (4) $$

where T is the temperature, r is radius and g is the gravitational acceleration. ρ, α and C_p are respectively the density, the coefficient of volumetric thermal expansion and the specific heat of the considered one phase material at a given depth.

In two-phase systems, the effects of phase changes need to be considered (Solomatov, 2007). For a zone of partial melting, the density ρ', the coefficient of volumetric thermal expansion α' and the specific heat C_p of the partially molten material are given as follows:

$$ \frac{1}{\rho'} = \frac{1-\phi}{\rho_s} + \frac{\phi}{\rho_m} \qquad (5) $$

<figure><img src="image_78.png" /><figcaption>Fig. 5. Adiabatic profiles (black solid lines computed from Thomas and Asimow, 2013), solidus (green) and liquidus (red) (from Andrault et al., 2011) considering a chondritic composition for the early mantle. In this context, as the adiabats are steeper than the liquidus in the whole mantle, an initially fully molten mantle will start to solidify from the core-mantle boundary toward the Earth's surface. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

$$ \alpha' = \alpha + \frac{\Delta\rho}{\rho'(T_{liq} - T_{sol})} \qquad (6) $$

$$ C_p = C_p + \frac{\Delta H}{T_{liq} - T_{sol}} \qquad (7) $$

where φ is the melt fraction at a given depth, ρ_s is the density of the solid phase, ρ_m is the density of the molten phase, Δρ is the density difference between solid and liquid. T_{liq} and T_{sol} are the liquidus and solidus temperature at a given depth. ΔH is the latent heat released during solidification. The equation for such adiabats is given by:

$$ \left(\frac{\partial T}{\partial r}\right)_S = -\frac{\alpha' g}{C_p} T \qquad (8) $$

The thermodynamical parameters for the molten magma ocean are closely related to its chemical composition. Volume and elastic parameters of silicate liquids have been recently characterized up to 140 GPa using shock compression experiments (Mosenfelder et al., 2007, 2009; Thomas et al., 2012; Thomas and Asimow, 2013). In Fig. 5 we also reported the adiabats computed from Thomas and Asimow (2013) considering a chondritic composition for the early mantle. The adiabatic temperature profiles are calculated by numerical integration of Eqs. (1) and (5) using a fourth-order Runge-Kutta method (Press et al., 1993). We assumed here a multi-component system with a chondritic composition (62% enstatite + 24% forsterite + 8% fayalite + 4% anorthite + 2% diopside). Using fourth-order Birch-Murnaghan/Mie-Grüneisen equation of state fits for molten silicate liquids from Thomas and Asimow (2013), we obtained the melt density ρ_m, the volumetric thermal expansion α as a function of pressure as well as the specific heat C_p of the molten material for this multicomponent assemblage (see Monteux et al., 2016 for more details). Fig. 5 illustrates that a chondritic two-phase adiabat is steeper than the purely liquid or solid one-phase chondritic adiabats (Solomatov, 2007).

Depending on the slope difference between the adiabats and the liquidus, the crystallization will occur either from the bottom toward the surface or from the top toward the bottom (e.g. Solomatov, 2007; Labrosse et al., 2007; Monteux et al., 2016). Hence, determining precisely the adiabatic profiles within the magma ocean is also a key aspect to constrain the thermo-chemical evolution of the early terrestrial mantle (Thomas and Asimow, 2013). When the melt fraction of the magma ocean decreases down to a critical value of 40% the cooling dynamics are not governed by the liquid phase but by the solid phase (Abe, 1997). Upon cooling, both the solid and liquid compositions will diverge from the early mantle composition. The density contrast between the liquid and solid phases may lead to a separation between the two phases as a consequence of crystal sedimentation and/or transcrystalline migration of the liquid (e.g. Schiano et al., 2006). This differentiation may lead to the separation between compatible elements that will easily reach the solid phase (e.g. Mg, Cr) while the incompatible elements will remain in the liquid phase (e.g. Al, Na, Fe). The density contrast will affect the chemical segregation and, hence, the geodynamical consequences will differ depending if the buoyancy between the two phases is positive or negative even if the vigor of the convection will potentially erase the chemical contrasts. The ascent or descent of the liquid phase within the partially molten mantle will either contribute to early volcanism or to the formation of a dense basal magma ocean. This aspect is currently under debate: Nomura et al. (2011) suggest an iron rich liquid phase denser than the mantle while Andrault et al. (2012) suggest a buoyant liquid phase that, as a consequence, leads to the formation of a deep residual layer enriched in compatible elements.
## 4.1. Implications for the magma ocean during accretion

The observations of a high liquidus temperature for a chondritic mantle as well as a large gap between its solidus and liquidus temperatures have important implications for the properties of the magma ocean during accretion. As reported in Fig. 5, a fully molten chondritic mantle that would extend to a pressure range of 40–60 GPa (pressure/depth of magma ocean range derived from siderophile elements) exhibits temperatures in the range of 3000–3500 K (e.g. Bouhifd and Jephcoat, 2011; Righter, 2011; Siebert et al., 2012; Fischer et al., 2015). Such conditions correspond to a surface potential temperature of about 2500 K (e.g. Andrault et al., 2011). Such a hot surface is not stable and is only compatible with a transient magma ocean. In fact, surface temperatures higher than 1800 K prevent formation of an (H₂O-CO₂)-rich atmosphere required to produce an efficient thermal blanket to the magma ocean (Abe and Matsui, 1988; Zahnle et al., 1988). Above 1800 K, it appears that significant silicate vaporization occurs, and such rock vapor atmosphere conducts heat easily so that the magma ocean would cool down in a few thousand years (Nagahara and Ozawa, 1996; Sleep et al., 2001). A conclusion of transient magma ocean during the first 100 Ma of Earth's history drawn from melting curves of chondritic mantle is consistent with geochemical evidence (Tucker and Mukhopadhyay, 2014). In this latter study, the authors present evidence for at least two separate magma ocean episodes on early Earth based on the ratio of primordial ³He and ²²Ne in the present-day mantle.

## 4.2. Role on mantle chemical heterogeneities

As noted above the abundance of moderately siderophile elements in the Earth's mantle are better explained with scenarios where the pressure and temperature of core-mantle equilibration increase with time (e.g. Wood et al., 2006; Fischer et al., 2015), hence where the magma ocean's floor deepened throughout Earth's accretion. In the majority of these core-segregation models, the magma ocean floor is located at half or 2/3 of the depth of the core-mantle boundary, which increases progressively on the course of the accretion (e.g. Wade and Wood, 2005; Rubie et al., 2015a; Fischer et al., 2015). This implies that as soon as the Earth achieved about ~30% of its present-day mass, the floor of the magma ocean reached a pressure sufficient to enter the stability domain of bridgmanite (Al-bearing (Mg,Fe)SiO₃ with perovskite-type structure, Bg). Therefore, the equilibration between magma-ocean, Bg-rich solid mantle and segregating metallic phases, may have had an important role during Earth's differentiation.

Archean rocks show anomalies of ¹⁴²Nd in comparison to modern samples, which suggests an early differentiation of the primitive mantle through early decay of the extinct radionuclide ¹⁴⁶Sm (Boyet and Carlson, 2006; Caro et al., 2006; Rizo et al., 2013). The strong fractionation between Sm and Nd required for explaining these isotopic anomalies support a scenario of a magma ocean at ~4.47 Ga (Rizo et al., 2011). The modern terrestrial samples show relatively homogeneous ¹⁴²Nd/¹⁴⁴Nd ratios, resulting from efficient mantle convection and mixing of the early-formed enriched and depleted reservoirs. This constant isotopic ratio recorded by samples from the surface is however higher than those measured in the majority of chondrites, with the exception of some enstatite chondrites (Gannoun et al., 2011). The question of the starting material for modeling the long-term evolution of the silicate Earth is essential. Although none of the chondrite group known so far seems to reproduce exactly the terrestrial composition, enstatite chondrites are closest to Earth's composition when looking at the abundance of ¹⁴²Nd and of many other refractory elements (Warren, 2011). The stable Nd isotopic composition of terrestrial samples is different from that of carbonaceous and ordinary chondrites (Carlson et al., 2007; Burkhardt et al., 2016). This suggests that the ¹⁴²Nd difference is not due to the ¹⁴⁶Sm decay but was caused by isotopic variation imparted during the formation of Solar System. For

enstatite chondrites, the difference is more elusive since their stable Nd isotopic composition are within the terrestrial field considering the analytical error and that they have 10 ppm deficit in ¹⁴²Nd relative to the Earth. If Earth is mainly made by enstatite chondrites, the calculated ¹⁴⁷Sm/¹⁴⁴Nd ratio of the silicate Earth would be 0.200, which is 2% higher than the average chondritic ¹⁴⁷Sm/¹⁴⁴Nd ratio and within the range of values previously measured for whole-rock chondrites. Bouvier and Boyet (2016) have also identified primitive Solar System material (calcium-aluminium-rich refractory inclusions) without isotopic anomalies in Nd compared to the terrestrial composition but this type of material cannot have formed Earth. Therefore, the preservation of an unsampled or hidden reservoir due to early differentiation from magma ocean crystallization is still a conceivable hypothesis but the volume is probably more reduced than estimates proposed by Boyet and Carlson (2005). As a matter of fact, seismic tomography showed the presence of anomalously high seismic velocities at the D" layer, that can be interpreted in several ways, including the occurrence of partial melts, enriched in incompatible elements and representing leftovers from an ancient basal magma ocean, in a scenario where the magma ocean is crystallizing from the mid-lower mantle (Labrosse et al., 2007). The recent discovery of large excesses of ¹⁸²W that seems to be correlated to high ³He/⁴He ratios in modern flood basalts indicates long-term preservation of W isotopic signatures formed in the first 50 Ma of Earth's history (Rizo et al., 2016). The notoriety of this undegassed hidden reservoir hypothesis is typically due to the fact that it can resolve several geochemical paradoxes such as the noble gas isotopic signatures (Coltice et al., 2011).

In fact, a remarkable feature of the core-mantle boundary is the existence of ULVZs (Ultra Low Velocity Zones). The ULVZs are characterized by S- and P- seismic velocity reductions and are usually associated with partial melting of silicates directly above the core (Williams and Garnero, 1996; Revenaugh and Meyer, 1997; Vidale and Hedlin, 1998; Lay et al., 2004; Rost et al., 2005) although alternative interpretations are possible: solid-state origin (Fe-enriched (Mg,Fe)O) (e.g. Wicks et al., 2010; Bower et al., 2011; Brown et al., 2015) or accumulated sediments on top of the Earth's core (Buffett et al., 2000). More details on the origin of the ULVZs can also be found in the review papers by Tackley (2012) or Carlson et al. (2014).

In order to trap and isolate enriched melts deep at the CMB, their density should stay higher than the solid mantle. Two recent experimental studies aimed at determining the density contrast between Bg and silicate melt by measuring the partition coefficients of Fe between Bg and melt (D<sub>Fe</sub><sup>Bg/melt</sup>) at various pressures, down to CMB conditions (Andrault et al., 2012; Nomura et al., 2011). These studies provided contrasting results where D<sub>Fe</sub><sup>Bg/melt</sup> was of 0.2–0.3 in Nomura et al. (2011) and Tateno et al. (2014) and 0.5–0.6 in Andrault et al. (2012), leading to contradictory interpretations for the fate of deep silicate melts, which are either denser or more buoyant than the solid silicate. A recent study conducted at fixed pressure (25 GPa) has proposed another explanation to this controversy, showing that the degree of melting (F) can significantly vary D<sub>Fe</sub><sup>Bg/melt</sup> (Boujibar et al., 2016). This study shows that at 25 GPa, D<sub>Fe</sub><sup>Bg/melt</sup> is close to 1 at extremely low degree of melting and becomes constant and equaling 0.3 ± 0.1 when F is higher than 10%. Therefore, the degree of melting is another parameter that should be taken into account in order to address mantle dynamics and the fate of a molten enriched reservoir at the base of the lower mantle.

## 4.3. Role on the redox state of the mantle

Several studies have previously proposed that Bg crystallization had major effects on the redox state of the mantle (Frost et al., 2004; Wood et al., 2006). Bg has the ability to integrate high levels of Fe³⁺ in its structure. In particular, it was previously shown that the presence of Al enhances the incorporation of Fe³⁺, owing to the coupled substitution of Al and Fe³⁺ for Mg and Si in their respective dodecahedral and octahedral sites (Andrault et al., 2007; Lauterbach et al., 2000). This
substitution mechanism is so stable that it was suggested to force $Fe^{2+}$ to disproportionate into $Fe^{3+}$ and $Fe^{0}$ through the reaction:

$$3 \text{Fe}^{2+} = \text{Fe}^{0} + 2 \text{Fe}^{3+} \quad (2)$$

Therefore, since the magma ocean's floor reached a depth sufficiently high to crystallize Bg, the magma ocean/solid mantle boundary was presumably acting as an oxygen pump, where $Fe^{3+}$ was stabilized in Bg and progressively released to the magma ocean, while $Fe^{0}$ was systematically trapped in the solid lower mantle and subsequently falling to the core. Therefore the crystallization of the magma ocean was considered as a plausible mechanism for raising the oxygen state of the Earth's mantle (Wade and Wood, 2005). This idea was based on the fact that Al-content increases the solubility of $Fe^{3+}$ in Bg independently from the oxygen fugacity (Frost et al., 2004; Frost and Langenhorst, 2002). However recent experimental studies showed that the oxygen fugacity (nature of the starting composition or capsule material) has an important effect on $Fe^{3+}$ incorporation in Bg (Nakajima et al., 2012; Boujibar et al., 2016). In particular, recent experimental data on equilibrium between Bg, silicate melt and liquid metal showed that at an $f_{O_2}$ of ~IW-2 (2 log units below IW buffer), the $Fe^{3+}$-content in Bg is relatively low (~20%) and weakly dependent on the Al-content. Aluminium starts having an important positive effect on $Fe^{3+}$ in Bg only when the oxygen fugacity is higher than ~IW-1 to IW (Fig. 6) (Boujibar et al., 2016). These results are challenging previous planetary differentiation models and imply that Bg crystallization during the cooling of the magma ocean can only have a substantial effect on the redox state of the mantle after the end of core segregation. Given that the final $f_{O_2}$ equilibration between core and mantle is about IW-2 (based on present-day FeO content of the mantle of ~8 wt% and on the 85 wt% Fe of the core, estimated from the Earth's core density deficit), one can conclude that any increase of the Earth's redox state before the completion of core segregation should originate from other mechanisms than Bg crystallization from the magma ocean.

## 5. Concluding remarks

Earth's core formation models that include a combination of a multistage core-mantle differentiation with N-body accretion simulations lead to coherent view of the early history of our planet. The actual models suggest that accretion of Earth was heterogeneous with the

<figure><img src="image_79.png" /><figcaption>**Fig. 6.** $Fe^{3+}/\sum Fe$ in bridgmanite as a function of Al-content for different $f_{O_2}$ conditions ranging from IW-2 to IW-1 (black circles) (Boujibar et al., 2016; Terasaki et al., 2007; Frost et al., 2004), IW-1 to IW (gray diamonds) (Frost et al., 2004) and close or higher than IW buffer (where capsule material was metallic Fe, Au, Re or Pt) (Frost and Langenhorst, 2002; Frost et al., 2004; Lauterbach et al., 2000; McCammon et al., 2004; Saikia et al., 2009).</figcaption></figure>

initial 60–70% of Earth's mass forming from reducing material and the final 30–40% from relatively oxidized material. In addition, the available observations point to the occurrence of magma oceans in the early evolution of terrestrial planets. By applying this scenario we found that the Earth's core is a plausible reservoir for primordial helium. It is also possible to match the observed subchondritic Nb/Ta ratio of the bulk silicate Earth, highlighting that planetary accretion of reduced materials played an important role in the chemical evolution of Earth, and more generally, that Nb/Ta ratio can be used to trace prevailing oxygen fugacities during the segregation of planetary cores.

The available experiments of a high liquidus temperature for a chondritic mantle, as well as observations based on the ratio of primordial helium and neon isotopes, lead to a conclusion of transient magma ocean during the first 100 Ma of Earth's history. At least two separate magma ocean episodes on early Earth are suggested from noble gases isotopes.

To solve the existing controversy about the partitioning of Fe between bridgmanite and silicate melt, it seems that the degree of melting should be taken into account in order to address mantle dynamics and the fate of a molten enriched reservoir at the base of the lower mantle.

## Acknowledgments

This paper was presented in the 10th Silicate Melt Workshop – La Pétite Pierre. We thank Don Dingwell, Kai Uwe Hess and Pascal Richet for the organization of the workshop and for being the guest editors of this special issue. Part of this work was funded by the CNRS-INSU (program: PNP), ANR (project: Oxy-Deep), and ERC (project: ISOREE). This research was also financed by the French Government Laboratory of Excellence initiative n°ANR-10-LABX-0006, the Région Auvergne and the European Regional Development Fund. This is Laboratory of Excellence ClerVolc contribution number 231.

## References

* Abe, Y., 1997. Thermal and chemical evolution of the terrestrial magma ocean. Phys. Earth Planet. Inter. 100, 27–39.

* Abe, Y., Matsui, T., 1988. Evolution of an impact-generated $H_2O-CO_2$ atmosphere and formation of a hot proto-ocean on Earth. J. Atmos. Sci. 45, 3081–3101.

* Albarède, F., 2009. Volatile accretion history of the terrestrial planets and dynamic implications. Nature 461, 1227–1233.

* Allègre, C.J., Staudacher, T., Sarda, P., Kurz, M., 1983. Constraints on evolution of Earth's mantle from rare gas systematics. Nature 303, 762–766.

* Amelin, Y., Krot, A.N., Hutcheon, I.D., Ulyanov, A.A., 2002. Lead isotopic ages of chondrules and calcium-aluminium-rich inclusions. Science 297, 1678–1683.

* Andrault, D., Bolfan-Casanova, N., Bouhifd, M.A., Guignot, N., Kawamoto, T., 2007. The role of Al-defects on the equation of state of Al-(Mg,Fe)SiO$_3$ perovskite. Earth Planet. Sci. Lett. 263, 167–179.

* Andrault, D., Bolfan-Casanova, N., Lo Nigro, G., Bouhifd, M.A., Garbarino, G., Mezouar, M., 2011. Solidus and liquidus profiles of chondritic mantle: implication for melting of the Earth across its history. Earth Planet. Sci. Lett. 304, 251–259.

* Andrault, D., Petitgirard, S., Lo Nigro, G., Devidal, J.-L., Veronesi, G., Garbarino, G., Mezouar, M., 2012. Solid-liquid iron partitioning in Earth's deep mantle. Nature 487, 354–357.

* Andrault, D., Pesce, G., Bouhifd, M.A., Bolfan-Casanova, N., Hénot, J.-M., Mezouar, M., 2014. Melting of subducted basalt at the core-mantle boundary. Science 344, 892–895.

* Armytage, R., Jephcoat, A.P., Bouhifd, M.A., Porcelli, D., 2013. Metal-silicate partitioning of iodine at high pressures and temperatures: implications for the Earth's core and $^{129}$Xe budgets. Earth Planet. Sci. Lett. 373, 140–149.

* Bajgain, S., Ghosh, D.B., Karki, B.B., 2015. Structure and density of basaltic melts at mantle conditions from first-principles simulations. Nat. Commun. 6:8578. http://dx.doi.org/10.1038/ncomms9578.

* Bergin, E.A., Blake, G.A., Ciesla, F., Hirschmann, M.M., Li, J., 2015. Tracing the ingredients for a habitable Earth from interstellar space through planet formation. Proc. Natl. Acad. Sci. 112, 8965–8970.

* Bouhifd, M.A., Jephcoat, A.P., 2003. The effect of pressure on partitioning of Ni and Co between silicate and iron-rich metal liquids: a diamond-anvil cell study. Earth Planet. Sci. Lett. 209, 245–255.

* Bouhifd, M.A., Jephcoat, A.P., 2006. Aluminium control of argon solubility in silicate melts under pressure. Nature 439, 961–964.

* Bouhifd, M.A., Jephcoat, A.P., 2011. Convergence of Ni and Co metal-silicate partition coefficients in the deep magma-ocean and light-elements core solubility at high pressures. Earth Planet. Sci. Lett. 307, 341–348.

* Bouhifd, M.A., Jephcoat, A.P., Kelley, S.P., 2008. Argon solubility drop in silicate melts at high pressures: a review of recent experiments. Chem. Geol. 256, 252–258.

* Bouhifd, M.A., Jephcoat, A.P., Heber, V.S., Kelley, S.P., 2013. Helium in an early Earth's core. Nat. Geosci. 6, 982–986.
Bouhifd, M.A., Boyet, M., Cartier, C., Hammouda, T., Bolfan-Casanova, N., Devidal, J.L.,
Andrault, D., 2015. Superchondritic Sm/Nd ratio of the Earth: impact of Earth's core
formation. Earth Planet. Sci. Lett. 413, 158–166.

Boujibar, A., Andrault, D., Bouhifd, M.A., Bolfan-Casanova, N., Devidal, J.L., Trcera, N., 2014. Metal-silicate partitioning of sulfur, new experimental and thermodynamic constraints on planetary accretion. Earth Planet. Sci. Lett. 391, 42–54.

Boujibar, A., Andrault, D., Bolfan-Casanova, N., Bouhifd, M.A., Monteux, J., 2015. Cosmochemical fractionation by collisional erosion during the Earth's accretion. Nat. Commun. 6:8295. http://dx.doi.org/10.1038/ncomms9295.

Boujibar, A., Bolfan-Casanova, N., Andrault, D., Bouhifd, M.A., Trcera, N., 2016. Incorporation of Fe²⁺ and Fe³⁺ in bridgmanite during magma ocean crystallization. Am. Mineral. 101, 1560–1570.

Bouvier, A., Boyet, M., 2016. Primitive solar system materials and Earth share a common initial $^{142}$Nd abundance. Nature 537, 399–402.

Bouvier, A., Wadhwa, M., 2010. The age of the solar system redefined by the oldest Pb-Pb age of a meteoritic inclusion. Nat. Geosci. 3, 637–641.

Bower, D.J., Wicks, J.K., Gurnis, M., Jackson, J.M., 2011. A geodynamic and mineral physics model of a solid-state ultralow-velocity zone. Earth Planet. Sci. Lett. 303, 193–202.

Boyet, M., Carlson, R.W., 2005. $^{142}$Nd evidence for early (>4.53 Ga) global differentiation of the silicate Earth. Science 309, 576–581.

Boyet, M., Carlson, R.W., 2006. A new geochemical model for the Earth's mantle inferred from $^{146}$Sm-$^{142}$Nd systematics. Earth Planet. Sci. Lett. 250, 254–268.

Boyet, M., Carlson, R.W., Borg, L.E., Horan, M., 2015. Sm-Nd systematics of lunar ferroan anorthositic suite rocks: constraints on lunar crust formation. Geochim. Cosmochim. Acta 148, 203–218.

Brennecka, G.A., Wadhwa, M., 2012. Uranium isotope compositions of the basaltic angrite meteorites and the chronological implications for the early Solar System. Proc. Natl. Acad. Sci. U. S. A. 109, 9299–9303.

Brooker, R.A., Du, Z., Blundy, J.D., Kelley, S.P., Allan, N.L., Wood, B.J., Chamorro, E.M., Wartho, J.A., Purton, J.A., 2003. The "zero charge" partitioning behaviour of noble gases during mantle melting. Nature 423, 738–741.

Brown, S.P., Thorne, M.S., Miyagi, L., Rost, S., 2015. A compositional origin to ultralow velocity zones. Geophys. Res. Lett. 42:1039–1045. http://dx.doi.org/10.1002/2014GL062097.

Buffett, B.A., Garnero, E.J., Jeanloz, R., 2000. Sediments at the top of Earth's core. Science 290, 1338–1342.

Burkhardt, C., Borg, L.E., Brennecka, G.A., Shollenberger, Q.R., Dauphas, N., Kleine, T., 2016. A nucleosynthetic origin for the Earth's anomalous $^{142}$Nd composition. Nature 537, 394–398.

Canup, R.M., 2008. Accretion of the Earth. Phil. Trans. R. Soc. A 366, 4061–4075.

Canup, R.M., 2012. Forming a Moon with an Earth-like composition via a giant impact. Science 338, 1052–1055.

Carlson, R.W., Boyet, M., Horan, M., 2007. Chondrite barium, neodymium, and samarium isotopic heterogeneity and early Earth differentiation. Science 316, 1175–1178.

Carlson, R.W., et al., 2014. How did early Earth become our modern world? Annu. Rev. Earth Planet. Sci. 42, 151–178.

Caro, G., Bourdon, B., Birck, J.L., Moorbath, S., 2006. High precision $^{142}$Nd/$^{144}$Nd measurements in terrestrial rocks: constraints on the early differentiation of the Earth's mantle. Geochim. Cosmochim. Acta 70, 164–191.

Cartier, C., Hammouda, T., Boyet, M., Bouhifd, M.A., Devidal, J.L., 2014. Redox control on Nb/Ta fractionation during planetary accretion. Nat. Geosci. 7, 573–576.

Cartier, C., Hammouda, T., Boyet, M., Mathon, O., Testemale, D., Moine, B.N., 2015. Evidence for Nb$^{2+}$ and Ta$^{3+}$ in silicate melts under highly reducing conditions: A XANES study. Am. Mineral. 100, 2152–2158.

Chabot, N.L., Draper, D.S., Agee, C.B., 2005. Conditions of core formation in the Earth: Constraints from nickel and cobalt partitioning. Geochim. Cosmochim. Acta 69, 2141–2151.

Chambers, J.E., 2004. Planetary accretion in the inner solar system. Earth Planet. Sci. Lett. 223, 241–252.

Chamorro-Perez, E., Gillet, P., Jambon, A., 1996. Argon solubility in silicate melts at very high pressures. Experimental set-up and preliminary results for silica and anorthite melts. Earth Planet. Sci. Lett. 145, 97–107.

Clesi, V., Bouhifd, M.A., Bolfan-Casanova, N., Manthilake, G., Fabbrizio, A., Andrault, D., 2016. Effect of H₂O on metal-silicate partitioning of Ni, Co, V, Cr, Mn and Fe: implications for the oxidation state of the Earth and Mars. Geochim. Cosmochim. Acta 192, 97–121.

Coltice, N., Moreira, M., Herlund, J., Labrosse, S., 2011. Crystallization of a basal magma ocean recorded by helium and neon. Earth Planet. Sci. Lett. 308, 193–199.

Connelly, J.N., Bizarro, M., Krot, A.N., Nordlund, A., Wielandt, D., Ivanova, M.A., 2012. The absolute chronology and thermal processing of solids in the Solar protoplanetary disk. Science 338, 651–655.

Dahl, T.W., Stevenson, D.J., 2010. Turbulent mixing of metal and silicate during planet accretion and interpretation of the Hf-W chronometer. Earth Planet. Sci. Lett. 295, 177–186.

Dasgupta, R., Chi, H., Shimizu, N., Buono, A.S., Walker, D., 2013. Carbon solution and partitioning between metallic and silicate melts in a shallow magma ocean: implications for the origin and distribution of terrestrial carbon. Geochim. Cosmochim. Acta 102, 191–212.

Dauphas, N., Chen, J.H., Zhang, J., Papanastassiou, D.A., Davis, A.M., Travaglio, C., 2014. Calcium-48 isotopic anomalies in bulk chondrites and achondrites: Evidence for a uniform isotopic reservoir in the inner protoplanetary disk. Earth Planet. Sci. Lett. 407, 96–108.

Deguen, R., Olson, P., Cardin, P., 2011. Experiments on turbulent metal-silicate mixing in a magma ocean. Earth Planet. Sci. Lett. 310, 303–313.

Deguen, R., Landeau, M., Olson, P., 2014. Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans. Earth Planet. Sci. Lett. 391, 274–287.

Doremus, R.H., 1966. Physical solubility of gases in fused silica. J. Am. Ceram. Soc. 49, 461–462.

Elkins-Tanton, L.T., 2012. Magma oceans in the inner solar system. Annu. Rev. Earth Plan- et. Sci. 40, 113–139.

Fiquet, G., Auzende, L., Siebert, J., Corgne, A., Bureau, H., Ozawa, H., Garbarino, G., 2010. Melting of peridotite to 140 GPa. Science 329, 1516–1518.

Fischer, R.A., Ciela, F.J., 2014. Dynamics of the terrestrial planets from a large number of N-body simulations. Earth Planet. Sci. Lett. 392, 28–38.

Fischer, R.A., Nakajima, Y., Campbell, A.J., Frost, D.J., Harries, D., Langenhorst, F., Miyajima, N., Pollock, K., Rubie, D.C., 2015. High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O. Geochim. Cosmochim. Acta 167, 177–194.

Frost, D.J., Langenhorst, F., 2002. The effect of Al₂O₃ on Fe-Mg partitioning between magnesiowüstite and magnesium silicate perovskite. Earth Planet. Sci. Lett. 199, 227–241.

Frost, D.J., Liebske, C., Langenhorst, F., McCammon, C.A., Tronnes, R.G., Rubie, D.C., 2004. Experimental evidence for the existence of iron-rich metal in the Earth's lower mantle. Nature 428, 409–412.

Gaillard, F., Scaillet, B., 2014. A theoretical framework for volcanic degassing chemistry in a comparative planetology perspective and implications for planetary atmospheres. Earth Planet. Sci. Lett. 403, 307–316.

Gaillard, F., Scaillet, B., Pichacant, M., Iacono-Marziano, G., 2015. The redox geodynamics linking basalts and their mantle sources through space and time. Chem. Geol. 418, 217–233.

Gannoun, A., Boyet, M., Rizo, H., El Goresy, A., 2011. $^{146}$Sm-$^{142}$Nd systematics measured in enstatite chondrites reveals a heterogeneous distribution of $^{142}$Nd in the solar nebula. PNAS 108, 7693–7697.

Gessmann, C.K., Rubie, D.C., 2000. The origin of the depletions of V, Cr and Mn in the mantles of the Earth and Moon. Earth Planet. Sci. Lett. 184, 95–107.

Goldschmidt, V.M., 1937. The principles of distribution of chemical elements in minerals and rocks. The Seventh Hugo Müller Lecture, Delivered before the Chemical Society on March 17th, 1937.

Greenwood, R.C., Franchi, I.A., Jambon, A., Buchanan, P.C., 2005. Widespread magma oceans on asteroidal bodies in the early solar system. Nature 435, 916–918.

Grossman, L., 1972. Condensation in the primitive solar nebula. Geochim. Cosmochim. Acta 105, 597–619.

Grossman, L., Beckett, J.R., Fedkin, A.V., Simon, S.B., Ciesla, F.J., 2008. Redox conditions in the solar nebula: observational, experimental, and theoretical constraints. Rev. Mineral. Geochem. 68, 93–140.

Guillot, B., Sator, N., 2012. Noble gases in high-pressure silicate liquids: a computer simulation study. Geochim. Cosmochim. Acta 80, 51–69.

Halliday, A.N., 2013. The origins of volatiles on the terrestrial planets. Geochim. Cosmochim. Acta 105, 146–171.

Halliday, A.N., Wood, B.J., 2009. How did Earth accrete? Science 325, 44–45.

Harper, C.L., Jacobsen, S.B., 1996. Noble gases and Earth's accretion. Science 273, 1814–1818.

Hartmann, W.K., Davis, D.R., 1975. Satellite-sized planetesimals and lunar origin. Icarus 24, 504–519.

Heber, V.S., Brooker, R.A., Kelley, S.P., Wood, B.J., 2007. Crystal-melt partitioning of noble gases (helium, neon, argon, krypton, and xenon) for olivine and clinopyroxene. Geochim. Cosmochim. Acta 71, 1041–1061.

Herzberg, C., et al., 2013. Nickel and helium evidence for melt above the core-mantle boundary. Nature 493, 393–397.

Hirschmann, M.M., 2016. Constraints on the early delivery and fractionation of Earth's major volatiles from C/H, C/N, and C/S ratios. Am. Mineral. 101, 540–553.

Hirschmann, M.M., Dasgupta, R., 2009. The H/C ratios of Earth's near-surface and deep reservoirs, and consequences for deep Earth volatile cycles. Chem. Geol. 262, 4–16.

Hirschmann, M.M., Withers, A.C., Ardia, P., Foley, N.T., 2012. Solubility of molecular hydrogen in silicate melts and consequences for volatile evolution of terrestrial planets. Earth Planet. Sci. Lett. 345–348, 38–48.

Hofmann, A.W., 1988. Chemical differentiation on the Earth: the relationship between mantle, continental crust, and oceanic crust. Earth Planet. Sci. Lett. 90, 297–314.

Honda, M., McDougall, I., Patterson, D.B., Doulgeris, A., Clague, D.A., 1991. Possible solar noble-gas component in Hawaiian basalts. Nature 349, 149–151.

Javoy, M., 1995. The integral enstatite chondrite model of the Earth. Geophys. Res. Lett. 22, 2219–2222.

Javoy, M., et al., 2010. The chemical composition of the Earth: enstatite chondrite models. Earth Planet. Sci. Lett. 293, 259–268.

Jephcoat, A.P., Bouhifd, M.A., Porcelli, D., 2008. Partitioning experiments in the laser-heated diamond anvil cell: volatile content in the Earth's core. Phil. Trans. R. Soc. A 366, 4295–4314.

Kadik, A.A., Koltashev, V.V., Kryukova, E.B., Plotnichenko, V.G., Tsekhonya, T.I., Kononkova, N.N., 2015. Solubility of nitrogen, carbon, and hydrogen in FeO-N₂O-Al₂O₃-SiO₂ melt and liquid iron alloy: influence of oxygen fugacity. Geochem. Int. 53, 849–868.

Kargel, J.S., Lewis, J.S., 1993. The composition and early evolution of Earth. Icarus 105, 1–25.

Ke, Y., Solomatov, V.S., 2009. Coupled core–mantle thermal evolution of early Mars. J. Geophys. Res. Planets 114 (13), 1–12.

Kleine, T., Touboul, M., Bourdon, B., Nimmo, F., Mezger, K., Palme, H., Jacobsen, S.B., Yin, Q.-Z., Halliday, A.N., 2009. Hf-W chronology of the accretion and early evolution of asteroids and terrestrial planets. Geochim. Cosmochim. Acta 73, 5150–5188.

Kruijer, T.S., Sprung, P., Kleine, T., Leya, I., Burkhardt, C., Wieler, R., 2012. Hf-W chronometry of core formation in planetesimals inferred from weakly irradiated iron meteorites. Geochim. Cosmochim. Acta 99, 287–304.

Kruijer, T.S., Fischer-Gödde, M., Kleine, T., Burkhardt, C., Wieler, R., 2014. Nucleosynthetic W isotope anomalies and the Hf-W chronometry of Ca-Al-rich inclusions. Earth Plan- et. Sci. Lett. 403, 317–327.
Labrosse, S., Hernlund, J.W., Coltice, N., 2007. A crystallizing dense magma ocean at the base of the Earth's mantle, Nature 450, 866–869.

Larimer, J.W., Buseck, P.R., 1974. Equilibration temperatures in enstatite chondrites. Geochim. Cosmochim. Acta 38, 471–477.

Lauterbach, S., McCammon, C.A., van Aken, P., Langenhorst, F., Seifert, F., 2000. Mössbauer and ELNES spectroscopy of (Mg,Fe)(Si,Al)O₃ perovskite: a highly oxidized component of the lower mantle. Contrib. Mineral. Petrol. 138, 17–26.

Lay, T., Garnero, E.J., Williams, Q., 2004. Partial melting in a thermo-chemical boundary layer at the base of the mantle. Phys. Earth Planet. Inter. 146, 441–467.

Li, J., Agee, C.B., 1996. Geochimistry of mantle–core differentiation at high pressure. Nature 381, 686–689.

Li, J., Agee, C.B., 2001. The effect of pressure, temperature, oxygen fugacity and composition on partitioning of nickel and cobalt between liquid Fe–Ni–S and liquid silicate: implications for the Earth's core formation. Geochim. Cosmochim. Acta 65, 1821–1832.

Li, Y., Dasgupta, R., Tsuno, K., Monteleone, B., Shimizu, N., 2016. Carbon and sulfur budget of the silicate Earth explained by accretion of differentiated planetary embryos. Nat. Geosci. 7, 781–785.

Lupton, J.E., Craig, H., 1975. Excess ³He in oceanic basalts: evidence for terrestrial primordial helium. Earth Planet. Sci. Lett. 26, 133–139.

Macpherson, C.G., Hilton, D.R., Sinton, J.M., Poreda, R.J., Craig, H., 1998. High ³He/⁴He ratios in the Manus backarc basin: Implications for mantle mixing and the origin of plumes in the western Pacific Ocean. Geology 26, 1007–1010.

Malavergne, V., Tarrida, M., Combes, R., Bureau, H., Jones, J., Schwandt, C., 2007. New high-pressure and high-temperature metal/silicate partitioning of U and Pb: implications for the cores of the Earth and Mars. Geochim. Cosmochim. Acta 71, 2637–2655.

Mann, U., Frost, D.J., Rubie, D.C., 2009. Evidence for high pressure core-mantle differentiation from the metal-silicate partitioning of lithophile and weakly-siderophile elements. Geochim. Cosmochim. Acta 73, 7360–7386.

Marty, B., 2012. The origins and concentrations of water, carbon, nitrogen and noble gases on Earth. Earth Planet. Sci. Lett. 313–314, 56–66.

Marty, B., Alexander, C.M.O.'D., Raymond, S.N., 2013. Primordial origins of Earth's carbon. Rev. Mineral. Gechem. 75, 149–181.

McCammon, C.A., Lauterbach, S., Seifert, F., Langenhorst, F., van Aken, P.A., 2004. Iron oxidation state in lower mantle mineral assemblages. I. Empirical relations derived from high-pressure experiments. Earth Planet. Sci. Lett. 222, 435–449.

McDonough, W.F., Sun, S.-S., 1995. The composition of the Earth. Chem. Geol. 120, 223–253.

Mizuno, H., Nakazawa, K., Hayashi, C., 1980. Dissolution of the primordial rare gases into the molten Earth's material. Earth Planet. Sci. Lett. 50, 202–210.

Monteux, J., Ricard, Y., Coltice, N., Dubuffet, F., Ulrova, M., 2009. A model of metal–silicate separation on growing planets. Earth Planet. Sci. Lett. 287, 353–362.

Monteux, J., Araudrault, D., Samuel, H., 2016. On the cooling of a deep terrestrial magma ocean, Earth and planet. Sci. Lett. 448, 140–149.

Morbidelli, A., Lunine, J.J., O'Brien, D.P., Raymond, S.N., Walsh, K.J., 2012. Building terrestrial planets. Annu. Rev. Earth Planet. Sci. 40, 251–275.

Moreira, M., 2013. Noble gas constraints on the origin and evolution of Earth's volatiles. Geochim. Perspect. 2, 229–403.

Moreira, M., Charnoz, S., 2016. The origin of the neon isotopes in chondrites and on Earth. Earth Planet. Sci. Lett. 433, 249–256.

Mosenfelder, J.L., Asimow, P.D., Ahrens, T.J., 2007. Thermodynamic properties of Mg₂SiO₄ liquid at ultra-high pressures from shock measurements to 200 GPa on forsterite and wadsleyite. J. Geophys. Res. Solid Earth 112:B06208. http://dx.doi.org/10.1029/2006JB004364.

Mosenfelder, J.L., Asimow, P.D., Frost, D.J., Rubie, D.C., Ahrens, T.J., 2009. The MgSiO₃ system at high pressure: thermodynamic properties of perovskite, postperovskite, and melt from global inversion of shock and static compression data. J. Geophys. Res. Solid Earth 114:B01203. http://dx.doi.org/10.1029/2008JB005900.

Mukhopadhyay, S., 2012. Early differentiation and volatile accretion recorded in deep-mantle neon and xenon. Nature 486, 101–104.

Münker, C., Pfänder, J.A., Weyer, S., Büchl, A., Kleine, T., Mezger, K., 2003. Evolution of planetary cores and the Earth-Moon system from Nb/Ta systematics. Science 301, 84–87.

Nagahara, H., Ozawa, K., 1996. Evaporation of forsterite in H-2 gas. Geochim. Cosmochim. Acta 60, 1445–1459.

Nakajima, M., Stevenson, D.J., 2015. Melting and mixing states of the Earth's mantle after the Moon-forming impact. Earth Planet. Sci. Lett. 427, 286–295.

Nakajima, Y., Frost, D.J., Rubie, D.C., 2012. Ferrous iron partitioning between magnesium silicate perovskite and ferropericlase and the composition of perovskite in the Earth's lower mantle. J. Geophys. Res. 117, 1–12.

Nebel, O., van Westrenen, W., Vroon, P.Z., Wille, M., Raith, M.M., 2010. Deep mantle storage of the Earth's missing niobium in late-stage residual melts from a magma ocean. Geochim. Cosmochim. Acta 74, 4392–4404.

Nevins, D., Spera, F.J., 1998. Molecular dynamics simulations of molten CaAl₂Si₂O₈: dependence of structure and properties on pressure. Am. Mineral. 83, 1220–1230.

Newsom, H.E., 1990. Accretion and core formation in the Earth: evidence from siderophile elements. In: Newson, H.E., Jones, J.H. (Eds.), Origin of the Earth. Oxford University Press, New York, pp. 273–288.

Niwa, K., Miyakawa, C., Yagi, T., Matsuda, J., 2013. Argon solubility in SiO₂ melt under high pressures: a new experimental result using laser-heated diamond anvil cell. Earth Planet. Sci. Lett. 363, 1–8.

Nomura, R., Ozawa, H., Tateno, S., Hirose, K., Hernlund, J.W., Muto, S., Ishii, H., Hiraoka, N., 2011. Spin crossover and iron-rich silicate melt in the Earth's deep mantle. Nature 473, 199–202.

Nomura, R., Hirose, K., Uesugi, K., Ohishi, Y., Tsuchiyama, A., Miyake, A., Ueno, Y., 2014. Low core–mantle boundary temperature inferred from the solidus of pyrolite. Science 343, 522–525.

O'Brien, D.P., Morbidelli, A., Levison, H.F., 2006. Terrestrial planet formation with strong dynamical friction. Icarus 184, 39–58.

Ohtani, E., Yurimoto, H., 1996. Element partitioning between metallic liquid, magnesiowüstite, and silicate liquid at 20 GPa and 2500 °C: a secondary ion mass spectrometric study. Geophys. Res. Lett. 23, 1993–1996.

O'Neill, H.S.C., Palme, H., 1998. Composition of the silicate Earth: implications for accretion and core formation. In: Jackson, I. (Ed.), The Earth's Mantle, Composition, Structure, and Evolution. Cambridge University Press, pp. 3–126.

O'Neill, H.S.C., Canil, D., Rubie, D.C., 1998. Oxide-metal equilibria to 2500 °C and 25 GPa: Implications for core formation and the light component in the Earth's core. J. Geophys. Res. 103, 12239–12260.

Pepin, R.O., Porcelli, D., 2006. Xenon isotope systematics, giant impacts, and mantle degassing on the early Earth. Earth Planet. Sci. Lett. 250, 470–485.

Porcelli, D., Ballentine, C.J., 2002. Models for the distribution of terrestrial noble gases and evolution of the atmosphere. Rev. Mineral. Geochem. 46, 411–480.

Porcelli, D., Halliday, A.N., 2001. The core as a possible source of mantle helium. Earth Planet. Sci. Lett. 192, 45–56.

Porcelli, D., Woolum, D., Cassen, P., 2001. Deep Earth rare gases: initial inventories, capture from the solar nebula, and losses during Moon formation. Earth Planet. Sci. Lett. 193, 237–251.

Pradhan, G.K., Fiquet, G., Siebert, J., Auzende, A.-L., Morard, G., Antonangeli, D., Garbarino, G., 2015. Melting of MORB at core-mantle boundary. Earth Planet. Sci. Lett. 431, 247–255.

Press, W.H., Teukolsky, S.A., Vetterling, W.T., Flannery, B.P., 1993. Numerical Recipes in FORTRAN the Art of Scientific Computing. second ed. Cambridge University Press, New York, NY, USA.

Raymond, S.N., Quinn, T., Lunine, J.J., 2006. High-resolution simulations of the final assembly of Earth-like planets I. Terrestrial accretion and dynamics. Icarus 183, 265–282.

Revenaugh, J., Meyer, R., 1997. Seismic evidence of partial melt within a possibly ubiquitous low-velocity layer at the base of the mantle. Science 277, 670–673.

Righter, K., 2011. Prediction of metal–silicate partition coefficients for siderophile elements: an update and assessment of PT conditions for metal–silicate equilibrium during accretion of the earth. Earth Planet. Sci. Lett. 304, 158–167.

Rizo, H., Boyet, M., Blichert-Toft, J., Rosing, M., 2011. Combined Nd and Hf isotope evidence for deep-placed source of Isua lavas. Earth Planet. Sci. Lett. 312, 267–279.

Rizo, H., Boyet, M., Blichert-Toft, J., Rosing, M.T., 2013. Early mantle dynamics inferred from Nd-142 variations in Archean rocks from southwest Greenland. Earth Planet. Sci. Lett. 377, 324–335.

Rizo, H., et al., 2016. Preservation of Earth-forming events in the tungsten isotopic composition of modern flood basalts. Science 352, 809–812.

Roskosz, M., Bouhifd, M.A., Jephcoat, A.P., Marty, B., Mysen, B.O., 2013. Nitrogen solubility in molten metal and silicate at high pressure and temperature. Geochim. Cosmochim. Acta 121, 15–28.

Rost, S., Garnero, E.J., Williams, Q., 2005. Seismological constraints on a possible plume root at the core-mantle boundary. Nature 435, 666–669.

Rubie, D.C., Melosh, H.J., Reid, J.E., Liebske, C., Righter, K., 2003. Mechanisms of metal–silicate equilibration in the terrestrial magma ocean. Earth Planet. Sci. Lett. 205, 239–255.

Rubie, D.C., Gessmann, C.K., Frost, D.J., 2004. Partitioning of oxygen during core formation on the Earth and Mars. Nature 429, 58–61.

Rubie, D.C., Frost, D.J., Mann, U., Asahara, Y., Nimmo, F., Tsuno, K., Kegler, P., Holzheid, A., Palme, H., 2011. Heterogeneous accretion, composition and core–mantle differentiation of the Earth. Earth Planet. Sci. Lett. 301, 31–42.

Rubie, D.C., et al., 2015a. Accretion and differentiation of the terrestrial planets with implications for the compositions of early-formed solar system bodies and accretion of water. Icarus 248, 89–108.

Rubie, D.C., Nimmo, F., Melosh, H.J., 2015b. Formation of the Earth's core. treatise on geophysics, second edition. In: Stevenson, D. (Ed.), Evolution of the Earth 9. Elsevier, Amsterdam, pp. 43–79.

Rudnick, R., Barth, M., Horn, I., McDonough, W., 2000. Rutile-bearing refractory eclogites: missing link between continents and depleted mantle. Science 287, 278–281.

Samuel, H., 2012. A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans. Earth Planet. Sci. Lett. 313, 105–114.

Samuel, H., Tackley, P.J., Evonuk, M., 2010. Heat partitioning in terrestrial planets during core formation by negative diapirism. Earth Planet. Sci. Lett. 290, 13–19.

Saikia, A., Boffa Ballaran, T., Frost, D.J., 2009. The effect of Fe and Al substitution on the compressibility of MgSiO₃-perovskite determined through single-crystal X ray diffraction. Phys. Earth Planet. Inter. 173, 153–161.

Sanloup, C., et al., 2013. Structural change in molten basalt at deep mantle conditions. Nature 503, 104–107.

Schiano, P., Provost, A., Clocchiatti, R., Faure, F., 2006. Transcrystalline melt migration and Earth's mantle. Science 314, 970–974.

Shackelford, J.F., Studt, P.L., Fulrath, R.M., 1972. Solubility of gases in glass II. He, Ne, and H₂ in fused silica. J. Appl. Phys. 43, 1619–1626.

Shannon, R.D., 1976. Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides. Acta Crystallogr. Sect. A: Cryst. Phys. Diffr. Theor. Gen. Crystallogr. 32, 751–767.

Shelby, J.E., 1976. Pressure dependence of helium and neon solubility in vitreous silica. J. Appl. Phys. 47, 135–139.

Siebert, J., Corgne, A., Ryerson, F.J., 2011. Systematics of metal–silicate partitioning for many siderophile elements applied to Earth's core formation. Geochim. Cosmochim. Acta 75, 1451–1489.

Siebert, J., Badro, J., Antonangeli, D., Ryerson, F.J., 2012. Metal-silicate partitioning of Ni and Co in a deep magma ocean. Earth Planet. Sci. Lett. 321–322, 189–197.

Siebert, J., Badro, J., Antonangeli, D., Ryerson, F.J., 2013. Terrestrial accretion under oxidizing conditions. Science 339, 1194–1197.
Sleep, N.H., Zahnle, K.J., Neuhoff, P.S., 2001. Initiation of clement surface conditions on the earliest Earth. Proc. Natl. Acad. Sci. 98, 3666–3672.

Solomatov, V.S., 2007. Magma oceans and primordial mantle differentiation. In: Schubert, G. (Ed.)first ed. Treatise on Geophysics v. 9. Elsevier, pp. 91–119.

Stevenson, D.J., 1990. In: Newsom, H.E., Jones, J.H. (Eds.), Fluid dynamics of core formation, in the *Origin of the Earth*. Oxford University Press.

Studt, P.L., Shackelford, J.F., Fulrath, R.M., 1970. Solubility of gases in glass: a monoatomic model. J. Appl. Phys. 44, 2777–2780.

Tackley, P.J., 2012. Dynamics and evolution of the deep mantle resulting from thermal, chemical, phase and melting effects. Earth Sci. Rev. 110, 1–25.

Tateno, S., Hirose, K., Ohishi, Y., 2014. Melting experiments on peridotite to lowermost mantle conditions. J. Geophys. Res. Solid Earth 119, 4684–4694.

Terasaki, H., Frost, D.J., Rubie, D.C., Langenhorst, F., 2007. Interconnectivity of Fe–O–S liquid in polycrystalline silicate perovskite at lower mantle conditions. Phys. Earth Planet. Inter. 161, 170–176.

Thibault, Y., Walter, M.J., 1995. 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, 991–1002.

Thomas, C.W., Asimow, P.D., 2013. Direct shock compression experiments on premolten forsterite and progress toward a consistent high-pressure equation of state for CaO–MgO–Al₂O₃–SiO₂–FeO liquids. J. Geophys. Res. Solid Earth 118, 5738–5752.

Thomas, C.W., Liu, Q., Agee, C.B., Asimow, P.D., Lange, R.A., 2012. Multi-technique equation of state for Fe₂SiO₄ melt and the density of Fe-bearing silicate melts from 0 to 161 GPa. J. Geophys. Res. Solid Earth 117:B10206. http://dx.doi.org/10.1029/2012JB009403.

Tonks, W.B., Melosh, H.J., 1993. Magma ocean formation due to giant impacts. J. Geophys. Res. 98, 5319–5333.

Touboul, M., Puchtel, I.S., Walker, R.J., 2012. ¹⁸²W evidence for long-term preservation of early mantle differentiation products. Science 335, 1065–1069.

Touboul, M., Sprung, P., Aciego, S.M., Bourdon, B., Kleine, T., 2015. Hf–W chronology of the eucrite parent body. Geochim. Cosmochim. Acta 156, 106–121.

Trieloff, M., Kunz, J., 2005. Isotope systematics of noble gases in the Earth's mantle: possible sources of primordial isotopes and implications for mantle structure. Phys. Earth Planet. Inter. 148, 13–38.

Tsuno, K., Frost, D.J., Rubie, D.C., 2013. Simultaneous partitioning of silicon and oxygen into the Earth's core during early Earth differentiation. Geophys. Res. Lett. 40, 66–71.

Tucker, J.M., Mukhopadhyay, S., 2014. Evidence for multiple magma ocean outgassing and atmospheric loss episodes from mantle noble geases. Earth Planet. Sci. Lett. 393, 254–265.

Valley, J.W., et al., 2014. Hadean age for post-magma-ocean zircon confirmed by atom-probe tomography. Nat. Geosci. 7, 219–223.

Vidale, J.E., Hedlin, M.A.H., 1998. Evidence for partial melt at the core–mantle boundary north of Tonga from the strong scattering of seismic waves. Nature 391, 682–685.

Wacheul, J.B., Le Bars, M., Monteux, J., Aurnou, J.M., 2014. Laboratory experiments on the breakup of liquid metal diapirs. Earth Planet. Sci. Lett. 403, 236–245.

Wade, J., Wood, B.J., 2001. The Earth's "missing" niobium may be in the core. Nature 409, 75–78.

Wade, J., Wood, B.J., 2005. Core formation and the oxidation state of the Earth. Earth Planet. Sci. Lett. 236, 78–95.

Warren, P.H., 2011. Stable-isotopic anomalies and the accretionary assemblage of the Earth and Mars: a subordinate role for carbonaceous chondrites. Earth Planet. Sci. Lett. 311, 93–100.

Wetherill, G.W., 1985. Occurrence of giant impacts during the growth of terrestrial planets. Science 228, 877–879.

Wheeler, K.T., Walker, D., Fei, Y., Minarik, W.G., McDonough, W.F., 2006. Experimental partitioning of uranium between liquid iron sulfide and liquid silicate: implications for radioactivity in the Earth's core. Geochim. Cosmochim. Acta 70, 1537–1547.

White, B.S., Brearly, M., Montana, A., 1989. Solubility of argon in silicate liquids at high pressures. Am. Mineral. 74, 513–529.

Wicks, J.K., Jackson, J.M., Sturhahn, W., 2010. Very low sound velocities in iron-rich (Mg,Fe)O: Implications for the core–mantle boundary region. Geophys. Res. Lett. 37, L15304. http://dx.doi.org/10.1029/2010GL043689.

Wilde, S.A., Valley, J.W., Peck, W.H., Graham, C.M., 2001. Evidence from detrital zircons for the existence of continental crust and oceans on the Earth 4.4 Gyr ago. Nature 409, 175–178.

Williams, Q., Garnero, E.J., 1996. Seismic evidence for partial melt at the base of Earth's mantle. Science 273, 1528–1530.

Wood, B.J., Walter, M.J., Wade, J., 2006. Accretion of the Earth and segregation of its core. Nature 441, 825–833.

Wood, B.J., Li, J., Shahar, A., 2013. Carbon in the core: influence on the properties of core and mantle. Rev. Mineral. Geochem. 75, 231–250.

Yokochi, R., Marty, B., 2005. Geochemical constraints on mantle dynamics in the Hadean. Earth Planet. Sci. Lett. 238, 17–30.

Yoshino, T., Walter, M., Katsura, T., 2003. Core formation in planetesimals triggered by permeable flow. Nature 422, 154–157.

Zahnle, K.J., Kastings, J.F., Pollack, J.B., 1988. Evolution of a steam atmosphere during Earth's accretion. Icarus 74, 62–97.

Zhang, Y., Yin, Q.Z., 2012. Carbon and other light element contents in the Earth's core based on first-principles molecular dynamics. Proc. Natl. Acad. Sci. U. S. A. 109, 16579–16583.
# Dynamics and stability of an iron drop falling in a magma ocean

B. Qaddah<sup>a,b,*</sup>, J. Montieux<sup>a</sup>, V. Clesi<sup>a</sup>, M.A. Bouhifd<sup>a</sup>, M. Le Bars<sup>b</sup>

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

Check for
updates

## ARTICLE INFO

**Keywords:**

Core formation

Fluid dynamics

Numerical modeling

Metal drop

Magma ocean

## ABSTRACT

The latest stages of planetary accretion involved large impacts between differentiated bodies, hence large scale melting events. Consequently, the iron brought by the impactors sank within a deep magma ocean, before reaching the proto-core. Yet the fluid dynamics of this process remains poorly known. Here, we report numerical simulations of the sinking dynamics of an initially spherical liquid iron drop within a molten silicate phase, up to its possible fragmentation. We consider a 2D cylindrical axisymmetric geometry. We vary the viscosity of the molten silicates in the range of 0.05 Pa.s–100 Pa.s and the initial radius of the iron drop in the range of 1 mm–350 mm. Hence, we investigate Reynolds number in the range of [0.027–85600] and Weber number in the range of [0.073–7480]. Our numerical model constrains the morphology, dynamics and stability of the iron drop as a function of the dimensionless Weber and Reynolds numbers as well as of the viscosity ratio between the molten silicates and the liquid iron drop. In particular, we show that the maximal stable drop radius and the critical Weber number are monotonically increasing functions of the magma ocean viscosity. The momentum boundary layer thickness depends mainly on the drop radius and slightly on the magma ocean viscosity. Increasing the viscosity of the silicate phase prevents oscillations of the iron phase and limits the exchange surface. Oppositely, increasing the initial radius of the iron drop enhances its deformation and increases its relative exchange surface. Above the critical Weber number, we confirm that the fragmentation of the liquid iron occurs within a falling distance equal to 3.5–8 times the drop initial radius in the explored range of moderate Weber number, and we describe a variety of fragmentation regimes. Consequences for Earth’s formation models are briefly assessed.

## 1. Introduction

The Earth core formation is a complex process which remains actively debated theoretically, experimentally, and numerically (e.g. Stevenson, 1990; Tonks and Melosh, 1992; Chambers, 2004; O'Brien et al., 2006; Wood et al., 2006; Morbidelli et al., 2012; Deguen et al., 2014; Wacheul et al., 2014; Bouhifd et al., 2017). The core formation is highly dependent on accretion as both processes are contemporaneous (Kleine et al., 2002; Touboul et al., 2007). In the latest stages of planetary accretion, giant impacts have occurred [Fig. 1], leading for instance to the formation of the Moon following the collision between a large differentiated body of Mars size and the proto-Earth (Hartmann and Davis, 1975). The enormous amount of kinetic energy brought in by these collisions (Tonks and Melosh, 1992; Montieux et al., 2007; Samuel, 2012), the decay of short-lived radio-elements ($^{26}$Al and $^{60}$Fe) causing radioactive heating (Walter and Tronnes, 2004), and the heat dissipation from the conversion of potential energy (Montieux et al., 2009; Samuel, 2012) probably led to the formation of deep magma

oceans (Tonks and Melosh, 1992). Following the impact, the liquid iron from the impactor core spread and then sank into this less dense magma ocean as an immiscible fluid, leading to thermo-chemical exchanges between the two phases, before merging with the proto-core. This dynamical process may involve deformation of the initial drops coming from the primary breakup at impact, and possibly their breakup into even smaller droplets [Fig. 1] (Samuel, 2012).

Three scenarios have been proposed to characterise the motion of liquid iron within the magma ocean and the thermo-chemical equilibration between the two phases. In the first scenario, i.e. the so-called iron rain model (Stevenson, 1990; Karato and Rama Murthy, 1997; Rubie et al., 2003), it is suggested that an impactor's core with a diameter of [10–100] km immediately fragments into small droplets with a single characteristic diameter of the order of 1 cm, corresponding to the capillary size. All droplets descend independently towards the bottom of the magma ocean at the same velocity and without any further change in shape, leading to an efficient chemical equilibration at a distance of less than 200 m within the magma ocean. In contrast, Dahl

<footnote><sup>a</sup> Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, F-63000 Clermont-Ferrand, France</footnote>

<footnote><sup>b</sup> CNRS, Aix Marseille Univ, Centrale Marseille, IRPHE, Marseille, France</footnote>

<footnote>* Corresponding author.</footnote>

<footnote>E-mail address: qaddah@irphe.univ-mrs.fr (B. Qaddah).</footnote>
<figure><img src="image_81.png" /><figcaption>Fig. 1. Schematic of the metal/silicate separation during an impact between a differentiated planetesimal and the early Earth with a schematic of our computational domain.</figcaption></figure>

and Stevenson (2010) proposed a theoretical model where the large liquid iron impactor's core does not break up into the magma ocean if its initial size is large enough; it is only eroded by Rayleigh-Taylor and Kelvin-Helmholtz instabilities. Such a core formation episode would be far from thermo-chemical equilibrium between the impactor's core and the magma ocean. More recently, Deguen et al. (2011) and Deguen et al. (2014) used fluid dynamics experiments to show that a large core always breaks and forms a cloud of droplets with different sizes. But before doing so, it first evolves from its initial state as a turbulent thermal with strong entrainment and mixing, leading to rapid equilibration even before fragmentation.

In these previous studies, the influence of the viscosity ratio between the molten silicates and the liquid iron is neglected. Wacheul et al. (2014) developed analog experiments using gallium to mimic the liquid iron and a mixture of water and glycerol to mimic molten silicates. They found that the value of the viscosity ratio is very important and changes significantly the flow morphology. Later on, Wacheul and Le Bars (2018) validated the turbulent thermal model proposed by Deguen et al. (2014) using global temperature measurements of the liquid metal equilibration. However, numerous questions remain regarding the local dynamics at the scale of one drop, including: What is the maximum stable size of a drop as a function of the viscosity ratio? How is the exchange surface between iron drops and molten silicates influenced by the viscosity ratio? What is the size of the dynamical boundary layer across which the thermal and chemical exchanges occur? Answering those questions is fundamental for a better estimate of the equilibration length in planets, hence for relevant models of planet initial thermochemical state.

Several works have already investigated the dynamics of a buoyant bubble or drop in a viscous environment. For instance, a seminal experimental study of rising bubbles of hydrogen in aqueous sugar solutions is reported by Bhaga and Weber (1981) to determine the shape of the bubble, its terminal velocity, the geometry of its wake, and the flow streamline around it. Several experimental investigations of drop breakup in gas-liquid and liquid-liquid systems are reviewed by Pilch and Erdman (1987). An axisymmetric numerical model of a secondary breakup of a spherical liquid drop falling from rest due to gravity in another immiscible liquid is presented in Han and Tryggvason (1999) for different density and viscosity ratios, showing large drop deformation and various breakup modes. Another axisymmetric numerical investigation of the evolution of a large bubble rising by buoyancy in the presence of both capillary and viscous effects was carried out by Bonometti and Magnaudet (2006), in particular to determine the transition from a spherical cap to a toroidal shape. The influence of viscosity and density ratios is investigated systematically by Ohta et al.

(2009), Ohta et al. (2010), Ohta et al. (2014), Ohta and Sussman (2012). Finally, the first systematics on the asymmetric motion and the fragmentation modes of a bubble of gas rising in a liquid using three-dimensional numerical simulations are reported in Tripathi et al. (2015). However, previous works on the effects of viscosity contrast were motivated by industrial applications mainly, apart from 2 cases in Ichikawa et al. (2010) exploring the relative dynamical changes for viscosity ratios 0.1 and 10 between iron and a magma ocean. In the context of the core formation, the viscosity of the magma ocean strongly depends on its evolving temperature and pressure (Karki and Stixrude, 2010). Hence, the viscosity contrast between the metallic and silicate phases can differ by several orders of magnitude. We thus propose here to constrain the influence of the viscosity contrast between the metallic phase and the magma ocean focusing on the sinking of a unique iron droplet. We also investigate the influence of the droplet radius and derive scaling laws to characterise this fluid dynamics.

The study is organised as follows. Section 2 introduces the physical and numerical model, including the governing equations, the relevant non-dimensional parameters, the numerical model description, and our computational domain and mesh. Section 3 then presents our main numerical results including systematic studies of the drag coefficient, fragmentation modes, characteristic time and distance before breakup, and maximum stable radius before breakup. In Section 4, we discuss possible planetary applications, focusing on the potential efficiency of thermo-chemical exchanges around the drop which depends on the boundary layer thickness and the exchange surface between iron drop and molten silicates. Conclusions and future works are detailed in the final Section 5.

## 2. Physical and numerical model

### 2.1. Governing equations

We consider the sinking dynamics of an initially spherical iron drop, falling in a less dense and more viscous fluid under the action of gravity. We assume that both the liquid iron drop and the surrounding molten silicates behave as Newtonian, incompressible, and immiscible fluids with uniform surface tension, and constant density and viscosity within each fluid. We note here that we do not consider the thermo-chemical exchanges that may occur between the two phases. This point will be the subject of a separated study. Hence, the fluid dynamics is governed by the Navier-Stokes equations to characterize:

1. the conservation of mass

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$
## 2. the conservation of momentum

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = \nabla \cdot \left[ -P\mathbf{I} + \mu(\nabla \mathbf{u} + (\nabla \mathbf{u})^T) + \rho \mathbf{g} + \mathbf{F}_{st} \right] \quad (2) $$

with **u** the fluid velocity vector (m/s), ρ the fluid density (kg/m³), μ the fluid viscosity (Pa.s), t the time (s), P the fluid pressure (Pa), g the gravitational acceleration (m/s²), Fst the surface tension force (N/m³) and I the identity matrix.

We monitor the interface between the liquid iron drop and the molten silicates using the Level Set method, an Eulerian and implicit method frequently used in multiphase flow problems (e.g. Luo et al., 2006). It consists in defining a level set function φ, equal to 1 in iron and 0 in the surrounding silicates, and rapidly changing through the interface, whose position is determined by the isocontour φ = 0.5. The equation governing the transport and reinitialization of φ is:

$$ \frac{\partial \phi}{\partial t} + \mathbf{u} \cdot \nabla \phi = \gamma \nabla \cdot \left[ \epsilon \nabla \phi - \phi(1 - \phi) \frac{\nabla \phi}{|\nabla \phi|} \right] \quad (3) $$

with γ (m/s) and ϵ (m) the reinitialization parameters. ϵ determines the thickness of the layer around the interface, and is typically chosen equal to half the size of the characteristic mesh in the region explored by the interface. γ determines the amount of reinitialization: a suitable value for γ is the maximum velocity magnitude experienced in the model.

The density and dynamical viscosity are evaluated using the level set function:

$$ \rho = \rho_s + (\rho_d - \rho_s)\phi \quad (4) $$

$$ \mu = \mu_s + (\mu_d - \mu_s)\phi \quad (5) $$

where subscripts “s” and “d” stand for the molten silicates and the liquid iron drop respectively.

The surface tension force is determined by:

$$ \mathbf{F}_{st} = \nabla \cdot \mathbf{T} = \nabla \cdot (\sigma[\mathbf{I} + (-\mathbf{n}\mathbf{n}^T)]\Delta) \quad (6) $$

with σ (N/m) the surface tension coefficient, I the identity matrix, n the interface normal unit vector, and Δ the Dirac delta function, nonzero only at the fluid interface. The interface normal unit vector is calculated as

<figure><img src="image_82.png" /><figcaption>Fig. 2. Global view of our geometry with the initial drop at the top (left), and a zoom illustrating our manual method for adaptive mesh (right).</figcaption></figure>

$$ \mathbf{n} = \frac{\nabla \phi}{|\nabla \phi|} \quad (7) $$

The level set parameter φ is also used to approximate the delta function by a smooth function (Hu et al., 2014) defined by

$$ \Delta = 6|\phi(1 - \phi)|\nabla\phi. \quad (8) $$

## 2.2. Physical and non-dimensional parameters

According to the Buckingham Π theorem, in our study, there are 6 dimensionless numbers based on the control parameters (see Figs. 1 and 2). One could for instance choose the two aspect ratios related to the drop vs. the computational domain sizes, the density ratio, the viscosity ratio, and then define a priori Weber and Reynolds numbers based on the Newton theoretical velocity (see e.g. Wacheul et al., 2014). A complete dynamical study would then necessitate to explore the influence of those 4 later parameters in a computational box sufficiently large such that boundary conditions do not influence the dynamics. Here, we are interested in determining the dynamics of metal drops sinking into a magma ocean: we thus vary only the drop radius and the magma ocean viscosity, which are the main variables in the geophysical situation of interest, while we keep all the other dimensioned parameters fixed at their expected geophysical values (see Table 1).

In each simulation, starting from rest, the drop accelerates until reaching a constant terminal velocity, possibly with small oscillations around it. We monitor this mean terminal velocity V for each case. We end up our simulations when either the drop reaches a stable regime before the bottom of the domain or the drop fragments. In the later case the axisymmetric approximation is not relevant anymore. Each run is thus characterised by the 4 dimensionless numbers:

* the Reynolds number compares inertia and viscous effects: $Re = \frac{\rho_s V D}{\mu_s}$, with D = 2R the drop initial diameter (R is the initial radius). The Reynolds number determines the falling drop dynamical regime: Re < 1 implies a Stokes regime where the viscous effects dominate; Re = 1–500 implies an intermediate regime where both viscous effects and inertial forces are important; and Re > 500 implies a Newtonian regime where the inertial forces are dominant
<caption>**Table 1**</caption>

Symbol definitions and values of the physical and non-dimensional parameters used in this study. Reynolds and Weber numbers are defined a posteriori using the relevant terminal velocity, measured in each run.

<table><thead><tr><th></th><th>Symbol</th><th>Value or range</th></tr></thead><tbody><tr><td>Silicates density</td><td>ρ<sub>s</sub></td><td>3500 kg/m<sup>3</sup></td></tr><tr><td>Iron drop density</td><td>ρ<sub>d</sub></td><td>7500 kg/m<sup>3</sup></td></tr><tr><td>Density ratio</td><td>R<sub>ρ</sub></td><td>2.14</td></tr><tr><td>Iron drop viscosity</td><td>μ<sub>d</sub></td><td>0.005 Pa.s</td></tr><tr><td>Magma ocean viscosity</td><td>μ<sub>s</sub></td><td>0.05–100 Pa.s</td></tr><tr><td>Viscosity ratio</td><td>R<sub>μ</sub></td><td>10–20,000</td></tr><tr><td>Initial drop radius</td><td>R</td><td>1–350 mm</td></tr><tr><td>Surface tension coefficient</td><td>σ</td><td>1 N/m</td></tr><tr><td>Reynolds number</td><td>Re</td><td>0.027–85,600</td></tr><tr><td>Weber number</td><td>We</td><td>0.073–7480</td></tr></tbody></table>

(see Samuel, 2012).

* the Weber number measures the relative importance of inertia over surface tension: $We = \frac{\rho_s V^2 D}{\sigma}$. It controls and governs the deformation and breakup of the drop (Pilch and Erdman, 1987). Indeed inertia forces aim at deforming and fragmenting the iron drop, while the surface tension force prevents it from deformation and disruption. When $We \gg 1$, the inertia forces dominate leading to strong deformation and rapid fragmentation (e.g. Wacheul et al., 2014).

* the viscosity ratio is defined as the ratio of the silicates viscosity over the iron drop viscosity: $R_μ = \frac{\mu_s}{\mu_d}$. The viscosity of liquid iron weakly depends on pressure and ~ 5 × 10⁻³ Pa.s at both inner core boundary pressure conditions (Poirier, 1988) and at the Earth surface (Assael et al., 2006). Oppositely, the magma ocean viscosity varies between 10⁻⁴ and 100 Pa.s depending on the thermal and dynamic conditions (Karki and Stixrude, 2010; Rubie et al., 2003; Samuel, 2012). In our model, we fix the viscosity of the metallic phase and the viscosity contrast ranges between 10 and 20000.

* the density ratio compares the iron drop density over the silicates density: $R_ρ = \frac{ρ_d}{ρ_s}$. The densities of liquid iron and molten silicates depend on the temperature and pressure conditions (Assael et al., 2006; de Wijs et al., 1998; Samuel, 2012). In our models, for simplicity, we fix $ρ_d = 7500$ kg.m⁻³ and $ρ_s = 3500$ kg.m⁻³. As a consequence, in this study, the density ratio is constant and equals to 2.14.

After a large collision involving two differentiated protoplanets, the core of the impactor is first fragmented in large scale iron structures (Kendall and Melosh, 2016). In the magma ocean, these structures overcome a second fragmentation resulting in the formation of a cloud of droplets with typical size ranging between a few millimetres to several centimetres (Deguen et al., 2014; Wacheul et al., 2014). In our study, we characterize the dynamics of a droplet in the iron cloud and determine the influence of the viscosity contrast on this droplet. For that, we vary the initial radius of the drop from 1 to 350 mm, exploring a large range of expected sizes in the iron rain cloud.

All the parameters used in this paper are listed in Table 1.

## 2.3. Numerical method

In this work, we solve Eqs. (1)–(3) using the COMSOL Multiphysics software, based on the finite element method. We study a 2D axisymmetric geometry with no-slip boundary conditions at the lateral boundary and open boundary conditions at the top and bottom boundaries. To avoid the wall effects, the computational domain must be large enough. For instance, Bonometti and Magnaudet (2006) considered a computational domain of 5.5D × 12.6D to avoid contamination of the results. Their computations are stopped before the bubble arrives too close to the domain boundary. Ohta et al. (2010) used a computational domain of 4D × 6D, and considered a levitating drop

with an entering inflow at the top boundary equal to its terminal velocity. Samuel (2012) also considered a levitating drop but with a square computational domain of 5D × 5D, and the inflow velocity is self-adjusting at the instantaneous velocity of the center of mass of the drop. Here, we have checked the convergence of our results as a function of the domain size, and we do not see any significant difference once the domain is larger than 4D and longer than 75D. As sketched in Figs. 1 and 2, our computational domain is an axisymmetric cylinder of size ($r \times z$) = (6D × 100D).

The drop dynamics requires a fine mesh to capture all its details while conserving the overall mass of iron. Especially, the level set method has a tendency to diffuse numerical artifacts on the interface and then gives incorrect results. Since our computational domain is very large, we need an adaptive mesh strongly refined in the drop vicinity, which we have implemented manually. To do so, we divide our domain in several regions, as shown in Fig. 2. We mesh the region that the drop crosses over a given time interval with a very fine mesh of size $h_x = 0.025R$, small enough to capture all the dynamics on the interface without any significant error or numerical diffusion. We mesh the region behind the drop with a thin mesh $h_x = 0.07R$ to correctly capture the physics of the wake, which also influences the dynamics. The mesh around the drop gradually widens in the radial direction until reaching $h_x = 0.2R$, from where it remains fixed until the wall. Finally, our mesh in front of the drop in the vertical direction first increases in an intermediate region to $h_x = 0.07R$, and then keeps increasing until reaching a very large value $h_x = 5000R$. When the drop approaches the bottom of the finer mesh region, the simulation is stopped, the whole pattern is translated, and the simulation is restarted on this new grid. We have tested this procedure and found that the dynamics of the falling drop doesn't change significantly when using smaller $h_x$ (see Appendix). Also the overall iron mass during the course of each complete computation does not change by more than 0.7%, hence showing good numerical convergence. This method is relatively inexpensive in term of calculation time, and allows for a systematic study, with runs taking from 7 h to one week on a bi-processor, 3.2–3.6 GHZ, computer.

## 3. Systematic numerical study

We performed 84 simulations, whose dimensionless parameters are given in Table 2. In the following sections, we discuss our obtained numerical results in terms of drag coefficient, fragmentation modes, time and distance before breakup, and maximum stable drop radius. A particular attention is paid to the influence of the viscosity contrast between the metallic and silicate phases. The influence of the initial drop shape is shortly addressed at the end of this section; everywhere else, we start from a spherical drop.

### 3.1. Drag coefficient

The drag coefficient is a dimensionless number that quantifies the drag or resistance of the sinking drop into the magma ocean. It strongly depends on the viscosity ratio, density ratio, and surface tension, all of which control the interface conditions. In our case where the driving force is buoyancy and the initial geometry is spherical, we classically define the drag coefficient as

$$C_D = \frac{8 (\rho_d - \rho_s) g R}{3 \rho_s V^2}, \qquad (9)$$

with V the measured terminal velocity. Below, we first rapidly review the different theoretical models for the drag coefficient, and then compare with our numerical measurements.

#### 3.1.1. For small Reynolds numbers

Several previous studies focused on the analytical determination of the drag coefficient for spherical fluid particles (e.g. Clift and Gauvin,
<caption>**Table 2**</caption>

Non-dimensional parameters for all performed simulations used in this study. In all our simulations, $R_ρ = 2.14$.

<table><caption>**Table 2 (continued)**</caption><thead><tr><th>Simulation</th><th>Re</th><th>We</th><th>R<sub>μ</sub></th><th>R(mm)</th></tr></thead><tbody><tr><td>#1</td><td>19.7</td><td>0.139</td><td>10</td><td>1</td></tr><tr><td>#2</td><td>315</td><td>7.09</td><td>10</td><td>5</td></tr><tr><td>#3</td><td>372</td><td>8.22</td><td>10</td><td>6</td></tr><tr><td>#4</td><td>426</td><td>9.27</td><td>10</td><td>7</td></tr><tr><td>#5</td><td>553</td><td>10.9</td><td>10</td><td>10</td></tr><tr><td>#6</td><td>1610</td><td>37.0</td><td>10</td><td>25</td></tr><tr><td>#7</td><td>3280</td><td>95.8</td><td>10</td><td>40</td></tr><tr><td>#8</td><td>73.9</td><td>8.13</td><td>50</td><td>6</td></tr><tr><td>#9</td><td>93.0</td><td>9.64</td><td>50</td><td>8</td></tr><tr><td>#10</td><td>106</td><td>11.1</td><td>50</td><td>9</td></tr><tr><td>#11</td><td>108</td><td>10.4</td><td>50</td><td>10</td></tr><tr><td>#12</td><td>202</td><td>24.2</td><td>50</td><td>15</td></tr><tr><td>#13</td><td>322</td><td>37.0</td><td>50</td><td>25</td></tr><tr><td>#14</td><td>638</td><td>91.0</td><td>50</td><td>40</td></tr><tr><td>#15</td><td>1510</td><td>290</td><td>50</td><td>70</td></tr><tr><td>#16</td><td>14.3</td><td>1.82</td><td>100</td><td>6</td></tr><tr><td>#17</td><td>35.3</td><td>7.43</td><td>100</td><td>8</td></tr><tr><td>#18</td><td>44.2</td><td>8.74</td><td>100</td><td>10</td></tr><tr><td>#19</td><td>51.8</td><td>9.58</td><td>100</td><td>11</td></tr><tr><td>#20</td><td>66.0</td><td>13.5</td><td>100</td><td>11.5</td></tr><tr><td>#21</td><td>71.4</td><td>15.2</td><td>100</td><td>12</td></tr><tr><td>#22</td><td>92.4</td><td>20.3</td><td>100</td><td>15</td></tr><tr><td>#23</td><td>154</td><td>33.9</td><td>100</td><td>25</td></tr><tr><td>#24</td><td>314</td><td>87.8</td><td>100</td><td>40</td></tr><tr><td>#25</td><td>750</td><td>287</td><td>100</td><td>70</td></tr><tr><td>#26</td><td>11.9</td><td>3.35</td><td>200</td><td>6</td></tr><tr><td>#27</td><td>19.7</td><td>6.90</td><td>200</td><td>8</td></tr><tr><td>#28</td><td>25.1</td><td>9.02</td><td>200</td><td>10</td></tr><tr><td>#29</td><td>30.4</td><td>11.0</td><td>200</td><td>12</td></tr><tr><td>#30</td><td>36.4</td><td>14.6</td><td>200</td><td>13</td></tr><tr><td>#31</td><td>38.7</td><td>15.3</td><td>200</td><td>14</td></tr><tr><td>#32</td><td>42.5</td><td>17.2</td><td>200</td><td>15</td></tr><tr><td>#33</td><td>76.1</td><td>33.1</td><td>200</td><td>25</td></tr><tr><td>#34</td><td>146</td><td>75.7</td><td>200</td><td>40</td></tr><tr><td>#35</td><td>363</td><td>268</td><td>200</td><td>70</td></tr><tr><td>#36</td><td>637</td><td>580</td><td>200</td><td>100</td></tr><tr><td>#37</td><td>1.53</td><td>1.05</td><td>1000</td><td>8</td></tr><tr><td>#38</td><td>2.70</td><td>2.61</td><td>1000</td><td>10</td></tr><tr><td>#39</td><td>6.49</td><td>10.0</td><td>1000</td><td>15</td></tr><tr><td>#40</td><td>11.2</td><td>22.4</td><td>1000</td><td>20</td></tr><tr><td>#41</td><td>15.4</td><td>33.9</td><td>1000</td><td>25</td></tr><tr><td>#42</td><td>19.7</td><td>46.4</td><td>1000</td><td>30</td></tr><tr><td>#43</td><td>26.0</td><td>68.8</td><td>1000</td><td>35</td></tr><tr><td>#44</td><td>28.8</td><td>74.3</td><td>1000</td><td>40</td></tr><tr><td>#45</td><td>38.5</td><td>106</td><td>1000</td><td>50</td></tr><tr><td>#46</td><td>52.1</td><td>161</td><td>1000</td><td>60</td></tr><tr><td>#47</td><td>68.6</td><td>240</td><td>1000</td><td>70</td></tr><tr><td>#48</td><td>122</td><td>530</td><td>1000</td><td>100</td></tr><tr><td>#49</td><td>0.027</td><td>0.073</td><td>20000</td><td>15</td></tr><tr><td>#50</td><td>0.127</td><td>0.927</td><td>20000</td><td>25</td></tr><tr><td>#51</td><td>0.509</td><td>9.27</td><td>20000</td><td>40</td></tr><tr><td>#52</td><td>0.952</td><td>25.9</td><td>20000</td><td>50</td></tr><tr><td>#53</td><td>1.22</td><td>38.7</td><td>20000</td><td>55</td></tr><tr><td>#54</td><td>1.54</td><td>56.6</td><td>20000</td><td>60</td></tr><tr><td>#55</td><td>2.27</td><td>105</td><td>20000</td><td>70</td></tr><tr><td>#56</td><td>5.15</td><td>379</td><td>20000</td><td>100</td></tr><tr><td>#57</td><td>24.4</td><td>26.5</td><td>500</td><td>20</td></tr><tr><td>#58</td><td>7.56</td><td>42.5</td><td>2500</td><td>30</td></tr><tr><td>#59</td><td>8.42</td><td>49.5</td><td>2500</td><td>35</td></tr><tr><td>#60</td><td>9.82</td><td>61.2</td><td>2500</td><td>40</td></tr><tr><td>#61</td><td>12.3</td><td>84.7</td><td>2500</td><td>45</td></tr><tr><td>#62</td><td>2670</td><td>636</td><td>50</td><td>100</td></tr><tr><td>#63</td><td>24100</td><td>1380</td><td>10</td><td>150</td></tr><tr><td>#64</td><td>37100</td><td>2460</td><td>10</td><td>200</td></tr><tr><td>#65</td><td>67800</td><td>5480</td><td>10</td><td>300</td></tr><tr><td>#66</td><td>85600</td><td>7480</td><td>10</td><td>350</td></tr><tr><td>#67</td><td>4770</td><td>1350</td><td>50</td><td>150</td></tr><tr><td>#68</td><td>7360</td><td>2420</td><td>50</td><td>200</td></tr><tr><td>#69</td><td>13400</td><td>5380</td><td>50</td><td>300</td></tr><tr><td>#70</td><td>16960</td><td>7330</td><td>50</td><td>350</td></tr><tr><td>#71</td><td>2380</td><td>1350</td><td>100</td><td>150</td></tr><tr><td>#72</td><td>3660</td><td>2400</td><td>100</td><td>200</td></tr><tr><td>#73</td><td>6700</td><td>5340</td><td>100</td><td>300</td></tr></tbody></table>

<table><thead><tr><th>Simulation</th><th>Re</th><th>We</th><th>R<sub>μ</sub></th><th>R(mm)</th></tr></thead><tbody><tr><td>#74</td><td>8430</td><td>7290</td><td>100</td><td>350</td></tr><tr><td>#75</td><td>1180</td><td>1330</td><td>200</td><td>150</td></tr><tr><td>#76</td><td>1820</td><td>2370</td><td>200</td><td>200</td></tr><tr><td>#77</td><td>3360</td><td>5360</td><td>200</td><td>300</td></tr><tr><td>#78</td><td>4221.35</td><td>7273.386</td><td>200</td><td>350</td></tr><tr><td>#79</td><td>231</td><td>1270.5</td><td>1000</td><td>150</td></tr><tr><td>#80</td><td>502.25</td><td>3603.643</td><td>1000</td><td>250</td></tr><tr><td>#81</td><td>667.8</td><td>5309.01</td><td>1000</td><td>300</td></tr><tr><td>#82</td><td>837.9</td><td>7164.045</td><td>1000</td><td>350</td></tr><tr><td>#83</td><td>28.35</td><td>17.22</td><td>300</td><td>15</td></tr><tr><td>#84</td><td>24.57</td><td>19.2</td><td>400</td><td>18</td></tr></tbody></table>

1970; Hadamard, 1911; Mei et al., 1994) and for deformed fluid particles (e.g. Moore, 1965; Darton and Harrison, 1974; Clift et al., 1978; Loth, 2008). For spherical drops, the internal circulation is supposed to prevent forming any wake separation of the external flow. This in turn helps to prevent any change in shape and keeps the droplet spherical. Hadamard (1911) proposed an analytical relation at finite but low Reynolds Re ≪ 1. An asymptotic solution was derived by Harper and Moore (1968) for intermediate Reynolds numbers but remains limited to Re < 100 at finite values of the viscosity ratio. To fill the gap between Harper and Moore (1968) and Hadamard (1911), an empirical mixed Stokes correction factor resulting from the internal circulation was proposed by Mei et al. (1994) in the limit $R_μ → 0$:

$$C_{D_{R_μ→0}} = \frac{16}{Re} \left[ 1 + \left[ \frac{8}{Re} + \frac{1}{2} \left( 1 + \frac{3.315}{\sqrt{Re}} \right) \right]^{-1} \right]. \quad (10)$$

### 3.1.2. For moderate and large Reynolds numbers

For intermediate but increasing Reynolds numbers ($0.1 < Re < 2000$), the drop begins to distort, leading to changes in the wake from an attached laminar wake (spherical drop), to a separated laminar wake (deformable drop), to an unsteady transitional wake (breakup into droplets) and finally to a turbulent wake (catastrophic breakup) (Loth, 2008). The wake actually depends on the interplay between surface tension and hydrodynamic pressure stresses, hence on We. Davies and Taylor (1950) found a converged value $C_D = \frac{8}{3}$ for gas bubbles in a liquid at infinite Reynolds and Weber numbers. In order to encompass the influence of viscosity and a wide range of small to large Reynolds numbers, an empirical correlation was proposed by Darton and Harrison (1974) and Clift et al., 1978 for infinitely large We:

$$C_{D_{We→∞}} = \frac{8}{3} + \frac{24}{Re} \left( \frac{2 + 3R_{\mu}}{3 + 3R_{\mu}} \right). \quad (11)$$

To capture the drag in the intermediate range of Re and We values (of interest here), Loth (2008) proposed the following expression by combining Eqs. (10) and (11) and adding a functional dependence ($\Delta(C_D)^*$) on We, where $\Delta(C_D)^*$ increases monotonically with Weber number:

$$C_D = C_{DWe-0} + \Delta(C_D)^* [C_{DWe-∞} - C_{DWe-0}] \quad (12)$$

$$\Delta(C_D)^* = \tanh(0.021We^{1.6}). \quad (13)$$

Here $C_{DWe-0}$ is given by Eq. (10), the imposed spherical shape being related to an infinite surface tension, hence to We → 0. Experimentally, several investigations partly tackled the relevant limit of a metal drop in a more viscous environment, including: the fragmentation of liquid mercury drops sinking within water reported by Patel (1978) and Patel and Theofanous (1981), who found a drag coefficient equal to 2.5–3; the fragmentation of gallium drops in water by Kim et al. (1983), who found that the drag coefficient depends on the Reynolds number; and the analog model of Wacheul and Le Bars (2018), who reported a convergence value of the drag coefficient at large Reynolds equal to
<figure><img src="image_83.png" /><figcaption>Fig. 3. Drag coefficient as a function of the Reynolds number for various values of the silicates viscosity. Filled symbols are the analytical results given by Eq. (12), while empty symbols show our numerical results from the first 62 simulations in Table 2. Filled black symbols show numerical results for an un-deformable sphere $R_μ = 0$, and black stars are the analytical results from Eq. (10).</figcaption></figure>

3.7 ± 1.

### 3.1.3. Results and discussion

We have monitored the velocity $V$ of the sinking metallic drop in our simulations and obtained the corresponding values for $C_D$ from Eq. (9). Fig. 3 presents our measured $C_D$ as a function of the Reynolds number for various viscosity ratios, and the comparison of our results with the analytical solutions of Mei et al. (1994) (Eq. (10)) and Loth (2008) (Eq. (12)). For $We < 1$ and low Reynolds numbers (up to $Re = 20$ depending on $R_μ$), the droplet remains spherical and our results agree with the analytical model proposed by Mei et al. (1994): $C_D$ is inversely proportional to $Re$ and barely depends on the viscosity ratio. When flow separation occurs (i.e. when the Weber number exceeds unity), this expression is not valid anymore because the drop deforms: it takes the shape of an ellipsoid, a disc or a cup. This deformation increases the frontal area and leads to a transient increase of $C_D$ with $Re$, depending on $R_μ$. Qualitative agreement between our simulations and Eq. (12) is satisfactory, but small quantitative differences exist. This is not surprising, since the empirical formula Eq. (12) was calibrated on gas bubbles rising in a liquid, which are more sensitive to deformation than the liquid drops studied here. But no equivalent to Eq. (12) is yet available for liquid-liquid systems, and one has thus to rely on numerical results. Finally, for large $Re$, $C_D$ converges towards a constant value $3.5 \pm 0.5$, independently of $R_μ$. This value is above the analytical estimate of 8/3, but is compatible with the results of Wacheul and Le Bars (2018) who found a mean value of $C_D$ equal to $3.7 \pm 1$.

### 3.2. Fragmentation modes

The breakup mechanism is very sensitive to the Weber number, as well as to the viscosity ratio between the metal and the silicates. Two main fragmentation modes have been documented in previous experimental and numerical studies, namely the “bag breakup” just above the critical Weber number, and the “shear breakup” at large Weber number (see the limit cases in Fig. 4) (Krzeczkowski, 1980; Pilch and Erdman, 1987; Dai and Faeth, 2001). Several authors have tackled the description of the different breakup mechanisms in the intermediate range of Weber number, leading to a complex situation. Examples include (Krzeczkowski, 1980), who developed series of experiments for liquid droplets of water, methanol, ethanol, butanol and glycerine in an external air stream, and documented two breakup modes independent of the viscosity ratio: the “bag-jet” and the “transition” modes. Pilch and

<figure><img src="image_84.png" /><figcaption>Fig. 4. Regime diagram of the fragmentation modes as a function of $Re$ and $We$ numbers for various viscosity ratios: $R_μ = 10$ (blue dotted line), $R_μ = 50$ (green dotted line), $R_μ = 100$ (black dotted line), $R_μ = 200$ (orange dotted line), $R_μ = 1000$ (red dotted line), and $R_μ = 2500$ (turquoise dotted line). Our numerical results come from the first 62 simulations in Table 2. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

Erdman (1987) introduced two breakup modes for a single liquid drop within an external gas flow: the “bag-stamen” mode and the “sheet stripping” mode. Dai and Faeth (2001) used droplets of water and ethanol and defined a “bag-plume” mode, close to the “bag-jet” and “bag-stamen”, and a “plume-shear” mode, close to the “shear breakup”. Cao et al. (2007) pursued this experimental investigation and presented a new “dual-bag” breakup mode for $We = 28–41$. Recently, 3D numerical simulations of liquid droplets levitating in a gas flow with a uniform velocity have been performed by Kékesi et al. (2014), changing the viscosity and density ratios for a fixed Weber number value ($We = 20$). They identified 5 breakup regimes, as shown in Fig. 4: “thick rim shear” and “thick rim bag” (both close to “bag breakup”), “rim shear”, “jellyfish shear”, and “shear” breakup respectively. In our study, we use their classification.

For each mode of fragmentation, we describe below the temporal evolution of a typical spherical drop starting from rest at $t = 0$, until its breakup at $t = t_{bk}$.

* Thick rim shear (Fig. 5): the initial sphere rapidly deforms into a spherical cap ($t = 0.13$ s), then the hydrodynamic pressure force concentrates on the middle of the drop leading to the formation of a small bag ($t = 0.2$ s). Surface tension then prevents breakup and the drop returns to form a half-sphere ($t = 0.355$ s), before a second oscillation starts. Then, as the inertial forces and associated pressure are reinforced, the vertical thickness decreases even more, and a thin film forms near the symmetry axis, with a thicker rim ($t = 0.55$ s). The rim extends radially and drains out the liquid film ($t = 0.6$ s), until a hole occurs at the symmetry axis ($t = 0.62$ s). At this stage, the metal phase actually forms a ring that might persist for a longer time before fragmentation of Rayleigh-Plateau type, to which we cannot have access with our axisymmetric simulation. Kékesi et al. (2014) showed in 3D simulations that this last stage is actually rapid, i.e. a few milliseconds. Note also that thick rim shear seems marginal in the parameter space explored here, and is replaced rapidly by thick rim bag, except for cases with the smaller viscosity ratio $R_μ = 10$ (see Fig. 4).
<figure><img src="image_85.png" /><figcaption>Fig. 5. Deformation and thick rim shear breakup of an iron drop within a magma ocean. In this model, Re = 36.4, We = 14.6, Rµ = 200 (simulation #30 in Table 2).</figcaption></figure>

<figure><img src="image_86.png" /><figcaption>Fig. 6. Deformation and thick rim bag breakup of an iron drop within a magma ocean. In this model, Re = 42.5, We = 17.2, Rµ = 200 (simulation #32 in Table 2).</figcaption></figure>

* Thick rim bag (Fig. 6): this breakup process is similar to the thick rim shear, until fragmentation. Then, breakup appears almost simultaneously in the center of the drop ($t = 0.675$ s) and at the connection between the film and the rim, forming satellite droplets ($t = 0.68$ s).

* Rim shear (Fig. 7): again the initiation of the breakup is similar to

the two processes described above, with the formation of a spherical cap ($t = 0.131$ s), followed by a bag ($t = 0.2$ s), and back ($t = 0.37$ s). Then, the drop extends in both radial and vertical directions, forming a sheet ($t = 0.48$ s and $0.51$ s) under the action of the strong vorticity in the wake. Finally, surface tension drains the iron both into the center and into the rim ($t = 0.57$ s), while the sheet

<figure><img src="image_87.png" /><figcaption>Fig. 7. Deformation and rim shear breakup of an iron drop within a magma ocean. In this model, Re = 92.4, We = 20.3, Rµ = 100 (simulation #22 in Table 2).</figcaption></figure>
<figure><img src="image_88.png" /><figcaption>Fig. 8. Deformation and jellyfish shear breakup of an iron drop within a magma ocean. In this model, Re = 76.1, We = 33.1, Rμ = 200 (simulation #33 in Table 2).</figcaption></figure>

connecting the two thins, and finally breaks ($t = 0.59$ s).

* Jellyfish shear (Fig. 8): in this mode, no rim forms. Instead, a thin skirt forms while iron is mostly localized along the symmetry axis ($t = 0.4$ s and $t = 0.45$ s). Oscillations develop in the thin membrane, whose ends fold in and out. The velocity is concentrated in the wake of the silicates which accelerates the drop center relative to its membranes, leading to a rapid fragmentation in the connecting region ($t = 0.51$ s).

* Shear breakup (Fig. 9): the drop shape deforms first into a spherical cap ($t = 0.24$ s) and second into a skirt. Then, the rim does not thicken, but folds inside and turns around itself through the action of vorticity, leading to its thinning ($t = 0.4$ s). Rapidly, this thin sheet breaks into small droplets ($t = 0.42$ s). No oscillations are observed here, and the rapid fragmentation is mainly due to the vorticity in the wake behind the drop that is two to three times larger than the vorticity in the drop.

The deformation of the droplet depends on the viscosity of the external flow, the density ratio, the surface tension and the drop size. When Re and We are both small, the drop remains spherical without any deformation, due to the internal circulation within the drop that prevents forming any separated wake. When the hydrodynamic pressure force increases, a separated laminar wake and an external circulation behind the drop occur, leading to shape deformation. Increasing the Reynolds number, the drop first deforms to an ellipsoid, and possibly

<figure><img src="image_89.png" /><figcaption>Fig. 9. Deformation and shear breakup of an iron drop within a magma ocean. In this model, Re = 2670, We = 636, Rμ = 50 (simulation #62 in Table 2).</figcaption></figure>

<figure><img src="image_90.png" /><figcaption>Fig. 10. Critical Weber number as a function of the viscosity ratio. The blue diamond is the critical Weber number when the viscosity ratio equal to 1, as obtained by Villermaux and Bossa (2009). The red and blue lines represent the scaling laws derived from our data. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

gives rise to oscillations from oblate to prolate shapes. For a viscosity ratio less than 100, the drop keeps oscillating until the end of the domain, while for a viscosity ratio greater than 100, the velocity of the drop exceeds the oscillation velocity, leading to rapid damping of the initially excited oscillation. Further increasing the Reynolds number, the drop takes the shape of a spherical cap. And once hydrodynamic pressure overcomes surface tension (i.e. large enough Weber number), one of the five fragmentation modes occurs, depending also on the viscosity ratio. From our systematic study for different viscosity ratios, we built up a fragmentation regime diagram in terms of Reynolds and Weber numbers shown in Fig. 4. For a viscosity ratio $R_μ = 1000$, we observe only two modes, and drop oscillations are prevented by the high silicates viscosity. On the contrary for $R_μ = 50$, we observe the five different modes and strong drop oscillations. One should also notice that the critical Weber number depends on the viscosity ratio: for $R_μ = 10$, the breakup begins at $We \approx 9$ while for $R_μ = 1000$, it starts at $We \approx 34$. Systematic estimate of the critical Weber number is shown in Fig. 10. We consider in this figure that fragmentation occurs as soon as the first volume of iron separates from the main drop. The increase in the viscosity ratio leads to an increase in the critical Weber number following two empirical scaling laws depending on first-breakup regime. When the viscosity ratio is rather low (i.e. <300 typically), the drop fragments quickly as a thick rim shear after only one or two oscillations along a short drop path. In this case, critical Weber is defined as a low power law of the viscosity ratio by the following scaling law

$$We_c = 5.7 R_{\mu}^{0.187} \qquad (14)$$

On the other hand, at a viscosity ratio of ~300, a Jellyfish behaviour takes place, and between $R_μ = 300$ and $R_μ = 500$, the drop first breakup becomes slow and transitions from thick rim shear to Jellyfish. In this region and above, where shear breakup takes place, the critical Weber number depends more significantly on the viscosity ratio and is given by the following scaling law

$$We_c = 1.137 R_{\mu}^{0.483} \qquad (15)$$

valid at least in the range $R_μ = [300;2500]$ explored here.

From our simulation #41, we also report a new fragmentation mode illustrated in Fig. 11. This mode is initially similar to the jellyfish shear breakup except that the membrane fold is so important that it finally closes at the back, encapsulating the silicates within a bubble of iron (Fig. 11, $t = 3.7$ s). Then, the next oscillation of the membrane tears the
<figure><img src="image_91.png" /><figcaption>Fig. 11. The path of a new fragmentation mode where the iron structure encloses the silicates.</figcaption></figure>

glued back which breaks from the center, forming one central droplet and two satellites droplets. Note that Wacheul et al. (2014) observed in their experiments similar bubbles of metal enclosing the viscous ambiant fluid and exhibiting a lower velocity that pure iron droplets. While being dynamically intriguing, this mode remains very marginal.

### 3.3. Time and distance before breakup

After an impact, large drops deform and fragment at a distance and time that depend on the sinking dynamics. We define these time and distance as the breakup time $t_{bk}$ and distance $d_{bk}$. Deguen et al. (2014) predicted that for large Weber numbers, $t_{bk}$ and $d_{bk}$ reach an asymptotic regime. In this section we test this prediction. From our models, we determine $t_{bk}$ and $d_{bk}$ as the time and location where the first droplet or ligament of iron separates from the main drop. We then define the dimensionless breakup time $t_{bk}^*$ as Pilch and Erdman (1987):

$$t_{bk}^* = \frac{t_{bk} V}{D} \sqrt{\frac{\rho_d}{\rho_s}} \qquad (16)$$

where the break up time is normalized by the characteristic time of drop breakup by Kelvin-Helmholtz instability. We define the dimensionless breakup distance as:

$$d_{bk}^* = \frac{d_{bk}}{R} \qquad (17)$$

Fig. 12 shows the variations of $t_{bk}^*$ as a function of the Weber number for various viscosity ratios. From Fig. 12, we see that $t_{bk}^*$ decreases with the increase of We. We note that, for large We, $t_{bk}^*$ converges towards a constant value ranging between 1 and 1.8. This result is consistent with

<figure><img src="image_92.png" /><figcaption>Fig. 12. Non-dimensional break-up time as a function of Weber number for viscosity ratios $R_μ$ ranging between 10 and 1000.</figcaption></figure>

<figure><img src="image_93.png" /><figcaption>Fig. 13. Non-dimensional break-up distance $d_{bk}^*$ as a function of the Weber number for viscosity ratios $R_μ$ ranging between 10 and 1000.</figcaption></figure>

the result of Pilch and Erdman (1987) who found $t_{bk}^* = 1.25$ for We > 350. A dependence on the viscosity ratio is also present, especially at We ≤ 80, with a tendency for large $R_μ$ to stabilise the drops.

The non-dimensionalised breakup distance $d_{bk}^*$ is reported in Fig. 13. It also shows a dependence on the viscosity ratio, especially at low We, and a possible convergence towards an asymptotic value ranging between 2 and 4 at large We. For comparison, Deguen et al. (2014) found a breakup distance equal to 6–8 times the drop radius for We up to 3000 and $R_μ = 0.5$. Landeau et al. (2014) found a value equal to 3.5–7.5 times the radius for $25 < We < 1000$ and $0.45 < R_μ < 1.25$. Hence our results are in correct agreement with those experimental observations. This is all the more noticeable that Landeau et al. (2014) observed 3D turbulence in lab experiments at We > 600, effects that are missing from our axisymmetric simulations. Note finally that Wacheul and Le Bars (2018) found a breakup distance ($13 \pm 2)R$ for $10 < We < 50$ and $0.4 < R_μ < 700$, but their breakup criteria actually detected a “significant” and well-advanced breakup, thus implicitly leading to an overestimation of $d_{bk}^*$ compared to other results. Fig. 13 illustrates the influence of the viscosity ratio on the breakup distance. For low We, $d_{bk}^*$ ranges between 10 and 80 for the whole range of viscosity ratios used in our study. For large We, $d_{bk}^*$ ranges between 2 and 4.5. Increasing the viscosity ratio generically increases the breakup distance and as a consequence, the potential depth of equilibrium between iron drops and silicates.

### 3.4. Maximum stable drop radius

After an impact, the impactor's iron core disrupts into large scale drops (Kendall and Melosh, 2016). After this first fragmentation, a
<figure><img src="image_94.png" /><figcaption>Fig. 14. Stable drop radius as a function of plausible magma ocean viscosities. The green diamonds present our numerical results, the red squares correspond to the results of Rubie et al. (2003), the black circles display the results of Samuel (2012), the blue stars correspond the analytical results of the Eq. (18) and the black dash line presents the scaling law proposed by this study. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

second fragmentation occurs within the magma ocean resulting in the formation of a cloud of droplets (Deguen et al., 2014; Wacheul et al., 2014). At the end of this second fragmentation, iron droplets reach a maximum stable radius $R_{max}$. Following Wacheul et al. (2014) and using the Newtonian velocity scale accounting for the drag coefficient (Eq. (9)), $R_{max}$ is simply related to the critical Weber number by

$$R_{max} = \sqrt{\frac{We_c \frac{3C_D}{16} \sigma}{\Delta \rho g}} \quad (18)$$

Numerical data from Fig. 14 show the last drop radius for which deformation may occur but without fragmentation. According to Fig. 14, the maximal stable drop radius increases with the magma ocean viscosity. From our numerical data, we propose a scaling law to determine the maximum stable radius as a power function of magma ocean viscosity, valid for magma viscosities larger than ~ 0.1 Pa.s up to 100 Pa.s:

$$R_{max} = 0.014\mu_s^{0.32}(m). \quad (19)$$

Our results are in agreement with the analytical Eq. (18), confirming the self-consistency of our numerical results. Fig. 14 also shows the comparison of our results with the model of Rubie et al. (2003) and the model of Samuel (2012). Differences come from different values of $C_D$, where we use our effective numerical values rather than any theoretical estimation. It should also be noted that Rubie et al. (2003) and Samuel (2012) do not consider the dynamical influence of the viscosity contrast between metal and silicates in their models.

### 3.5. The influence of initial conditions

After an impact, the impactor's core is likely fragmented in metallic diapirs with a large range of shapes (spherical, oblate or prolate) (Kendall and Melosh, 2016). This deformation process is likely to affect the post-impact sinking dynamics of the metallic phase as we envision in our study by modifying our initial shape conditions. Bonometti and Magnaudet (2006) found that their final bubble topology and features vary with the initial shape of the bubble. Without trying to be exhaustive on this subject that would deserve a dedicated study on its own, we have thus investigated the behavior of a drop at the critical Weber number for two viscosity ratios ($R_μ$ = 200 and $R_μ$ = 1000),

<figure><img src="image_95.png" /><figcaption>Fig. 15. The final shape of the iron droplets for different initial conditions. Spherical case (left), prolate form (middle) and oblate form (right).</figcaption></figure>

starting from a spherical, oblate or prolate form with an ellipticity ranging between 0.59 and 0.82, maintaining the total volume of iron constant. Results for $R_μ$ = 1000 are shown in Fig. 15, reporting a significant sensitivity to initial conditions. For the initial oblate form, the drop does not fragment and converges to a skirt shape. On the contrary, the initial prolate shape disrupts more rapidly than the spherical case, the time and distance before breakup being respectively 27.4% and 28.3% smaller than the corresponding values for the spherical case. The shape is also different, with an increase of the surface in the prolate case of 10.9% compared to the spherical case. For the other viscosity ratio $R_μ$ = 200, we find that the fragmentation modes are identical for the three initial conditions, but differences on the order of 10% exist in the breakup time, distance before breakup, and surface. Among the three initial cases, the spherical initial condition always leads to a less rapid fragmentation and to a smaller surface.

## 4. Implications for the metal/silicate exchanges

In the present study, we focus on the dynamics of the two phase flow only, without considering neither the thermal evolution of the metallic droplet nor chemical exchanges that could occur between the metallic and silicate phases. Our results nevertheless enable to constrain two fundamental parameters for thermochemical equilibration, namely the thickness of the boundary layer at the iron/silicates interface, and the surface of the metallic droplet across which thermochemical exchanges are likely to occur. Those two points are addressed below. We then propose some first order estimates of the equilibration of an initially 10 km metallic core sinking and fragmenting within a magma ocean and compare our results with the classical iron-rain model (Rubie et al., 2003).

### 4.1. Boundary layer thickness

During the sinking of an iron drop, three boundary layers can form at the interface between the metalic and silicate phases: the dynamical, thermal and chemical boundary layers. The thickness of each boundary layer may have a strong influence on the thermo-chemical exchanges. In this study, we focus on the dynamical boundary layer $\Delta$. We thus consider a velocity profile along the normal to the drop interface at the front of the drop, as shown in Fig. 16: $\Delta$ is defined as the distance from the drop interface to the point where the velocity is equal to 10% the drop terminal velocity. In the case of a rigid sinking sphere, the thickness of the boundary layer scales as $Re^{-1/2}$ (Ulvrová et al., 2011). In our models, the metallic droplet is deformable and three parameters a priori control the flow, hence $\Delta$: $R_μ$, Re and We. Here, we seek a power law for the boundary layer thickness in terms of the two varied parameters in our study, namely the drop radius and the magma ocean viscosity.

Our numerical results are shown in Fig. 17. A good fit is found with the empirical scaling law

$$\Delta = 1.442R^{0.916}\mu_s^{0.062}(m). \quad (20)$$

Hence, the boundary layer thickness depends mainly on the drop radius and only slightly on the magma ocean viscosity. For a standard viscous boundary layer in the vicinity of a solid sphere, one would expect a scaling $\propto \mu_s^{0.5}$. This surprisingly low dependence on the ambient viscosity is probably due to the drop deformation and to the
<figure><img src="image_96.png" /><figcaption>**Fig. 16.** Schematic representation of the geometry used for the determination of the boundary layer thickness $\Delta$ at the drop surface: spherical case (top) and deformed drop (bottom).</figcaption></figure>

<figure><img src="image_97.png" /><figcaption>**Fig. 17.** Boundary layer thickness as a function of the varied parameters in this study (drop radius and magma ocean viscosity) for viscosity ratios $R_\mu$ ranging between 10 and 1000. Our numerical results come from the first 62 simulations in Table 2).</figcaption></figure>

<figure><img src="image_98.png" /><figcaption>**Fig. 18.** The streamlines around a large drop. In this model, $Re = 2670$, $We = 636$, $R_\mu = 50$ (simulation #62 in Table 2). The black lines are the streamlines, the red region represents the molten silicates and the blue region represents the deformable iron drop. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

associated internal circulations and wake. It presumably affects the equilibration. Note however that Eq. (20) is valid in the frontal section of the drop, but the boundary layer size is not constant along the drop surface. Also, the generated wake might encapsulate silicates that are entrained with the falling iron, as shown in Fig. 18 and observed numerically (Monteux and Arkani-Hamed, 2014) and experimentally (Fleck et al., 2018), for an even more viscous surrounding fluid. This also influences equilibration. Clearly, the equilibration process is complex, global, and deserves a study on its own.

## 4.2. Time evolution of the droplet surface

The exchange surface between two fluids qualifies the surface through which thermo-chemical transfers occur. As shown in the Figs. 5–9 and 11, the spherical drop deforms and oscillates during its sinking. This leads to variations in the potential exchange surface between the liquid iron and the molten silicates.

We monitored the surface $S$ of each drop from its initial position until reaching a stationary regime or until breakup. We normalized $S$ by the corresponding spherical surface ($S^* = S/(4\pi R^2)$). Fig. 19 shows the evolution of $S^*$ as a function of time for an initial radius $R = 10$ mm. Fig. 19 illustrates the influence of the viscosity contrast $R_\mu$. For $R_\mu = 1000$, the drop remains spherical and its normalised surface is constant ($S^* = 1$). When the magma ocean viscosity decreases, the deformation of the drop becomes significant and the oscillations of the drop lead to oscillations of its surface. For $R_\mu = 100$, the drop deforms and oscillates several times, but the surface tension is large enough to prevent the drop from breaking up, and a steady state is reached after 1.5 s. For $R_\mu = 10$ and $R_\mu = 50$, the drop surface oscillations are followed by a fragmentation (at 0.393 s and 0.6 s respectively). These large deformation processes significantly increase the droplet surface, which can increase by up to 200% before the fragmentation. Our results show that the viscosity ratio between the liquid iron and silicate phases strongly influences the shape of the metallic drops. As a consequence, considering a purely spherical drop when calculating the transfers across the exchange surface provides only a lower bound estimate which might be significantly off, especially for $R_\mu < 100$.

## 4.3. Application to a 10 km metal diapir

In the previous sections, we studied the dynamical properties of
<figure><img src="image_99.png" /><figcaption>Fig. 19. Non-dimensional surface exchange as a function of time for viscosity ratios $R_μ$ ranging between 10 and 1000. Here, we consider a drop of initial radius $R = 10$ mm. For viscosity ratio = 10 (black line), $Re = 553$ and $We = 10.9$. For viscosity ratio = 50 (red line), $Re = 108$ and $We = 10.4$. For viscosity ratio = 100 (blue line), $Re = 44.2$ and $We = 8.74$. For viscosity ratio = 200 (purple line), $Re = 25.1$ and $We = 9.02$. For viscosity ratio = 1000 (green line), $Re = 2.70$ and $We = 2.61$. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

sinking droplets as a function of the viscosity ratio between the iron droplet and the magma ocean. We showed that, within an iron cloud, each droplet is associated to a boundary layer thickness and an exchange surface that both depend on the droplet size and viscosity ratio. In our numerical models we do not solve the equations governing the thermo-chemical transfers between the metallic droplets and the magma ocean. However, we can estimate the efficiency of the potential thermo-chemical equilibration by determining the total potential exchange surface within a metallic cloud composed of small droplets and resulting from the fragmentation of a 10 km radius metallic diapir, which is the typical lengthscale of iron fragments dispersed after a large impact (Kendall and Melosh, 2016). We consider in these estimations that the 10 km metallic diapir is fragmented in a population of small droplets that compose a larger scale metallic cloud sinking through the magma ocean (Deguen et al., 2014).

To characterize the size of the droplets within the cloud, we used the analytical law for the distribution of droplets sizes obtained from the fragmentation of a large analog diapir with a viscosity ratio $R_μ = 50$ derived experimentally by Wacheul et al. (2014). We assume here that this analytical law is valid for any system with the same viscosity ratio once it is normalised by the characteristic radius of the system, i.e. the maximal stable radius $R_{max}$. The number of drops $n(R)$ for a given radius $R$ resulting from the breakup of a 10 km iron core is then associated to a gamma distribution:

$$n(R) = NR^{k-1}e^{-\frac{R}{\theta}} \qquad (21)$$

where $N$ is a multiplicative constant determined by volume conservation of iron, $k$ is the shape of the gamma distribution, and $\theta$ is the scale of the gamma distribution normalised by $R_{max}$. Wacheul et al. (2014) explicitly mentioned that the shape and scale might depend on the viscosity ratio. In the absence of any further data, we assume constant values with $\theta = 1.9$ and $k = 2.2$ for viscosity ratios ranging between 10 and 1000. Note that $R_{max}$ is nearly constant in the study of Wacheul et al. (2014) and Eq. (21) is obtained from a snapshot at a given time relatively shortly after the primary breakup, where transient large radius drops are still present. In our study, we have shown that $R_{max}$ is a function of the magma ocean viscosity (Fig. 14). Hence, for each magma ocean viscosity, we use $R_{max}$ from our numerical study (Fig. 14) and determine $N$ by volume conservation from a initially 10 km iron

<figure><img src="image_100.png" /><figcaption>Fig. 20. Normalised exchange surface $S^*$ as a function of magma ocean viscosity after the breakup of a 10 km metal diapir. In the sinking droplet population, no cut-off is considered at $R = R_{max}$.</figcaption></figure>

diapir radius. We then compute the number of drops of each size, and their relative surface of exchange normalised by their corresponding spherical surface.

We have shown in our numerical models that shape oscillations might influence the surface of the drops (Fig. 19), especially for the largest radii and smallest viscosity ratio. We consider here three models of droplet population within the cloud:

1. a non oscillating population of droplets with a constant radius for a fixed viscosity corresponding to the classical iron rain scenario (Rubie et al., 2003),

2. a non oscillating droplet population derived from our models with variable radii where the droplets remain spherical,

3. an oscillating population derived from our models with variable radii where the droplets deform.

Fig. 20 shows the total dimensionless surface within the fragmenting cloud for the three different cases. In all cases, the exchange surface decreases when increasing the magma ocean viscosity because the maximum stable drop radius also increases when increasing the magma ocean viscosity (see Fig. 14). In the iron rain scenario, the population is composed of droplets smaller than in the two cases derived from our models (oscillating and non oscillating) where large drops are allowed in the distribution. Hence the potential exchange surface is larger (up to 60%) in the iron rain scenario and the equilibration should be more efficient. For the population composed of oscillating droplets, the deformation leads to an increase of the surface of the sinking droplets especially for low magma ocean viscosities where $R_{max}$ is smaller than for large magma ocean viscosities. The total exchange surface in the oscillating case is larger than in the non-oscillating case (by 13–43%) but remains smaller than in the iron rain scenario (by 7–37%). For large magma ocean viscosity (i.e. large viscosity contrasts), the droplets do not oscillate and the exchange surface value is close for both oscillating and non oscillating populations. For low magma ocean viscosities (i.e. small viscosity contrasts), the deformations are so important in the oscillating population that the total surface is close to the surface obtained in the iron rain scenario.

It is to be noted that large metallic drops (with $R > R_{max}$) are allowed in the populations used to obtain Fig. 20. Such large drops are not stable but are likely to exist during the second phase of the fragmentation following an impact that occurs in the shallow magma ocean (Wacheul et al., 2014). Within a deep magma ocean, transient drops with radii $R > R_{max}$ will fragment in smaller droplets and disappear from the populations used in Fig. 20. We can thus redo our calculation
<figure><img src="image_101.png" /><figcaption>**Fig. 21.** Normalised exchange surface $S^*$ as a function of magma ocean viscosity after the breakup of a 10 km metal diapir. In the sinking droplet population, a cut-off is considered at $R = R_{max}$.</figcaption></figure>

by imposing a cut-off at $R_{max}$ using the values from our numerical models (Fig. 14). Fig. 21 illustrates the total potential exchange surface within a sinking cloud of metallic droplets as a function of the magma ocean viscosity considering the same initial volume as in Fig. 20. Imposing a cut-off at $R_{max}$ generates a population of smaller droplets than in Fig. 20. In the oscillating and non oscillating populations used in Fig. 21 many droplets are even smaller than the maximum stable drop radius of the iron rain model. Hence Fig. 21 shows that the total surface of the metallic droplets for both the oscillating and non oscillating populations is larger than the total surface derived from the iron rain scenario. Still, the non oscillating model will be less efficient for thermo-chemical equilibrium than the oscillating model and large magma ocean viscosities (i.e. large viscosity contrasts) significantly reduce (by a factor 10) the potential exchange surface between the silicate and metallic phases.

Results from Figs. 20 and 21 both underline the competition between (1) the deformation processes at the scale of the droplet that enhance the equilibration between the iron and silicate phases and (2) the viscosity contrast that limits the potential exchange surface by allowing large drops. In the shallow part of the magma ocean where large transient metallic drops are likely to exist, the thermo-chemical equilibration will be less efficient than in the deepest part of the magma ocean. However, this conclusion needs to be constrained by implementing thermo-chemical transfers in our dynamical models and using realistic partition coefficients and conductivities that should vary with pressure and temperature.

## 5. Conclusions and future works

During the late stages of planetary accretion, large impacts between differentiated protoplanets have strongly influenced the thermo-chemical state of the future terrestrial planets. Following the impact and the formation of a deep magma ocean, the metallic phase from the impactor has overcome strong deformation and fragmentation processes before reaching the deepest part of the magma ocean. The dynamics of this fragmentation probably played a key role on the equilibration efficiency between the metallic phase from the impactor and the impacted proto-mantle.

We have performed axisymmetric numerical simulations to model the sinking dynamics of an initially spherical liquid iron drop within a molten silicate phase. We have explored a large range of relevant parameters, considering initial drop radii in the range of [1–350] mm resulting from the breakup of an initially larger metallic diapir, and magma ocean viscosities in the range [0.05–100] Pa.s. This large range

of viscosity is meant to encompass all relevant geophysical situations, depending on pressure, depth, temperature, composition and impact history.

For large Reynolds numbers, we have showed that the drag coefficient for all viscosity ratios converges towards a unique value on the order of $3.5 \pm 0.5$. We have found that the maximum stable drop radius increases as a function of the magma ocean viscosity following a power law that scales with $\mu_s^{0.32}$. The corresponding critical Weber number increases as a function of the viscosity ratio following two power laws depending on the breakup regime, that scale with $R_\mu^{0.187}$ and $R_\mu^{0.483}$ for the range of low and large magma ocean viscosity respectively. We have identified five breakup mechanisms depending on the Weber and Reynolds numbers and on the viscosity ratio between silicates and iron phases. We have also shown that the initial shape of the metallic drop (spherical, prolate or oblate) strongly influences its subsequent dynamics by modifying its final shape before breakup, its potential exchange surface, and the time and distance before breakup. We have emphasized the effect of the viscosity contrast on the potential exchange surface between the iron phase and the molten surrounding silicates. Indeed, an increase of the viscosity ratio decreases the potential exchange surface between the iron drop and the molten silicates while it increases the time and distance before breakup. We have shown that the dynamical boundary layer thickness increases as a function of the drop radius and the magma ocean viscosity following power law that scales with $R_\mu^{0.916}$ and $\mu_s^{0.062}$. Finally, we have implemented our numerical results in a first order model to characterise the ability of an initially 10 km metallic diapir to exchange with its surrounding environment while fragmenting. Our dynamical model predicts that potential thermo-chemical equilibration within a fragmenting cloud depends on the depth at which the cloud is fragmenting: thermo-chemical equilibration should be less efficient in the shallowest part of the magma ocean than in the deepest part.

The next step is now to implement in our models the resolution of the equations governing the thermo-chemical exchanges between the metallic phase and the magma ocean. The chemical and thermal diffusion of a sinking undeformable sphere has been extensively studied in the chemical/heat transfer literature (Levich, 1962; Clift et al., 1978). In their numerical models, Ulvrová et al. (2011) evaluated time scales of chemical equilibration within an undeformable metallic droplet sinking through a deformable medium. If the exchange dynamics in the case of a sinking sphere is strongly constrained, the dynamics of thermo-chemical equilibration between a deformable droplet and its environment and the influence of the viscosity contrast between the two phases still deserve extensive studies.

As shown in our study, the depth of the magma ocean could also influence the fragmentation dynamics. As large transient droplets are more likely in the early fragmentation regime, full thermo-chemical equilibration could be achieved deeper in the magma ocean. The transfer parameters governing the equilibration rate are also pressure/temperature dependent. The viscosity of the molten silicates (Karki and Stixrude, 2010), the partition coefficient of lithophile/siderophile elements (Bouhifd and Jephcoat, 2003) and the thermal conductivity of iron (de Koker et al., 2012) are all affected by pressure changes with depth. An exhaustive study accounting for realistic parameters in the context of a deep magma ocean will surely help to constrain the thermo-chemical signature at the end of the core-mantle separation.

## Acknowledgements

The authors acknowledge funding by the ERC under the European Union's Horizon 2020 research and innovation program through Grant No. 681835- FLUDYCO-ERC-2015-CoG. This research also received funding from the French PNP program (INSU-CNRS), the French Government Laboratory of Excellence initiative No. ANR-10-LABX-0006 and the Région Auvergne. This paper is Laboratory of Excellence ClerVolc contribution No. 338.
## Appendix A. Computational resolution and convergence

In order to confirm that our mesh correctly captures the dynamics of the falling drop within the magma ocean and to determine the quality and limits of our numerical model, we performed systematics tests on two representative cases.

In the first case, we focus on the dynamics of a strongly oscillating drop that does not converge to a single shape and remains in oscillation (simulation #18 in Table 2). This corresponds to a worst case scenario from a numerical point of view. Fig. 22(left) shows the normalised boundary layer thickness on the drop radius as a function of different grid sizes $h_x = R/15, R/25, R/33.4, R/40, R/50, R/66.7$. We observe a reasonable convergence of the numerical results from $h_x = R/40$, with changes limited to 1.4%. To quantify the dynamic difference between $h_x = R/40$ and the finest mesh $h_x = R/66.7$, we calculate the normalised exchange surface as a function of time in Fig. 22(Right). We note that the first two drop oscillations have the same dynamics; then, differences occur. Yet the dynamical times remain close, and the evolution of the exchange surface is almost the same between the two mesh sizes, with a relative maximum error ≈ 3.5%.

In this paper, we are also interested in the fragmentation modes. To confirm that these fragmentation modes don't depend on the mesh size and don't come from any numerical artifact, we calculate an extremely distorted drop which breaks up after a few oscillations in a second test case (simulation #83 in Table 2). Fig. 23 shows the final drop shape just after the breakup for the grid sizes $h_x = R/40$ and $h_x = R/66.7$. This figure confirms that from $h_x = R/40$, the fragmentation mode and the final drop shape do not change significantly with the grid resolution.

We finally compare the cost (CPU) for one second of simulation of the first test case with different mesh sizes in Table 3. The finest mesh is 4 times more expensive than $h_x = R/40$. The purpose of this article being to perform a systematic study to capture the drop dynamics for a large range of radius and ambient viscosity, we have chosen the mesh size of $h_x = R/40$, which is a good compromise between simulation costs and global dynamical results.

<figure><img src="image_102.png" /><figcaption>Fig. 22. Comparison of different mesh sizes of a strongly oscillating drop (simulation #18 in Table 2). Left: Normalised boundary layer thickness as a function of grid sizes. Right: Normalised exchange surface as a function of time.</figcaption></figure>

<figure><img src="image_103.png" /><figcaption>Fig. 23. Comparison of the final shape of the fragmented drop for two different mesh sizes. (simulation #83 in Table 2).</figcaption></figure>

<table><caption>**Table 3**

Comparison of the cost of one simulation for three different grid sizes.</caption><thead><tr><td>Grid size ($h_x$)</td><td>CPU</td></tr></thead><tbody><tr><td>R/40</td><td>51h15</td></tr><tr><td>R/50</td><td>121h5</td></tr><tr><td>R/66.7</td><td>205h</td></tr></tbody></table>

## References

Assael, M.J., Kakosimos, K., Banish, R.M., Brillo, J., Egry, I., Brooks, R., Quested, P.N., Mills, K.C., Nagashima, A., Sato, Y., Wakeham, W.A., 2006. Reference data for the density and viscosity of liquid aluminum and liquid iron. J. Phys. Chem. Ref. Data 35, 285–300.

Bhaga, D., Weber, M.E., 1981. Bubbles in viscous liquids: shapes, wakes and velocities. Fluid Mech. 105, 61–85.

Bonometti, T., Magnaudet, J., 2006. Transition from spherical cap to toroidal bubbles. Phys. Fluids 18, 052102.

Bouhifd, A., Clési, V., Boujibar, A., Bolfan-Casanova, N., Cartier, C., Hammouda, T., Boyet, M., Manthilake, G., Monteux, J., Andrault, D., 2017. Silicate melts during earth's core formation. Chem. Geol. 461, 128–139.

Bouhifd, M., Jephcoat, A., 2003. The effect of pressure on partitioning of ni and co between silicate and iron-rich metal liquids: a diamond-anvil cell study. Earth Planet. Sci. Lett. 209, 245–255.

Cao, X.K., Sun, Z.G., Li, W.F., Liu, H.F., Yu, Z.H., 2007. A new breakup regime of liquid
drops identified in a continuous and uniform air jet flow. Phys. Fluids 19, 057103.

Chambers, J.E., 2004. Planetary accretion in the inner solar system. Earth Planet. Sci. Lett. 223, 241–251.

Clift, R., Gauvin, W.H., 1970. The motion of particles in turbulent gas streams. 1, 14–28.

Clift, R., Grace, J.R., Weber, M.E., 1978. Bubbles, Drops and Particles. Academic Press, New York.

Dahl, T.W., Stevenson, J., 2010. Turbulent mixing of metal and silicate during planet accretion and interpretation of the hf-w chronometer. Earth Planet. Sci. Lett. 295, 177–186.

Dai, Z., Faeth, G.M., 2001. Temporal properties of secondary drop breakup in the multimode breakup regime. Multiph. Flow 27, 217–236.

Darton, R.C., Harrison, D., 1974. The rise of single gas bubbles in liquid fluidized beds. Trans. Inst. Chem. Eng. 52, 301–306.

Davies, R., Taylor, G.I., 1950. The mechanics of large bubbles rising through extended liquids and through liquids in tubes. Proc. R. Soc. London 200.

de Koker, N., Steinle-Neumann, G., Vlček, V., 2012. Electrical resistivity and thermal conductivity of liquid Fe alloys at high p and t, and heat flux in earth's core. Proc. Natl. Acad. Sci. U.S.A. 109, 4070–4073.

de Wijs, G.A., Kresse, G., Vocadlo, L., Dobson, D., Alfè, D., Gillan, M.J., Price, G.D., 1998. The viscosity of liquid iron at the physical conditions of the earth's core. Nature 392, 805–807.

Deguen, R., Landeau, M., Olson, P., 2014. Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans. Earth Planet. Sci. Lett. 391, 274–287.

Deguen, Renaud, Olson, Peter, Cardin, Philippe, 2011. Experiments on turbulent metal-silicate mixing in a magma ocean. Earth Planet. Sci. Lett. 310 (3–4), 303–313.

Fleck, J.R., Rains, C.L., Weeraratne, D.S., Nguyen, C.T., Brand, D.M., Klein, S.M., McGehee, J.M., Rincon, J.M., Martinez, C., Olson, P.L., 2018. Iron diapirs entrain silicates to the core and initiate thermochemical plumes. Nat. Commun.

Hadamard, J.S., 1911. Mouvement permanent lent d'une sphere liquide et visqueuse dans un liquide visqueux. C. R. Acad. Sci. 152, 1735–1738.

Han, J., Tryggvason, G., 1999. Secondary breakup of axisymmetric liquid drop. Phys. Fluids 11.

Harper, J.F., Moore, D.W., 1968. The motion of a spherical liquid drop at high reynolds number. Fluid Mech. 32, 367–391.

Hartmann, W.K., Davis, D.R., 1975. Satellite-sized planetesimals and lunar origin. Icarus 24, 504–514.

Hu, J., Jia, R., Wan, K.T., Xiong, X., 2014. Simulation of droplet impingement on a solid surface by the level set method. COMSOL Conference in Boston.

Ichikawa, H., Labrosse, S., Kurita, K., 2010. Direct numerical simulation of an iron rain in the magma ocean. Geophys. Res. 115, B01404.

Karato, S.I., Rama Murthy, V., 1997. Core formation and chemical equilibrium in the earth physical considerations. Phys. Earth Planet. Int. 100, 61–79.

Karki, B.B., Stixrud, L.P., 2010. Viscosity of mgsio3 liquid at earth's mantle conditions: implications for an early magma ocean. Science 328, 740–742.

Kékesi, T., Amberg, G., Wittberg, L.P., 2014. Drop deformation and breakup. Int. J. Multiph. Flow 66, 1–10.

Kendall, J.D., Melosh, H., 2016. Differentiated planetesimal impacts in a terrestrial magma ocean: fate of the iron core. Earth Planet. Sci. Lett. 448, 24–33.

Kim, D.S., Burger, M., Frohlich, G., Unger, H., 1983. Experimental investigation of hydrodynamic fragmentation of gatliumn drops in water flows. Int. Mtg on Light Water Reactor Severe Accident Evaluation, Cambridge, Mass.

Kleine, T., Münker, C., Mezger, K., Palme, H., 2002. Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry. Nature 418, 952–955.

Krzeczkowski, S.A., 1980. Measurements of liquid droplet disintegration mechanism. Multiph. Flow 6, 227–239.

Landeau, M., Deguen, R., Olson, P., 2014. Experiments on the fragmentation of a buoyant liquid volume in another liquid. Fluid Mech. 749, 478–518.

Levich, V., 1962. Physicochem. Hydrodyn.

Loth, E., 2008. Quasi-steady shape and drag of deformable bubbles and drops. Multiph. flow 34, 523–546.

Luo, X.-Y., Ni, M.-J., Ying, A., Abdou, M., 2006. Application of the level set method for multi-phase flow computation in fusion engineering. Fusion Eng. Des. 81,

1521–1526.

Mei, R., Klausner, J.F., Lawrence, C.J., 1994. A note on the history force on a spherical bubble at finite reynolds number. Phys. Fluid 6, 418–420.

Monteux, J., Arkani-Hamed, J., 2014. Consequences of giant impacts in early Mars: core merging and Martian dynamo evolution. J. Geophys. Res. (Planets) 119, 480–505. https://doi.org/10.1016/j.icarus.2013.05.008.

Monteux, J., Coltice, N., Dubuffet, F., Ricard, Y., 2007. Thermo-mechanical adjustment after impacts during planetary growth. Geophys. Res. Lett. 34, 24201–24205.

Monteux, J., Ricard, Y., Coltice, N., Dubuffet, F., Ulvrova, M., 2009. A model of metal-silicate separation on growing planets. Earth Planet. Sci. Lett. 287, 353–362. https://doi.org/10.1016/j.epsl.2009.08.020.

Moore, D.W., 1965. The velocity of rise of distorted gas bubbles in a liquid of small viscosity. Fluid Mech. 23, 749–765.

Morbidelli, A., Lunine, J.I., O'Brien, D.P., Raymond, S.N., Walsh, K.J., 2012. Building terrestrial planets. Earth Planet. Sci. 40, 451–475.

O'Brien, David P., Morbidelli, Alessandro, Levison, Harold F., 2006. Terrestrial planet formation with strong dynamical friction. Icarus 184, 39–58. https://doi.org/10.1016/j.icarus.2006.04.005.

Ohta, M., Sussman, M., 2012. The buoyancy-driven motion of a single skirted bubble or drop rising through a viscous liquid. Phys. Fluids 24, 112101.

Ohta, M., Akama, Y., Yoshida, Y., Sussman, M., 2009. Three-dimensional simulations of vortex ring formation from falling drops in an immiscible viscous liquid. Chem. Eng. Jpn. 42, 648–655.

Ohta, M., Yamaguchi, S., Yoshida, Y., Sussman, M., 2010. The sensitivity of drop motion due to the density and viscosity ratio. Phys. Fluids 22, 072102.

Ohta, M., Akama, Y., Yoshida, Y., Sussman, M., 2014. Influence of the viscosity ratio on drop dynamics and breakup for a drop rising in an immiscible low-viscosity liquid. Fluid Mech. 752, 383–409.

Patel, P.D., 1978. Hydrodynamic fragmentation of drops (Ph.D. Dissertation). Purdue Univ. West Lafayette, Ind.

Patel, P.D., Theofanous, T.G., 1981. Hydrodynamic fragmentation of drops. Fluid Mech. 103, 207–223.

Pilch, M., Erdman, C.A., 1987. Use of breakup time data and velocity history data to predict the maximum size of stable fragments for acceleration-induced breakup of a liquid drop. Multiph. Flow 13, 741–757.

Poirier, J.P., 1988. Transport properties of liquid metals and viscosity of the earth's core. Geophys. J. 92, 99–105.

Rubie, D.C., Melosh, H.J., Reid, J.E., Liebske, C., Righter, K., January 2003. Mechanisms of metal-silicate equilibration in the terrestrial magma ocean. Earth Planet. Sci. Lett. 205, 239–255.

Samuel, H., 2012. A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans. Earth Planet. Sci. Lett. 313–314, 105–114.

Stevenson, D.J., 1990. Fluid dynamics of core formation. Origin of the Earth 231–249.

Tonks, W.B., Melosh, H.J., 1992. Core formation by giant impacts. Icarus 100, 326–346.

Touboul, M., Kleine, T., Bourdon, B., Palme, H., Wieler, R., 2007. Late formation and prolonged differentiation of the Moon inferred from W isotopes in lunar metals. Nature 450, 1206–1209.

Tripathi, M.K., Sahu, K.C., Govindarajan, R., 2015. Dynamics of an initially spherical bubble rising in quiescent liquid. Nat. Commun.

Ulvrová, M., Coltice, N., Ricard, Y., Labrosse, S., Dubuffet Velimski, F., Šràmek, M., 2011. Compositional and thermal equilibration of particles, drops, and diapirs in geophysical flows. Geochem. Geophys. Geosyst. 12.

Villermaux, E., Bossa, B., 2009. Single-drop fragmentation determines size distribution of rain drops. Nat. Phys. 5, 697–702.

Wacheul, J.B., Le Bars, M., 2018. Experiments on fragmentation and thermo-chemical exchanges during planetary core formation. Phys. Earth Planet. Int. 276, 134–144.

Wacheul, J.B., LeBars, M., Monteux, J., Aurnou, J.M., 2014. Laboratory experiments on the breakup of liquid metal diapirs. Earth Planet. Sci. Lett. 403, 236–245.

Walter, M.J., Tronnes, R.G., 2004. Early earth differentiation. Earth Planet. Sci. Lett. 225, 253–269.

Wood, B.J., Walter, M.J., Wade, J., 2006. Accretion of the earth and segregation of its core. Nature 441, 825–833.
# Dynamics of core-mantle separation: Influence of viscosity contrast and metal/silicate partition coefficients on the chemical equilibrium

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

Check for
updates

V. Clesi<sup>a,∗</sup>, J. Monteux<sup>a</sup>, B. Qaddah<sup>a,b</sup>, M. Le Bars<sup>b</sup>, J.-B. Wacheul<sup>b</sup>, M.A. Bouhifd<sup>a</sup>

## ABSTRACT

The composition of the Earth's core and mantle is set by the chemical equilibrium between metals and silicates during core/mantle segregation. The metallic core separated from the mantle by gravitational descent in the form of diapirs in a magma ocean, and therefore the dynamics of the diapir's downward movement has an influence on the chemical equilibrium. In this study, we characterize the descent of metallic droplets into a molten silicate using numerical models. By varying the silicate and metal viscosities (between 0.1 and 1000 Pa·s for each phase) as well as the partition coefficient between metal and silicate (D<sup>met/sil</sup>, varying between 1 and 1000), we obtained quantifying parametrizing equations for the degree of equilibrium between molten metal and molten silicate, in a regime characterized by low We (We < 10) and low Re (10<sup>-3</sup> < Re < 10<sup>2</sup>). We showed that the main parameters controlling the equilibrium for a siderophile element are the viscosity of the silicate and the partition coefficient. We applied our parameterization for Ni and Co in the context of late accretion on Earth so as to quantify the variation of the Ni/Co ratio after a large impact as a function of the magma ocean viscosity, for an iron-rain scenario of metal/silicate segregation. Using previous models (Canup, 2004) of the Moon-forming impact, we showed that the Moon formation had an effect on the current Ni/Co ratio. Depending on the radius of Theia's core and the viscosity of the magma ocean produced after the impact between the proto-Earth and Theia, the Moon formation could account for 0.45% to 3% of the current Ni/Co ratio for magma ocean viscosities of 0.1 to 100 Pa·s, respectively.

## 1. Introduction

The main process occurring on Earth during the first 100–150 million years was the segregation of its metallic Fe-rich core (e.g. Bouhifd et al., 2017; Rubie et al., 2015; Wood et al., 2006; and references therein). One major consequence of this segregation is the depletion of the Earth's mantle in the siderophile elements relative to primitive solar system abundances (e.g. Jones and Drake, 1986).

The most recent Earth core formation models include a combination of multistage core-mantle differentiation with N-body accretion simulations, combined with continuous Earth core formation where pressure, temperature, oxygen fugacity and chemical composition of both the mantle and core vary during accretion of the Earth (e.g. Rubie et al., 2015, and references therein). These models consider that the accretion of the Earth occurred during a series of large impact events (e.g. Canup, 2008; Wetherill, 1985; O'Brien et al., 2006; Monteux et al., 2009; Nakajima and Stevenson, 2015). More importantly, these models provide results for the metal-silicate partitioning of a significant number of siderophile elements, which are consistent with observed mantle abundances given that the oxygen fugacity increases from about IW-5 to IW-2 (5 to 2 log units lower than the Iron-Wüstite buffer) during core

segregation (e.g. Cartier et al., 2014; Rubie et al., 2015; Wade and Wood, 2005).

In this context, the behavior of Ni and Co (two refractory and moderately siderophile elements) has been considered to provide an important indication of the conditions of Earth core formation. The consensus is now that metal-silicate equilibration at high pressure, in the range of 40–60 GPa (corresponding to depths of 1000–1500 km), was required to produce the observed Ni and Co depletions in the mantle (Bouhifd and Jephcoat, 2011; Siebert et al., 2012; Righter, 2011; Fischer et al., 2015; Clesi et al., 2016, for some of the most recent studies). Similar conclusions were reached based on the metal-silicate partitioning of lithophile and weakly-siderophile elements (e.g. Mann et al., 2009). One can note here that the conditions of Earth core formation as derived from metal-silicate partitioning of several elements cannot be used as an argument for single-stage core formation. This is highly unlikely given that Earth core formation occurred over a series of large impact events.

Most core formation models are based on metal-silicate partitioning experiments where both phases are fully molten and both thermodynamic and kinetic equilibrium are reached. This type of scenario favors rapid core - mantle segregation, which is confirmed by isotopic

<footnote><sup>a</sup> Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magnas et Volcans, F-63000 Clermont-Ferrand, France</footnote>

<footnote><sup>b</sup> CNRS, Aix Marseille Université, Centrale Marseille, Institut de Recherche sur les Phénomènes Hors Equilibre, UMR 7342, 49, Rue F. Joliot-Curie - B.P. 146, 13384 Marseille Cedex 13, France</footnote>

<footnote>* Corresponding author at: Laboratoire de Géologie de Lyon- Terre, Planètes, Environnement, CNRS UMR 5276, ENS de Lyon, Université Lyon 1, Campus de la Douà, Batiment Geode, 2 Rue Raphaël Dubois, 69622, Villeurbanne Cedex, France.
E-mail address: vincent.clesi@univ-lyon1.fr (V. Clesi).</footnote>
studies on Hf/W and short lived isotopes (Kleine et al., 2002; Jacobsen et al., 2008; Kleine et al., 2009; Kleine and Rudge, 2011). If the kinetic equilibrium is well constrained, the extent to which the mantle and core equilibrate with each other is model dependent. So we know that if the metal were to segregate by forming large diapirs descending through the silicate magma ocean, the timescale for the equilibration would be too long to explain the isotopic ratios of Hf/W (Rubie et al., 2003). Therefore, the metal could have segregated by forming smaller diapirs, or droplets, whose radius varied between 0.1 and 10 cm, forming a cloud of iron in the magma ocean (Deguen et al., 2014; Wacheul et al., 2014). An intermediate scenario combining iron-rain mingling for the first part of segregation in the magma ocean, followed by creation of large descending diapirs also exists, which can explain the equilibrium in a deep magma ocean. In such a scenario, the metal falls through the molten silicate in the form of droplets, and then forms a metallic pond at the base of the magma ocean. Gravitational instability then allows metal to form large diapirs which descend through the solid layer of the mantle into the proto-core (e.g. Stevenson, 1981; Monteux et al., 2015; Fleck et al., 2018). This study focuses on the iron-rain phase of core-mantle segregation, during which chemical equilibrium is reached (Rubie et al., 2003, 2011), as large diapirs have less chemical interaction with the surrounding fluid (Ulvrová et al., 2011; Wacheul et al., 2014). In the scenario of segregation by iron rain (Fig. 1), the extent to which the magma ocean is in equilibrium with the core depends on the properties of the metallic flow. The physics of the fluid flow is dependent on various parameters, one of the most important being the viscosity of the surrounding fluid (Ke and Solomatov, 2009; Monteux et al., 2009; Ricard et al., 2009). Numerical and analog simulations of a metallic diapir descending through a surrounding fluid with different viscosities have shown that the diapir is more stable when surrounding viscosity increases from low values (~$10^{-1}$ Pa·s) to high values (~100 Pa·s), (Samuel, 2012; Wacheul et al., 2014). This has an impact on chemical equilibrium, since break-up of the diapir increases the surface of exchange between metal and silicate and small droplets sink more slowly than large diapirs (Ulvrová et al., 2011; Wacheul et al., 2014). The efficiency of equilibrium is also increased when internal convection favors homogenization of concentrations in the metallic droplets (Ulvrová et al., 2011).

The goal of the present study is to characterize the influence of viscosity (metal and silicate phases) and chemical equilibrium on the global equilibrium rate in a magma ocean. To do that, the most relevant methods are analog and numerical modeling (e.g. Samuel and Tackley, 2008; Šrámek et al., 2010; Deguen et al., 2011; Wacheul et al., 2014). Numerical simulations are more practical than analog ones for studying the effects of different phenomena happening simultaneously. In addition, it is difficult to find analogs for both metal and silicate which will

react exactly like these two phases, and even more complicated to find analogs that behave chemically and physically in the same way as molten metal and molten silicate. Given the limitations of analog experiments for multiple physical modeling, we chose to perform numerical simulations of a falling metallic sphere into a molten silicate using COMSOL Multiphysics 5.0 modeling software. This approach allows the chemical equilibrium for different types of siderophile elements to be determined. In our simulations we combined chemical diffusion calculations with different types of flow. In particular, these simulations led us to explore the contamination of a silicate column by an excess of a siderophile element initially present in the diapir. Since viscosity influences the type of flow, which in turn is crucial in core-mantle equilibrium modeling (Rubie et al., 2003, 2011; Samuel, 2012), viscosity must have an influence on the extent of core-mantle equilibrium. Furthermore, the equilibrium depends on the partition coefficients of the elements between metal and silicate ($D^{\text{met/sil}}$), which depend on pressure, temperature, oxygen fugacity and composition of metal and silicate. Therefore, we tested the contamination of silicate by generic siderophile elements using four different values of $D^{\text{met/sil}}$ (1, 10, 100, 1000).

Since major chemical transfers likely occurred during the iron rain process (Rubie et al., 2003) our simulations looked at the equilibrium of droplets with the surrounding silicate. In general, large diapirs in a turbulent environment tend to break into small droplets, forming a cloud in which each droplet moves at the same speed (Wacheul and Le Bars, 2017, 2018). In this study, we look at how one droplet of a given radius ($R_{\text{Fe}}$ = 1 cm), which can be deformed but cannot break (see Wacheul et al., 2014), interacts with the surrounding silicate. Then this behavior is extrapolated to the entire droplet cloud, and we can estimate the extent to which equilibrium is reached between the magma ocean and the falling metal during planetary accretion.

As a result of these simulations, we constrain the chemical exchanges between a small metallic diapir and a silicate volume as a function of silicate viscosity (varying from 0.1 Pa·s to 1000 Pa·s), the viscosity ratio between metal and silicate, and time. We thus parameterize the effect of varying viscosity on metal-silicate equilibration, and we use this parameterization to re-evaluate models of core-mantle segregation. Moreover, we characterize the effect of changing the magma ocean viscosity on the chemical composition of the mantle. Finally, we discuss the potential consequences of giant impacts such as the Moon-forming impact on the chemical composition of the Earth's mantle.

## 2. Physical model

To model the sinking of a metallic droplet into a silicate liquid, we

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

**Fig. 1.** Schematics of core-mantle segregation. First, the impactor (metal + silicate) enters a turbulent magma ocean. The impactor is fragmented and melted on impacting the surface of the proto-planet. The metallic phase then forms a diapir which, if the magma ocean is turbulent enough, fragments into a cloud of droplets ranging in size from a few millimeters to 10 cm maximum (Rubie et al., 2011; Deguen et al., 2014; Wacheul et al., 2014; Wacheul and Le Bars, 2017, 2018). Each droplet interacts with the surrounding silicate, exchanging heat and chemical elements, until it reaches thermo-chemical equilibrium with its environment.
<caption>**Table 1**</caption>

Parameters and variables used in this study. When values are not given, the variable is calculated later in the text (according to the specific case).

<table><thead><tr><th>Parameter</th><th>Name</th><th>Values in this study</th></tr></thead><tbody><tr><td colspan="3">All models</td></tr><tr><td>ρ</td><td>Density of the fluid</td><td>Silicate: 3300 - metal: 8000 kg·m<sup>-3</sup></td></tr><tr><td>t</td><td>Time</td><td>0.6-125 s</td></tr><tr><td>u</td><td>Velocity vector</td><td>m·s<sup>-1</sup></td></tr><tr><td>R<sub>Fe</sub></td><td>Radius of the diapir</td><td>1 cm</td></tr><tr><td>h<sub>max</sub></td><td>Maximum unit cell size</td><td>0.1 × R<sub>Fe</sub></td></tr><tr><td colspan="3">Two-phase flow</td></tr><tr><td>P</td><td>Dynamic Pressure of the fluid</td><td>1 atm</td></tr><tr><td>μ</td><td>Dynamic viscosity of the fluid</td><td>10<sup>-1</sup>-10<sup>3</sup> Pa·s</td></tr><tr><td>g</td><td>Acceleration of gravity</td><td>9.81 m·s<sup>-2</sup></td></tr><tr><td>l</td><td>Distance fluid - initial interface</td><td>m</td></tr><tr><td>F<sub>st</sub></td><td>Surface tension force</td><td>Calculated for σ = 1 N·m<sup>-1</sup>, in N·m<sup>-3</sup></td></tr><tr><td>Φ</td><td>Volume fraction of the fluid</td><td>Silicate: Φ = 0 metal: Φ = 1 no unit</td></tr><tr><td>ε<sub>ls</sub></td><td>Parameter controlling interface reinitialization</td><td>h<sub>max</sub>/2</td></tr><tr><td>γ</td><td>Parameter controlling velocity reinitialization</td><td>0.001-0.3 m·s<sup>-1</sup></td></tr><tr><td colspan="3">Transport of diluted species</td></tr><tr><td>c<sub>i</sub></td><td>Concentration of element i</td><td>mol·m<sup>-3</sup></td></tr><tr><td>k<sub>c</sub></td><td>Diffusion coefficient of element i</td><td>10<sup>-6</sup> m<sup>2</sup>·s<sup>-1</sup></td></tr><tr><td>R<sub>i</sub></td><td>Reaction rate for element i</td><td>mol·m<sup>-3</sup>·s<sup>-1</sup></td></tr><tr><td>N<sub>i</sub></td><td>Molar flux of element i</td><td>0 mol·m<sup>-2</sup>·s<sup>-1</sup></td></tr><tr><td colspan="3">Scaling of the study</td></tr><tr><td>Re</td><td>Reynolds number</td><td>2.5 × 10<sup>-5</sup> to 110</td></tr><tr><td>We</td><td>Weber number</td><td>2.3 × 10<sup>-5</sup> to 3.5</td></tr><tr><td>R<sub>μ</sub></td><td>Viscosity ratio (=μ<sub>met</sub>/μ<sub>sil</sub>)</td><td>10<sup>-4</sup> to 10<sup>4</sup></td></tr><tr><td>R<sub>C</sub></td><td>Ratio of the silicate volume reacted</td><td>See Section 3</td></tr></tbody></table>

used the software COMSOL Multiphysics, and in particular the modules “Computational Fluids Dynamics” and “Chemical Reaction Engineering”. In each module, a set of equations is defined and solved at each time step. These equations are described in the following subsections. The parameters we used for these equations are listed in Table 1.

## 2.1. Two-phase flow model

In our study, the fluid dynamics are governed by the Navier-Stokes equations that characterize:

* Conservation of mass

$$ \nabla \cdot \mathbf{u} = 0 \qquad (1) $$

* Conservation of momentum

$$ \rho \frac{\partial \mathbf{u}}{\partial t} + \rho (\mathbf{u} \cdot \nabla) \mathbf{u} = - \nabla P I + \nabla \cdot [\mu (\nabla \mathbf{u} + (\nabla \mathbf{u})^T)] + \rho g + F_{st} \qquad (2) $$

Eq. (1) is the conservation of mass for an incompressible fluid with **u** the flow velocity vector. Eq. (2) describes the conservation of momentum, with acceleration of the fluid ($\frac{\partial \mathbf{u}}{\partial t}$) and an inertia term ($\rho(\mathbf{u} \cdot \nabla) \mathbf{u}$) where $\rho$ is fluid density. The first term on the right-hand side of Eq. (2) ($-\nabla P$) is the effect of the dynamic pressure $P$ on the fluid. The second term ($\nabla \cdot [\mu(\nabla \mathbf{u} + (\nabla \mathbf{u})^T)]$) is the component representing the effect of viscous forces on the fluid with $\mu$ the fluid viscosity. The third term ($\rho g$) is the gravitational force that applies to the whole domain with **g** the gravitational acceleration vector. The fourth term $F_{st} = \nabla \cdot (\sigma(I - nn^T)\delta)$ is the surface tension force with $\sigma$ the surface tension coefficient, $I$ the identity matrix, $n$ a unit vector normal to the surface of interest and $\delta$ a smeared out Dirac function located at the interface. Since the volume of metal is small, we ignored the Coriolis

forces in Eq. (2).

We monitored the interface between the liquid iron droplet and the molten silicates using the Level Set method, a Eulerian and implicit method used in multiphase flow problems (e.g. Qaddah et al., 2019). For that, we defined a function $\Phi$, characterizing the silicate when $\Phi = 0$ and the metal when $\Phi = 1$. The boundary between the two phases was then set at $\Phi = 0.5$. These values were obtained by modifying the level set field using a smeared out Heaviside function (Olsson and Kreiss, 2005), which allows direct calculation of the volume fraction of a given phase by integration of the variable $\Phi$ on the volume considered. The equation governing the transport of $\Phi$ is:

$$ \frac{\partial \Phi}{\partial t} + \mathbf{u} \cdot \nabla \Phi = \gamma \nabla \left( \epsilon_{ls} \nabla \Phi - \Phi(1 - \Phi) \frac{\nabla \Phi}{|\nabla \Phi|} \right) \qquad (3) $$

with $\gamma$ (m/s) and $\epsilon_{ls}$ (m) the reinitialization parameters. The $\gamma$ parameter limits the numerical diffusion of the metallic phase during the simulation and needs to be adjusted empirically when the viscosity varies. Indeed, $\gamma$ is a parameter that determines the amount of reinitialization or stabilization and must be carefully adjusted for each specific problem. If $\gamma$ is too low, the thickness of the interface might not remain constant, and oscillations in $\Phi$ could appear because of numerical instabilities. On the other hand, if $\gamma$ is too high, the interface moves incorrectly. The ideal value of $\gamma$ was given by the maximum velocity reached by the diapir during its descent, as detailed in Section 3. The range of values for $\gamma$ is given in Table 1. $\epsilon_{ls}$ is the parameter controlling the width of the interface between fluids during reinitialization. This equation (Eq. (3)) is a combination of the equations presented in Olsson and Kreiss (2005), modified for the COMSOL inbuilt solver. This solver allows the advection (on the left-hand side of the equation) and the reinitialization parameter (right-hand side of the equation) to be solved at the same time. The dynamic viscosity and the density are evaluated using the level set function and are defined by $\mu = \mu_{sil} + \Phi(\mu_{met} - \mu_{sil})$ and $\rho = \rho_{sil} + \Phi(\rho_{met} - \rho_{sil})$, respectively (subscript met refers to metal and sil to silicates). The resolution of the level set equation without using the COMSOL Multiphysics module is detailed in Olsson and Kreiss, 2005. Furthermore, this method has been already validated and applied to a similar problem in Qaddah et al., 2019.

The temperature of both the metallic droplets and the surrounding molten silicates is difficult to constrain. It is related to the history of both phases before the merging and involves radioactive heating, viscous dissipation and accretionary heating. In our study we consider that the impact that precedes the iron droplet sinking homogenizes the temperature of both phases. Hence, we do not consider heat conservation in the equations solved in our models, and focus instead on the chemical equilibration processes alone.

## 2.2. Chemical model

The chemical diffusion and reactions of elements are described by transport equations including an advective term to account for the movement of the phase and conservation of mass, a diffusion term (Fick's law) to account for diffusive transfer, and a reaction term to account for the chemical transfer based on the partition coefficient between metal and silicate. In our models we solve the following equation:

$$ \frac{\partial c_i}{\partial t} + \nabla \cdot (-k_c \nabla c_i) + \mathbf{u} \cdot \nabla c_i = R_i \qquad (4) $$

Eq. (4) is the transport equation for diluted species including Fick's law for the diffusion of an element ($\nabla \cdot (-k_c \nabla c_i)$) with the addition of an advective term ($\mathbf{u} \cdot \nabla c_i$) to account for motions within the fluids. $c_i$ is the concentration of element i in the phase of interest, $k_c$ the diffusion coefficient and $R_i$ the reaction rate for element i.

In our case, there are no external sources of elements ($N_i = 0$, Fig. 2).
<figure><img src="image_106.png" /><figcaption>Fig. 2. Schematic representation of the initial conditions and geometry of our models. The characteristics of the metal and silicate phases are listed respectively in red and black. The boundary conditions are detailed in blue (top), solid red line (side), green (bottom and global conditions). The red arrow illustrates the symmetry axis (red dotted line) used to define the 3D output. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

We assumed the diffusion coefficient ($k_c$) to be the same for both phases, and independent of the partition coefficient ($D_i$). Therefore we defined two functions to describe the behavior of element $i$ in each phase. One function ($X_i^{sil}$) describes the behavior of element $i$ in the silicate, and therefore tends toward 0 when $\Phi = 1$, the other ($X_i^{met}$) described the behavior of element $i$ in the metal and tends toward 0 when $\Phi = 0$. The combination of the two functions give us the global concentration of element $i$ in our computational domain. The link between the two functions, the partitioning coefficient and the phase field is made by defining a reaction at the interface between the metal and silicate (i.e. when $\Phi = 0.5$): $R_i$ is the reaction rate necessary to reach the chemical equilibrium defined by the metal-silicate partition coefficient of element $i$ ($D_i$). $D_i$ is the Nernst partition coefficient given by $D_i = \frac{X_i^{met}}{X_i^{sil}}$, where $X_i^{met}$ is the mass fraction of element $i$ in the metallic phase, and $X_i^{sil}$ the mass fraction of the same element in the silicate phase at equilibrium. $R_i$ is given by Eq. (5):

$$ \frac{dc_i}{dt} = R_i \qquad (5) $$

where $c_i$ is the concentration in either silicate or metal in mol·m⁻³, obtained from the density of the phase where the element is $\rho^{phase}$ (for either a metal or silicate phase), the molar mass of the element ($M_i$) and the mass fraction of the element $i$ in the phase considered is represented by $c_i = \frac{X_i^{phase} \rho^{phase}}{M_i}$. $R_i$ is automatically calculated in a stationary state from the initial conditions and a given partition coefficient so as to reach equilibrium using COMSOL's in-built solver before the simulation. The reaction rate is high enough to maintain the equilibrium between two calculation steps: the typical time for equilibration is $10^{-5}$ s, which is less than the time step solved by solver. $R_i$ is negative when removing siderophile elements from the metallic phase, and positive when creating the same element in the silicate phase. The reaction function is only defined for the metal/silicate boundary, given by the condition $\Phi = 0.5$. The higher the value of $D_i^{met/sil}$, the closer the conditions are to equilibrium, and the lower the reaction rate ($R_i$) will be. For the domain where $\Phi = 0$ or 1, the equilibrium is set to 0, and therefore the reaction rate is 0. For a given time step, the in-built COMSOL Multiphysics 5.0 solver calculate the Navier-Stokes equations (Eqs. (1) and (2)) using finite element method, yielding the velocity field. Then the level-set

equation is computed at the same time the advection-diffusion (Eq. (4)) equation is solved by finite-element methods on the entire computational domain, with the reaction $R_i$ only defined where $\Phi = 0.5$. The results of the calculation are then used as starting condition for the next time step.

## 2.3. Geometry, mesh and initial conditions

First, we defined a geometry for the simulation, using one of the defined boundaries for the eulerian solver available with COMSOL Multiphysics 5.0 software geometries. The best choice for our type of simulation would have been a 3-D model with a falling sphere inside it, so as to compare our simulation to analog simulations (Deguen et al., 2011, 2014; Wacheul et al., 2014). Since this type of simulation is very costly in terms of computational time, especially for fluid flow simulations, we used a 2D-axisymmetric geometry, which solves the equation on a 2D surface, and uses the symmetry conditions to give a 3D output. This is a good compromise between a 2D and 3D simulation when the diapir is not fragmented. Our simulations focused on a weakly deformable falling sphere (i.e. no break up) with a radius of 1 cm in a cylinder whose diameter was set during the calibration of the simulation.

The geometry and boundary conditions used in this study are presented in Fig. 2. In this geometry, a rotational symmetry axis is defined. For the fluid flow, the boundaries of our calculation domain are characterized as follows:

- on the top boundary, the condition is given by $P = 1$ atm,

- on the bottom boundary, the condition is set to be a null velocity field

- on the side boundary, the condition is a no slip condition

- on all the boundary there is a no flow condition for the level set equation.

The conditions for chemical diffusion is a constant budget in chemical elements throughout the numerical simulation, i.e. no elemental flux ($N_i = 0$), inward or outward, is allowed on all the external boundaries.
For our calculation domain we used a fixed mesh of triangles, with three main parameters: $h_{max}$, $h_{min}$ and a curving parameter. The $h_{max}$ parameter gives the maximum size of the triangle side, $h_{min}$ gives the minimum size and the curving parameter gives the arc interception for curving triangles at spherical boundaries. The parameter controlling the accuracy and the computation time for the calculation is $h_{max}$. The best value for this parameter was determined in the resolution study, see Section 2.4.

The initial conditions are described in Fig. 2. The viscosities were varied for each simulation, but the geometry of the set-up, the initial concentrations, the densities and the droplet radius were kept constant. To avoid numerical instabilities inherent in $X_{sil}^0 = 0\%$wt, the initial concentration was set to $X_{sil}^0 = 10^{-4}\%$wt for the silicate with a density $\rho_{sil} = 3300 \text{ kg.m}^{-3}$; and $X_{met}^0 = 20\%$ wt for the metal with a density $\rho_{met} = 8000 \text{ kg.m}^{-3}$. The radius of the droplet was $R_{Fe} = 1$ cm, falling in a cylinder with a 5 cm radius and a 28 cm height, with the center of the droplet placed initially at 26 cm height.

## 2.4. Tests of the COMSOL simulation

To correctly set up the computational domain in terms of size and mesh, two effects need to be accounted for:

(1) Grid resolution, which has to be sufficient to capture the physics of small droplet deformation.

(2) Boundary effects, to avoid a major influence of the outer border of the domain on the flow and the generation of artifacts in the velocity field.

To minimize the numerical diffusion of iron during the diapir's descent, we tested several maximum sizes for the mesh grid component $h_{max}$. Fig. 3 illustrates the time evolution of the ratio between the calculated volume over the theoretical volume of one droplet (given by $V_{theoretical} = \frac{4}{3}\pi R_{Fe}^3$ with $R_{Fe} = 1$ cm) for $h_{max}$ values ranging between

<figure><img src="image_107.png" /><figcaption>Fig. 3. Evolution of the metallic volume ratio ($V_{calculated}/V_{theoretical}$) during a simulation for different $h_{max}$ values, for a silicate viscosity of 100 Pa·s. The colored and dashed line represent the evolution of $V_{calculated}/V_{theoretical}$ for different mesh sizes with time. The scale at the top shows the distance traveled by the diapir over a given time.</figcaption></figure>

0.06 and 0.2 cm. The simulations were performed using the same initial geometry as in Fig. 2, with a silicate and a metal viscosity of 100 Pa·s (corresponding to a Reynolds number of ~2. $10^{-3}$, see Eq. (6)) and no chemical component to solve. The precision of the calculation increased as the mesh size decreased, especially when there was deformation. On the other hand, the smaller the mesh grid components were, the longer it took to complete the calculation. As can be seen in Fig. 3, the precision is better for a mesh component size of 0.1 cm than for smaller (0.06 and 0.08 cm) or larger (0.2 cm) mesh component sizes. In addition, a maximum component size allowed a calculation time of a few hours (~6 h), whereas a mesh grid size of 0.06 cm yielded calculation times from a dozen hours to 2 days for the less viscous flows. Therefore, for the rest of the study, we used a maximum size of 0.1 cm for the mesh grid (i.e. $R_{Fe}/10$). The mesh size could be adapted to the type of flow for each simulation: for example, for a higher degree of diapir deformation during the simulation (typically for less viscous flows), the mesh size had to be decreased to avoid numerical artifacts. To increase the precision, COMSOL also made it possible to set a minimum size for the mesh grid component so as to get an adapted mesh, with smaller elements at the interfaces between the two phases, and larger elements as the distance from the interface increases. Varying the parameter $h_{min}$ did not affect the simulation for a given value of $h_{max}$, so it was kept at a constant value of $h_{min} = h_{max}/40$.

The second point was to avoid boundary effects on the fluid flow during the simulation. We tested several domain sizes and checked the different velocity field components (r- and z-components). We chose a height equal to $28 \times R_{Fe}$ and a width equal to $5 \times R_{Fe}$ which was sufficient to avoid the recirculating flow associated to the descending diapir having a large effect on the velocity field (Fig. 4), while keeping computational time within reasonable bounds. This size of domain still affected the velocity of the diapir (Chang, 1961 and Fig. 4). Fig. 4 shows that increasing the width to $6 \times R_{Fe}$ changed the velocity field, albeit very slightly, which indicates a diminution of the wall effect: the mean velocity of the diapir increased by 8% compared to the mean velocity for a width of $5 \times R_{Fe}$. However, increasing the width of the domain from $5 \times R_{Fe}$ to $6 \times R_{Fe}$ led to a significant increase in computational time (by at least 46%). Given the limited effect of increasing the domain width and the significant effects on the computational time, we considered this smaller width to be suitable for the study. We applied the same reasoning to choose the height of our domain. To test the resolution of the diffusion in Eq. (4), we benchmarked the validity of the COMSOL solver against steady and unsteady analytical solutions for diffusion problems (Crank, 1975). To test the advective part, we compared the flow velocities from our models with those velocities obtained by Samuel (2012) (see Section 3.2). Our results showed an acceptable agreement, despite a large discrepancy between theoretical results and observed velocities (45% for the largest one) for flows with the largest velocities, i.e. in the intermediate regime. To explain the observed differences, we calculated the mesh Peclet number $Pe_{mesh}$ after each simulation with $Pe_{mesh} = \frac{v h_{max}}{k_c}$, where v is the velocity of the diapir, $h_{max}$ is the mesh size and $k_c$ is the diffusion coefficient ($10^{-6} \text{ m}^2\text{s}^{-1}$, which is higher than a typical chemical diffusion, but high enough to limit numerical diffusion, see Qaddah et al., accepted manuscript). In the reference case and for the other cases within the Stokes regime, the value of $Pe_{mesh}$ was sufficiently low ($Pe_{mesh} < 10$) meaning that numerical diffusion was negligible (see Mittal and Jain, 2012). For the cases within the intermediate regime, $Pe_{mesh}$ was higher ($Pe_{mesh} > 50$), which indicates possible numerical diffusion that likely affected our results. To correct this effect a much smaller mesh size would have been necessary. However, the mesh size necessary to prevent numerical diffusion in our models would have led to unreasonable computational times. Therefore, for the lowest viscosity used in our models some error calculations might have occurred in our results, but are negligible in the final model (see Sections 5 and 6).
<figure><img src="image_108.png" /><figcaption>Fig. 4. Evolution of the velocity field, r-component (on the left) and z-component (on the right), for different widths of computational domain. Widths vary from 2 × RFe to 6 × RFe. The chosen width for the simulations is 5xRFe. Note that the results presented here are calibration runs, with a height of 30xRFe. The same simulations have been performed at varying heights, leading to a final height of 28 × RFe (see Fig. 2 for the calibrated geometry).</figcaption></figure>

## 3. Results

In this section we present the results of our simulations as well as a definition of the non-dimensional numbers and parameters used to interpret the results of the simulation in terms of degree of equilibrium. A reference run was arbitrarily defined to set a reference time evolution for the non-dimensional parameters. Then we varied the viscosity of the silicate and metallic phase for a given partition coefficient value, and finally we varied the value of the partition coefficient for a given viscosity.

### 3.1. Reference case and characteristic non-dimensional numbers

In this section we present the results for a reference case, where μsil = μmet = 100 Pa·s. The temporal evolution of the reference case flow is presented in Fig. 5. This figure shows that for this case both the diapir deformation and the numerical diffusion are weak. The metallic droplet moves a distance of 10 × RFe (10 cm) in 12.1 s.

<figure><img src="image_109.png" /><figcaption>Fig. 5. Evolution of the volume fraction of the silicate during the simulation. When equal to 1 (red), the fluid is a silicate, when equal to 0 (blue) the fluid is a metal. The limit between metal and silicate is the black line between the red and blue parts (Φts = 0.5). The time necessary for the diapir to reach the distance of 10RFe is 12.1 s. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

To compare this flow to other flows with different viscosities, we need to use non-dimensional numbers. Since this study focuses on the viscosity variations and the characteristics of the flow, the Reynolds number (Re) is the more adequate non-dimensional number, which is given by Eq. (6):

$$Re = \frac{\rho_{sil} R_{Fe} v_{diapir}}{\mu_{sil}} \qquad (6)$$

where ρsil is the density of the surrounding silicate, RFe is the initial radius of the diapir, vdiapir is the velocity of the diapir and μsil is the silicate viscosity (see Table 1). This number expresses the ratio of inertial forces to viscous forces. When Re < 1, the flow is considered to be a Stokes flow, which means that diapir deformation is limited during the descent. When 1 < Re < 500, the flow is an intermediate flow between Newtonian and Stokes flow (Samuel, 2012). In our study, we did not investigate the case of Re > 500, which is a highly turbulent Newtonian flow.

The velocity of the diapir is strongly governed by the viscosity of the surrounding silicate, which affects the value of Re (e.g. Samuel, 2012, and references therein). In order to compare different values of the viscosity ratio Rμ = μmet/μsil, we defined the reference time for each calculation as the time necessary for a metallic droplet with radius RFe to move through a distance d = 10 × RFe.

The other non-dimensional number which controls the flow is the Weber number (We), which is the ratio of the inertial forces on the surface tension:

$$We = \frac{\rho_{sil} v_{diapir}^2 R_{Fe}}{\sigma} \qquad (7)$$

with σ the surface tension (see Table 1 for value). The break-up of the diapir is controlled by the Weber number (Wacheul et al., 2014). When We < 6 the diapir is stable in Stokes and intermediate flows, and when We > 6 the surface tension is no longer high enough to prevent the diapir breaking up (Wacheul and Le Bars, 2017). In this case and in the rest of the present study, the diapirs are always stable with We < 6. For the reference case, the diapir velocity is about vdiapir ≈ 8 × 10⁻³ m·s⁻¹. This velocity yields low values of Re and We (2.2 × 10⁻³ and 1.2 × 10⁻³, respectively), which explains the diapir's stability during its descent (see Fig. 5).

To study the chemical equilibrium between metal and silicate, we calculated a non-dimensional concentration, or more precisely a non-dimensional mass fraction termed X'. This non dimensional mass fraction was calculated so as to satisfy three conditions:
(i) If there is no change in the initial concentration in the silicate, X' is equal to 0 (no equilibrium at all)

(ii) If equilibrium is reached, X' is equal to 1.

(iii) X' is proportional to the concentration in the silicate ($X' \propto X$)

The first condition is reached if $X = X_{sil}^0$. Therefore $X' = 0$ implies a null numerator for $X'$, i.e. $X' \propto X - X_{sil}^0$. For the second condition, it is reached if $X = X_{met}^{eq}/D^{met/sil}$, by definition of the partition coefficient, where $X_{met}^{eq}$ is the concentration of the element in the metallic phase at equilibrium. In our case, the element is a siderophile ($D^{met/sil} \ge 1$) and the initial concentration is high ($X_{met}^0 = 0.2$). Therefore, at the end of the simulation we could approximate the equilibrium concentration by $X_{met}^{eq} \sim X_{met}^0$. Satisfying the condition $X' = 1$, with $X' \propto X - X_{sil}^0$ and the non-dimensionality of $X'$ yield to the following formula given in Eq. (8):

$$X' = \frac{X - X_{sil}^0}{X_{met}^0 / D^{met/sil} - X_{sil}^0} \quad (8)$$

$X'$ is a useful parameter for tracking the evolution of concentration in a particular case, but it is not enough to quantify how much of the silicate is affected by a change in concentration. Moreover, this parameter is not suitable for comparing different studies, because it is strongly affected by the duration of the sinking. The evolution of the non-dimensional concentration, $X'$, in the reference case is shown in Fig. 6. Equilibrium in the silicate is reached when $X' = 1$ (area in red). This area forms a very narrow zone around the diapir, and most of the silicate is not chemically equilibrated with the metal, because the simulation run time is not long enough.

The parameter $R_C$ is the ratio of the silicate volume that has reacted with the metal to the total volume of silicate. The definition is $R_C = V_{reacted}^{sil}/V_{total}^{sil}$, where $V_{reacted}^{sil}$ is the volume of silicate where $X > X_{sil}^0$, meaning $R_C$ is the volume fraction of silicate that has reacted with the metal. To visualize how the ratio is calculated, a representation of the integrated volume is shown in Fig. 7. In this figure, the areas colored in red represent the volume of silicate that has evolved in terms of chemical composition, while the areas in blue represent that which has not been contaminated by the diapir's descent.

For each simulation, the volume chemically affected by the diapir's descent increases with time. The time evolution of the $R_C$ value for the reference case is shown in Fig. 8. $R_C$ is increasing with time following

<figure><img src="image_110.png" /><figcaption>Fig. 6. Evolution of adimensional concentration of the siderophile element in the silicate, for a partition coefficient $D^{met/sil} = 1$. For clarity purposes the metal is shown by the grey area surrounded by a black line. The black line is the boundary between metal and silicate and defines the surface where the exchange reaction happens. The adimensional concentration is given by the following formula given in Eq. (9). In this case $X_{sil}^0 = 10^{-5}$, $X_{met}^0 = 0.2$ and $D^{met/sil} = 1$, and the chemical Peclet number is $Pe \sim 800$.</figcaption></figure>

<figure><img src="image_111.png" /><figcaption>Fig. 7. Evolution of the volume of silicate undergoing reaction during the simulation. In red, the volume of silicate affected chemically. In blue, the volume of silicate unaffected by the passage of the diapir. The metal is colored in grey. The ratio $R_C$ defined in the study corresponds to the volume of the red areas divided by the total volume of silicate (red + blue areas). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_112.png" /><figcaption>Fig. 8. Evolution of $R_C$ as a function of time in the reference case defined in the text ($\mu_{sil} = \mu_{met} = 100$ Pa·s). The four lines illustrate four different values of $D^{met/sil}$ ranging from 1 (black) to $10^3$ (blue). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

the same trend in each case. The difference is in the absolute values of $D^{met/sil}$; the more elevated the value is, the closer the initial conditions are to equilibrium conditions, which leads to a lower reaction rate.

### 3.2. Range of Re and We values

The aim of our models was to constrain the influence of the viscosity contrast between the silicate and metallic phases. This led to a wide