V Clesi

J Monteux

B Qaddah

M Le Bars

J.-B Wacheul

M A Bouhifd

Dynamics of core-mantle separation: Influence of viscosity contrast and metal/silicate partition coefficients on the chemical equilibrium

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 met/sil , 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 -3 <Re<10 2 ).

 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.

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. [START_REF] Bouhifd | Silicate melts during Earth's core formation[END_REF][START_REF] Rubie | Accretion and differentiation of the terrestrial planets with implications for the compositions of early-formed Solar System bodies and accretion of water[END_REF][START_REF] Wood | Accretion of the Earth and segregation of its core[END_REF]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. [START_REF] Jones | Geochemical constraints on core formation in the Earth[END_REF].

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. [START_REF] Canup | Accretion of the Earth[END_REF][START_REF] Wetherill | Occurrence of giant impacts during the growth of the terrestrial planets[END_REF][START_REF] O'brien | Terrestrial planet formation with strong dynamical friction[END_REF]Monteux et al., 2009;[START_REF] Nakajima | Melting and mixing states of the Earth's mantle after the Moon-forming impact[END_REF]. 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. [START_REF] Cartier | Redox control of the fractionation of niobium and tantalum during planetary accretion and core formation[END_REF]Rubie et al., 2015;[START_REF] Wade | Core formation and the oxidation state of the Earth[END_REF].

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 [START_REF] Bouhifd | Convergence of Ni and Co metal-silicate partition coefficients in the deep magma-ocean and coupled silicon-oxygen solubility in iron melts at high pressures[END_REF][START_REF] Siebert | Metal-silicate partitioning of Ni and Co in a deep magma ocean[END_REF]Righter, 2011;[START_REF] Fischer | High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O[END_REF][START_REF] Clesi | Effect of H2O on metal-silicate partitioning of Ni, Co, V, Cr, Mn and Fe: Implications for the oxidation state of the Earth and Mars[END_REF] for some of the most recent studies). Similar conclusions were reached based on the metalsilicate partitioning of lithophile and weakly-siderophile elements (e.g. [START_REF] Mann | Evidence for high-pressure core-mantle differentiation from the metal-silicate partitioning of lithophile and weakly-siderophile elements[END_REF].

One can note here that the conditions of Earth core formation as derived from metalsilicate 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 studies on Hf/W and short lived isotopes [START_REF] Kleine | Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry[END_REF][START_REF] Jacobsen | 26Al-26Mg and 207Pb-206Pb systematics of Allende CAIs: Canonical solar initial 26Al/27Al ratio reinstated[END_REF]. 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 [START_REF] Rubie | Mechanisms of metalsilicate equilibration in the terrestrial magma ocean[END_REF]. Therefore, the metal could have segregated by forming smaller diapirs, or droplets, whose radius varied between 0.1 to 10 cm, forming a cloud of iron in the magma ocean [START_REF] Deguen | Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans[END_REF][START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. 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. [START_REF] Stevenson | Models of the Earth's core[END_REF][START_REF] Monteux | Giant impacts, heterogeneous mantle heating and a past hemispheric dynamo on Mars[END_REF][START_REF] Fleck | Iron diapirs entrain silicates to the core and initiate thermochemical plumes[END_REF]. This study focuses on the iron-rain phase of core-mantle segregation, during which chemical equilibrium is reached [START_REF] Rubie | Mechanisms of metalsilicate equilibration in the terrestrial magma ocean[END_REF][START_REF] Rubie | Heterogeneous accretion, composition and core-mantle differentiation of the Earth[END_REF], as large diapirs have less chemical interaction with the surrounding fluid (Ulvrova et al., 2011;[START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. 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 [START_REF] Ke | Coupled core-mantle thermal evolution of early Mars[END_REF]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), [START_REF] Samuel | A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans[END_REF][START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. 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 (Ulvrova et al., 2011;[START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. The efficiency of equilibrium is also increased when internal convection favors homogenization of concentrations in the metallic droplets (Ulvrova 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. [START_REF] Samuel | Dynamics of core formation and equilibration by negative diapirism[END_REF]Sramek et al., 2010;[START_REF] Deguen | Experiments on turbulent metal-silicate mixing in a magma ocean[END_REF][START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF].

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 [START_REF] Rubie | Mechanisms of metalsilicate equilibration in the terrestrial magma ocean[END_REF]2011;[START_REF] Samuel | A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans[END_REF], 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 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 met/sil (1, 10, 100, 1000).

Since major chemical transfers likely occurred during the iron rain process [START_REF] Rubie | Mechanisms of metalsilicate equilibration in the terrestrial magma ocean[END_REF] 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 andLe Bars, 2017, 2018). In this study, we look at how one droplet of a given radius (RFe = 1 cm), which can be deformed but cannot break (see [START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF], 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.

Physical Model

To model the sinking of a metallic droplet into a silicate liquid, we 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.

Two-phase flow model

In our study, the fluid dynamics are governed by the Navier-Stokes equations that characterize: 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 Φ, characterizing the silicate when Φ = 0 and the metal when Φ = 1. The boundary between the two phases was then set at Φ = 0.5. These values were obtained by modifying the level set field using a smeared out Heaviside function [START_REF] Olsson | A conservative level set method for two phase flow[END_REF], which allows direct calculation of the volume fraction of a given phase by integration of the variable Φ on the volume considered. The equation governing the transport of Φ is:

 Conservation of mass ∇ • 𝒖 = 0 (1)  Conservation of momentum 𝜌 𝜕𝒖 𝜕𝑡 + 𝜌(𝒖 • ∇)𝒖 = -∇𝑃𝑰 + ∇ • [𝜇(∇𝒖 + (∇𝒖) 𝑇 )] + 𝜌𝒈 + 𝑭 𝒔𝒕 (
∂Φ 𝜕𝑡 + 𝒖 • ∇Φ = 𝛾∇ (𝜖 𝑙𝑠 ∇Φ -Φ(1 -Φ) ∇Φ |∇Φ| ) (3)
with γ (m/s) and 𝜖 𝑙𝑠 (m) the reinitialization parameters. The γ parameter limits the numerical diffusion of the metallic phase during the simulation and needs to be adjusted empirically when the viscosity varies. Indeed, γ is a parameter that determines the amount of reinitialization or stabilization and must be carefully adjusted for each specific problem.

If γ is too low, the thickness of the interface might not remain constant, and oscillations in Φ could appear because of numerical instabilities. On the other hand, if γ is too high, the interface moves incorrectly. The ideal value of γ was given by the maximum velocity reached by the diapir during its descent, as detailed in section 3. The range of values for γ is given in Table 1. 𝜖 𝑙𝑠 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 [START_REF] Olsson | A conservative level set method for two phase flow[END_REF], modified for the COMSOL in-built solver. This 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.

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:

𝜕𝑐 𝑖 𝑑𝑡 + ∇ • (-𝑘 𝑐 ∇𝑐 𝑖 ) + 𝒖 • ∇𝑐 𝑖 = 𝑅 𝑖 (4)
Eq. 4 is the transport equation for diluted species including Fick's law for the diffusion of an element (∇ • (-𝑘 𝑐 ∇𝑐 𝑖 )) with the addition of an advective term (𝒖 • ∇𝑐 𝑖 ) to account for motions within the fluids. ci is the concentration of element i in the phase of interest, kc the diffusion coefficient and Ri the reaction rate for element i.

In our case, there are no external sources of elements (Ni=0, Fig. 2).

We assumed the diffusion coefficient (kc) to be the same for both phases, and independent of the partition coefficient (Di ). Therefore we defined two functions to describe the behavior of element i in each phase. One function (𝑋 𝑖 𝑠𝑖𝑙 ) describes the behavior of element i in the silicate, and therefore tends toward 0 when Φ = 1, the other (𝑋 𝑖 𝑚𝑒𝑡 ) described the behavior of element i in the metal and tends toward 0 when Φ =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 Φ = 0.5):

Ri is the reaction rate necessary to reach the chemical equilibrium defined by the metalsilicate partition coefficient of element i (Di ). Di is the Nernst partition coefficient given by 𝐷 𝑖 = 𝑋 𝑖 𝑚𝑒𝑡 𝑋 𝑖 𝑠𝑖𝑙 , where 𝑋 𝑖 𝑚𝑒𝑡 is the mass fraction of element i in the metallic phase, and 𝑋 𝑖 𝑠𝑖𝑙 the mass fraction of the same element in the silicate phase at equilibrium. Ri is given by Eq. 5:

𝑑𝑐 𝑖 𝑑𝑡 = 𝑅 𝑖 (5)
Where ci is the concentration in either silicate or metal in mol.m -3 , obtained from the density of the phase where the element is 𝜌 𝑝ℎ𝑎𝑠𝑒 (for either a metal or silicate phase), the molar mass of the element (Mi) and the mass fraction of the element i in the phase considered is represented by 𝑐 𝑖 = 𝑋 𝑖 𝑝ℎ𝑎𝑠𝑒 𝜌 𝑝ℎ𝑎𝑠𝑒 𝑀 𝑖 ⁄ . Ri 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. Ri 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 Φ = 0.5. The higher the value of Di met/sil , the closer the conditions are to equilibrium, and the lower the reaction rate (Ri) will be.

For the domain where Φ = 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 advectiondiffusion (Eq. 4) equation is solved by finite-element methods on the entire computational domain, with the reaction Ri only defined where Φ = 0.5. The results of the calculation are then used as starting condition for the next time step.

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 be a 3-D model with a falling sphere inside it, so as to compare our simulation to analog simulations [START_REF] Deguen | Experiments on turbulent metal-silicate mixing in a magma ocean[END_REF][START_REF] Deguen | Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans[END_REF][START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. Since this type of simulation is very costly in terms of computational time, especially for fluid flow simulations, we used a 2D-axisymetric 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 (Ni = 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: hmax, hmin and a curving parameter. The hmax parameter gives the maximum size of the triangle side, hmin 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 hmax. 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 𝑋 𝑠𝑖𝑙 0 = 0 %wt , the initial concentration was set to 𝑋 𝑠𝑖𝑙 0 = 10 -4 %wt for the silicate with a density ρsil = 3300 kg.m -3 ; and 𝑋 𝑚𝑒𝑡 0 = 20 % wt for the metal with a density ρmet = 8000 kg.m -3 . The radius of the droplet was RFe = 1cm, 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.

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 hmax. 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 hours), 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. RFe/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 hmin did not affect the simulation for a given value of hmax, so it was kept at a constant value of hmin = hmax/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 (rand z-components). We chose a height equal to 28xRFe and a width equal to 5xRFe 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 6x RFe 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 5xRFe.

However, increasing the width of the domain from 5xRFe to 6xRFe 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 [START_REF] Crank | The Mathematics of Diffusion[END_REF]. 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 𝑃𝑒 𝑚𝑒𝑠ℎ after each simulation with 𝑃𝑒 𝑚𝑒𝑠ℎ = 𝑣ℎ 𝑚𝑎𝑥

𝑘 𝑐

, where v is the velocity of the diapir, hmax is the mesh size and kc is the diffusion coefficient (10 -6 m 2 .s -1 , which is higher than a typical chemical diffusion, but high enough to limit numerical diffusion (see Qaddah et al., 2020). In the reference case and for the other cases within the Stokes regime, the value of Pemesh was sufficiently low (Pemesh < 10) meaning that numerical diffusion was negligible (see [START_REF] Mittal | Redefined cubic B-splines collocation method for solving convection-diffusion equations[END_REF]. For the cases within the intermediate regime, Pemesh was higher (Pemesh > 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 Section 5 and 6).

Results

In this section we present the results of our simulations as well as a definition of the nondimensional 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.

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

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

𝑅𝑒 = 𝜌 𝑠𝑖𝑙 𝑅 𝐹𝑒 𝑣 𝑑𝑖𝑎𝑝𝑖𝑟 𝜇 𝑠𝑖𝑙 (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 [START_REF] Samuel | A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans[END_REF]. 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:

𝑊𝑒 = 𝜌 𝑠𝑖𝑙 𝑣 𝑑𝑖𝑎𝑝𝑖𝑟 2 𝑅 𝐹𝑒 𝜎 (7)
with σ the surface tension (see Tab. 1 for value). The break-up of the diapir is controlled by the Weber number [START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. 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 [START_REF] Wacheul | Fall and fragmentation of liquid metal in a viscous fluid[END_REF]. 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 -3 m.s -1 . This velocity yields low values of Re and We

(2.2×10 -3 and 1.2×10 -3 , 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 nondimensional 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 (𝑋 ′ ∝ 𝑋)
The first condition is reached if 𝑋 = 𝑋 𝑠𝑖𝑙 0 . Therefore X' = 0 implies a null numerator for X', i.e. 𝑋 ′ ∝ 𝑋 -𝑋 𝑠𝑖𝑙 0 . For the second condition, it is reached if 𝑋 = 𝑋 𝑚𝑒𝑡 𝑒𝑞 𝐷 𝑚𝑒𝑡/𝑠𝑖𝑙 ⁄ , by definition of the partition coefficient, where 𝑋 𝑚𝑒𝑡 𝑒𝑞 is the concentration of the element in the metallic phase at equilibrium. In our case, the element is a siderophile (𝐷 𝑚𝑒𝑡/𝑠𝑖𝑙 ≥ 1)

and the initial concentration is high (𝑋 𝑚𝑒𝑡 0 = 0.2). Therefore, at the end of the simulation we could approximate the equilibrium concentration by 𝑋 𝑚𝑒𝑡 𝑒𝑞 ~ 𝑋 𝑚𝑒𝑡 0 . Satisfying the condition X' = 1, with 𝑋 ′ ∝ 𝑋 -𝑋 𝑠𝑖𝑙 0 and the non-dimensionality of X' yield to the following formula given in Eq. 8:

𝑋 ′ = 𝑋-𝑋 𝑠𝑖𝑙 0 𝑋 𝑚𝑒𝑡 0 /𝐷 𝑚𝑒𝑡 𝑠𝑖𝑙 ⁄ -𝑋 𝑠𝑖𝑙 0 (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 nondimensional 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. 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 RC value for the reference case is shown in Fig. 8. RC is increasing with time following 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.

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 range of Re and We values in our numerical simulations. The viscosity of the silicates surrounding the metal droplet controls the flow regime and the sinking velocity. The evolution of the time necessary to reach the distance of 10×RFe, which depends mostly on the silicate viscosity [START_REF] Samuel | A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans[END_REF] references therein), is presented in Fig. 9. The higher the viscosity of the silicate, the longer it takes for the diapir to reach a set distance. The changes in flow regime result in different expressions for the diapir velocity. In our study, the range of Re values only allows Stokes regime flows (Re < 1) or intermediate regime flows (1<Re<500), with two different expressions for the diapir's velocity. In a Stokes flow the velocity is given by:

𝑣 𝑆𝑡𝑜𝑘𝑒𝑠 = 2 9 (𝜌 𝑚𝑒𝑡 -𝜌 𝑠𝑖𝑙 )𝑔𝑅 𝐹𝑒 2 𝜇 𝑠𝑖𝑙 (9)
For intermediate flows, the sinking velocity is given by:

𝑣 𝑖𝑛𝑡𝑒𝑟𝑚𝑒𝑑𝑖𝑎𝑡𝑒 = √ (𝜌 𝑚𝑒𝑡 -𝜌 𝑠𝑖𝑙 )𝑔𝑅 𝐹𝑒 𝜌 𝑠𝑖𝑙 𝐶 𝐷 2 (10)
with CD the drag coefficient. In our case, the drag coefficient is approximated by 𝐶 𝐷 = 12 𝑅𝑒 + 0.3, following [START_REF] Samuel | A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans[END_REF]. From Eq. 10, we can also calculate the velocity in each drop flow regime: Stokes, intermediate, and Newton (e.g. Qaddah et al., 2019). The theoretical times associated with the velocity are compared to the actual times observed in simulations in Fig. 9. It shows a good agreement between theoretical and observed diapir velocity, despite some scatter. This scatter is consistent with the wall effect of our simulations: the variation in velocity observed in the Stokes flow is of the order of 10 %, which is consistent with the variation expected if the wall effect is negligible (see Section 2.4). Comparison with the typical equilibrium time from [START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF] in Fig. 9 shows that, except for the intermediate flows, the fall time chosen in our study is long enough to reach full equilibrium (or more than 90 % of equilibrium in some cases) inside the metal, leaving the silicate as the limiting phase for equilibrium. The dispersion of the obtained reference time (Eqs. 8 and 9) in Fig. 9 for a chosen value of µsil is considered therefore as a negligible discrepancy in the rest of the study.

In Fig. 10 we plot the Re values as a function of the We values for each simulation (i.e.

for each viscosity contrast used in our calculations). The velocity of the diapir (vdiapir) used to calculate Re and We is determined from our numerical models by the time necessary for the diapir to move through a distance of 10×RFe. This time is dependent on the viscosity and illustrated in Fig. 9. Fig. 10 (left) shows that in all our numerical simulations, the metallic droplet is stable (i.e. no break up occurs) even though some models are in the our study, which focuses on small droplets in the iron rain scenario.

Influence of the viscosity ratio

The viscosity ratio is an important parameter controlling the dynamics of diapir descent, both its shape (Qaddah et al., 2019) and stability [START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. In this section, we focus on the influence of the viscosity ratio Rµ (=µmet/µsil) by varying the value of the silicate viscosity for a given metal viscosity. Changing the viscosity of metal affects its ability to deform and, hence, its terminal velocity, but this effect is not as strong as the the sinking velocity scales with µsil -1 (Eq. ( 9)) while in the intermediate regime, the sinking velocity scales with CD -1 . Given the expression of the drag coefficient (derived from [START_REF] Samuel | A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans[END_REF]; see Section 3.2), and the expression of Re (Eq. ( 6)), the sinking velocity in the intermediate regime scales with (ρsilµsil) -0.5 , as given by Eq. ( 10). Hence, Rµ has a greater effect on Stokes flows than intermediate flows, because the changes in diapir velocity are more pronounced for Stokes flows, as shown in Fig. 5.

Although the evolutions of RC as a function of Rµ are quite similar in the five panels of Fig. 12, some differences are noticeable. For the lowest values of µmet, the RC ratio increases at high Rµ values while for high values of µmet, the RC ratio is nearly constant or decreases with high Rµ (except for D met/sil =1000). This could be explained by the fact that a higher viscosity ratio tends to stabilize the diapir and slow down its descent [START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. In a Stokes regime, it leads to an increase in the time necessary to reach the same distance, in other words a longer reaction time; but in intermediate regimes this time does not change significantly, while the diapir is less deformed, leading to a smaller surface of exchange between metal and silicate. We illustrate this point in Fig. 13, for a silicate viscosity of 0.1 Pa.s: the normalized surface of the diapir varies from 1 to 2.3 with decreasing values of Rµ, which means that low Rµ for high Re (~100) values lead to an increase of exchange surface by 130% relatively to the initial diapir surface. The same effect is observed for µsil =1 Pa.s (Re~10), to a lesser degree: for low metallic viscosity (Rµ = 0.1) the diapir surface increases by 20%. This leads to an increase in Rc value for all partition coefficients except for D met/sil = 1000, for which an increase in diapir surface has a negligible effect, the deformation occurring over a short timescale at a low reaction rate.

This effect is weak in our study compared to the effect of Re and D met/sil (see Section 3.4), because the We range of our study does not allow major diapir deformation. For the deformation rate to have a strong impact on chemical equilibrium, the We values have to be much higher (We> 10, see [START_REF] Lherm | Small-Scale Metal/Silicate Equilibration During Core Formation: The Influence of Stretching Enhanced Diffusion on Mixing[END_REF].

Influence of the partition coefficient

We have shown in Figs. 11 and 12 that chemical equilibration was less efficient when the value of the viscosity ratio between the iron and silicate phases was higher. This conclusion also stands for the whole range of partition coefficients used in our models. As detailed in Section 2.2, the higher the value of D met/sil , the closer the conditions are to equilibrium conditions, and the lower the reaction rate (Ri) will be. However, the thermodynamical equilibrium is reached faster, which limits the amount of exchange between the metallic droplet and the surrounding silicate. Therefore, the contaminated volume represented by RC decreases when D met/sil increases (see Figs. 11 and 12). This effect is relatively small: an increase of 3 orders of magnitude decreases Rc by less than one order of magnitude.

Parameterization of the chemical contamination

As shown in the previous sections, the contamination of silicate by a siderophile element depends on the Reynolds number, viscosity ratio, and the values of metal-silicate partition coefficients. To use our results in an iron rain scenario following an impact during planetary formation, it is necessary to quantify the relative effect of each parameter on chemical contamination. In section 3.1 and 3.2 we chose a distance of 10×RFe and the relevant time (Fig. 5) at which RC is considered to have reached a stationary value. With this assumption the error is limited, since RC tends to reach a plateau in all our simulations (see Fig. 8). Therefore, we consider the RC values obtained from our models at a distance of 10×RFe in our parameterization.

We show in Sections 3.2 and 3.3 that increasing the Reynolds number and viscosity ratio both tend to decrease RC. In Figs. 11 and 12, log(Rc) decreases linearly with both log(Re) and log(Rµ) independently of the value of D met/sil with a small error in the prediction (R 2 > 0.85 -0.9). Concerning the effect of partition coefficients, (Section 3.4), all the curves in Figs. 11 and 12 are parallel, showing that there is a linear relationship between log(Rc) and log(D met/sil ). We parameterize the evolution of log(RC) as a function of logRe, logRµ and logD met/sil . The method used here is a multi-linear regression fit on all the data retrieved from the simulations:

logRC = a×logD met/sil + b×logRe + c×logRµ + d (11)
where a, b, c and d are constants fitted to the data by least-square regressions. The values of the parameters are presented in Table 2. As shown in Fig. 14, the parameterization of Eq. ( 11) shows a relatively good fit to the data retrieved from calculations. This kind of fit is not perfect, and could certainly be improved by including non-linear dependencies. It is also worth noting that this parameterization is valid for a narrow range of Re and We values, and its precision is lower at high Re values (see section 2.4). For instance, different results can be found in turbulent cases (high Re, [START_REF] Deguen | Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans[END_REF] and/or with deformed diapirs (high We, [START_REF] Lherm | Small-Scale Metal/Silicate Equilibration During Core Formation: The Influence of Stretching Enhanced Diffusion on Mixing[END_REF]. However, this form of equation (Eq. ( 11)) is practical to use in models of planetary formation with limited error on the value of RC, therefore we chose this form to simplify the calculations in sections 5 and 6.

This parameterization shows that the main parameters controlling the contamination of a liquid silicate reservoir by a siderophile element during the sinking of a metallic droplet are the Reynolds number and the metal-silicate partition coefficients.

The higher the D met/sil is, the less the silicate will be contaminated, which reflects the fact that we consider here a siderophile element. The Reynolds number has the same effect as D met/sil : the higher Re is, the lower RC will be. This high absolute value of b reflects the importance of the silicate viscosity, which controls the diapir velocity, as well as the type of flow and the deformation of the diapir. The fact that the parameter b is negative shows that increasing the value of the Reynolds number will decrease RC. Indeed, decreasing the viscosity increases the Reynolds number (Eq. 6) as well as the diapir velocity (e.g. [START_REF] Samuel | A re-evaluation of metal diapir breakup and equilibration in terrestrial magma oceans[END_REF]Wacheul et al., 2014, and Fig. 5).

From a single droplet to the large impact context

In the previous section, we derived a parameterization (Eq. 11) to estimate the volume fraction of a magma ocean contaminated by a unique 1cm-radius metallic droplet. In this section we apply this parameterization to large impacts in the context of late accretion.

The parameterization was calculated for a pertinent range of metallic and silicate viscosities. Numerical simulations and laboratory experiments estimated the viscosity of silicate melts at high pressure within the range of 0.01 Pa.s to 1 Pa.s. [START_REF] Liebske | Viscosity of peridotite liquid up to 13 GPa: Implications for magma ocean viscosities[END_REF][START_REF] Karki | Viscosity of MgSiO3 liquid at Earth's mantle conditions: Implications for an early magma ocean[END_REF][START_REF] Karki | Simulation of Silicate Melts Under Pressure[END_REF]. This range is stable because the viscosity tends to increase with pressure, but decrease with increasing temperature. As for the viscosity of metallic liquid, few experiments have been made, but the experimental results also show that the range of viscosity for liquid iron alloy is stable, at around 0.01-0.1 Pa.s [START_REF] Kono | Highpressure viscosity of liquid Fe and FeS revisited by falling sphere viscometry using ultrafast X-ray imaging[END_REF]. Therefore, the effect of viscosity on the equilibrium rate calculated in our simulation can be extrapolated for a high pressure context in an entirely liquid magma ocean. The higher values (10 -1000 Pa.s) of µsil tested in our simulation can also be applied at the end of the magma ocean stage when the silicate liquid is starting to crystallize (but before the mushy stage where the flow dynamics is governed by the viscosity of the solid silicate phase).

A more limiting factor to the extrapolation of our parameterization to planetary formation is the size of our computational domain. However, it is possible under certain conditions, which are listed below.

When a differentiated impactor collides with a growing planet, two fragmentation mechanisms of its metallic core occur. The first fragmentation is induced by the impact itself that stretches and disperses the impactor's core [START_REF] Kendall | Differentiated planetesimal impacts into a terrestrial magma ocean: fate of the iron core[END_REF]. The lunar crater observations showed that the most probable impact angle is α = 45° [START_REF] Shoemaker | Interpretation of lunar craters[END_REF]. When considering an impactor with a metallic core radius of Rcore impacting a magma ocean at an angle α = 45°, it is possible to apply our parameterization:

the dispersion of the impactor's core material becomes the same as our study computational domain. The mantle fraction affected by the impact (i.e. the volume of mantle into which the metallic droplets are initially spread out) can be approximated by the volume of a cylindrical portion with an angle , a thickness h and a radius L. The volume of mantle affected is then 𝑉 𝑎𝑓𝑓𝑒𝑐𝑡𝑒𝑑 = 𝛼ℎ𝐿 2 /2. For Rimp=100 km, Rcore=50 km, and =/4, [START_REF] Kendall | Differentiated planetesimal impacts into a terrestrial magma ocean: fate of the iron core[END_REF] give h=200km, L=2000 km: we can then estimate Vaffected = 3.14 x 10 17 m 3 . After the first fragmentation induced by the impact, a second rapid fragmentation occurs leading to the formation of much smaller droplets of different sizes. [START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF] showed that metallic diapirs will fragment into droplets with a mean radius of between 4 mm and 20 mm. Assuming an average droplet radius of 1 cm [START_REF] Rubie | Mechanisms of metalsilicate equilibration in the terrestrial magma ocean[END_REF] and that all the impactor's core diapir is fragmented into these cmradius droplets, an impactor's core with radius Rcore=50 km will fragment into ≈10 20 cmradius droplets. This yields a concentration of 400 droplets/m 3 , or 1 droplet per 2.5 liters, which yields a metal/silicate volume ratio of 0.0016 within Vaffected. Currently, no scaling law describes the relation between (L, h, 𝛼) and Rimp (and as a consequence Rcore).

According to [START_REF] Kendall | Differentiated planetesimal impacts into a terrestrial magma ocean: fate of the iron core[END_REF], for Rcore=50km, L≈40Rcore and h ≈ 4Rcore and by making the assumption that these two relations are also valid for different Rcore values, we can infer that the ratio of the impactor's core volume over the affected volume is constant for any impactor size.

In the numerical models described in the previous sections, we consider a metallic droplet sinking into a large volume of molten silicates to avoid boundary effects. In our study the ratio between the volume of the iron droplet and the volume of the computational domain is 0.0019, which is relatively close to the volume ratio of 0.0016 expected after an impact.

After a large impact, the impactor's core is dispersed within a volume 𝑉 𝑎𝑓𝑓𝑒𝑐𝑡𝑒𝑑 = 𝛼ℎ𝐿 2 /2. Assuming the formation of a cloud of cm-radius droplets with an iron fraction 𝜒 𝐹𝑒 = 0.0016 within 𝑉 𝑎𝑓𝑓𝑒𝑐𝑡𝑒𝑑 , the volume where chemical exchange between the droplets and the magma ocean will occur is:

𝑉 𝑒𝑥𝑐ℎ = 𝑅 𝑐 𝑉 𝑎𝑓𝑓𝑒𝑐𝑡𝑒𝑑 (12)
Using Eq. ( 11), Tab. 2 and considering L=40Rcore and h = 4Rcore [START_REF] Kendall | Differentiated planetesimal impacts into a terrestrial magma ocean: fate of the iron core[END_REF], we can calculate the fraction of the volume of the Earth's mantle Vexch/Vmantle in which chemical exchange is likely to occur. Fig. 15 shows Vexch/Vmantle as a function of the impactor radius and illustrates the influence of the partition coefficient and the viscosity contrast between the magma ocean and the liquid iron on this ratio.

As the impactor size increases, the mantle volume chemically affected by the impact increases by 𝑅 𝑖𝑚𝑝 3 . This relation is constrained by the simple formula we have used to relate Vaffected to Rcore. Given the current knowledge of the impact-induced fragmentation mechanisms of an iron core on a planetary scale, this simplification is a first step toward a global understanding of the chemical equilibration occurring after a large impact. Fig.

(left) illustrates that an increase in partition coefficient leads to a decrease in the ratio

Vexch/Vmantle in agreement with Eq. ( 11) and the values from Tab. 2. Fig. 15 (left) shows that for Rcore=1000 km, increasing the partition coefficient by 2 orders of magnitude decreases the value of the ratio Vexch/Vmantle by a factor of 3. This effect is significant but less important than the influence of the viscosity ratio. Fig. 15 (right) shows that increasing the magma ocean viscosity relative to a constant liquid iron viscosity of 0.1 Pa.s from 0.1 Pa.s to 100 Pa.s leads to an increase in Vexch/Vmantle by a factor of more than 5. According to Fig. 15, a very large impact such as that which led to the formation of the Moon-Earth system, involving a highly viscous magma ocean, would strongly enhance the chemical equilibration between the impactor's core and the Earth's mantle. This result needs to be nuanced, since a large impact is likely to increase temperature and therefore lower the viscosity of the magma ocean.

Signature of a large impact on the mantle composition

This section aims at improving previous accretion models (e.g. [START_REF] Fischer | High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O[END_REF][START_REF] Siebert | Metal-silicate partitioning of Ni and Co in a deep magma ocean[END_REF][START_REF] Wood | Core formation and the oxidation state of the Earth: Additional constraints from Nb, V and Cr partitioning[END_REF] by accounting for a possible chemical disequilibrium at a given depth. Previous models assumed that at a given step in accretion, equilibrium is reached throughout the entire mantle [START_REF] Wood | Core formation and the oxidation state of the Earth: Additional constraints from Nb, V and Cr partitioning[END_REF][START_REF] Siebert | Metal-silicate partitioning of Ni and Co in a deep magma ocean[END_REF][START_REF] Boujibar | Metal-silicate partitioning of sulphur, new experimental and thermodynamic constraints on planetary accretion[END_REF][START_REF] Clesi | Effect of H2O on metal-silicate partitioning of Ni, Co, V, Cr, Mn and Fe: Implications for the oxidation state of the Earth and Mars[END_REF]. This equilibrium is set by the final pressure on reaching equilibrium, usually at a depth corresponding to 50 % of the core mantle boundary depth (or close to this depth, for instance in [START_REF] Rubie | Accretion and differentiation of the terrestrial planets with implications for the compositions of early-formed Solar System bodies and accretion of water[END_REF], in order to fit the Ni and Co partitioning behavior [START_REF] Bouhifd | The effect of pressure on partitioning of Ni and Co between silicate and iron-rich metal liquids: a diamond-anvil cell study[END_REF]2011;[START_REF] Fischer | High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O[END_REF].

In this section and the following one, the hypothesis is the same as in previous models of accretion, but with the notable exception that not all of the mantle reaches equilibrium.

The amount of mantle equilibrated is defined using Eqs. ( 11) and ( 12), i.e. our parameterization is applied directly to a classical model of equilibration in a magma ocean, thus adding one step of complexity to previous models.

5.1.

Geochemical elements of interest

To estimate the chemical equilibrium between the impactor's core and the impacted mantle, we focus our calculations on the chemical behavior of moderately siderophile elements such as Ni and Co for which D met/sil is strongly dependent on pressure [START_REF] Bouhifd | The effect of pressure on partitioning of Ni and Co between silicate and iron-rich metal liquids: a diamond-anvil cell study[END_REF], 2011;[START_REF] Siebert | Metal-silicate partitioning of Ni and Co in a deep magma ocean[END_REF][START_REF] Fischer | High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O[END_REF]. These two elements are important, since the models of deep magma ocean are designed to explain their relatively high abundances in the BSE [START_REF] Drake | Determining the composition of the Earth[END_REF]. The maximum depth of the magma ocean (around 50% the CMB depth throughout accretion) is derived from the pressure for which the partitioning behavior of Ni and Co yields a 19.05±2 ratio in the BSE [START_REF] Bouhifd | The effect of pressure on partitioning of Ni and Co between silicate and iron-rich metal liquids: a diamond-anvil cell study[END_REF], 2011). In the previous sections, we considered a generic, moderately siderophile element with a metal/silicate partition coefficient ranging from 1 to 1000. This range is large enough to extrapolate the trend derived from our simulations to the behavior of Ni and Co at high pressure. For Ni, this partition coefficient ranges between ≈400 for P =1 GPa, and ≈ 20 for 40<P<60 GPa. For Co, the partition coefficient ranges between ≈ 100 for P = 1 GPa and ≈ 20, for 40<P<60 GPa [START_REF] Bouhifd | Convergence of Ni and Co metal-silicate partition coefficients in the deep magma-ocean and coupled silicon-oxygen solubility in iron melts at high pressures[END_REF] and references therein). In the models presented below, we assume that the pressure of equilibrium after an impact corresponds to the last stages of equilibrium in a deep magma ocean, at between 50 and 60 GPa (see Appendix B in [START_REF] Clesi | Effect of H2O on metal-silicate partitioning of Ni, Co, V, Cr, Mn and Fe: Implications for the oxidation state of the Earth and Mars[END_REF] for the corresponding calculations and the full evolution of pressure, temperature and fO2 during the accretion) with values of Ni and Co partition coefficients of between 20 and 80 for the corresponding fO2. In the following section, we consider that the iron content of the mantle is close to the BSE content (8% wt) and is not affected by the impact. The fO2 is then considered constant relatively to the iron-wüstite buffer and its value is logfO2 ≈ -2 ΔIW.

In this case, the Ni and Co partition coefficients have converged toward the same value:

D met/sil ≈ 20, which is used in the following calculations.

5.2.

Chemical signature of a large impact on the Ni/Co ratio

Here we characterize the effect of a single impactor on the Ni/Co ratio in the mantle. We We consider that the silicate fraction of the impactor is lost during impact and that only its metallic core is trapped within the impacted mantle. The impactor has a radius of between 10 km and 1000 km and the same composition as the Earth's core [START_REF] Mcdonough | 3.16-Compositional model for the Earth's core[END_REF]. Given the range of impactor size, the mass accreted is negligible compared to the Earth's mass: with a mean density of 9400 kg.m -3 , a metallic impactor of 1000 km radius adds a total mass of ~ 3.10 22 kg, which is 0.5% of the Earth's mass. This choice of size and composition for the impactor allows us to focus on the change in mantle concentration, without changing the final core concentration, especially the light element (Si and S) content, in the core.

Furthermore, this impactor composition, and its relatively low mass, would not significantly change the oxygen fugacity and thus does not affect the partitioning behavior of Ni and Co.

This section aims at characterizing the partitioning behavior of Ni and Co as a function of the impactor radius and the viscosity contrast Rµ computing Eqs. ( 11) and ( 12). The first term of Eq. ( 11) illustrates the influence of the partition coefficients for Ni and Co. These partition coefficients are controlled mainly by pressure and temperature, as well as the composition of the metallic phase, and oxygen fugacity [START_REF] Bouhifd | The effect of pressure on partitioning of Ni and Co between silicate and iron-rich metal liquids: a diamond-anvil cell study[END_REF], 2011[START_REF] Siebert | Metal-silicate partitioning of Ni and Co in a deep magma ocean[END_REF][START_REF] Fischer | High pressure metal-silicate partitioning of Ni, Co, V, Cr, Si, and O[END_REF]. To obtain partitioning behavior independent of oxygen fugacity, we use the exchange partition coefficient Kd, which is the thermodynamical constant of the reaction:

𝑀𝑂 𝑛/2 + 𝑛 2 𝐹𝑒 → 𝑛 2 𝐹𝑒𝑂 + 𝑀 ( 13
)
Where M is the element considered (in our case Ni or Co), and n is the valence of the element M (in our case, n= 2). In reaction ( 13), the oxides are in the silicate phase, while the reduced element is in the metallic phase. The exchange partition coefficient is given by

𝐾 𝑑 = 𝑐 𝑀 𝑐 𝑀𝑂 𝑛/2 × ( 𝑐 𝐹𝑒𝑂 𝑐 𝐹𝑒 ) 𝑛 2 ⁄
where 𝑐 𝑀 and 𝑐 𝐹𝑒 are the molar fractions of M and Fe in the metallic phase, and 𝑐 𝐹𝑒𝑂 and 𝑐 𝑀𝑂 𝑛/2 are the molar fractions of FeO and 𝑀𝑂 𝑛/2 in the silicate phase. To estimate the values of the partition coefficients, we use the parameterization of Kd from [START_REF] Clesi | Effect of H2O on metal-silicate partitioning of Ni, Co, V, Cr, Mn and Fe: Implications for the oxidation state of the Earth and Mars[END_REF], where we ignore the effect of water, as well as the carbon content in the metallic phase. The exchange partition coefficient, Kd, is given by:

𝑙𝑜𝑔𝐾 𝑑 𝑒𝑙𝑒𝑚𝑒𝑛𝑡 = 𝑏 𝑒𝑙 𝑇 + 𝑐 𝑒𝑙 𝑃 𝑇 + 𝑒 𝑒𝑙 𝑙𝑜𝑔(1 -𝜒 𝑆𝑖 𝑚𝑒𝑡𝑎𝑙 ) + 𝑔 𝑒𝑙 𝑙𝑜𝑔(1 -𝜒 𝑆 𝑚𝑒𝑡𝑎𝑙 ) + ℎ 𝑒𝑙 (14)
The parameters bel, cel, eel, gel and hel are detailed in [START_REF] Clesi | Effect of H2O on metal-silicate partitioning of Ni, Co, V, Cr, Mn and Fe: Implications for the oxidation state of the Earth and Mars[END_REF], and are different for each element studied. P and T are the pressure (in GPa) and temperature (in Kelvin) of equilibrium, respectively, corresponding to the pressure at 50% of CMB depth (i.e., 1450 km), and the liquidus temperature associated with it [START_REF] Andrault | Solidus and liquidus profiles of chondritic mantle: Implication for melting of the Earth across its history[END_REF]. We assume here an equilibrium pressure of 63 GPa, and an equilibrium temperature of 3450 K. 𝜒 𝑆𝑖 𝑚𝑒𝑡𝑎𝑙 and 𝜒 𝑆 𝑚𝑒𝑡𝑎𝑙 are the mass fractions of silicon and sulfur in the impactor. We assume here that 𝜒 𝑆𝑖 𝑚𝑒𝑡𝑎𝑙 =0.06 and 𝜒 𝑆 𝑚𝑒𝑡𝑎𝑙 =0.019 [START_REF] Mcdonough | 3.16-Compositional model for the Earth's core[END_REF]. Kd is then converted into a Nernst partition coefficient (mass ratio) following:

𝐷 𝑚𝑒𝑡/𝑠𝑖𝑙 = 𝑀 𝑜𝑥𝑖𝑑𝑒 ×𝑀 𝐹𝑒 𝑀 𝐹𝑒𝑂 ×𝑀 𝑒𝑙𝑒𝑚𝑒𝑛𝑡 × 𝐾 𝑑 𝑒𝑙𝑒𝑚𝑒𝑛𝑡 × 𝐷 𝐹𝑒 𝑚𝑒𝑡/𝑠𝑖𝑙 (15)
where Melement is the molar mass of Ni or Co, Moxide is the molar mass of NiO or CoO, and

MFe and MFeO are the molar masses of Fe and FeO respectively. 𝐷 𝐹𝑒 𝑚𝑒𝑡/𝑠𝑖𝑙 is the partition coefficient of iron and equals 13.65 in the BSE model. We consider that the impactor does not significantly change the oxygen fugacity and therefore does not influence the 𝐷 𝐹𝑒 𝑚𝑒𝑡/𝑠𝑖𝑙 value.

The second term in Eq. ( 11) characterizes the influence of the Reynolds number, which is given by Eq. ( 6). As stated in Section 4, we consider that the metal fragments into droplets with a radius of 1 cm. The sinking velocity of these droplets is given by Eq. ( 9). We assume here a silicate density of ρmantle =4500 kg.m -3 (mean density of the Earth's mantle), a metallic density of ρmetal= 9400 kg.m -3 , which is consistent with iron rich liquid alloys at these pressures [START_REF] Morard | The Earth's core composition from high pressure density measurements of liquid iron alloys[END_REF] and the current Earth's value for the acceleration of gravity (i.e., 9.81 m.s -2 ). The viscosity of the silicate phase in our calculations ranges from 0.1 to 100 Pa.s. The third term in Eq. ( 11) measures the influence of the viscosity ratio Rµ =µmet/µsil. In our calculations, we consider a constant value for the metallic viscosity (= 0.1 Pa.s), consistent with iron liquid viscosity at high pressure [START_REF] Kono | Highpressure viscosity of liquid Fe and FeS revisited by falling sphere viscometry using ultrafast X-ray imaging[END_REF].

Once the parameter Rc is calculated for Ni and Co, the impactor equilibrates with a volume of mantle defined by Eq. ( 12) in section 4, considering that Vimpactor/Vaffected = 0.0016, independent of the size of the impactor. As RC expresses the volume fraction of the silicate that has reacted with the metal, it may overestimate the volume of the mantle that has equilibrated with the metal phase. Hence the parameter RC should be considered as a firstorder approximation of the volume equilibrated. Therefore, for Eq. ( 12), the mass of mantle that is equilibrated is obtained by simple multiplication by ρmantle. for different viscosities as a function of the impactor radius are presented in Fig. 16.

Fig. 16 shows that a single impactor, even a large one, does not significantly affect the Ni/Co ratio in the mantle. The maximum variation is 0.25 %, and is obtained for a large impactor (1000 km radius) for a relatively high viscosity of the magma ocean (100 Pa.s).

For low viscosities (0.1 and 1 Pa.s), a metallic impactor weakly affects the Ni/Co ratio in the mantle (less than 0.1 % change for a 1000 km radius impactor). Such a small variation confirms that high spatial resolution models are not necessary even for large Pemesh values obtained for low viscosities. For a relative error of 100% in the values of Rc , at a viscosity of 0.1 Pa.s, it would only induce a variation in Ni/Co of between 0 and 0.2% within the final model.

Eq. ( 16) shows that the final concentration is proportional to the volume equilibrated if mantle with a BSE composition. Therefore, when the mass of mantle equilibrated is higher, the variation in Ni or Co concentration compared to the BSE in the mantle is higher. If the disequilibrium is high, the mequilibrated is low, and therefore there is less variation in the Ni/Co ratio. As shown in Fig. 16, the amount of disequilibrium is higher for a low viscosity magma ocean, as well as for small metallic impactors, and therefore there is no change in the moderately siderophile element concentrations in the mantle under these conditions. This means that, to have a significant effect on the moderately siderophile element concentrations in the BSE, late accretion impactors need to be large and impact a viscous magma ocean in order to increase the degree of equilibrium. Even in this case, if the amount of Ni and Co in the impactor is not significantly higher than that in the Earth's core, the change in BSE composition will not be significant, as illustrated in Fig. 16.

Discussion

6.1.

Earth accretion models and magma ocean properties

In the previous section we characterized the effects of a late impact, negligible in mass, composed only of a metallic phase. However, most of the disequilibrium happens during accretion. Therefore, Eq. ( 16) is not relevant, at least not in this form, in the context of a growing planet. Assuming that the mantle of the impactor and of the growing Earth is well mixed before equilibration with the metallic fraction, we can define the concentration of the post-impact mantle, 𝜒 𝑝𝑜𝑠𝑡-𝑖𝑚𝑝𝑎𝑐𝑡

𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒

. By modifying Eq. ( 16), and knowing the composition of the impactor's metallic phase it is possible to calculate the composition of a growing planetary mantle before it equilibrates with the metallic phase of the impactor.

It is given by:

𝜒 𝑝𝑜𝑠𝑡-𝑖𝑚𝑝𝑎𝑐𝑡 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 = 𝑚 𝑚𝑎𝑛𝑡𝑙𝑒 𝑚 𝑚𝑎𝑛𝑡𝑙𝑒 -𝑚 𝑒𝑞𝑢𝑖𝑙𝑖𝑏𝑟𝑎𝑡𝑒𝑑 (𝜒 𝑓𝑖𝑛𝑎𝑙 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 - 𝑚 𝑒𝑞𝑢𝑖𝑙𝑖𝑏𝑟𝑎𝑡𝑒𝑑 𝑚 𝑚𝑎𝑛𝑡𝑙𝑒 𝜒 𝑒𝑞𝑢𝑖𝑙𝑏𝑟𝑎𝑡𝑒𝑑 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 ) (17)
Since mequilibrated is dependent on the size of the impactor core and magma ocean viscosity, it is possible to calculate the composition of the mantle for a 90% accreted Earth hit by a Mars-sized impactor (Moon-forming impact) as a function of the impactor's core radius and viscosity of the magma ocean. Fig. 17 shows the corresponding results for comparison with Fig. 16: in this case, the greater the degree of equilibrium, the lower the initial Ni/Co ratio.

The results displayed in Fig. 17 show that it is possible to reach the same final Ni/Co ratio for the Earth's mantle with drastically different properties of the impactor and the magma ocean. From Fig. 17, two post-impact models can be derived:

-A model with a low Rcore/Rimpactor ratio and low magma ocean viscosity, which yields a lower degree of equilibrium between metal and silicate. In this kind of model, the composition of the Earth's mantle is more dependent on the equilibrium conditions between the metallic phase and the silicate phase than on the composition of the impactor. For instance, the models presented in e.g. [START_REF] Burbine | Determining the possible building blocks of the Earth and Mars[END_REF]), Rai & van Westrenen (2013) or [START_REF] Dauphas | Calcium-48 isotopic anomalies in bulk chondrites and achondrites: Evidence for a uniform isotopic reservoir in the inner protoplanetary disk[END_REF] are based on mixing different chondritic compositions where the main parameter fitted is the isotopic consistency, and they are therefore predominantly compositionallyderived models. For these models to be consistent, the impactor's core needs to be small (Rcore/Rimpactor < 0.4, left-hand side of Fig. 17), and the magma ocean viscosity low (µmagma ocean ≈ 0.1-1 Pa.s), so as to achieve a low equilibrium rate between metal and silicate.

-A model with a high Rcore/Rimpactor ratio, and high viscosity of the magma ocean, which yields greater equilibrium between metal and silicate. In this kind of model, the main changes in Earth mantle composition are due more to the composition of the impactor and less to the chemical reactions between metal and silicate. For instance, the models presented in e.g. [START_REF] Wood | Core formation and the oxidation state of the Earth: Additional constraints from Nb, V and Cr partitioning[END_REF][START_REF] Siebert | Metal-silicate partitioning of Ni and Co in a deep magma ocean[END_REF][START_REF] Boujibar | Metal-silicate partitioning of sulphur, new experimental and thermodynamic constraints on planetary accretion[END_REF] are based on the metal-silicate partitioning behavior, where the discriminating parameter is the relative core/mantle abundances of siderophile elements, and they are therefore predominantly equilibrium-derived models. In this kind of model, a minimal equilibrium needs to be reached, and therefore the impactor's core needs to be large (Rcore/Rimpactor > 0.4, right-hand side of Fig. 17), and/or the magma ocean viscosity needs to be high (µmagma ocean ≈ 10-100 Pa.s).

showed that the volume of silicate equilibrated with a small volume of metal is mostly controlled by the partition coefficients and the viscosity of the silicate. The viscosity ratio between metal and silicate is a second order parameter but still has a significant effect.

The scope of this study does not include other elements that have different partitioning behavior (V, Cr, Mn etc.) and future models of equilibrium between metal and silicate will have to integrate these elements to get a clearer view of chemical exchanges within the Early Earth. Furthermore, the range of Re and We covered in this study does not cover the whole range of Re and We numbers during the accretion event (Fig. 10). Therefore, in the future, we will need to model more complex flows (such as shown in Samuel, 2012 or [START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF], to widen the scope of our interpretations (see Section 4 for the extrapolation conditions). Future models should also account for the thermal equilibration during metal/silicate separation [START_REF] Wacheul | Experiments on fragmentation and thermo-chemical exchanges during planetary core formation[END_REF] to characterize the effects of temperature on the partition coefficient values. Other phenomena, while important, are ignored due to being beyond the scope of the study, for instance the effect of the droplet size and shape (Qaddah et al., 2019), or stretch-enhancing diffusion [START_REF] Lherm | Small-Scale Metal/Silicate Equilibration During Core Formation: The Influence of Stretching Enhanced Diffusion on Mixing[END_REF], or the possibility of large core-merging events [START_REF] Landeau | Core merging and stratification following giant impact[END_REF]. All these phenomena are beyond the scope of this study and should be taken into account in further work and models of accretion.

Nevertheless, we proposed a parameterization of the disequilibrium between a silicate magma ocean and a metallic droplet for moderately siderophile element behavior. Using

Captions of Figures and Tables

Figure 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 [START_REF] Rubie | Heterogeneous accretion, composition and core-mantle differentiation of the Earth[END_REF][START_REF] Deguen | Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans[END_REF][START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]Wacheul andLe Bars, 2017,2018). Each droplet interacts with the surrounding silicate, exchanging heat and chemical elements, until it reaches thermo-chemical equilibrium with its environment. 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 28x RFe (see Fig. 2 for the calibrated geometry). 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 (Φls = 0.5). The time necessary for the diapir to reach the distance of 10RFe is 12.1 s.

  the conservation of mass for an incompressible fluid with 𝒖 the flow velocity vector. Equation 2 describes the conservation of momentum, with acceleration of the fluid ( 𝜕𝒖 𝜕𝑡 ) and an inertia term (𝜌(𝒖 • ∇)𝒖) where ρ is fluid density. The first term on the right-hand side of Eq. 2 (-∇𝑃) is the effect of the dynamic pressure P on the fluid. The second term (∇ • [𝜇(∇𝒖 + (∇𝒖) 𝑇 )]) is the component representing the effect of viscous forces on the fluid with μ the fluid viscosity. The third term (𝜌𝒈) is the gravitational force that applies to the whole domain with g the gravitational acceleration vector. The fourth term 𝑭 𝒔𝒕 = ∇ • (𝜎(𝑰 -𝒏𝒏 𝑇 )𝛿) is the surface tension force with σ the surface tension coefficient, I the identity matrix, 𝒏 a unit vector normal to the surface of interest and δ a smeared out Dirac function located at the interface. Since the volume of metal is small, we ignored the Coriolis forces in Eq. 2.

  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 𝜇 = 𝜇 𝑠𝑖𝑙 + Φ(𝜇 𝑚𝑒𝑡 -𝜇 𝑠𝑖𝑙 ) and 𝜌 = 𝜌 𝑠𝑖𝑙 + Φ(𝜌 𝑚𝑒𝑡 -𝜌 𝑠𝑖𝑙 ), 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[START_REF] Olsson | A conservative level set method for two phase flow[END_REF]. Furthermore, this method has been already validated and applied to a similar problem inQaddah et al., 2019.

  Fig.3illustrates the time evolution of the ratio between the calculated volume over the theoretical volume of one droplet () for hmax values ranging between 0.06 and 0.2 cm.

  Stokes regime and others are in the Intermediate flow regime. In Fig.10(right) we compare the range of values with the range of We and Re values relevant to the geological context of a liquid metallic droplet sinking through a terrestrial magma ocean[START_REF] Bars | Laboratory experiments on the breakup of liquid metal diapirs[END_REF]. As shown in Fig.10(right), our range of values is limited compared to possible geophysical values. Larger diapirs leading to high Re and We values are not considered in

  Fig. 11 illustrates the influence of the Re number on the silicate volume chemically

  consider an initially fully accreted Earth. The reference composition of the impacted mantle is set to the BSE composition defined in McDonough and Sun, (1995) (Ni/Co = 19.05, 2000 ppm of Ni and 105 ppm of Co in the mantle).

Figure 2 :

 2 Figure 2: Schematic representation of the initial conditions and geometry of our models.

Figure 3 :

 3 Figure 3: Evolution of the metallic volume ratio (Vcalculated/Vtheoretical) during a simulation

Figure 4 :

 4 Figure 4: Evolution of the velocity field, r-component (on the left) and z-component (on

Figure 5 :

 5 Figure 5: Evolution of the volume fraction of the silicate during the simulation. When

Figure 6 :

 6 Figure 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

Figure 7 :

 7 Figure 7: Evolution of the volume of silicate undergoing reaction during the simulation.

Figure 8 :

 8 Figure 8: Evolution of RC as a function of time in the reference case defined in the text (µsil = µmet = 100 Pa.s). The four lines illustrate four different values of D met/sil ranging from 1 (black) to 10 3 (blue).

Figure 9 :

 9 Figure 9: Evolution of the time necessary for the diapir to move through 10RFe as a

Figure 10 :

 10 Figure 10: On the left: Re as a function of We obtained from our numerical simulations

Figure 11 :

 11 Figure 11: Evolution of RC as a function of Re during the sinking of a 1cm-radius droplet

Figure 12 :

 12 Figure 12: Evolution of RC as a function of Rµ for 5 different metallic viscosities (shown in the 5 boxes), different values of D met/sil and obtained at a distance of 10RFe. From top

Figure 13 :

 13 Figure 13: Evolution of log RC as a function of the normalized surface (left panel) and

Figure 14 :

 14 Figure 14: logRc calculated from the parameterization of Eq (11) vs logRc obtained from the simulations. The red line is the 1:1 line (perfect fit), and the dashed lines bound the 5σ confidence interval.

Figure 15 :

 15 Figure 15: Mantle fraction chemically contaminated by a moderately siderophile element after metallic core fragmentation as a function of the impactor size. In these figures we compute Eqs. 6, 11 and 12 using RFe=1cm, vdiapir=vstokes and μmet=0.1 Pa.s. The left panel shows the influence of the partition coefficient for a fixed viscosity ratio Rμ=1. The right panel shows the influence of the viscosity ratio for a fixed partition coefficient D met/sil =100.

Figure 16 :.

 16 Figure 16: Evolution of Ni/Co ratio in the mantle after an impact as a function of the

  The parameter RC is the ratio of the silicate volume that has reacted with the metal to the total volume of silicate. The definition is 𝑅 𝑐 = 𝑉 𝑟𝑒𝑎𝑐𝑡𝑒𝑑

	volume of silicate where X > 𝑿 𝒔𝒊𝒍
	𝑠𝑖𝑙	⁄	𝑉 𝑡𝑜𝑡𝑎𝑙 𝑠𝑖𝑙	, where 𝑉 𝑟𝑒𝑎𝑐𝑡𝑒𝑑 𝑠𝑖𝑙	is the

𝟎

; meaning RC is the volume fraction of silicate that has

  In this mass of mantle, the concentration of Ni or Co is in equilibrium and is given by 𝜒 𝑒𝑞 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 = 𝜒 𝑓𝑖𝑛𝑎𝑙 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 is the concentration of Ni or Co in the final mantle, 𝜒 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 is the initial concentration of Ni or Co in the mantle, 𝜒 𝑒𝑞 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 is the concentration at equilibrium, mequilibrated is the mass that is equilibrated and mmantle is the mass of the Earth's mantle. The final Ni/Co ratio is then derived from 𝜒 𝑓𝑖𝑛𝑎𝑙 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 for Ni and Co. The results of calculations

	where,
	𝑚𝑒𝑡𝑎𝑙 𝜒 𝑖𝑚𝑝𝑎𝑐𝑡𝑜𝑟 𝐷 𝑚𝑒𝑡/𝑠𝑖𝑙 . The final mass fraction of Ni or Co is given by:
	𝜒 𝑓𝑖𝑛𝑎𝑙 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 = 𝜒 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 + (𝜒 𝑒𝑞 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 -𝜒 𝑖𝑛𝑖𝑡𝑖𝑎𝑙 𝑠𝑖𝑙𝑖𝑐𝑎𝑡𝑒 )	𝑚 𝑒𝑞𝑢𝑖𝑙𝑖𝑏𝑟𝑎𝑡𝑒𝑑 𝑚 𝑚𝑎𝑛𝑡𝑙𝑒	(16)

Acknowledgments

This research 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. xxx. This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 681835 FLUDYCO-ERC-2015-CoG). This project has received funding from the European Research Council (ERC) under the European Union's Horizon 2020 research and innovation programme (grant agreement 716429). We thank Renaud Deguen for his fruitful comments.

To summarize, our results cannot constrain the viscosity of the magma ocean (especially given the error in our calculations for low viscosity calculations, see Section 2.4), but for a given model, our results can constrain the range of magma ocean viscosities in order to reach the BSE concentrations for moderately siderophile elements.

However, this work needs to be done for each accretion step or impact, and for each element, which would increase the number of plausible scenarios to accrete the Earth, but could also provide compositional constraints on irreconcilable scenarios, and therefore allow some impactor compositions to be excluded from a given model.

6.2.

Effect of Moon formation on the Ni/Co ratio in the Earth's mantle

In the previous section we showed that the properties of the magma ocean are model dependent: they have to be inferred from the impactor properties and the final BSE content chosen as a control of the model output. In this section, we infer the viscosity of a magma ocean for one particular model: the Moon formation by an impact with a Mars-sized impactor [START_REF] Canup | Simulations of a late lunar-forming impact[END_REF]. The Moon forming impactor (named Theia) is believed to be a large Mars-sized impactor which hit the proto-Earth at the end of accretion with an oblique trajectory [START_REF] Canup | Origin of the Moon in a giant impact near the end of the Earth's formation[END_REF][START_REF] Canup | Simulations of a late lunar-forming impact[END_REF]. In a simulation of this kind of impact, most of the metal and silicate of the impactor merges with the Earth [START_REF] Canup | Simulations of a late lunar-forming impact[END_REF]). If we consider mass of the Moon to be negligible compared to that of Theia, and ignore the loss of particles due to the impact, we can determine that the final 10% mass accreted to the Earth was due to the impact with Theia. The core of this impactor is 30% of the total mass of the impactor [START_REF] Canup | Simulations of a late lunar-forming impact[END_REF], with the same mean density as the Earth, so it is not improbable that the ratio Rcore/Rimpactor is comparable to that of the Earth or Mars today (given by the red and blue vertical lines in Figs. 17). In that case, Eq. ( 17) and Fig. 17 shows that some equilibrium occurred between the magma ocean formed by the impact and the impactor's core sinking into the Earth's core. For every magma ocean viscosity tested here, it is necessary to have interaction between Theia's core and the proto-Earth's mantle to get the final BSE concentration. For a magma ocean viscosity of 0.1 Pa.s, the proto-Earth's Ni/Co ratio is between 0.25 % and 0.45 % lower than the BSE values for a Rcore/Rimpactor comparable to those of Mars and the Earth respectively. For a magma ocean viscosity of 100 Pa.s, the proto-Earth Ni/Co ratio is between 1.45% and 3.00 % lower than the BSE values for a Rcore/Rimpactor comparable to those of Mars and the Earth, respectively. This means that the Moon forming impact had a significant effect on the BSE Ni/Co ratio, which can be high if the magma ocean formed after the impact was viscous (100 Pa.s) and the core radius of Theia was around 1850 km.

Conclusion

We developed numerical models to characterize the effects of viscosity and partition coefficients on the metal/silicate equilibrium for moderately siderophile elements. We this parameterization on Ni and Co behavior at the end of accretion, it is possible to define some constraints on the viscosity of the magma ocean. For a given model of accretion to fit the current BSE concentrations, the viscosity of the magma ocean and the impactor composition need to be changed accordingly: large impactors with large metallic cores (Rcore > 1100 km) and high magma ocean viscosity (10 -100 Pa.s) favor a high degree of metal/silicate equilibrium; while metal/silicate disequilibrium models can be achieved with moderately small metallic cores (Rcore ≈ 300-1000 km) and low viscosity (0.1 -1 Pa.s) magma ocean.

When applying our results to the Moon-forming impact, we showed that the current Ni/Co ratio in the BSE is affected by the Moon's formation, depending on the viscosity of the magma ocean formed after the impact. The Moon's formation can account for 0. The Rcore/REarth current ratio is shown by the blue dashed line, and the Rcore/RMars ratio by a red dashed line. In the calculation, the impactor is the same size as Mars (Rimpactor = 3390 km), and the impactor core radius ranges from 340 km to 2040 km. We represent here the variation in Ni/Co ratio before the impact relative to the BSE concentration but for a wellmixed mantle after impact: 𝛥 ( Table 1: 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 2. Values of fitted parameters for Eq. 11 using the least-square methods. The errors associated with the value are given by the 1σ error value. The data used for parameterization are the same as the data presented in sections 3.3 and 3.4. The comparison between calculated values from Eq. 11 and the values obtained after simulation is presented in Fig. 13.