<figure><img src="image_113.png" /><figcaption>Fig. 9. Evolution of the time necessary for the diapir to move through 10RFe as a function of the silicate viscosity. Black points are the values corresponding to each study. The red area corresponds to Stokes flows, and the blue area to intermediate flows. The clusters of black points correspond to the variation in velocity due to the variation in metal viscosity. The theoretical values are plotted in red and derived from Eqs. (9) and (10). Typical times of diapir equilibrium, from Wacheul et al. (2014), are shown in blue. The size of the marker is scaled logarithmically on the chemical Peclet number, ranging from 8 (smallest) to 3000 (largest). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

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 (Samuel, 2012 and 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:

$$v_{\text{Stokes}} = \frac{2 (\rho_{\text{met}} - \rho_{\text{sil}}) g R_{\text{Fe}}^2}{9 \mu_{\text{sil}}} \qquad (9)$$

For intermediate flows, the sinking velocity is given by:

$$v_{\text{intermediate}} = \sqrt[2]{\frac{(\rho_{\text{met}} - \rho_{\text{sil}}) g R_{\text{Fe}}}{\rho_{\text{sil}} C_D}} \qquad (10)$$

with C_D the drag coefficient. In our case, the drag coefficient is approximated by C_D = 12/Re + 0.3, following Samuel (2012). 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 Wacheul et al. (2014) 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

<figure><img src="image_114.png" /><figcaption>Fig. 10. On the left: Re as a function of We obtained from our numerical simulations (one red circle for each calculation). The blue area shows the Intermediate regime flow, and the red area the Stokes regime flow. The dashed area illustrates the domain where diapirs are stable. On the right: Comparison between the geologically plausible values of We and Re for the Earth's magma ocean and the values covered by our study. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
a chosen value of $\mu_{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 ($v_{diapir}$) 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 \times R_{Fe}$. 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 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 (Wacheul et al., 2014). 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 our study, which focuses on small droplets in the iron rain scenario.

### 3.3. 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 (Wacheul et al., 2014). In this section, we focus on the influence of the viscosity ratio $R_\mu (=\mu_{met}/\mu_{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 change in velocity implied by a change in ambient viscosity. This behavior is consistent with the equation of diapir velocity in Stokes and intermediate flows (Eqs. (9) and (10)), and is strongly dependent on the silicate viscosity. It is visible in Fig. 9, where a decrease in silicate viscosity from 1000 Pa·s to 10 Pa·s (Stokes flows) leads to an increase in sinking velocity from $10^{-3}$ m·s$^{-1}$ to $0.1$ m·s$^{-1}$ for a metal viscosity of 1 Pa·s. For a constant silicate viscosity of 10 Pa·s, a decrease of the metallic viscosity from 1000 Pa·s to 0.1 Pa·s only increases the velocity from 0.05 to $0.1$ m·s$^{-1}$. The viscosity ratio influences the sinking velocity and, hence, affects the *Re* number and the flow regime, but not as strongly as the silicate viscosity.

Fig. 11 illustrates the influence of the *Re* number on the silicate volume chemically contaminated during the sinking of a metallic droplet with $R_{Fe} = 1$ cm. For each case, the evolution of $R_C$ is shown as a function of *Re*, and for a distance of 10 cm ($10 \times R_{Fe}$). The time corresponding to this distance is shown in Fig. 5 as a function of viscosity. Except for the case with $D^{met/sil} = 10^3$, there is no linear correlation between $\log(R_C)$ and $\log(Re)$. However, we clearly show a transition between the Stokes regime and the Intermediate regime. In the Stokes regime flow, an increasing *Re* (i.e. a decreasing value of $\mu_{sil}$) leads to a significant decrease of $R_C$, while $R_C$ values seem to reach a plateau when the flow reaches the intermediate flow regime. Fig. 11 also shows that the approximation made by Rubie et al. (2003) to achieve an analytical solution for equilibrium is only valid for low *Re* flows (high silicate viscosity). As the density of the silicate is not significantly altered by the reaction, the parameter $R_c$ (volume ratio of silicate equilibrated) can be compared to their parameter $R_s$ (mass ratio of silicate equilibrated). For the lowest *Re* value, there is good agreement between our numerical data and the analytical calculation. However, extrapolating the analytical solution to higher *Re* flows leads to an overestimation of the equilibrated volume.

The $R_C$ ratio is presented in Fig. 12 as a function of $R_\mu$. A decrease in $R_C$ is observed with an increase in $R_\mu$ (i.e. a decrease of $\mu_{sil}$) for a given viscosity of the metallic phase. For the five parts of Fig. 12, which correspond to the five values of $\mu_{met}$, the trend and the magnitude of the values of $R_c$ are quite similar, meaning that the effect of the viscosity ratio is less important than the ambient viscosity, which governs the flow regime. For high values of $R_\mu$, the same plateau phenomenon as seen in Fig. 11 is observed, which corresponds to lower silicate viscosity values (intermediate flow). In the Stokes regime, the sinking velocity

<figure><img src="image_115.png" /><figcaption>Fig. 11. Evolution of $R_C$ as a function of *Re* during the sinking of a 1 cm-radius droplet for an adimensional time of 1. The Stokes flow regime (red area) and intermediate flow regime (blue area) are separated by a black vertical line at $Re = 1$. The dot color characterizes different values for $D^{met/sil}$ ranging from 1 to $10^3$. The size of the marker is scaled logarithmically on the chemical Peclet number, ranging from 8 (smallest) to 3000 (largest). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

scales with $\mu_{sil}^{-1}$ (Eq. (9)) while in the intermediate regime, the sinking velocity scales with $C_D^{-1}$. Given the expression of the drag coefficient (derived from Samuel, 2012; see Section 3.2), and the expression of *Re* (Eq. (6)), the sinking velocity in the intermediate regime scales with $(\rho_{sil}\mu_{sil})^{-0.5}$, as given by Eq. (10). Hence, $R_\mu$ 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 $R_C$ as a function of $R_\mu$ are quite similar in the five panels of Fig. 12, some differences are noticeable. For the lowest values of $\mu_{met}$, the $R_C$ ratio increases at high $R_\mu$ values while for high values of $\mu_{met}$, the $R_C$ ratio is nearly constant or decreases with high $R_\mu$ (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 (Wacheul et al., 2014). 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_\mu$, which means that low $R_\mu$ 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 $\mu_{sil} = 1$ Pa·s (*Re* ~ 10), to a lesser degree: for low metallic viscosity ($R_\mu = 0.1$) the diapir surface increases by 20%. This leads to an increase in $R_C$ 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 Lherm and Deguen, 2018).
<figure><img src="image_116.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_117.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_118.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_119.png" /><figcaption>e</figcaption></figure>

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

(caption on next page)
Fig. 12. Evolution of R<sub>C</sub> as a function of R<sub>μ</sub> for 5 different metallic viscosities (shown in the 5 boxes), different values of D<sup>met/sil</sup> and obtained at a distance of 10R<sub>Fe</sub>. From top left to bottom panel the metallic viscosity increases from 0.1 Pa·s to 1000 Pa·s. The silicate viscosities corresponding to R<sub>μ</sub> are indicated on the top x-axis. The size of the marker is scaled logarithmically with the chemical Peclet number, ranging from 8 (smallest) to 3000 (largest).

### 3.4. 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<sup>met/sil</sup>, the closer the conditions are to equilibrium conditions, and the lower the reaction rate (R<sub>i</sub>) 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 R<sub>C</sub> decreases when D<sup>met/sil</sup> increases (see Figs. 11 and 12). This effect is relatively small: an increase of 3 orders of magnitude decreases R<sub>c</sub> by less than one order of magnitude.

### 3.5. 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 × R<sub>Fe</sub> and the relevant time (Fig. 5) at which R<sub>C</sub> is considered to have reached a stationary value. With this assumption the error is limited, since R<sub>C</sub> tends to reach a plateau in all our simulations (see Fig. 8). Therefore, we consider the R<sub>C</sub> values obtained from our models at a distance of 10 × R<sub>Fe</sub> in our parameterization.

We show in Sections 3.2 and 3.3 that increasing the Reynolds number and viscosity ratio both tend to decrease R<sub>C</sub>. In Figs. 11 and 12, log(R<sub>c</sub>) decreases linearly with both log(Re) and log(R<sub>μ</sub>) independently of the value of D<sup>met/sil</sup> with a small error in the prediction (R<sup>2</sup> > 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(R<sub>c</sub>) and log(D<sup>met/sil</sup>). We parameterize the evolution of log(R<sub>C</sub>) as a function of logRe, logR<sub>μ</sub> and logD<sup>met/sil</sup>. The method used here is a multi-linear regression fit on all

the data retrieved from the simulations:

$$ \log R_C = a \times \log D^{\text{met/sil}} + b \times \log Re + c \times \log R_{\mu} + d \quad (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, Deguen et al., 2014) and/or with deformed diapirs (high We, Lherm and Deguen, 2018). However, this form of equation (Eq. (11)) is practical to use in models of planetary formation with limited error on the value of R<sub>C</sub>, 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<sup>met/sil</sup> 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<sup>met/sil</sup>: the higher Re is, the lower R<sub>C</sub> 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 R<sub>C</sub>. Indeed, decreasing the viscosity increases the Reynolds number (Eq. (6)) as well as the diapir velocity (e.g. Samuel, 2012; Wacheul et al., 2014, and Fig. 5).

### 4. 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 1 cm-radius metallic droplet. In this section we apply this parameterization to large impacts in the context of late accretion.

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

<figure><img src="image_122.png" /><figcaption>Fig. 13. Evolution of log R<sub>C</sub> as a function of the normalized surface (left panel) and evolution of the normalized surface as a function of R<sub>μ</sub> (right panel) for models with μ<sub>silicate</sub> = 0.1 Pa·s for the reference time defined in Fig. 9. The normalized surface of the diapir is the surface of the diapir in the simulation divided by S<sub>theoretical</sub> = 4πR<sub>Fe</sub><sup>2</sup> (initial surface of the sphere).</figcaption></figure>
<caption>**Table 2**</caption>

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.

<table><thead><tr><th>Parameters</th><th>a (Dmet/sil)</th><th>b (Re)</th><th>c (R<sub>μ</sub>)</th><th>d</th></tr></thead><tbody><tr><td>Values</td><td>-0.235</td><td>-0.283</td><td>0.011</td><td>-1.686</td></tr><tr><td>1σ</td><td>0.02</td><td>0.03</td><td>0.001</td><td>0.07</td></tr></tbody></table>

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. (Liebske et al., 2005; Karki and Stixrude, 2010; Karki et al., 2018). 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 (Kono et al., 2015). 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 μ<sub>sil</sub> 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 (Kendall and Melosh, 2016). The lunar crater observations showed that the most probable impact angle is α = 45° (Shoemaker, 1962). When considering an impactor with a metallic core radius of R<sub>core</sub> 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 V<sub>affected</sub> = αhL<sup>2</sup>/2. For R<sub>imp</sub> = 100 km, R<sub>core</sub> = 50 km, and α = π/4, Kendall and Melosh (2016) give h = 200 km, L = 2000 km: we can then estimate V<sub>affected</sub> = 3.14 × 10<sup>17</sup> m<sup>3</sup>. After the first fragmentation induced by the impact, a second rapid fragmentation occurs leading to the formation of much smaller droplets of different sizes. Wacheul et al. (2014) 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 (Rubie et al., 2003) and that all the impactor's core diapir is fragmented into these cm-radius droplets, an impactor's core with radius R<sub>core</sub> = 50 km will fragment into ≈ 10<sup>20</sup> cm-radius droplets. This yields a concentration of 400 droplets/m<sup>3</sup>, or 1 droplet per 2.5 l, which yields a metal/silicate volume ratio of 0.0016 within V<sub>affected</sub>. Currently, no scaling law describes the relation between (L, h, α) and R<sub>imp</sub> (and as a consequence R<sub>core</sub>). According to Kendall and Melosh (2016), for R<sub>core</sub> = 50 km, L ≈ 40R<sub>core</sub> and h ≈ 4R<sub>core</sub> and by making the assumption that these two relations are also valid for different R<sub>core</sub> 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

<figure><img src="image_123.png" /><figcaption>**Fig. 14.** logR<sub>c</sub> calculated from the parameterization of Eq. (11) vs logR<sub>c</sub> obtained from the simulations. The red line is the 1:1 line (perfect fit), and the dashed lines bound the 5σ confidence interval. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

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 V<sub>affected</sub> = αhL<sup>2</sup>/2. Assuming the formation of a cloud of cm-radius droplets with an iron fraction χ<sub>Fe</sub> = 0.0016 within V<sub>affected</sub>, the volume where chemical exchange between the droplets and the magma ocean will occur is:

$$V_{\text{exch}} = R_c V_{\text{affected}} \qquad (12)$$

Using Eq. (11), Table 2 and considering L = 40R<sub>core</sub> and h = 4R<sub>core</sub> (Kendall and Melosh, 2016), we can calculate the fraction of the volume of the Earth's mantle V<sub>exch</sub>/V<sub>mantle</sub> in which chemical exchange is likely to occur. Fig. 15 shows V<sub>exch</sub>/V<sub>mantle</sub> 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 R<sub>imp</sub><sup>3</sup>. This relation is constrained by the simple formula we have used to relate V<sub>affected</sub> to R<sub>core</sub>. 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. 15 (left) illustrates that an increase in partition coefficient leads to a decrease in the ratio V<sub>exch</sub>/V<sub>mantle</sub> in agreement with Eq. (11) and the values from Table 2. Fig. 15 (left) shows that for R<sub>core</sub> = 1000 km, increasing the partition coefficient by 2 orders of magnitude decreases the value of the ratio V<sub>exch</sub>/V<sub>mantle</sub> 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 V<sub>exch</sub>/V<sub>mantle</sub> 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
<figure><img src="image_124.png" /><figcaption>Fig. 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 $R_{Fe} = 1$ cm, $v_{diapir} = v_{stokes}$ and $\mu_{met} = 0.1$ Pa.s. The left panel shows the influence of the partition coefficient for a fixed viscosity ratio $R_{\mu} = 1$. The right panel shows the influence of the viscosity ratio for a fixed partition coefficient $D^{met/sil} = 100$.</figcaption></figure>

temperature and therefore lower the viscosity of the magma ocean.

## 5. Signature of a large impact on the mantle composition

This section aims at improving previous accretion models (e.g. Fischer et al., 2015; Siebert et al., 2012; Wood et al., 2008) 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 (Wood et al., 2008; Siebert et al., 2012; Boujibar et al., 2014; Clesi et al., 2016). 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 Rubie et al., 2015), in order to fit the Ni and Co partitioning behavior (Bouhifd and Jephcoat, 2003, 2011; Fischer et al., 2015).

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 (Bouhifd and Jephcoat, 2003, 2011; Siebert et al., 2012; Fischer et al., 2015). These two elements are important, since the models of deep magma ocean are designed to explain their relatively high abundances in the BSE (Drake and Righter, 2002). 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 (Bouhifd and Jephcoat, 2003, 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 (Bouhifd and Jephcoat, 2011 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 Clesi et al., 2016 for the corresponding calculations and the full evolution of pressure, temperature and $f_{O_2}$ during the accretion) with values of Ni and Co partition coefficients of between 20 and 80 for the corresponding $f_{O_2}$. 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 $f_{O_2}$ is then considered constant relatively to the iron-wüstite buffer and its value is $\log f_{O_2} \approx -2 \Delta \text{IW}$. In this case, the Ni and Co partition coefficients have converged toward the same value: $D^{met/sil} \approx 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 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). 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 (McDonough, 2003). 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⁻³, a metallic impactor of 1000 km radius adds a total mass of ~3.10²² 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_{\mu}$ 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 (Bouhifd and Jephcoat, 2003, 2011; Siebert et al., 2012; Fischer et al., 2015). To obtain partitioning behavior independent of oxygen fugacity, we use the exchange partition coefficient $K_d$, which is the thermodynamical constant of the reaction:
<figure><img src="image_125.png" /><figcaption>Fig. 16. Evolution of Ni/Co ratio in the mantle after an impact as a function of the impactor size for different silicate viscosities (0.1 to 100 Pa·s). We represent here the variation of the Ni/Co ratio after the impact relative to the initial BSE concentration: $A(\frac{Ni}{Co}) = 100 \times \frac{\frac{Ni}{Co}Final - \frac{Ni}{Co}BSE}{\frac{Ni}{Co}BSE}$. Since the impact happens on the BSE, the more the Ni/Co ratio is different to that of the BSE, the more the impactor has reached equilibrium.</figcaption></figure>

$$MO_{n/2} + \frac{n}{2}Fe \rightarrow \frac{n}{2}FeO + M \quad (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 $K_d = \frac{c_M}{c_{MO_{n/2}}} \times (\frac{c_{FeO}}{c_{Fe}})^{n/2}$ where $c_M$ and $c_{Fe}$ are the molar fractions of M and Fe in the metallic phase, and $c_{FeO}$ and $c_{MO_{n/2}}$ are the molar fractions of FeO and $MO_{n/2}$ in the silicate phase. To estimate the values of the partition coefficients, we use the parameterization of $K_d$ from Clesi et al. (2016), where we ignore the effect of water, as well as the carbon content in the metallic phase. The exchange partition coefficient, $K_d$, is given by:

$$\log K_d^{element} = \frac{b_{el}}{T} + c_{el} \frac{P}{T} + e_{el} \log(1 - \chi_{Si}^{metal}) + g_{el} \log(1 - \chi_{S}^{metal}) + h_{el} \quad (14)$$

The parameters $b_{el}$, $c_{el}$, $e_{el}$, $g_{el}$ and $h_{el}$ are detailed in Clesi et al. (2016), 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 (Andrault et al., 2011). We assume here an equilibrium pressure of 63 GPa, and an equilibrium temperature of 3450 K. $\chi_{Si}^{metal}$ and $\chi_{S}^{metal}$ are the mass fractions of silicon and sulfur in the impactor. We assume here that $\chi_{Si}^{metal} = 0.06$ and $\chi_{S}^{metal} = 0.019$ (McDonough, 2003). $K_d$ is then converted into a Nernst partition coefficient (mass ratio) following:

$$D^{met/sil} = \frac{M_{oxide} \times M_{Fe}}{M_{FeO} \times M_{element}} \times K_d^{element} \times D_{Fe}^{met/sil} \quad (15)$$

where $M_{element}$ is the molar mass of Ni or Co, $M_{oxide}$ is the molar mass of NiO or CoO, and $M_{Fe}$ and $M_{FeO}$ are the molar masses of Fe and FeO respectively. $D_{Fe}^{met/sil}$ 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 $D_{Fe}^{met/sil}$ 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 $\rho_{mantle} = 4500 \text{ kg.m}^{-3}$ (mean density of the Earth's mantle), a metallic density of $\rho_{metal} = 9400 \text{ kg.m}^{-3}$, which is consistent with iron rich liquid alloys at these pressures (Morard et al., 2013) and the current Earth's value for the acceleration of gravity (i.e., 9.81 m·s⁻²). 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_\mu = \mu_{met}/\mu_{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 (Kono et al., 2015).

Once the parameter $R_c$ is calculated for Ni and Co, the impactor equilibrates with a volume of mantle defined by Eq. (12) in Section 4, considering that $V_{impactor}/V_{affected} = 0.0016$, independent of the size of the impactor. As $R_c$ 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 $R_c$ should be considered as a first-order approximation of the volume equilibrated. Therefore, for Eq. (12), the mass of mantle that is equilibrated is obtained by simple multiplication by $\rho_{mantle}$. In this mass of mantle, the concentration of Ni or Co is in equilibrium and is given by

$$\chi_{eq}^{silicate} = \frac{\chi_{impactor}^{metal}}{D^{met/sil}} \quad \text{The final mass fraction of Ni or Co is given by:}$$

$$\chi_{final}^{silicate} = \chi_{initial}^{silicate} + (\chi_{eq}^{silicate} - \chi_{initial}^{silicate}) \frac{m_{equilibrated}}{m_{mantle}} \quad (16)$$

where, $\chi_{final}^{silicate}$ is the concentration of Ni or Co in the final mantle, $\chi_{initial}^{silicate}$ is the initial concentration of Ni or Co in the mantle, $\chi_{eq}^{silicate}$ is the concentration at equilibrium, $m_{equilibrated}$ is the mass that is equilibrated and $m_{mantle}$ is the mass of the Earth's mantle. The final Ni/Co ratio is then derived from $\chi_{final}^{silicate}$ for Ni and Co. The results of calculations 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 $Pe_{mesh}$ values obtained for low viscosities. For a relative error of 100% in the values of $R_c$, 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 $\chi_{eq}^{silicate}$ is higher than $\chi_{initial}^{silicate}$. In this case, late accretion of a metallic impactor, which yields a high metallic concentration in Ni and Co (similar to the Earth core concentration given in McDonough, 2003); with a high equilibrium pressure, which yields low partition coefficient values for Ni and Co, $(\chi_{eq}^{silicate} - \chi_{initial}^{silicate})$ is positive. The impactor is metallic and rich in Ni and Co (5.2%wt and 0.24%wt, respectively), and impacts a silicate 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 $m_{equilibrated}$ 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.

# 6. 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, $\chi_{post-impact}^{silicate}$. 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:

$$ \chi_{post-impact}^{silicate} = \frac{m_{mantle}}{m_{mantle} - m_{equilibrated}} \left( \chi_{final}^{silicate} - \frac{m_{equilibrated}}{m_{mantle}} \chi_{equilibrated}^{silicate} \right) \quad (17) $$

Since $m_{equilibrated}$ 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-

<figure><img src="image_126.png" /><figcaption>Fig. 17. Ni/Co ratio in the Earth's mantle before equilibration with the core of an impactor corresponding to 10% of the Earth's mass as a function of the radius of the impactor's core normalized to the impactor size and for different magma ocean viscosities. The $R_{core}/R_{Earth}$ current ratio is shown by the blue dashed line, and the $R_{core}/R_{Mars}$ ratio by a red dashed line. In the calculation, the impactor is the same size as Mars ($R_{impactor} = 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 well-mixed mantle after impact: $\Delta(\frac{Ni}{Co}) = 100 \times \frac{(\frac{Ni}{Co})_{post-impact} - (\frac{Ni}{Co})_{BSE}}{(\frac{Ni}{Co})_{BSE}}$. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

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 $R_{core}/R_{impactor}$ 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. Burbine and O'Brien (2004), Rai and van Westrenen (2013) or Dauphas et al. (2014) are based on mixing different chondritic compositions where the main parameter fitted is the isotopic consistency, and they are therefore predominantly compositionally-derived models. For these models to be consistent, the impactor's core needs to be small ($R_{core}/R_{impactor} < 0.4$, left-hand side of Fig. 17), and the magma ocean viscosity low ($\mu_{magma ocean} \approx 0.1-1$ Pa·s), so as to achieve a low equilibrium rate between metal and silicate.

- A model with a high $R_{core}/R_{impactor}$ 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. Wood et al., 2008, Siebert et al., 2012 and Boujibar et al., 2014 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 ($R_{core}/R_{impactor} > 0.4$, right-hand side of Fig. 17), and/or the magma ocean viscosity needs to be high ($\mu_{magma ocean} \approx 10-100$ Pa·s).

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 (Canup, 2004). 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 (Canup and Asphaug, 2001; Canup, 2004). In a simulation of this kind of impact, most of the metal and silicate of the impactor merges with the Earth (Canup, 2004). 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 (Canup, 2004), with the same mean density as the Earth,
so it is not improbable that the ratio $R_{core}/R_{impactor}$ is comparable to that of the Earth or Mars today (given by the red and blue vertical lines in Fig. 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 $R_{core}/R_{impactor}$ 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 $R_{core}/R_{impactor}$ 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.

## 7. Conclusion

We developed numerical models to characterize the effects of viscosity and partition coefficients on the metal/silicate equilibrium for moderately siderophile elements. We 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 Wacheul et al., 2014), 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 (Wacheul and Le Bars, 2018) 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 (Lherm and Deguen, 2018), or the possibility of large core-merging events (Landeau et al., 2016). 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 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 ($R_{core} > 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 ($R_{core} \approx 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.25 to 0.45% of the final BSE Ni/Co ratio for a low viscosity (0.1 Pa·s) magma ocean, and for 1.45% to 3% of the final BSE Ni/Co ratio for a high viscosity (100 Pa·s) magma ocean.

## CRediT authorship contribution statement

V. Clesi:Conceptualization, Methodology, Investigation, Formal analysis, Data curation, Visualization, Writing - original draft, Writing - review & editing.J. Monteux:Conceptualization, Writing - review & editing, Resources, Visualization, Supervision, Funding acquisition.B. Qaddah:Methodology, Investigation, Visualization, Writing - review & editing.M. Le Bars:Supervision, Validation, Writing - review & editing.J.-B. Wacheul:Validation, Writing - review & editing.M.A. Bouhifd:Conceptualization, Supervision, Funding acquisition, Writing - review & editing.

## Declaration of competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

## 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. 421. 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.

## References

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 (1-2), 251-259.

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 (1-2), 245-255.

Bouhifd, M.A., Jephcoat, A.P., 2011. 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. Earth Planet. Sci. Lett. 307 (3-4), 341-348.

Bouhifd, M.A., Clesi, 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.

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

Burbine, T.H., O'Brien, K.M., 2004. Determining the possible building blocks of the Earth and Mars. Meteorit. Planet. Sci. 39 (5), 667-681.

Canup, R.M., 2004. Simulations of a late lunar-forming impact. Icarus 168 (2), 433-456.

Canup, R.M., 2008. Accretion of the Earth. Philos. Trans. R. Soc. A Math. Phys. Eng. Sci. 366 (1883), 4061-4075.

Canup, R.M., Asphaug, E., 2001. Origin of the Moon in a giant impact near the end of the Earth's formation. Nature 412 (6848), 708.

Cartier, C., Hammouda, T., Boyet, M., Bouhifd, M.A., Devidal, J.-L., 2014. Redox control of the fractionation of niobium and tantalum during planetary accretion and core formation. Nat. Geosci. 7 (8), 573.

Chang, I.D., 1961. On the wall effect correction of the Stokes drag formula for axially symmetric bodies moving inside a cylindrical tube. Z. Angew. Math. Phys. 12 (1), 6-14.

Clesi, V., Bouhifd, M.A., Bolfan-Casanova, N., Manthilake, G., Fabbrizio, A., Andrault, D., 2016. 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. Geochim. Cosmochim. Acta 192, 97-121.

Crank, J., 1975. The Mathematics of Diffusion, 2nd ed. Clarendon Press, Oxford viii, 414 p.

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 (3-4), 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.

Drake, M.J., Righter, K., 2002. Determining the composition of the Earth. Nature 416 (6876), 39.

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

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. 9 (1), 71.

Jacobsen, B., Yin, Q.-Z., Moynier, F., Amelin, Y., Krot, A.N., Nagashima, K., Hutcheon, I.D., Palme, H., 2008. 26Al-26Mg and 207Pb-206Pb systematics of Allende CAIs: Canonical solar initial 26Al/27Al ratio reinstated. Earth Planet. Sci. Lett. 272 (1-2), 353-364.

Jones, J.H., Drake, M.J., 1986. Geochemical constraints on core formation in the Earth. Nature 322 (6076), 221.

Karki, B.B., Stixrude, L.P., 2010. Viscosity of MgSiO3 liquid at Earth's mantle conditions: implications for an early magma ocean. Science 328 (5979), 740-742.

Karki, B.B., Ghosh, D.B., Bajgain, S.K., 2018. Simulation of silicate melts under pressure. In: Magmas Under Pressure, pp. 419-453.

Ke, Y., Solomatov, V.S., 2009. Coupled core-mantle thermal evolution of early Mars. J. Geophys. Eng. Planets 114 (E7).

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

Kleine, T., Rudge, J.F., 2011. Chronometry of meteorites and the formation of the Earth and Moon. Elements 7 (1), 41-46.

Kleine, T., Munker, C., Mezger, K., Palme, H., 2002. Rapid accretion and early core formation on asteroids and the terrestrial planets from Hf-W chronometry. Nature 418 (6901), 952-955.

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 (17), 5150-5188.

Kono, Y., Kenney-Benson, C., Shibazaki, Y., Park, C., Shen, G., Wang, Y., 2015. High-pressure viscosity of liquid Fe and FeS revisited by falling sphere viscometry using ultrafast X-ray imaging. Phys. Earth Planet. Inter. 241, 57-64.

Landeau, M., Olson, P., Deguen, R., Hirsh, B.H., 2016. Core merging and stratification following giant impact. Nat. Geosci. 9 (10), 786.

Lherm, V., Deguen, R., 2018. Small-scale metal/silicate equilibration during core formation: the influence of stretching enhanced diffusion on mixing. J. Geophys. Res. Solid Earth 123 (12), 10,496-10,516.

Liebske, C., Schmickler, B., Terasaki, H., Poe, B.T., Suzuki, A., Funakoshi, K.I., Funakochi, K.I., Ando, R., Rubie, D.C., 2005. Viscosity of peridotite liquid up to 13 GPa: implications for magma ocean viscosities. Earth Planet. Sci. Lett. 240 (3-4), 589-604.

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 (24), 7360-7386.

McDonough, W.F., 2003. 3.16-Compositional model for the Earth's core. In: Treatise on Geochemistry, pp. 547-568.

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

Mittal, R.C., Jain, R.K., 2012. Redefined cubic B-splines collocation method for solving convection-diffusion equations. Appl. Math. Model. 36 (11), 5555-5573.

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 (3-4), 353-362.

Monteux, J., Amit, H., Choblet, G., Langlais, B., Tobie, G., 2015. Giant impacts, heterogeneous mantle heating and a past hemispheric dynamo on Mars. Phys. Earth Planet. Inter. 240, 114-124.

Morard, G., Siebert, J., Andrault, D., Guignot, N., Garbarino, G., Guyot, F., 2013. F. and D.

Antonangeli. The Earth's core composition from high pressure density measurements of liquid iron alloys. Earth Planet. Sci. Lett. 373, 169-178.

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.

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

Olsson, E., Kreiss, G., 2005. A conservative level set method for two phase flow. J. Comput. Phys. 210 (1), 225-246.

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

Qaddah, B., Monteux, J., Le Bars, M., April 2020. Thermal Evolution of a Metal Drop Falling in a less Dense, more Viscous Fluid. Physical Review Fluids Accepted manuscript.

Rai, N., van Westrenen, W., 2013. Core-mantle differentiation in Mars. J. Geophys. Eng. Planets 118 (6), 1195-1203.

Ricard, Y., Šránek, O., Dubuffet, F., 2009. A multi-phase model of runaway core-mantle's segregation in planetary embryos. Earth Planet. Sci. Lett. 284 (1-2), 144-150.

Righter, K., 2011. Prediction of metal-silicate partition coefficients for siderophile elements: an update and assessment of P-T conditions for metal-silicate equilibrium during accretion of the Earth. Earth Planet. Sci. Lett. 304 (1-2), 158-167.

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 (3-4), 239-255.

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 (1-2), 31-42.

Rubie, D.C., Jacobson, S.A., Morbidelli, A., O'Brien, D.P., Young, E.D., de Vries, J., Nimmo, F., Palme, H., Frost, D.J., 2015. 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.

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

Samuel, H., Tackley, P.J., 2008. Dynamics of core formation and equilibration by negative diapirism. Geochem. Geophys. Geosyst. 9 (6).

Shoemaker, E.M., 1962. Interpretation of lunar craters. In: Physics and Astronomy of the Moon, pp. 283-359.

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.

Šránek, O., Ricard, Y., Dubuffet, F., 2010. A multiphase model of core formation. Geophys. J. Int. 181 (1), 198-220.

Stevenson, D.J., 1981. Models of the Earth's core. Science 214 (4521), 611-619.

Ulrová, M., Coltice, N., Ricard, Y., Labrosse, S., Dubuffet, F., Velišský, J., Šránek, O., 2011. Compositional and thermal equilibration of particles, drops, and diapirs in geophysical flows. Geochem. Geophys. Geosyst. 12 (10).

Wacheul, J.-B., Le Bars, M., 2017. Fall and fragmentation of liquid metal in a viscous fluid. Physical Review Fluids 2 (9), 090507.

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

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., 2005. Core formation and the oxidation state of the Earth. Earth Planet. Sci. Lett. 236 (1-2), 78-95.

Wetherill, G.W., 1985. Occurrence of giant impacts during the growth of the terrestrial planets. Science 228 (4701), 877-879.

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

Wood, B.J., Wade, J., Kilburn, M.R., 2008. Core formation and the oxidation state of the Earth: additional constraints from Nb, V and Cr partitioning. Geochim. Cosmochim. Acta 72 (5), 1415-1426.
# Thermal evolution of a metal drop falling in a less dense, more viscous fluid

B. Qaddah<sup>1,2,*</sup> J. Monteux,<sup>2</sup> and M. Le Bars<sup>1</sup>

<sup>1</sup>CNRS, Aix Marseille Université, Centrale Marseille, IRPHE, 13013 Marseille, France

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

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

(Received 20 December 2019; accepted 16 April 2020; published 15 May 2020)

The initial state of terrestrial planets was partly determined, during accretion, by the fall of metal drops in a liquid magma ocean. Here we perform systematic numerical simulations in two-dimensional cylindrical axisymmetric geometry of these falling dynamics and associated heat exchanges at the scale of one single drop for various initial sizes and ambient viscosities. We explore Reynolds number in the range 0.05–48, viscosity ratios in the range 50–4000, Weber number in the range 0.04–5, and Peclet number in the range 70–850. We show that heat exchange between the two phases occurs predominantly at the front section of the drop. Our systematic, parametric study shows that the thermal boundary layer thickness, the depth and time for equilibration, the Nusselt number, and the magma ocean volume affected by thermal exchanges all scale as power laws of the Peclet number. Because of drop distortions, these scaling laws deviate from the classical balances considering only heat diffusion through a laminar thermal boundary layer. Finally, when considering a temperature-dependent viscosity of the ambient fluid, we show that a low-viscosity layer surrounds the drop, which influences the thermal evolution of nondeformable, low-Reynolds-number drops only and decreases the breakup distance for some limited breakup modes.

DOI: 10.1103/PhysRevFluids.5.053801

## I. INTRODUCTION

Core formation of terrestrial planets is a complex process contemporaneous with planetary accretion [1,2]. Its fluid dynamics and thermodynamics have been addressed in numerous studies (e.g., Refs. [3–8]). During the last stages of, e.g., Earth accretion, giant impacts likely occurred between the proto-Earth and up to Mars-sized differentiated bodies [9]. The kinetic energy released during such collisions [10,11], the radioactive heating caused by the disintegration of short-lived radio elements [12], and the heat dissipation resulting from the conversion of potential energy during core formation and core/mantle separation [13], melt part or all of the Earth mantle [11]. Following each impact, the iron core of the impactor thus spread and sank into a deep magma ocean. There, the metal further fragmented into blobs of different sizes, ranging from millimeter drops up to kilometer diapirs before assembling with the Earth protocore ([4,8,14,15]). Thermochemical exchanges occurred between the fragmented metal drops and the liquid magma ocean during their sinking, determining the initial thermal and chemical state of the planet ([4,16–18]). Past studies have provided many scenarios to characterize and quantify the thermochemical exchanges. References [17] and [19] modelled the diffusive equilibration through a laminar thermal boundary layer of, respectively, a cloud of uniform drops and a large diapir of iron. Reference [18] further evaluated the influence of drop deformations in Ref. [17]

<footnote>*Corresponding author: qaddah@irphe.univ-mrs.fr</footnote>
scenario. Reference [20] solved the fully coupled dynamical and thermal/chemical equations, but for a fixed spherical geometry only. Reference [4] used experiments where a large volume of immiscible fluid falls into a less dense ambient to show that the smallest scale of turbulence—rather than diffusion through a laminar boundary layer—leads to rapid thermochemical equilibration, even before fragmentation. Reference [7] confirmed this conclusion in their analog model, measuring the global cooling of a large volume of hot Galinstan after its fall through a deep tank of viscous oil. Yet a systematic temporal description of heat exchanges at the scale of one falling, freely evolving drop, is still missing.

Importantly, the magma ocean viscosity highly depends on its evolving temperature and pressure [21]. Therefore, the viscosity ratio between the magma ocean and iron drops can vary by several orders of magnitude as a function of depth, of time after impact, etc. Following and extending an abundant literature in different contexts (e.g., Refs. [22–26]), analog experiments in Refs. [7,8] and numerical simulations at the scale of one metal drop in Ref. [27] showed that the viscosity contrast indeed plays an important role in iron drops shape, velocity, and fragmentation. Reference [27] predicted that thermochemical exchanges should increase with drop deformation and oscillations; but they did not explicitly solve for the fully coupled dynamical and thermal equations. This is the purpose of the present paper. Open questions include the following: How and where do heat exchanges occur? Do the drop deformation/oscillations indeed favor heat exchanges? What are the characteristics time and depth needed to reach equilibration between the two phases? And what is the influence of a temperature-dependent viscosity of the magma?

The paper is organized as follows. Section II introduces the physical and numerical models, with the governing equations, the nondimensional parameters, and the numerical method. Section III presents in detail a reference case, describing its mechanical and thermal behavior, average temperature evolution, heat transfer at the drop interface, and the magma ocean volume heated during the drop sinking. In Sec. IV, we present the main numerical results from our systematic parametric study and derive generic scaling laws for the above detailed parameters. Section V then focuses on changes induced by a temperature-dependent viscosity in the magma ocean. Conclusions and future works are outlined in Sec. VI.

## II. PHYSICAL AND NUMERICAL MODELS

### A. Governing equations

We consider an initially spherical, liquid metal drop of radius $R$, falling in an initially motionless, less dense and more viscous surrounding fluid (i.e., a magma ocean) under the action of gravity. The initial temperature of the liquid drop and of the magma ocean strongly depends on the growth history of the protoplanet before the impact and on its initial heating caused by short lived elements [13]. Here we consider that the liquid metal drop is hotter than the magma ocean, with uniform initial temperatures in both phases. Both phases behave as Newtonian, incompressible, and immiscible fluids with uniform surface tension and constant density and viscosity within each fluid at first. In Sec. V, we also consider a temperature-dependent magma ocean viscosity. The dynamical and thermal evolution of the falling drop and ambient liquid is governed by the Navier-Stokes and heat transfer equations, describing

(i) The mass conservation:

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

with $\mathbf{u}$ the fluid velocity vector (m s$^{-1}$).

(ii) The momentum conservation:

$$ \rho \left( \frac{\partial \mathbf{u}}{\partial t} + \mathbf{u} \cdot \nabla \mathbf{u} \right) = \nabla \cdot \left\{ -PI + \mu \left[ \nabla \mathbf{u} + (\nabla \mathbf{u})^T \right] \right\} + \rho g + F_{\text{st}}, \qquad (2) $$
with $\rho$ the fluid density (kg m$^{-3}$), $\mu$ the fluid dynamic viscosity (Pa s), $t$ the time (s), $P$ the fluid pressure (Pa), $g$ the gravitational acceleration (m s$^{-2}$), $F_{st}$ the surface tension force (N m$^{-3}$) and $I$ the identity matrix.

(iii) The heat conservation:

$$ \rho C_p \left[ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right] + \nabla \cdot [-k \nabla T] = 0, \quad (3) $$

with $C_p$ the heat capacity at constant pressure (J kg$^{-1}$ K$^{-1}$), $T$ the fluid temperature (K), and $k$ the fluid thermal conductivity (W m$^{-1}$ K$^{-1}$). No heat source is considered in our model. In this work, we ignore the effect of viscous heating. Viscous heating during km-scale metal diapir sinking can significantly increase the temperature of the diapir and of the surrounding material depending on the viscosity contrast between the metal and silicate phases [13,28]. When considering the full mantle/core separation within terrestrial planets, viscous heating should be accounted for in the heat budget [29]. In this study we focus on small droplets (millimeters to centimeters in radius) falling over short distances (up to 200 times their radius). Our results (not shown here) confirm that thermal contribution from viscous heating is negligible at such small scales.

To monitor the interface between the falling drop and the magma ocean, we use the level set method, a Eulerian and implicit method frequently used in multiphase flow problems (e.g., Ref. [30]). The level set function $\phi$ equals to 1 in the metal drop and 0 in the ambient liquid and rapidly changes through the interface, whose position is determined by the isocontour $\phi = 0.5$. The transport and reinitialization of the level set function $\phi$ are governed by:

$$ \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 (4) $$

with $\gamma$ (m/s) and $\epsilon$ (m) the reinitialization parameters. $\gamma$ determines the reinitialization amount: A suitable value for $\gamma$ is the maximum velocity magnitude experienced in the model. $\epsilon$ determines the layer thickness around the interface and is equal to half the size of the characteristic mesh in the region explored by the interface. The density and dynamical viscosity are evaluated using the level set function:

$$ \rho = \rho_m + (\rho_d - \rho_m)\phi, \quad (5) $$

$$ \mu = \mu_m + (\mu_d - \mu_m)\phi, \quad (6) $$

where subscripts “m” and “d” denote the magma ocean and the liquid metal drop, respectively. The surface tension force is determined by:

$$ F_{st} = \nabla \cdot T = \nabla \cdot \{\sigma[I + (-nn^T)]\delta\}, \quad (7) $$

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

$$ n = \frac{\nabla \phi}{|\nabla \phi|}. \quad (8) $$

The level set parameter $\phi$ is also used to approximate the $\delta$ function by a smooth function [31] defined by

$$ \delta = 6 |\phi(1 - \phi)| |\nabla \phi|. \quad (9) $$

Note that in this work, we focus on thermal exchanges only and do not calculate the chemical exchanges between the two phases, which are also extremely interesting from a geological point of view (see, e.g., Ref. [32]). Indeed, the main difficulties in solving for chemical exchanges are (i) the small value of the chemical diffusivities and (ii) the presence of a partition coefficient at the
<table><caption>TABLE I. Symbol definitions and values of the physical and nondimensional parameters used in this study.</caption><thead><tr><th></th><th>Symbol</th><th>Value or range</th></tr></thead><tbody><tr><td>Magma ocean density</td><td>ρ<sub>m</sub></td><td>3500 kg m<sup>-3</sup></td></tr><tr><td>Metal drop density</td><td>ρ<sub>d</sub></td><td>7500 kg m<sup>-3</sup></td></tr><tr><td>Metal drop viscosity</td><td>μ<sub>d</sub></td><td>0.005 Pa s</td></tr><tr><td>Magma ocean viscosity</td><td>μ<sub>m</sub></td><td>0.25–20 Pa s</td></tr><tr><td>Initial drop radius</td><td>R</td><td>4–25 mm</td></tr><tr><td>Surface tension coefficient</td><td>σ</td><td>1 N m<sup>-1</sup></td></tr><tr><td>Magma ocean heat capacity</td><td>Cp<sub>m</sub></td><td>667 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Metal heat capacity</td><td>Cp<sub>d</sub></td><td>800 J kg<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Magma ocean conductivity</td><td>k<sub>m</sub></td><td>10 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Metal conductivity</td><td>k<sub>d</sub></td><td>100 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Viscosity ratio</td><td>R<sub>μ</sub></td><td>50–4000</td></tr><tr><td>Density ratio</td><td>R<sub>ρ</sub></td><td>2.14</td></tr><tr><td>Reynolds number</td><td>Re</td><td>0.05–48</td></tr><tr><td>Weber number</td><td>We</td><td>0.04–5</td></tr><tr><td>Peclet number</td><td>Pe</td><td>70–850</td></tr><tr><td>Nusselt number</td><td>Nu</td><td>1–6</td></tr></tbody></table>

metal/magma interface, meaning that species concentration at the moving interface is discontinuous [32]. Such challenges are at present beyond the scope of our numerical study. We nevertheless argue that the main conclusions shown here for heat exchanges also give some clues for chemical exchanges and hence should be accounted for in geochemical models of planet building [33].

## B. Physical and nondimensional parameters

The main parameters that characterize the dynamical and thermal evolution of a falling drop in a more viscous medium are the viscosity, density, thermal conductivity, heat capacity, and initial temperature of the two fluids; the initial drop size; the gravity; and the surface tension between the two phases. In the geophysical problem of interest (i.e., core formation), the magma ocean viscosity and the metal drop initial radius vary over a wide range of values, while the other parameters are roughly constant (even if rigorously, the thermal conductivity and heat capacity of a magma ocean moderately depend on its composition (e.g., Refs. [34,35]), and the thermal conductivity and heat capacity of metal drops moderately depend on temperature and pressure (e.g., Refs. [36,37]). Hence in this study, we vary these two parameters *R* and *μ*<sub>m</sub>, in the accessible, relevant ranges 4–25 mm and 0.25–20 Pa s, respectively, while we keep all the other parameters fixed at their representative geophysical values (see Table I).

In our simulations, the drop falls from rest, accelerates until reaching a constant terminal velocity, possibly with small oscillations around it, and exchanges heat with the ambient liquid. We continue our simulations until the drop reaches a stable dynamical regime and its temperature contrast with the ambient reaches less than 20% of its initial value. During the fall, we monitor the average velocity *V* and the average temperature *T* of the drop (minus the initial ambient temperature). The dynamical and thermal evolution of each drop is then characterized by the following output dimensionless numbers, where we use the terminal average velocity (note that if it has not yet reached a steady value by the end of the computation, we use estimate from our previous study [27]):

(i) The Reynolds number ($Re = \frac{ρ_m V R}{μ_m}$) is the ratio of inertial to viscous forces. Three different regimes are possible: The Stokes regime corresponds to $Re < 1$ where the viscous effects dominate; the intermediate regime corresponds to $Re = 1–500$ where both viscous and inertial forces are important; and the Newtonian regime corresponds to $Re > 500$ where the inertial forces are
<table><caption>TABLE II. Dimensional and nondimensional parameters for all performed simulations used in this study. Pe<sub>h</sub> is the grid Peclet number (see Sec. II C).</caption><thead><tr><th>Simulation</th><th>R (mm)</th><th>μ<sub>m</sub> (Pa s)</th><th>R<sub>μ</sub></th><th>Re</th><th>We</th><th>Pe</th><th>Pe<sub>h</sub></th><th>ΔT (K)</th></tr></thead><tbody><tr><td>no. 1</td><td>8</td><td>20</td><td>4000</td><td>0.05</td><td>0.04</td><td>70</td><td>0.43</td><td>100</td></tr><tr><td>no. 2</td><td>8</td><td>10</td><td>2000</td><td>0.2</td><td>0.15</td><td>137.26</td><td>0.84</td><td>100</td></tr><tr><td>no. 3</td><td>8</td><td>5</td><td>1000</td><td>0.76</td><td>0.52</td><td>254</td><td>1.55</td><td>100</td></tr><tr><td>no. 3a (T-dependent μ<sub>m</sub>)</td><td>8</td><td>5</td><td>1000</td><td>0.81</td><td>0.59</td><td>271</td><td>1.65</td><td>100</td></tr><tr><td>no. 3b (T-dependent μ<sub>m</sub>)</td><td>8</td><td>5</td><td>1000</td><td>0.98</td><td>0.85</td><td>327</td><td>2</td><td>1000</td></tr><tr><td>no. 4</td><td>8</td><td>1</td><td>200</td><td>9.5</td><td>3.23</td><td>635</td><td>3.89</td><td>100</td></tr><tr><td>no. 5</td><td>8</td><td>0.5</td><td>100</td><td>21.8</td><td>4.26</td><td>728.36</td><td>4.46</td><td>100</td></tr><tr><td>no. 5a (T-dependent μ<sub>m</sub>)</td><td>8</td><td>0.5</td><td>100</td><td>22.4</td><td>4.48</td><td>747.6</td><td>4.57</td><td>100</td></tr><tr><td>no. 5b (T-dependent μ<sub>m</sub>)</td><td>8</td><td>0.5</td><td>100</td><td>23.52</td><td>4.94</td><td>785</td><td>4.8</td><td>1000</td></tr><tr><td>no. 6</td><td>8</td><td>0.25</td><td>50</td><td>47.6</td><td>5</td><td>793.73</td><td>4.86</td><td>100</td></tr><tr><td>no. 6a (T-dependent μ<sub>m</sub>)</td><td>8</td><td>0.25</td><td>50</td><td>47.6</td><td>5</td><td>793.73</td><td>4.86</td><td>100</td></tr><tr><td>no. 6b (T-dependent μ<sub>m</sub>)</td><td>8</td><td>0.25</td><td>50</td><td>47.6</td><td>5</td><td>793.73</td><td>4.86</td><td>1000</td></tr><tr><td>no. 7</td><td>4</td><td>1</td><td>200</td><td>2.2</td><td>0.34</td><td>145.67</td><td>0.89</td><td>100</td></tr><tr><td>no. 8</td><td>6</td><td>1</td><td>200</td><td>5.67</td><td>1.53</td><td>378.19</td><td>2.32</td><td>100</td></tr><tr><td>no. 9</td><td>10</td><td>1</td><td>200</td><td>12.7</td><td>4.61</td><td>847.42</td><td>5.2</td><td>100</td></tr><tr><td>no. 10</td><td>25</td><td>1</td><td>200</td><td>45</td><td>23.2</td><td>3008</td><td>18.45</td><td>100</td></tr><tr><td>no. 10a (T-dependent μ<sub>m</sub>)</td><td>25</td><td>1</td><td>200</td><td>41.13</td><td>19.5</td><td>2745</td><td>16.8</td><td>1000</td></tr></tbody></table>

dominant. Here the Reynolds number ranges from 0.05 to 48, and hence our drops are in the Stokes to intermediate regimes.

(ii) The Weber number ($We = \frac{\rho_m V^2 R}{\sigma}$) compares the inertial and surface tension forces. It governs the deformation, breakup, and terminal shape of a drop (see, e.g., Refs. [27,38]). When $We < O(1)$, the drop remains spherical without any change of its morphology, while increasing Weber number leads to stronger and stronger deformation, then to fragmentation above a threshold which increases with the viscosity ratio, starting from ~3 for viscosity ratio $\le 1$ (see, e.g., Refs. [27,39]). Here the Weber number ranges from 0.04 to 5, considering stable, potentially deformable drops only.

(iii) The Peclet number ($Pe = \frac{\rho_m C_p V R}{k_m}$) compares the rate of heat advection to diffusion at the drop scale. Here the Peclet number ranges from 70 to 850, so heat transfer is strongly affected by advection.

(iv) The Nusselt number [$Nu = \frac{R\overline{T} \cdot n}{(T_{int} - T_m)}$] compares the measured, averaged heat transfer at the drop interface to a purely conductive case, with $\overline{T}_{int}$ the mean temperature at the interface and $T_m$ the magma ocean temperature far from the drop. Here the Nusselt number ranges from 1 to 6, hence confirming the important role of advection in heat transfer.

All relevant parameter values are given in Table I.

### C. Numerical model

We solve Eqs. (1)–(4) using axisymmetric simulations with COMSOL Multiphysics software, based on the finite-element method. The details of our 17 runs for this study are listed in Table II. Each run represents 2 to 4 weeks’ computation time on a biprocessor, eight-core, 3.2- to 3.6-GHz workstation. The axisymmetric geometry assumption is validated in Ref. [40] for a Weber number up to 120. For the dynamics, we use open conditions at the top and bottom boundaries and no-slip conditions at the lateral boundary. For the temperature, we consider no flux conditions at all boundaries. The computational domain must be large enough to allow for convergence without any wall effects. Here we chose an axisymmetric cylinder of size ($r \times z$) = (12R × 200R), which is sufficiently large to reach a statistically steady motion (see our previous study [27]) and to follow equilibration up to a 80% decrease of the initial temperature anomaly.
<figure><img src="image_128.png" /><figcaption>FIG. 1. A zoom illustrating our adaptive mesh (left) and the method for mesh evolution over time when the drop reaches the bottom of the finest mesh region (right).</figcaption></figure>

To capture precisely the dynamical and thermal evolution of the drop, a fine mesh is required. For that, we use an adaptive mesh with a high resolution in the drop vicinity. As shown in Fig. 1, we divide our domain into several regions where the cell sizes vary between $h = 0.015R$ and $h = 1.5R$. As detailed in Ref. [27], the simulation is programed to stop when the drop reaches the bottom of the finest mesh region: The whole mesh pattern is then translated and the simulation is restarted on this new grid [see Fig. 1 (right)].

Mesh convergence has been checked, always using standard Lagrange mesh elements of type P2-P3 (quadratic for the pressure and temperature fields and cubic for the velocity field). In Ref. [27], we showed that a mesh size $h = 0.025R$ (or smaller) allows us to capture the falling drop dynamics. Here we further performed two tests for the resolution of thermal transfers. In the first one, we compare the numerical and analytical heat transfers by thermal diffusion from a motionless spherical drop. The analytical solution for the radial temperature profile from the drop center to a given distance ($r$) is given by [41]

$$T = 0.5T_0 \left[ \text{erf} \frac{R+r}{2\sqrt{Dt}} + \text{erf} \frac{R-r}{2\sqrt{Dt}} \right] - \frac{T_0}{r} \sqrt{\frac{Dt}{\pi}} \left[ e^{-\frac{(R-r)^2}{4Dt}} - e^{-\frac{(R+r)^2}{4Dt}} \right] \quad (10)$$

with $T$ the temperature anomaly ($T_0$ its initial value) and $D$ the thermal diffusivity ($m^2 s^{-1}$). Figure 2 (left) shows the excellent agreement of our numerical results.

For the second test, we calculated, for different minimum grid sizes $h = 0.05R, 0.035R, 0.025R, 0.015R, 0.01R$, the normalized average drop temperature as a function of normalized time for our reference case no. 4 in Table II (see details in next section). Figure 2 (right) shows a reasonable convergence of the numerical results from $h = 0.015R$, with a relative maximum error $\simeq 3.5\%$, while the drop mass during the course of this simulation does not change by more than 0.4% from its initial mass. Therefore, we confirm that our mesh $h = 0.015R$ captures correctly the thermal evolution of the metal drop.

Finally, the grid Peclet number $Pe_h = \frac{Vh}{D_{\text{int}}}$, with $D_{\text{int}} = \frac{D_d+D_m}{2}$ the mean thermal diffusivity, is an appropriate parameter in the convection-diffusion equation to determine whether the heat transfer in the system corresponds to numerical artifacts or not. Reference [42] found very accurate solutions for grid Peclet number up to 10. As shown in Table II, our runs in this study remain below this
<figure><img src="image_129.png" /><figcaption>FIG. 2. Left: Comparison of our numerical results with the analytical results in Ref. [41] for the radial profile of the normalized temperature anomaly for a nondeformable and immobile drop. Right: Comparison of the time evolution of the normalized mean temperature anomaly determined with different mesh sizes for a deformable and mobile drop (simulation no. 4 in Table II).</figcaption></figure>

criterion for stability, except for run no. 10, which we hence do not consider for heat transfer studies in the following sections.

### III. REFERENCE CASE

In this section, we present our reference case, from which, in the next section, we then change the drop size keeping the viscosity ratio constant, and the magma ocean viscosity keeping the drop size constant. This reference case corresponds to simulation no. 4 in Table II: The drop initial radius is 8 mm, the magma ocean viscosity is 1 Pa s, the metal viscosity is 0.005 Pa s (viscosity contrast 200), and the initial temperature difference between the metal drop and the magma ocean is 100 K.

The drop motion from rest and its thermal evolution are shown in Fig. 3. The spherical drop accelerates due to gravity and rapidly deforms into a spherical cap. Then, surface tension equilibrates inertia at the drop interface and prevents any further change in morphology. Higher temperatures are concentrated at the front of the drop, while the thermal wake behind the drop expands as a function of time. To characterize this dynamics, we define and compute the following quantities.

#### A. Average drop velocity

During its fall, the drop mean velocity varies with time depending on the drop morphology. In our simulations, we compute it as:

$$V(t) = \frac{\int_{v} U(r, z, t)[\phi(r, z, t) \ge 0.5]dv}{\int_{v} [\phi(r, z, t) \ge 0.5]dv}, \qquad (11)$$

with $U(r, z, t)$ the local velocity magnitude (m/s), $dv = 2\pi rdrdz$ accounting for axisymmetric cylindrical geometry, and $[\phi(r, z, t) \ge 0.5]$ the Boolean operator allowing to only capture the iron drop volume. Note that we use the mean magnitude velocity of the drop, which is actually very close to its average falling velocity (change of less than 3% in our reference case) because vertical motions largely predominate. We normalize $V(t)$ by the free fall Newton velocity $U_N = \sqrt{\frac{\Delta\rho gR}{\rho_m}}$ and time by the diffusion time $t_{dif} = \frac{R^2 \rho_m C \rho_m}{k_m}$. Figure 4 shows the result for our reference case. The drop rapidly accelerates from rest up to $t^* = 0.02$, then several small oscillations occur before converging toward its asymptotic terminal fall velocity, which we use to compute the output dimensionless numbers.
<figure><img src="image_130.png" /><figcaption>FIG. 3. Dynamical and thermal evolutions of the metal drop as a function of time for simulation no. 4 in Table II. The color presents the temperature anomaly normalized by its initial value $T^*$. The black solid line separates the metallic material from the magma ocean. From left to right, the time normalized by the thermal diffusion time is $t^* = 0, 0.01, 0.02$, and 0.05.</figcaption></figure>

## B. Average temperature evolution of the metal drop

As shown in Fig. 3, the liquid drop exchanges heat with the liquid magma ocean and loses its heat as a function of time and depth. We determine the mean average temperature anomaly of the drop compared to the magma ocean as:

$$T(t) = \frac{\int_{v} T(r, z, t)[\phi \ge 0.5]dv}{\int_{v} [\phi \ge 0.5]dv} \quad (12)$$

<figure><img src="image_131.png" /><figcaption>FIG. 4. Normalized average velocity of the metal drop as a function of normalized time for simulation no. 4 in Table II.</figcaption></figure>
<figure><img src="image_132.png" /><figcaption>FIG. 5. Normalized average temperature anomaly of the metal drop as a function of normalized time (left) and depth (right) for simulation no. 4 in Table II. The black lines represent our numerical results and the red dashed lines with stars represent the exponential fits from Eq. (13). The correlation coefficient for the two fits is 0.9999.</figcaption></figure>

We normalize $T(t)$ by its initial value, giving $T^*(t)$. Its evolution for the reference case as a function of time or as a function of depth of its center of mass (normalized by the initial radius $R$) is shown in Fig. 5. After a rapid adjustment of the initial temperature jump at the interface, $T^*$ decreases exponentially toward equilibrium. We determine the best exponential fits according to

$$T^* = e^{(-t/t_c)} \quad \text{and} \quad T^* = C_z e^{(-z/l_c)}, \qquad (13)$$

where $C_z$ is a constant (see Fig. 5 in red with stars). In geophysical science, the characteristic time ($t_c$) and length ($l_c$) for equilibration are very important parameters [7], used to set the degree of equilibration between iron and silicate of magma ocean in planet building models [33]. Note that chemical equilibrium is more difficult to reach than thermal equilibrium because of the lower value of the chemical diffusivity. The chemical equilibration between the metallic phase and the magma ocean is a key process and constraints on this dynamics contribute to our understanding of the core formation in terms of timing. However, it is also a process that is complicated to compute because of the presence of a partition coefficient at the interface. The results shown below for thermal equilibrium may be considered as first-order trends for chemical equilibrium that will nevertheless deserve a dedicated study.

### C. Volume of heated magma during the drop sinking

The magma ocean temperature increases during the fall of the drop, especially in its wake. It is important in geophysics to quantify how much of the magma is contaminated (here, thermally affected) by the formation of the iron core: This would for instance affect the initial structure and heat budget of the Earth's mantle. In our simulations, we calculate the volume of magma ocean affected by thermal exchange as:

$$V_o = \int_{v} [T > T_c][\phi \le 0.5] dv, \qquad (14)$$

with $T_c$ a chosen temperature anomaly ($K$) and $[\phi \le 0.5]$ the Boolean operator allowing to only capture the magma ocean volume. An example of the captured volume of heated magma ocean is shown in Figure 6 (left). We performed tests with several $T_c$ corresponding to 0.1, 1, and 10 K. Results of the volume normalized by the initial drop volume as a function of the drop depth are shown in Figure 6 (right). The three curves corresponding to different equilibrium degrees are parallel, until a depth equals to 50R where $T_c = 10$ K is strongly affected by thermal diffusion in the magma ocean. The same behavior would of course take place at longer time/depth for the
<figure><img src="image_133.png" /><figcaption>FIG. 6. Left: Volume of heated magma ocean ($T_c = 1$ K) once the drop reaches a time $t^* = 0.067$ and a depth $z = 37R$ in our simulation no. 4 in Table II. Right: Normalized heated magma volume as a function of normalized depth for different values of $T_c$.</figcaption></figure>

other $T_c$. But our purpose here is to look at the relative initial behavior before the diffusion effect dominates, depending on the drop radius and magma viscosity. In the following, for our parametric study, we thus compare heated volumes considering $T_c = 1$ K at three different depths: 10R, 20R, and 30R.

### D. Heat transfer at the drop interface

To evaluate the Nusselt number, we compute the mean temperature anomaly and the mean temperature gradient at the interface, using, respectively,

$$ \overline{T}_{\text{int}} = \frac{\int_{v} T \phi(1-\phi) dv}{\int_{v} \phi(1-\phi) dv}, \qquad (15) $$

$$ \overline{\nabla T}_{\text{int}} = \frac{\int_{v} \nabla T \phi(1-\phi) dv}{\int_{v} \phi(1-\phi) dv}, \qquad (16) $$

where $\phi(1-\phi)$ allows us to capture only the interface region between the two phases. Figure 7 shows the temperature gradient in and around the drop at a given time for our reference simulation. Most heat transfer between the metal drop and the magma ocean occurs at the drop front side, while the back half only accounts for about 15% of the total. The existence of a hot thermal wake and the external recirculation behind the drop that encapsulates and entrains magma with the sinking drop are the main reasons limiting back thermal exchanges [27,43].

A close view of the temperature field and of the thermal boundary layer is shown in Fig. 8. We define the thermal boundary layer thickness as the distance from the drop interface to a point where the temperature anomaly reaches 1% of its interfacial value. The boundary layer is very thin at the drop front section, where most thermal exchanges are carried out. It significantly increases behind the drop. An example of thickness measurement is shown in Fig. 8 (right) at the thinnest position, i.e., at the front of the drop. In order to give a global estimate around the drop, we also define the average boundary layer thickness as

$$ \delta_{T_{\text{av}}} = \frac{\overline{T}_{\text{int}}}{\overline{\nabla T}_{\text{int}}}. \qquad (17) $$
<figure><img src="image_134.png" /><figcaption>FIG. 7. Temperature gradient ($K/m$) in and around the drop for our reference simulation no. 4 in Table II at time $t^* = 0.055$.</figcaption></figure>

## IV. PARAMETRIC STUDY AND SCALING ANALYSIS

We now present our systematic exploration of the parameter space, changing the magma ocean viscosity or the drop initial radius as shown in Table II, from our reference simulation no. 4. We then analyze our results in terms of scaling laws. In particular we determine the influence of the Peclet number on the dimensionless parameters defined in the previous section: the time and length of equilibration, the normalized thermal boundary layer thickness, the Nusselt number, and the dimensionless magma ocean volume affected by thermal exchanges.

### A. Time and length of equilibration

As expected from previous studies of drop dynamics, but largely neglected in geophysical applications (see, e.g., discussions in Refs. [8,27]), the characteristic time and length of equilibrium depend on the viscosity contrast between the metal drop and the magma ocean. As presented in Fig. 9, increasing the viscosity ratio increases the required time for thermal equilibration (Fig. 9,

<figure><img src="image_135.png" /><figcaption>FIG. 8. Left: Temperature anomaly (color scale) and contour of the thermal boundary layer (red line) around the drop (black line) for our reference simulation no. 4 in Table II at time $t^* = 0.067$. Middle: zoom at the drop front. Right: Temperature profile of the thermal boundary layer in front of the drop and determination of the local boundary layer thickness.</figcaption></figure>
<figure><img src="image_136.png" /><figcaption>FIG. 9. Characteristic time (left) and length (right) of equilibration as a function of the viscosity contrast [see Eq. (13)].</figcaption></figure>

left), mostly because a larger ambient viscosity limits advective heat exchanges. We can further observe in Fig. 9 (left) two different slopes corresponding to two different dynamical regimes of the drop: For $R_μ > 1000$, the drop remains spherical and the equilibration characteristic time is highly dependent on the viscosity ratio, while for $R_μ < 200$, the drop deforms which limits the influence of $R_μ$. The evolution of the equilibration length is nonmonotonic (Fig. 9, right), because the magma ocean viscosity also influences the falling velocity. Globally, with a more viscous ambient fluid, the thermal equilibrium between the two phases occurs less deeply in the magma ocean once the stable drop is formed, as we predicted in our previous work [27].

The drop size at a given viscosity ratio also influences the thermochemical equilibration (for a fixed spherical shape, see Ref. [20]). Increasing the initial radius of the liquid metal drop decreases the surface of exchange over volume initial ratio and increases the drop falling velocity: We thus expect an increase of both the characteristic time and length, as confirmed in Fig. 10. In Ref. [27], we predicted that increasing the drop initial radius also leads to surface extension due to drop distortion, hence to faster equilibration compared to a purely spherical drop: This effect is, however, limited, because drop deformation mostly occurs at the back of the drop, while heat exchange takes place mostly at the front.

Those two series of results can be rationalized by considering dimensionless properties as a function of the Peclet number. Here we normalize the equilibrium time by the thermal diffusion time and the equilibrium length by the initial drop radius. Results are shown in Fig. 11. The thermal equilibrium time for a high Peclet, spherical drop theoretically scales as $Pe^{-0.5}$ (see, e.g., Ref. [32]).

<figure><img src="image_137.png" /><figcaption>FIG. 10. Characteristic time (left) and length (right) of equilibration as a function of the drop initial radius [see Eq. (13)].</figcaption></figure>
<figure><img src="image_138.png" /><figcaption>FIG. 11. Dimensionless equilibrium time (left) and length (right) as a function of the Peclet number.</figcaption></figure>

Here we find

$$t_c^* = 2.45\text{Pe}^{-0.59\pm0.01}, \qquad (18)$$

where the shown uncertainty in the exponent comes from considering separately numerical simulations performed while changing the viscosity ratio and those performed while changing the initial sphere radius. This acceleration of the equilibration compared to the theoretical model may be a signature of the drop deformation at large Peclet. Indeed, performing the same fit only on the subset of data points where the drop remains spherical gives an exponent $0.51 \pm 0.01$.

Correspondingly, the length needed to reach the equilibration increases monotonically with Peclet number following

$$l_c^* = 2.37\text{Pe}^{0.41\pm0.01}, \qquad (19)$$

in agreement with a quasiconstant falling velocity, i.e., $l_c = V \times t_c$. In conclusion, when the flow advection measured by the Peclet number increases, thermal equilibration occurs faster but deeper in the magma ocean (after drop formation), which is of importance for Earth's building models [33].

## B. Heat exchanges

Simply equilibrating large-scale heat advection along the drop surface with heat diffusion through the thermal boundary layer perpendicular to it, the theoretical size of the thermal boundary layer normalized by the drop radius classically scales as $Pe^{-0.5}$: This was verified numerically in Ref. [20] for a rigid falling sphere. As shown in Fig. 12 (left), we also recover this scaling at the drop front, where most heat exchange occurs: The best fit gives

$$\delta_T^* = 2.53\text{Pe}^{-0.5\pm0.01}. \qquad (20)$$

With the deformable drops considered here, however, the boundary layer thickness strongly varies around the metal drop: We thus compute a characteristic average value $\delta_{T_{av}}$ from Eq. (17), normalized by the initial radius. Results are shown in Fig. 12 (right), with the best-fit scaling law

$$\delta_{T_{av}}^* = 10.7\text{Pe}^{-0.6\pm0.01}. \qquad (21)$$

Over the explored range, the average thickness is at least three times larger than the front one. The larger-than-expected 0.6 exponent might again be a signature of drop distortion at large Peclet number.

Then, by definition, the Nusselt number averaged over the drop surface should scale as $1/\delta_{T_{av}}^*$. This is indeed recovered, as shown in Fig. 13 with the best-fit scaling law

$$\text{Nu} = 0.08\text{Pe}^{0.63\pm0.02}. \qquad (22)$$

Note that Nu varies over time in many cases in this study: We consider here its time-averaged value.
<figure><img src="image_139.png" /><figcaption>FIG. 12. Dimensionless thickness of the thermal boundary layer at the drop front as a function of Peclet number (left). Dimensionless average thickness of the thermal boundary layer around the drop as a function of Peclet number (right).</figcaption></figure>

Comparing this result with the one obtained in the previous section, we also recover that the typical dimensionless time for equilibration $t_c^*$ scales like 1/Nu, as expected. This proves the self-consistency of our measurements.

### C. Volume of heated magma ocean

In the wake of the drop, the magma ocean is thermally affected by its passage. Because of heat diffusion within the magma ocean, the affected volume widens with time. But on the short times, at a given depth of the drop $z$, it can be estimated by simply balancing the heat that has passed through the drop interface with the heat accumulated in this volume of magma ocean, whose temperature has increased by a given amount larger or equal to $T_c$. We then predict

$$Vo^* = \frac{\text{Volume of heated magma}}{\text{Initial drop volume}} \sim 3 \frac{\text{Nu}}{\text{Pe}} \frac{z}{R} \frac{\Delta T_0}{T_c}, \quad (23)$$

with $\Delta T_0$ the initial temperature difference between the two phases. Results in Fig. 14 show a good agreement with a small relative error for various depths $z = 10R, 20R, 30R$ and temperature

<figure><img src="image_140.png" /><figcaption>FIG. 13. Nusselt number as a function of Peclet number</figcaption></figure>
<figure><img src="image_141.png" /><figcaption>FIG. 14. Dimensionless heated magma ocean volume as a function of $X$, and comparison with our scaling law (23). Our numerical results are represented with different symbols for different given depths.</figcaption></figure>

contrast $T_c = 1K$, using the previously determined scaling for $Nu$. This volume can be compared to the volume crossed by the falling drop, which writes in its dimensionless form

$$V_{\text{crossed}}^* = \frac{\pi R^2 z}{4/3\pi R^3} = \frac{3z}{4R}, \qquad (24)$$

which also scales linearly with depth. Our result does not exhibit any strong increase of the volume effectively contaminated vs. this simple geometrical estimate: this is because, even in our most turbulent cases, the falling drop encapsulates in its wake a small volume of silicate and entrains it toward the deep, while the rest of the ambient fluid remains largely motionless or laminar. No additional mixing is generated. Our conclusion should nevertheless be re-evaluated in the presence of global magma motions, coming, e.g., from turbulent convection, depending on the typical excited time- and length scales. But this is beyond the scope of the present paper. Our conclusion should also be re-evaluated while considering large-scale viscous heating which, for a turbulent wake of a large diapir or a large cloud of drops, might generate strong additional heating long after the iron fall.

## V. INFLUENCE OF A TEMPERATURE-DEPENDENT VISCOSITY

The magma ocean viscosity increases from the surface to the base of the magma ocean because of the combined effects of temperature and pressure [21]. Those large-scale variations are irrelevant in the context of our local study; nevertheless, during the drop fall over the typical length of ~200R considered here, the temperature of the liquid magma ocean increases because of heat exchanges with the drop, therefore its viscosity decreases at a given pressure (i.e., depth) following the equation in Ref. [21] for anhydrous liquid

$$\mu_T = 0.00033e^{[6400/(T-1000)]} \text{Pa s.} \qquad (25)$$

Here we investigate the influence of such a temperature-dependent viscosity on the dynamics, stability, and thermal evolution of a drop. For a given initial viscosity ratio, the initial magma ocean temperature is computed using (25), and we then consider an initial temperature difference with the hot metal drop of $\Delta T = 100$ K or $\Delta T = 1000$ K. We focus on three reference cases, only changing
<figure><img src="image_142.png" /></figure>

FIG. 15. Thermal evolution of the nondeformable metal drop as a function of normalized time (left). The color scale presents the temperature anomaly normalized by its initial value $T^*$, and the black solid line shows the surface of the metal drop. From left to right, $t^* = 0, 0.0167, 0.033, 0.05$, and $0.067$. The right figure shows the temperature gradient (K/m) in and around the spherical drop. In this model, $Re = 0.76$, $Pe = 254$, and $R_\mu = 1000$ (simulation no. 3 in Table II).

the viscosity ratio (see Table II). Those cases are first presented in the absence of temperature-dependent viscosity:

(i) The nondeformable drop case (simulation no. 3 in Table II, $Re = 0.76$ and $Pe = 254$) is presented in Fig. 15. Drop's shape remains mostly constant. The temperature is close to uniform inside the drop because of internal recirculation and progressively decreases in the wake. Most heat transfer takes place at the front.

(ii) The weakly deformable drop case (simulation no. 5 in Table II, $Re = 21.8$ and $Pe = 728.36$) is presented in Fig. 16. The drop quickly deforms toward a spherical cap, associated with a strong release of heat in its wake. The system then reaches a quasi-steady state, with the most surprising feature being a more rapid temperature decrease in the drop than in the wake, leading to a positive heat transfer from the silicate to the iron in the drop's back. Nevertheless, most heat transfer still takes place at the front.

(iii) Finally, the strongly deformable drop case (simulation no. 6 in Table II, $Re = 47.6$ and $Pe = 793.73$) is presented in Fig. 17. Here the drop shape keeps oscillating over the whole depth because of competing surface tension and inertial forces. This induces oscillatory temperature changes in the wake, associated to periodic thermal plumes emitted from the drop sides and strong inhomogeneous fluxes within the magma ocean. Nevertheless, heat transfers between iron and silicate are still largely focused at the drop's front, explaining that this case is not associated with any specific signature in the previous section.

We now rerun those three cases using the previously defined temperature-dependent viscosity.

## A. Drop dynamics

Figure 18 shows the normalized average velocity of each drop as a function of the normalized time for constant magma ocean viscosity in black lines and for temperature-dependent viscosity with $\Delta T = 100$ K ($\Delta T = 1000$ K) in red (blue) dashed lines. The normalized average velocity for the nondeformable drop (Fig. 18, left) with $\Delta T = 100$ K increases moderately by about 4.8% compared to the constant viscosity case, keeping the spherical shape without change in the morphology. For $\Delta T = 1000$ K, the drop deforms from the back side and its velocity strongly increases by about 26%; concomitantly, the average viscosity around the drop decreases by 93%, explaining this significant dynamical change. Such significant changes of the drag in the Stokes regime for an underformable drop were already addressed in Refs. [44,45].
<figure><img src="image_143.png" /><figcaption>FIG. 16. Thermal evolution of the weakly deformable metal drop as a function of normalized time (left). The color scale presents the temperature anomaly normalized by its initial value $T^*$, and the black solid line shows the surface of the metal drop. From left to right, $t^* = 0.015, 0.024, 0.033, 0.047, \text{ and } 0.06$. The right figure shows the temperature gradient (K/m) in and around the deformable drop. In this model, $Re = 21.8$, $Pe = 728.36$, and $R_\mu = 100$ (simulation no. 5 in Table II).</figcaption></figure>

<figure><img src="image_144.png" /><figcaption>FIG. 17. Thermal evolution of the strongly deformable metal drop as a function of normalized time (left). The color scale presents the temperature anomaly normalized by its initial value $T^*$, and the black solid line shows the surface of the metal drop. From left to right, $t^* = 0.02, 0.033, 0.047, 0.06$, and $0.077$. The right figure shows the temperature gradient (K/m) in and around the deformable drop. In this model, $Re = 47.6$, $Pe = 793.73$, and $R_\mu = 50$ (simulation no. 6 in Table II).</figcaption></figure>
<figure><img src="image_145.png" /></figure>

FIG. 18. Comparison of the normalized average velocity field as a function of normalized time to show the influence of a temperature-dependent viscosity for a nondeformable spherical drop (left), a weakly deformable drop (middle) and a strongly deformable drop (right). The black lines present the constant magma ocean viscosity. The red (blue) dashed lines present the temperature-dependent viscosity with $\Delta T = 100$ K ($\Delta T = 1000$ K).

On the contrary, for the weakly deformable drop (Fig. 18, middle) and the strongly deformable drop (Fig. 18, right), no significant dynamical change is observed, despite a large viscosity decrease in the associated thermal boundary layer (by 76% and 92%, respectively, for $\Delta T = 1000$ K; see also Fig. 19). This indicates that the drop dynamics is already mostly inviscid in the absence of temperature-dependent viscosity, as can be guessed from the corresponding values of the Reynolds number (see caption of Figs. 16 and 17).

## B. Drop fragmentation

In the most extreme previous case (the strongly deformable drop with $\Delta T = 1000$ K), Fig. 19 indicates that the hot silicate layer surrounding the drop has a viscosity value of the order 0.05 Pa s: Assuming such a viscosity uniformly in the ambient magma ocean, our previous study [27] predicts a rapid fragmentation, which is not observed here. Actually, the viscosity constrast at the interface does not influence the drop stability, because the thermal boundary layer where the low viscosity is localized is very thin in front of the drop, compared to any relevant dynamical length scale; besides, the wake has no influence on the drop. Temperature-dependent viscosity only affects the drop's

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

FIG. 19. View of the viscosity field (in Pa s) around the strongly deformable drop when a temperature-dependent viscosity is considered.
<figure><img src="image_147.png" /><figcaption>FIG. 20. Final form at breaking of the drop for a constant viscosity magma ocean (left) and a temperature-dependent viscosity magma ocean with $\Delta T = 1000$ K (middle). On the right, view of the viscosity field (in Pa.s) around the drop in the later case. In this model, $Re = 50$, $Pe = 3330$, and $R_\mu = 200$ (simulation no. 10 in Table II).</figcaption></figure>

dynamics and fragmentation when flows and temperature change on comparable scales in upfront or sides locations.

To further prove this, we performed a simulation with a drop radius of 25 mm and a viscosity ratio of 200 for constant magma ocean viscosity first and then for $\Delta T = 1000$ K-dependent magma ocean viscosity (simulation no. 10 in Table II). Figure 20 (left) shows the jellyfish fragmentation mode when constant magma ocean viscosity is considered (more details can be found in Ref. [27]). The normalized time ($t_{bk}^* = \frac{t_{bk}V}{2R}\sqrt{\frac{\rho_d}{\rho_m}}$) and distance ($d_{bk}^* = \frac{d_{bk}}{R}$) of breakup equal 3.8 and 10.7, respectively. For temperature-dependent viscosity (Fig. 20, middle), another jellyfish fragmentation mode is observed, with time and distance of breakup equal to 2.3 and 5.6, respectively. Actually, here the extended jellyfish membranes form filaments, where thermal effects are very important. These membranes sink into a less-viscous medium whose thickness is comparable to the membranes thickness (see Fig. 20, right). Thus, the filaments move freely in this zone with low viscous constraints and finally separate from the drop volume, leading to fragmentation.

Note finally that the local dimensionless numbers, taken at the mean temperature (and associated viscosity) of the thermal boundary layer, become $Re = 1000$, $We = 19$, and $R_\mu = 10$. From Figs. 4, 12, and 13 of our previous study [27], we find consistent results for the fragmentation mode of the drop, as well as for its time and distance of breakup, accounting for these dimensionless numbers (to compare with our previous study, Re and We should be multiplied in this study by a factor 2 because in Ref. [27], Re and We were calculated using the drop diameter).

Hence, temperature-dependent magma ocean viscosity influences the drop fragmentation in the filament forming regimes only. In such regimes, the results in Ref. [27] as a function of Re, We, and $R_\mu$ already allow predicting fragmentation mode, breaking time and distance by considering local values of the dimensionless numbers.

### C. Thermal evolution

In Fig. 21, we compare the normalized mean iron temperature for each case as a function of normalized time. For a nondeformable drop (Fig. 21, left), the average temperature for $\Delta T = 100$ K (red dashed lines) and $\Delta T = 1000$ K (blue dashed lines) temperature-dependent viscosity decreases by about 6.1% and 20.8% with respect to the constant viscosity case (black line). As a consequence, it significantly changes the characteristic time and length of equilibration (see Table III). For a weakly deformable drop (Fig. 21, middle), the average temperature evolution with $\Delta T = 1000$ K changes only by about 13.6%; no difference is measured for a strongly deformable drop (Fig. 21, right). As for the fall velocity, temperature-dependent viscosity influences significantly the thermal
<figure><img src="image_148.png" /><figcaption>FIG. 21. Comparison of the dimensionless temperature as a function of normalized time for a nondeformable drop (left), a weakly deformable drop (middle) and a strongly deformable drop (right). The black lines present the constant magma ocean viscosity. The red (blue) dashed lines present the temperature-dependent viscosity with $\Delta T = 100$ K ($\Delta T = 1000$ K).</figcaption></figure>

evolution of nondeformable drops only (i.e., the drops in the Stokes regime). Note, finally, that in the intermediate regime, the characteristics time and length of equilibration decrease slightly with a temperature-dependent viscosity compared to the constant magma ocean viscosity case, whereas in the Stokes regime, the equilibrium characteristic length increases with a temperature-dependent viscosity, in contrast to the characteristic time of equilibration (see Table III).

Regarding the volume of heated magma ocean shown in Fig. 22, no significant change is observed for all cases with $\Delta T = 100$ K (red dashed lines). The volume of heated magma ocean increases by about 1.5 to 2 times compared to the constant viscosity case (black line) for $\Delta T = 1000$ K (blue dashed lines) due to the initial larger temperature difference; but the curves remain parallel, hence exhibiting no specific dynamical signature.

Finally, we have measured the heat transfer between the two phases for all drops with a temperature-dependent viscosity. As shown in Table IV, in agreement with our previous observations, we do not notice any significant change compared to our scaling law for constant viscosity, except maybe for the less deformable drop: Indeed, heat transfer is mostly due to diffusive exchanges through a thin thermal boundary layer at the drop front, which are not affected by temperature-dependent viscosity. The same conclusion was reached experimentally in Ref. [45] in his study of diapiric magma transport through the lithosphere, hence in the limit of small Reynolds number, nondeformable drop, and using an exponential form of the temperature-dependent viscosity law. Associated numerical simulations in Ref. [44] in the limit of infinite Prandtl number exhibited a very weak dependence on the total viscosity variation. Additional systematic studies, both experimental and numerical, are necessary to explore the parameter range relevant to

<table><caption>TABLE III. The change in characteristics time and length for cases also studied with a temperature-dependent viscosity.</caption><thead><tr><th>Simulation</th><th>ΔT</th><th>t<sub>c</sub><sup>*</sup></th><th>I<sub>c</sub><sup>*</sup></th></tr></thead><tbody><tr><td>no. 3 (constant magma ocean viscosity)</td><td>100</td><td>0.1</td><td>24.25</td></tr><tr><td>no. 3a (T-dependent μ<sub>m</sub>)</td><td>100</td><td>0.0975</td><td>25.6</td></tr><tr><td>no. 3b (T-dependent μ<sub>m</sub>)</td><td>1000</td><td>0.082</td><td>27</td></tr><tr><td>no. 5 (constant magma ocean viscosity)</td><td>100</td><td>0.05</td><td>35.14</td></tr><tr><td>no. 5a (T-dependent μ<sub>m</sub>)</td><td>100</td><td>0.0496</td><td>33.33</td></tr><tr><td>no. 5b (T-dependent μ<sub>m</sub>)</td><td>1000</td><td>0.045</td><td>31.25</td></tr><tr><td>no. 6 (constant magma ocean viscosity)</td><td>100</td><td>0.045</td><td>32</td></tr><tr><td>no. 6a (T-dependent μ<sub>m</sub>)</td><td>100</td><td>0.044</td><td>29.4</td></tr><tr><td>no. 6b (T-dependent μ<sub>m</sub>)</td><td>1000</td><td>0.042</td><td>27.7</td></tr></tbody></table>
<figure><img src="image_149.png" /><figcaption>FIG. 22. Comparison of normalized heated magma ocean volume as a function of normalized time for a nondeformable drop (left), a weakly deformable drop (middle) and a strongly deformable drop (right). The black lines present the constant magma ocean viscosity. The red (blue) dashed lines present the temperature-dependent viscosity with $\Delta T = 100$ K ($\Delta T = 1000$ K).</figcaption></figure>

our configuration. In the meantime, we suggest that using our scaling law [Eq. (22)] is the most acceptable solution.

## VI. CONCLUSIONS AND FUTURE WORK

We have carried out series of numerical simulations to characterize the thermal exchanges between a falling drop and a viscous ambient fluid, exploring a new parameter range relevant for the geophysical application of a hot liquid iron drop falling in a magma ocean. We have shown that because of drop distortions, thermal equilibration properties slightly change from the theoretical predictions based on diffusive heat exchanges through a laminar thermal boundary layer. We have also tested that accounting for a temperature-dependent viscosity in the magma ocean barely influences the obtained results, except for limited cases like nondeformable, low-Reynolds-number drops or jellyfish fragmentation mode. Our most relevant results for geophysical application are the scaling laws for the normalized length of equilibration and for the Nusselt number, which both increase monotonically with the Peclet number as $Pe^{0.41\pm0.01}$ and $Pe^{0.63\pm0.02}$, respectively.

Future work should now consider (i) the presence of convection in the ambient magma ocean, which could affect our conclusions providing the associated velocity at the drop scale is at least of similar order as the falling drop velocity, hence requiring extremely turbulent regimes, and (ii) chemical exchanges between iron and silicates, which will determine the initial chemical state of the considered planet (see, e.g., Ref. [46]). Chemical and thermal constitutive equations being similar, similar equilibration scaling laws and dynamics are nevertheless expected.

<table><caption>TABLE IV. Comparison of the Nusselt numbers for cases of temperature-dependent viscosity between our numerical results and our scaling law for constant viscosity Eq. (22).</caption><thead><tr><th>Simulation</th><th>Nu: measured</th><th>Nu: Eq. (22)</th><th>Relative error (%)</th></tr></thead><tbody><tr><td>no. 3a</td><td>2.65</td><td>2.72</td><td>2.6</td></tr><tr><td>no. 3b</td><td>2.95</td><td>3.07</td><td>4.1</td></tr><tr><td>no. 5a</td><td>5.15</td><td>5.17</td><td>0.4</td></tr><tr><td>no. 5b</td><td>5.21</td><td>5.33</td><td>2.3</td></tr><tr><td>no. 6a</td><td>5.28</td><td>5.36</td><td>1.5</td></tr><tr><td>no. 6b</td><td>5.35</td><td>5.36</td><td>0.2</td></tr></tbody></table>
ACKNOWLEDGMENTS

The authors acknowledge funding by the European Research Council 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. 406. We are grateful to the Mésocentre Clermont Auvergne University for providing help as well as computing and storage resources.

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

[2] M. Touboul, T. Kleine, B. Bourdon, H. Palme, and R. Wieler, Late formation and prolonged differentiation of the Moon inferred from W isotopes in lunar metals, *Nature* **450**, 1206 (2007).

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

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

[5] M. Landeau, R. Deguen, and P. Olson, Experiments on the fragmentation of a buoyant liquid volume in another liquid, *Fluid Mech.* **749**, 478 (2014).

[6] D. J. Stevenson, Fluid dynamics of core formation, *Origin of the Earth* (Oxford University Press, 1990), pp. 231–249.

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

[8] J. B. Wacheul, M. LeBars, J. Monteux, and J. M. Aurnou, Laboratory experiments on the breakup of liquid metal diapirs, *Earth Planet. Sci. Lett.* **403**, 236 (2014).

[9] W. K. Hartmann and D. R. Davis, Satellite-sized planetesimals and lunar origin, *Icarus* **24**, 504 (1975).

[10] J. Monteux, N. Coltice, F. Dubuffet, and Y. Ricard, Thermo-mechanical adjustment after impacts during planetary growth, *Geophys. Res. Lett.* **34**, 24201 (2007).

[11] W. B. Tonks and H. J. Melosh, Core formation by giant impacts, *Icarus* **100**, 326 (1992).

[12] M. J. Walter and R. G. Tronnes, Early earth differentiation, *Earth Planet. Sci. Lett.* **225**, 253 (2004).

[13] J. Monteux, Y. Ricard, N. Coltice, F. Dubuffet, and M. Ulvrova, A model of metal-silicate separation on growing planets, *Earth Planet. Sci. Lett.* **287**, 353 (2009).

[14] H. Genda, R. Brasser, and S. J. Mojzsis, The terrestrial late veneer from core disruption of a lunar-sized impactor, *Earth Planet. Sci. Lett.* **480**, 25 (2017).

[15] J. D. Kendall and H. Melosh, Differentiated planetesimal impacts in a terrestrial magma ocean: Fate of the iron core, *Earth Planet. Sci. Lett.* **448**, 24 (2016).

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

[17] D. C. Rubie, H. J. Melosh, J. E. Reid, C. Liebske, and K. Righter, Mechanisms of metal-silicate equilibration in the terrestrial magma ocean, *Earth Planet. Sci. Lett.* **205**, 239 (2003).

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

[19] T. W. Dahl and J. Stevenson, Turbulent mixing of metal and silicate during planet accretion and interpretation of the hf-w chronometer, *Earth Planet. Sci. Lett.* **295**, 177 (2010).

[20] M. Ulvrová, N. Coltice, Y. Ricard, S. Labrosse, F. Dubuffet, J. Velimský, and M. Šrámek, Compositional and thermal equilibration of particles, drops, and diapirs in geophysical flows, *Geochem. Geophys. Geosyst* **12**, Q10014 (2011).

[21] B. B. Karki and L. P. Stixrude, Viscosity of mgsio3 liquid at earth's mantle conditions: Implications for an early magma ocean, *Science* **328**, 740 (2010).
[22] T. Bonometti and J. Magnaudet, Transition from spherical cap to toroidal bubbles, Phys. Fluids **18**, 052102 (2006).

[23] M. Ohta and M. Sussman, The buoyancy-driven motion of a single skirted bubble or drop rising through a viscous liquid, Phys. Fluids **24**, 112101 (2012).

[24] M. Ohta, Y. Akama, Y. Yoshida, and M. Sussman, Three-dimensional simulations of vortex ring formation from falling drops in an immiscible viscous liquid, Chem. Eng. Jpn. **42**, 648 (2009).

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

[26] M. Ohta, Y. Akama, Y. Yoshida, and M. Sussman, Influence of the viscosity ratio on drop dynamics and breakup for a drop rising in an immiscible low-viscosity liquid, Fluid Mech. **752**, 38 (2014).

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

[28] H. Samuel, P. J. Tackley, and M. Evonuk, Heat partitioning in terrestrial planets during core formation by negative diapirism, Earth Planet. Sci. Lett. **290**, 13 (2010).

[29] F. M. Flasar and F. Birch, Energetics of core formation: A correction, J. Geophys. Res. **78**, 6101 (1973).

[30] X-Y. Luo, M.-J. Ni, A. Ying, and M. Abdou, Application of the level set method for multi-phase flow computation in fusion engineering, Fus. Eng. Des. **81**, 1521 (2006).

[31] J. Hu, R. Jia, K. T. Wan, and X. Xiong, Simulation of droplet impingement on a solid surface by the level set method, in *Proceedings of the COMSOL Conference in Boston* (COMSOL, Boston, 2014).

[32] V. Lherm and R. Deguen, Small-scale metal/silicate equilibration during core formation: The influence of stretching enhanced diffusion on mixing, J. Geophys. Res.: Solid Earth **123**, 10,496 (2018).

[33] D. C. Rubie, F. Nimmo, and H. J. Melosh, *Formation of the Earth's Core, Treatise on Geophysics – Evolution of the Earth* (2nd ed.) (Elsevier, Amsterdam, 2015), pp. 43–79.

[34] A. M. Hofmeister, Mantle values of thermal conductivity and the geotherm from phonon lifetimes, Science **283**, 1699 (1999).

[35] C. W. Thomas and P. D. Asimow, Direct shock compression experiments on premolten forsterite and progress toward a consistent high-pressure equation of state for cao-mgo-al2o3-sio2-feo liquids, J. Geophys. Res.: Solid Earth **118**, 5738 (2013).

[36] N. de Koker, G. Steinle-Neumann, and V. Vlček, 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. USA **109**, 4070 (2012).

[37] H. Ichikawa, T. Tsuchiya, and Y. Tange, The p-v-t equation of state and thermodynamic properties of liquid iron, J. Geophys. Res.: Solid Earth **119**, 240 (2014).

[38] M. Pilch and C. A. Erdman, 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, Multiphase Flow **13**, 741 (1987).

[39] E. Villermaux and B. Bossa, Single-drop fragmentation determines size distribution of rain drops, Nat. Phys. **5**, 697 (2009).

[40] S. S. Jain, N. Tyagi, R. S. Prakash, R. V. Ravikrishna, and G. Tomar, Secondary breakup of drops at moderate Weber numbers: Effect of Density ratio and Reynolds number, Int. J. Multiphase Flow **117**, 25 (2019).

[41] J. Crank, *The Mathematics of Diffusion* (2nd ed.) (Clarendon Press, Oxford, 1975).

[42] R. C. Mittal and R. K. Jain, Redefined cubic b-splines collocation method for solving convection-diffusion equations, Appl. Math. Model. **36**, 5555 (2012).

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

[44] S. F. Daly and A. Raefsky, On the penetration of a hot diapir through a strongly temperature-dependent viscosity medium, Geophys. J. Int. **83**, 657 (1985).

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

[46] M. Bouhifd and A. Jephcoat, 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 (2003).
## 4 Impacts géants et évolution primitive des dynamos planétaires

Le champ magnétique terrestre est la conséquence de mouvements rapides du fer dans la partie liquide du noyau terrestre (voir ex. Jones (2015)). L'état du noyau (liquide ou solide) est gouverné par l'efficacité du transport de chaleur dans le manteau et par l'énergie thermique accumulée dans le noyau au cours de son histoire. Comprendre l'état thermique du noyau terrestre et son évolution au cours du temps permet donc aussi de mieux comprendre l'évolution du champ magnétique des objets du système solaire et ainsi comprendre pourquoi certaines planètes en ont un et d'autres non (Buffett, 2002).

Parmi les planètes telluriques et les lunes du système solaire, Ganymède, Mercure et la Terre possèdent actuellement un champ magnétique. Les données paléo-magnétiques les plus anciennes retrouvées à la surface de la Terre étant datées à 3.45 Ga, l'initiation de la géodynamo a potentiellement eu lieu peu de temps après la formation du noyau (Tarduno et al., 2010). Vénus possède un noyau liquide mais ne montre pas de signe de dynamo (Langlais et al., 2010). Mars et probablement la Lune ont enregistré dans leur croûte des preuves d'un champ magnétique interne transitoire (Dyal et al., 1974; Acuña et al., 1999). Les dynamos lunaire et martienne auraient été initiées tôt dans l'histoire de ces objets pour s'arrêter ensuite entre 4.1 et 3.75 Ga (Acuña et al., 1999; Lillis et al., 2008). Bien qu'il y ait aujourd'hui une recrudescence des preuves de ces dynamos transitoires, leur chronologie et les mécanismes de leur initiation et de leur arrêt restent des énigmes majeures.

Pour savoir si une dynamo planétaire peut être générée dans le noyau, il est nécessaire de connaitre l'état thermique interne de la planète et notamment la différence de température entre le noyau et le manteau environnant (i.e. le flux de chaleur). En effet, la génération d'un champ magnétique requiert trois ingrédients : un fluide métallique liquide avec une grande conductivité électrique, une source d'énergie apportée par la convection thermo-chimique et une dynamique favorable au sein du noyau (principalement de l'hélicité induite par de la rotation) (Moffatt, 1978; Olson, 2007). Pour des corps où le noyau est plus petit (ex. la Lune), d'autres mécanismes que la convection doivent être invoqués (Dwyer et al., 2011; Le Bars et al., 2011). Généralement, si nous considérons que le noyau est composé essentiellement de fer et que la planète tourne suffisamment vite sur elle-même, la deuxième condition peut être caractérisée en suivant l'évolution du flux de chaleur à la frontière manteau/noyau (CMB). Plus la convection au sein du noyau sera vigoureuse et complexe, plus la probabilité d'avoir une dynamo sera grande (Kutzner and Christensen, 2002; Christensen and Aubert, 2006). Une convection vigoureuse se traduira aussi par un flux de chaleur important à la CMB.

De nombreux modèles se sont intéressés à l'influence d'impacts géants à la fin de l'accélération des planètes sur les dynamos primitives (Roberts et al., 2009; Arkani-Hamed and Olson, 2010; Monteux et al., 2011). En effet, ces impacts géants déposent non seulement une quantité importante de chaleur dans le manteau susceptible de réduire considérablement l'efficacité de refroidissement du noyau et ainsi d'arrêter la dynamo, mais ils apportent aussi potentiellement
une quantité importante de fer via le noyau de l'impacteur qui peut se réchauffer pendant sa chute et modifier le bilan thermique du noyau.

Les cadres temporels des activités magnétiques des satellites et planètes du Système Solaire, et notamment le fait que ces activités aient pu exister durant les premiers millions d'années de leur histoire, suggèrent que ces dynamos sont peut-être liées aux étapes d'accrétion et de différenciation primitives lors de la formation des planètes (Stevenson, 2003). Grâce aux modèles développés au cours de mon parcours, j'ai proposé de nouveaux mécanismes d'initiation de dynamos planétaires, j'ai relié ces mécanismes aux processus de différenciation et j'ai apporté des contraintes quant aux temps, tailles et rhéologies nécessaires pour la génération de tels champs magnétiques.

## 4.1 Initiation des dynamos au sein des proto-planètes indifférenciées

Il est difficile de contraindre comment le champ magnétique d'une planète se met en marche. Dans la section suivante, je propose de détailler les conditions d'initiation de ce champ magnétique selon deux modalités : (1) au cours de la chute d'un important volume de fer et (2) une fois cet important volume de fer au centre de la planète

### 4.1.1 Au cours de la chute d'un diapir de fer

Qu'ils aient pour origine des phénomènes de séparation locale suite à un impact ou l'assemblage de deux noyaux après impact entre deux protoplanètes différenciées ou d'autres processus de différenciation (percolation, gouttelettes de fer au sein d'un océan magmatique, instabilités Rayleigh-Taylor ...), les diapirs de fer ont fortement influencé la formation du noyau des planètes telluriques et ce pendant les premiers millions d'années de leur histoire (Stevenson, 1981). La variété des processus de différenciation implique une large gamme de tailles caractéristiques de ces diapirs de fer qui peut varier entre quelques centimètres et plusieurs centaines de kilomètres. La vitesse de chute de la phase métallique dépend essentiellement de la taille des diapirs et de la rhéologie du manteau dans lequel ils plongent. La large gamme possible de ces deux paramètres (Karato and Murthy, 1997) implique donc des vitesses de chute variées qui peuvent atteindre dans les cas les plus favorables plusieurs kilomètres par seconde.

Lors de la chute du diapir de fer fondu, les interactions avec le manteau environnant plus visqueux vont engendrer un champ de vitesse axisymétrique au sein du diapir. J'ai déterminé si l'écoulement au sein d'un diapir de fer pouvait générer un champ magnétique lors de sa chute vers le centre d'une planète en rotation (Monteux et al., 2012). Les résultats des modèles numériques que j'ai développés en collaboration avec H. Amit (CR LPGNantes), N. Schaeffer et P. Cardin (ISTerre Grenoble) montrent que l'écoulement produit dans le diapir peut générer un champ magnétique seulement s'il plonge dans une planète en rotation et si les forces inertielles dans le diapir sont suffisamment importantes (Monteux et al., 2012). Les effets de la
rotation et des forces inertielles produisent de l'hélicité dans l'écoulement ce qui est nécessaire pour générer une dynamo (Moffatt, 1978). L'initiation d'une dynamo est plus facile pour un diapir qui plonge à partir de l'équateur que pour un diapir qui plonge au pôle en suivant l'axe de rotation de la planète. Nos modèles ont permis d'apporter des contraintes sur les tailles de diapir et les viscosités du manteau nécessaires pour initier de telles dynamos mais aussi sur l'importance de la rotation et de la turbulence dans le diapir (Figure 13). Nous avons aussi caractérisé l'intensité et la morphologie de ce champ magnétique.

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

FIGURE 13 – Caractéristiques des diapirs pouvant engendrer un champ magnétique pendant leur chute en fonction de leur taille $R$ et de la viscosité du manteau dans lequel ils plongent $\eta_s$. Les symboles verts représentent les dynamos cinématiques produites numériquement (les cercles pour des diapirs plongeant le long de l'axe de rotation de la planète et les carrés pour les diapirs plongeant à l'équateur). La ligne bleue représente la plus petite valeur du nombre de Reynolds magnétique pour laquelle une dynamo est obtenue dans nos modèles. La ligne en trait tiretés rouge $Ro = 1$ représente la limite du nombre de Rossby étudiée dans nos modèles (régime rotationnel). Pour des valeurs de $Ro$ plus grandes (non étudiées), une dynamo pourrait être générée. Pour de gros diapirs ($R > 10$ km) et de faibles viscosités du manteau ($\eta_s < 10^{10}$ Pa.s), le nombre de Reynolds du manteau $Re_s > 10$ et le diapir est susceptible de se dissocier (Monteux et al., 2012).

Nos modèles apportent une explication alternative (et/ou complémentaire) aux dynamos classiques lunaire et martienne engendrées par de la convection dans le noyau (Stegman et al., 2003; Stanley et al., 2008; Amit et al., 2011) qui étaient considérées jusqu'ici comme le seul modèle possible de génération de champ magnétique. Le mécanisme que nous proposons
sons s'ajoute aux modèles de dynamos par forçage mécanique provoquée par des mouvements différentiels entre le manteau et le noyau lunaire (Le Bars et al., 2011; Dwyer et al., 2011). Si nous plaçons nos résultats dans un contexte à n-diapirs lors de la différenciation d'une planète, nos modèles constituent une explication complémentaire des champs magnétiques fossiles observés à la surface de Mars et de la Lune (Hood et al., 2001; Langlais et al., 2004).

### 4.1.2 Une fois au centre de la protoplanète

Les mouvements convectifs et la génération d'une dynamo dans un noyau planétaire métallique dépendent de son état thermique et de la différence de température avec le manteau environnant. L'état thermique de la phase métallique séparée de la phase silicatée après impact météoritique une fois qu'elle a atteint le centre d'une planète en croissance dépend fortement du volume de diapir de fer et des contrastes de viscosité entre les phases métalliques et silicatées. Comme discuté précédemment (voir section 3), pour des rayons de diapir de fer suffisamment grands (supérieurs à 70 km), la dissipation thermique résultant de la déformation visqueuse du manteau environnant engendre une augmentation de la température du diapir métallique par rapport à sa température juste après l'impact (Monteux et al., 2009). Un protonoyau initialement plus chaud que le manteau environnant peut expliquer la présence d'une dynamo primitive au sein des planètes telluriques (Williams and Nimmo, 2004). Savoir si la différence de température entre le protonoyau et le manteau est suffisante pour générer une dynamo est fonction de la rhéologie du manteau et du volume du diapir métallique.

J'ai utilisé le modèle numérique développé durant ma thèse (Monteux et al., 2009) pour identifier les tailles caractéristiques de planètes, de diapir métallique et les contrastes de viscosité nécessaires à une dynamo primitive générée au sein d'un protonoyau différencié après un impact météoritique. Mes résultats montrent notamment que ce sont les conditions d'accrétion qui contrôlent l'existence du champ magnétique primitif (Monteux et al., 2011). En effet, j'ai montré que si une planète peut atteindre une taille significative sans subir de différenciation majeure, un impact géant à la fin de son accrétion peut différencier un protonoyau suffisamment chaud pour initier une dynamo interne pendant plusieurs millions d'années (Figure 14). Ainsi, si un corps avec un rayon d'environ 500 km impacte une planète indifférenciée de la taille de Mars ($R = 3400$ km), une dynamo sera initiée (Figure 14). La durée de cette dynamo évoluera entre quelques milliers d'années et plusieurs millions d'années selon les modes de transfert thermique à la frontière protonoyau/manteau.

## 4.2 Influence des impacts géants sur les dynamos planétaires

En plus d'être des hypothèses plausibles pour expliquer de nombreuses caractéristiques du système solaire comme le système Terre-Lune (Canup, 2004), les dichotomies magnétique et
<figure><img src="image_151.png" /></figure>

FIGURE 14 – Flux de chaleur à la frontière du protonoyau en fonction de la taille du proto-
noyau (symboles noirs) (Monteux et al., 2011). Ce protonoyau est formé par un épisode de
différenciation suite à un impact géant. La forme des symboles est liée à la taille de la planète
impactée. La ligne en traits tiretés représente la valeur seuil à dépasser pour générer une dy-
namo dans un cas où le protonoyau et le manteau ont la même viscosité. La ligne pointillée
correspond au même seuil mais en considérant que le protonoyau est beaucoup moins visqueux
que le manteau. Si le flux de chaleur est supérieur à la valeur seuil pendant plus de 50 Ma, un
deuxième symbole blanc représente la valeur du flux de chaleur correspondant.

topographique martiennes (Marinova et al., 2008; Nimmo et al., 2008) et le faible volume du
manteau de Mercure par rapport à son noyau (Benz et al., 1988), les impacts météoritiques
géants sont de plus en plus évoqués pour expliquer l'initiation ou l'arrêt de la dynamo des
planètes telluriques (Roberts et al., 2009; Arkani-Hamed and Olson, 2010; Reese et al., 2010).
Dans les modèles proposés, les impacts font intervenir des corps avec des rayons typiquement
compris entre 100 et 1000 km. Au moment de ces impacts, le bolide impacteur et la cible
impactée étaient probablement différenciés. Si la vitesse d'impact est faible (de l'ordre de
la vitesse d'échappement) l'efficacité de l'accrétion est proche de 100 % et tout le matériel
de l'impacteur se combine à la planète impactée (Asphaug, 2010). Ainsi, de grands volumes
de fer ont potentiellement pu être apportés au sein de la planète impactée et influencer la
génération d'une dynamo. Au vu des différences de viscosité et ainsi de temps caractéristique
de réajustement thermique qu'il existe entre le manteau et le noyau, le chauffage par impact
géant a pu jouer un rôle à plus long terme sur les dynamos planétaires que le chauffage du
4.2.1 Impacts géants et assemblages de noyaux

Lors d'un impact géant, la quasi totalité du noyau de l'impacteur s'est assemblée avec le noyau déjà existant (Canup, 2004). Une partie, voire la totalité du plongement de la phase mé-tallique, a pu se faire par diapirisme d'importants volumes de fer. Les études montrent que les impacts géants et l'assemblage des deux noyaux ont pu avoir des conséquences importantes d'un point de vue thermique (Arkani-Hamed and Olson, 2010), d'un point de vue chimique (Dahl and Stevenson, 2010) et du point de vue de l'histoire magnétique de l'objet impacté (Reese et al., 2010). Ces études se sont intéressées aux conséquences thermo-mécaniques d'un impact géant sur une planète différenciée. Cependant, aucune d'entre elles ne considère cor-rectement l'assemblage entre le noyau de l'impacteur et le noyau du corps impacté. De tels phénomènes ont pu avoir des conséquences importantes sur l'activité magnétique de la planète impactée.

J'ai développé avec M. Jellinek, C. Johnson (UBC, Vancouver) (Monteux et al., 2013) et avec J. Arkani-Hamed (University of Toronto) (Monteux and Arkani-Hamed, 2014)) des modèles numériques afin de caractériser la dynamique d'assemblage des deux noyaux et no-tamment l'accommodation par déformation du noyau préexistant ainsi que la modification du flux de chaleur au niveau de la CMB (Figure 15). Après un impact géant, une quantité im-portante de l'énergie cinétique de l'impact est convertie en chaleur au sein du manteau de la planète impactée. En augmentant la température dans le manteau, cet effet a tendance à dimi-nuer le flux de chaleur à la CMB mais aussi à créer une stratification thermique sous la CMB qui arrête la dynamo (Arkani-Hamed and Olson, 2010). Cependant, lors de la chute du noyau de l'impacteur, le chauffage visqueux convertissant l'énergie potentielle du diapir en chaleur peut potentiellement augmenter la température du noyau plongeant. Ce deuxième effet aug-menterait ainsi la température dans le noyau de la planète et modifierait la capacité du noyau à se refroidir par convection et ainsi son activité magnétique. Si la planète impactée ne possédait pas de dynamo avant l'impact géant, l'assemblage des noyaux pourrait en initier une.

Nos résultats montrent que les paramètres tels que le volume du noyau de l'impacteur re-tenu au sein de la planète impactée ainsi que la profondeur à laquelle est enseveli ce matériel gouvernent le temps caractéristique d'assemblage des deux noyaux mais aussi l'intensité de la dissipation visqueuse dans le noyau plongeant. Nos résultats montrent aussi qu'augmenter significativement le contraste de viscosité entre le noyau de l'impacteur et le manteau peut réduire la durée de l'assemblage des noyaux à moins de 1000 ans mais n'affecte pas la tempé-rature d'assemblage (Monteux and Arkani-Hamed, 2014). De plus, quelque soit l'efficacité du mélange thermique à l'intérieur du noyau martien, l'assemblage des noyaux ne fait que retarder le redémarrage de la dynamo préalablement arrêtée par la stratification thermique post-impact (Monteux et al., 2013; Monteux and Arkani-Hamed, 2014). Le temps nécessaire pour redé-marrer la dynamo après impact géant et assemblage des noyaux est de l'ordre de 150 – 200 Ma (Monteux and Arkani-Hamed, 2014).
<figure><img src="image_152.png" /></figure>

FIGURE 15 – Température normalisée (gauche) et composition (droite) (pour le fer $f = 1$, pour les silicates $f = 0$) aux temps $t = 0$ (a), $t = 0.6$ Ma (b), $t = 6$ Ma (c) and $t = 38.8$ Ma (d) (résultats obtenus avec une viscosité uniforme, $R_{planete} = 3400$ km et $R_{impacteur} = 600$ km) (Monteux et al., 2013).

### 4.2.2 Influence de l'état thermique post-impact du manteau sur la dynamo

Les propriétés thermiques du manteau jouent un rôle majeur sur la capacité du noyau sous-jacent à se refroidir et donc sur sa propension à générer un champ magnétique. Ainsi, la présence d'anomalies thermiques à la base du manteau peuvent engendrer des variations locales du flux de chaleur à la frontière manteau noyau. Si localement le manteau inférieur est chaud, le flux de chaleur à la CMB sera faible alors que des zones mantelliques profondes plus froides peuvent favoriser des flux de chaleur importants. Ces variations locales de flux de chaleur plus ou moins étendues peuvent avoir des conséquences importantes sur les mouvements convectifs à l'intérieur du noyau et sur la structure du champ magnétique généré.

Le contrôle du manteau sur les dynamos planétaires est souvent étudié en imposant des flux de chaleur hétérogènes à la limite noyau-manteau sur la frontière extérieure des simulations numériques de dynamo. Ces modèles entrent généralement dans deux catégories principales : Soit ils sont proportionnels aux modèles de tomographie sismique du manteau inférieur de la Terre pour simuler des conditions réalistes, soit ils sont représentés par des harmoniques sphériques simples pour une compréhension physique fondamentale. Cependant, en réalité,
<figure><img src="image_153.png" /></figure>

FIGURE 16 – Representation schématique illustrant la complexité dynamique possible dans le manteau inférieur de Mars et pouvant motiver des modèles alternatifs de flux de chaleur à la CMB (À partir de Amit et al. (2015)).

la dynamique du manteau inférieur est beaucoup plus compliquée et ces modèles de flux de chaleur à la CMB peuvent être un peu trop simplistes. Pour Mars, des modèles alternatifs ont été développés en considérant un chauffage localisé du manteau induit par un panache mantellique (Sreenivasan and Mark Jellinek, 2012) ou un impact géant (Kuang et al., 2014). En effet, lors d'un impact géant, le manteau se réchauffe significativement ce qui entraîne une modification du flux de chaleur à la frontière manteau/noyau (Figure 16). Ceci a pour effet de modifier les modalités de refroidissement du noyau et donc sa capacité à générer un champ magnétique (Buffett, 2002).

J'ai modélisé l'effet d'un impact géant sur le champ magnétique de Mars en considérant une anomalie thermique-post impact dans le manteau martien. J'ai déterminé les variations latérales de flux de chaleur ainsi que la réduction du flux de chaleur moyen à la CMB en fonction de la position géographique et de la taille de l'impact. J'ai ainsi pu imposer un flux de chaleur hétérogène à la frontière manteau/noyau comme condition aux limites dans un code de dynamo numérique (code MAGIC) et suivre l'évolution du champ magnétique martien après impact géant (Monteux et al., 2015). Un impact polaire (i.e. sur l'axe de rotation) entraîne une dichotomie magnétique Nord/Sud qui est plus forte que la dichotomie Est/Ouest générée par un impact à l'équateur. L'amplitude de la dichotomie est contrôlée par la vigueur de la convection elle-même fonction des variations latérales du flux de chaleur à la CMB. Nos résultats montrent qu'une hétérogénéité de flux de chaleur générée par un impact géant est plus efficace qu'une anomalie synthétique de degré 1 du flux de chaleur pour générer des dichotomies magnétiques hémisphériques. Nos résultats montrent qu'un impact entre Mars et un objet avec un rayon de ~ 1000 km a même pu enregistrer la dichotomie magnétique observée dans le champ crustal
martien si le refroidissement de l'océan magmatique post-impact a été suffisamment rapide
(avec des temps inférieurs à 100 kyr).

Références personnelles

Monteux, J., Jellinek, A. M., Johnson, C. L., 2011, Why might planets and moons have early dynamos ?, Earth Planet. Sci. Lett., 310, 349–359. LIEN HAL

Monteux, J., Schaeffer, N., Amit, H., Cardin, P., 2012, Can a sinking metallic diapir gene-
rate a dynamo ?, J. Geophys. Res. (Planets), 117 (16), 10005. LIEN HAL

Monteux, J., Jellinek, A. M., Johnson, C. L., 2013. Dynamics of core merging after a mega-impact with applications to Mars' early dynamo ? Icarus 226, 20–32. LIEN HAL

Monteux, J., J., Arkani-Hamed, 2014. Consequences of giant impacts in early Mars : Core merging and Martian dynamo evolution, J. Geophys. Res. (Planets), Volume 119, Issue 3, pages 480–505. LIEN HAL

Monteux, J., Amit, H., Choblet, G., Langlais, B., Tobie, G., 2015, Giant impacts, heterogeneous mantle heating and a past hemispheric dynamo on Mars, Phys. Earth Planet. Int., 240, 114-124. LIEN HAL

Amit, H., Choblet, G., Olson, P., **Monteux, J.**, Deschamps, F., Langlais, B., Tobie, G., 2015 Towards more realistic core-mantle boundary heat flux patterns : a source of diversity in planetary dynamos, Prog. Earth. Planet. Sci., 2-26. LIEN HAL
<figure><img src="image_154.png" /></figure>

ELSEVIER

Volume 301, Issues 1–2

3 January 2011

ISSN 0012-821X

EARTH
& PLANETARY
SCIENCE LETTERS

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

Editors

Rick Carlson, Washington, DC

Peter B. deMenocal, Palauisdes, NY

T. Mark Harrison, Los Angeles, CA

Yanick Ricard, Lyon

P. Shearer, La Jolla, CA

Tilman Spohr, Berlin

Lars Stixrud, London

(This is a sample cover image for this issue. The actual cover is not yet available at this time.)

This article appeared in a journal published by Elsevier. The attached copy is furnished to the author for internal non-commercial research and education use, including for instruction at the authors institution and sharing with colleagues.

Other uses, including reproduction and distribution, or selling or licensing copies, or posting to personal, institutional or third party websites are prohibited.

In most cases authors are permitted to post their version of the article (e.g. in Word or Tex form) to their personal website or institutional repository. Authors requiring further information regarding Elsevier's archiving and manuscript policies are encouraged to visit:

http://www.elsevier.com/copyright
# Why might planets and moons have early dynamos?

J. Monteux<sup>a,*</sup>, A.M. Jellinek<sup>a</sup>, C.L. Johnson<sup>a,b</sup>

## ARTICLE INFO

**Article history:**

Received 17 September 2010

Received in revised form 5 August 2011

Accepted 9 August 2011

Available online xxxx

Editor: T. Spohn

**Keywords:**

early dynamos

growing planets

meteoritic impacts

numerical modeling

## ABSTRACT

Remanent crustal magnetization of martian and lunar crustal rocks plausibly records magnetic fields related to core dynamos active during the first few hundred Myr of evolution of these bodies. Such early fields suggest that core dynamos may occur during or as a result of accretion. We investigate whether the processes governing the segregation and sinking of metallic Fe after a large impact can yield thermal conditions in the core that favor dynamo action on growing planets. Depending on the sizes of the impactor and planet, as well as the temperature-dependence of the viscosity, we identify conditions in which an early transient core dynamo is possible. We also consider the effect of a molten layer surrounding the protocore on the duration of this early dynamo. Our results suggest that dynamos can be initiated in bodies with a radius of 3500 km radius or greater under Earth-like conditions for ohmic dissipation in the core, and in smaller bodies if a less restrictive critical magnetic Reynolds number condition is applied. These dynamos may persist for several kyr to several Myr depending on the heat transfer regime at the protocore-mantle boundary.

© 2011 Elsevier B.V. All rights reserved.

## 1. Introduction

Among the terrestrial planets and moons of the inner Solar System, Mars and possibly the Moon record evidence for early and short-lived (a few hundred Myr) internally generated magnetic fields. Mars Global Surveyor and Lunar Prospector magnetometer data (Acuña et al., 1999; Hood et al., 2003) confirm the presence of ancient magnetized crust on these bodies. Although later dynamo fields (beginning around or after 3.9 Ga) have been suggested for both bodies (Cisowski et al., 1983; Schubert et al., 2000), most evidence points to earlier dynamos, that ceased by ~3.9 Ga (Garrick-Bethell et al., 2009; Johnson and Phillips, 2005; Lawrence et al., 2008; Lillis et al., 2008). Moreover, recent analyses of paleomagnetism in angrite meteorites suggest that an early and internally generated dynamo occurred within the angrite parent body (Weiss et al., 2008). Although there is increasing evidence for the occurrence of these transient dynamos, their origin and timing remain a major enigma. Here we investigate a plausible temporal link between early dynamos and the large impacts that characterize the late accretion histories of terrestrial bodies in our solar system. Assuming chondritic bulk compositions, these impacts can lead to extensive mantle melting, iron-silicate separation, iron segregation and ultimately core formation (Monteux et al., 2009; Ricard et al., 2009; Srámek et al., 2010; Tonks and Melosh, 1992). If dynamo action depends primarily on the rate of core heat loss (Buffett, 2002), then a crucial issue for explaining

early dynamos is to understand the thermal states of both the impacted planet and the nascent core: under what conditions will the resulting core cooling be sufficiently large to favor dynamo action?

Large meteorite impacts can govern the early stages of planetary formation. Fig. 1 shows a schematic temporal evolution of our impact-induced dynamo model from time $t_1$ (impact time) to time $t_5$ (dynamo is initiated). A collision between an impactor and an undifferentiated protoplanet leads to an increase in mass of the impacted body, the formation of an impact basin, the propagation of shock waves (Fig. 1, $t_1$) and a localized temperature increase, $\Delta T_0$, in a spherical region below the impact site (Fig. 1, $t_2$). All these processes are very fast compared to the following thermo-chemical re-equilibration. The temperature increase is constant within a spherical volume $V_{ic}$ called isobaric core with radius $R_{ic}$ and rapidly decreases away from it (Monteux et al., 2007; Pierazzo et al., 1997; Senshu et al., 2002). Assuming that meteorites impact at the escape velocity of the impacted body, this temperature increase scales with the size of the impacted planet and is superimposed on the initial mean temperature, $T_0$, of the planet before impact. For sufficiently large planets, the temperature increase can overcome the melting temperature of the metallic phase leading to metal-silicate separation (Monteux et al., 2009) (Fig. 1, $t_3$). The dense metallic phase sinks toward the center of the impacted body. How this occurs is debated. Proposed mechanisms include percolation through a solid mantle (Shannon and Agee, 1996; Stevenson, 1990), settling of cm-sized droplets in a magma ocean [e.g. (Rubie et al., 2003)] and to large mantle diapirs [e.g. (Ricard et al., 2009; Senshu et al., 2002; Tonks and Melosh, 1992)]. In addition, it has been proposed that metallic diapirs can sink toward the center of the planet via fracturation and diking (Solomatov, 2000; Stevenson, 2003b; Tonks and Melosh, 1992) or

<footnote><sup>a</sup> Department of Earth and Ocean Sciences, University of British Columbia, Canada</footnote>

<footnote><sup>b</sup> Planetary Science Institute, Tucson, AZ, United States</footnote>

<footnote>* Corresponding author.
E-mail address: jmonteux@eos.ubc.ca (J. Monteux).</footnote>
<figure><img src="image_156.png" /><figcaption>**Fig. 1.** Schematic view of the chemical equilibration following a large impact on an undifferentiated protoplanet. Within the isobaric region that results from the dissipation of the shock wave (time step t₁), the temperature increase (orange circle, t₂) melts the metal (red) that segregates rapidly from silicates (green) (t₃). The initial position of the inertial center of the iron diapir is indicated with a black dot (i.c.). Then the metallic phase sinks by a diapiric instability (t₄). Subsequent convection within the protocore can initiate an early magnetic field (t₅). The top left panel shows where viscous heating (H<sub>v</sub>) occurs for a spherical density anomaly (i.e. a metallic protocore) outlined by the white dashed line within a Rybczynski–Hadamard analytical velocity field (see text for details). H<sub>v</sub> ~ ε² where ε is the dimensionless strain rate tensor. The viscous heating is shown for 2 viscosity ratios, λ, between the hot protocore and the cold undifferentiated material. For λ = 1 (uniform viscosity), viscous heating can occur within the sphere, while for λ ≪ 1, viscous heating only occurs outside the sphere. Here, H<sub>v</sub> is normalized by the maximum value of H<sub>v</sub> when λ = 1. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

viscous deformation (Honda et al., 1993; Monteux et al., 2009) depending on the thermal state of the planet. Here, we consider a simplified scenario in which a single large impact into an initially undifferentiated planet leads to a single Fe-diapir that sinks through an effectively viscous mantle at a Stokes velocity (Fig. 1, t₄).

Whether convection and dynamo generation can occur in a liquid iron protocore depends on the temperature difference between the protocore and the surrounding mantle. A protocore sufficiently hotter than the surrounding silicate mantle is one way to explain the presence of an early dynamo on terrestrial planets and moons. Previous studies have examined the following possible causes of an early hot protocore and dynamo: radioactive heating (Nimmo, 2009; Williams and Nimmo, 2004; Yoshino et al., 2003), strong cooling by mantle convection (Breuer and Spohn, 2003; Elkins-Tanton et al., 2005; Nimmo and Stevenson, 2000; Stegman et al., 2003) or impact heating (Monteux et al., 2009; Reese and Solomatov, 2010). Here we focus on whether the thermal conditions resulting from a single large impact and the subsequent differentiation event are also favorable to dynamo initiation. In this case, the temperature difference between the protocore and the overlying mantle depends on the thermal evolution of the metallic diapir as it sinks through a colder mantle toward the center of the planet (Fig. 1, t₅). Because the heat content and buoyancy of the diapir depend on its volume, one critical parameter will be the radius of the metallic diapir R<sub>Fe</sub>. An additional potentially important contribution to the initial core–mantle temperature difference will be the extent to which the gravitational potential energy released on sinking is converted into heat in the metallic diapir

and in the surrounding mantle via viscous dissipation. The importance of this effect, and where it occurs (i.e., primarily in the diapir or surrounding mantle), will depend critically on the strongly temperature-dependent viscosity η(T) of the rocks composing the protoplanetary mantle. For impacted planets of given radii R we use the numerical model of Monteux et al. (2009) to identify R<sub>Fe</sub> – R – η(T) conditions favoring dynamo action. In addition, we investigate the dependence of the dynamo duration on the protocore size, and on conditions that lead to melting of the overlying mantle material.

## 2. Thermo-chemical model

### 2.1. Initial state

The initial structure and thermal state of a growing planet depend particularly on the characteristics of its accretion (Kaula, 1979; Safronov, 1978). Indeed, differentiated terrestrial planets probably formed rapidly from accumulation of chondritic material (Agee, 1997). In more detail, during accretion, heating driven by a combination of the dissipation of impact energy and the decay of short lived radionuclides such as <sup>26</sup>Al and/or <sup>60</sup>Fe (Monteux et al., 2007; Senshu et al., 2002; Yoshino et al., 2003) increases the mean internal temperature and gives rise to a radial temperature gradient that depends on the accretion rate relative to the rate of radiative cooling to space. Where this growth rate is very high in comparison to surface cooling, this heating can ultimately cause partial or complete melting (Yoshino et al., 2003) and extensive internal differentiation. In contrast, where the growth rate is very low in comparison to surface cooling, little additional heating occurs (Wood et al., 2006) and the growing planet remains an undifferentiated mechanical mixture of the accreted impactors. Hence, depending on its growth rate, the frequency of impacts and the size of the impactors, a protoplanet can potentially grow to a Mars size radius in ≈ 1 Myr without undergoing any significant episodes of melting and any subsequent global differentiation processes (Senshu et al., 2002). Short-lived isotopic constraints restrict the timing of core formation on terrestrial planets to be within the first tens of Myr following accretion (Kleine et al., 2002; Touboul et al., 2007; Yin et al., 2002). From these models, the martian core formation occurred faster (within 30 Myr) than the Earth's core formation (within 50 Myr). Consequently, between 1 Myr (the growing planet reaches a Mars size radius with no or little differentiation) and 50 Myr (the core is fully formed), large impacts can lead to episodic large differentiation events. Depending on the thermal states of the resulting protocore and surrounding mantle, subsequent cooling may drive a core dynamo that generates a magnetic field observable at the surface.

Whether an impact onto an undifferentiated planet leads to melting and core formation depends critically on the fraction γ of the initial kinetic energy E<sub>kin</sub> of the impactor that is dissipated within the planet itself as a result of work done by shock waves (Monteux et al., 2009; Tonks and Melosh, 1992). E<sub>kin</sub> is proportional to the impact velocity v<sub>imp</sub><sup>2</sup>, which can be decomposed into two contributions:

$$v_{\text{imp}}^2 = v_{\text{esc}}^2 + v_{\infty}^2. \qquad (1)$$

Here, v<sub>esc</sub> is the escape velocity of the impacted planet and v<sub>∞</sub> is the velocity of the impactor at a distance much greater than that over which the gravitational attraction of the impacted planet is important. A conservative and well-constrained estimate of E<sub>kin</sub> is possible if we take v<sub>∞</sub> = 0 and assume v<sub>imp</sub> = v<sub>esc</sub>. We note that in the case of an impact on a moon orbiting around a larger central body, the gravitational attraction of the central body can increase this impact velocity and we address this issue in Section 5. Laboratory experiments and modeling studies suggest that the energy released within the planet is mainly dissipated by frictional heating and melting within a spherical region (an “isobaric core”) with a volume V<sub>ic</sub> that is at most 3 times larger than the volume of the impactor V<sub>imp</sub> (Croft, 1982; O’Keefe and Ahrens, 1977; Pierazzo et al.,
1997; Senshu et al., 2002). Within this region the temperature increases uniformly from $T_0$ to $T_0 + \Delta T_0$ where $T_0$ is the temperature of the impacted body before the impact and $\Delta T_0$ is the temperature increase due to impact heating. Experiments suggest also that there is an additional thin spherical shell of damaged and heated material outside the isobaric core. Assuming that the work done by shock waves in this shell causes no melting, the excess temperature $\Delta T_0$ decreases rapidly and smoothly with distance $r$ across the shell as approximately $\Delta T_0(R_{ic}/r)^m$. Following Senshu et al. (2002), and fitting the decay of peak pressure with distance away from the edge of the isobaric core $m \approx 4.4$ (Monteux et al., 2007). Hence, the amount of energy $E_{ext}$ dissipated as heat outside the isobaric core is:

$$E_{ext} = 2\pi \overline{\rho_0 C_p} \int_0^\pi \int_{R_{ic}}^{+\infty} \Delta T_0 \left(\frac{R_{ic}}{r}\right)^m r^2 \sin \theta d\theta dr \quad (2)$$

with $\theta$ as the zenith angle from the impact site and $\overline{\rho_0 C_p}$ is the average specific heat of the impacted body. For the axisymmetric geometry of the isobaric core, the previous equation can be reduced to:

$$E_{ext} = \overline{\rho_0 C_p} V_{ic} \Delta T_0 \left( \frac{3(2m-5)}{2(m-3)(m-2)} \right). \quad (3)$$

Assuming that $v_{imp} = v_{esc}$ (Kokubo and Ida, 1996) and that $v_{ic} = 3v_{imp}$ (Senshu et al., 2002), a balance between the kinetic energy delivered to heat the growing planet ($\gamma E_{kin}$) and the energy used to heat up and melt the isobaric core and to heat up the surrounding material without melting is

$$\frac{4\pi}{9} \gamma G \rho_0^2 V_{ic} R^2 = \overline{\rho_0 C_p} \Delta T_0 V_{ic} + E_{ext} + V_{ic} (\rho_{Fe} f_0 L_{Fe} + \rho_{Si} (1-f_0) L_{Si}). \quad (4)$$

After some algebra we obtain the temperature increase at the impact site as a function of the impacted planet radius (Monteux et al., 2007)

$$\Delta T_0 = \frac{1}{h_m \overline{\rho_0 C_p}} \left[ \frac{4\pi}{9} \gamma \rho_0^2 G R^2 - (\rho_{Fe} f_0 L_{Fe} + \rho_{Si} (1-f_0) L_{Si}) \right]. \quad (5)$$

Here, G is the gravitational constant, $\rho_0$ is the density of the undifferentiated planetary material, R is the radius of the impacted planet. $L_{Fe}$ and $L_{Si}$ are the latent heat required for melting respectively the metallic and the silicate phases. $f_0$ is the volume fraction of metal in the impacted planet. From Eqs. (3), (4), and (5), $h_m$ represents the amount of heat that is used to increase the temperature inside and outside the isobaric core relative to the amount of heat used to increase the temperature by $\Delta T_0$ within the core and (Monteux et al., 2007)

$$h_m = 1 + \frac{3(2m-5)}{2(m-3)(m-2)} \approx 2.7. \quad (6)$$

Because $\Delta T_0$ depends on $\gamma$, the value for $\gamma$ requires some discussion. O'Keefe and Ahrens (1977) show that $\gamma$ depends on the impact velocity and that the fraction of energy converted into heat increases with impact velocity. Our assumption that $v_{imp} = v_{esc}$ implies that our modeled impactors are in a “slow impactor regime”. From the measurements of O'Keefe and Ahrens (1977) a reasonable value for $\gamma$ is 0.3. We apply this value in our calculations and discuss the influence of the uncertainty of this choice in Section 5. An additional issue that requires mention is the assumption that there is no melting in the spherical shell around the isobaric core. For $\gamma = 0.3$, calculation of $\Delta T_0$ as a function of planet radius shows that this approximation is strictly appropriate for planetary radii less than about 4000–4500 km. For larger planets, latent heat may be consumed outside the isobaric core and so $\Delta T_0$ calculated with Eq. 5 may be overestimated. The magnitude of the overestimate is, however, unclear because the effect of melting on the presence and structure of the spherical shell region

bounding the isobaric core is unknown. Nevertheless, assuming that the geometry of the region is similar to the case without melting, calculations show that the uncertainty in the value for $\gamma$ will have a greater influence on our estimates for $\Delta T_0$.

The magnitude of the temperature increase in Eq. 5 depends on the two terms in square brackets: the left-hand term describes the kinetic energy released to heat the mantle on impact and the right-hand term is the latent heat consumed to melt the silicate and iron components of the mantle. For iron–silicate segregation to occur in the heated volume, $\Delta T_0$ must be positive and a sketch of the thermal state following impact is shown in Fig. 1 ($t_2$). Thus, equating the two terms in parentheses in Eq. 5 gives a critical impacted planet radius above which segregation is possible:

$$R_{crit} > \left[ \frac{\rho_{Fe} f_0 L_{Fe} + \rho_{Si} (1-f_0) L_{Si}}{4\pi / 9\gamma \rho_0^2 G} \right]^{1/2}. \quad (7)$$

For parameter values listed in Table 1, $R_{crit} = 1620$ km. Hence, impact heating on a Moon to a Mars-sized body will result in the separation of the dense iron phase from the silicate phase as a metallic diapir (Fig. 1). We note that whereas our assumption that $v_{imp} = v_{esc}$ is appropriate for most planets, for moons this simplification underestimates the impact velocity, which will be governed by the gravitational potential of the central body and leads to a lower bound on $\Delta T_0$.

This temperature increase is superimposed on the pre-impact thermal regime of the protoplanet. The background thermal state depends on the accretionary process, which can lead to a temperature profile that increases with radius toward the surface (Senshu et al., 2002), and on the strong early rate of radiogenic heating, which augments the mean internal temperature of the planet (Yoshino et al., 2003). For simplicity we assume that $T_0$ is uniform and close to the eutectic temperature of the Fe–FeS system at 1 bar (Fei et al., 1997) and that the impacted body before the impact is a homogenous mixture of silicate and metal.

For a given initial bulk composition, the radius of the diapir separated after a large impact $R_{Fe}$ is a function of the impactor size $R_{imp}$ and $f_0$ in the impacted planet and scales as $R_{Fe} = (3f_0)^{1/3} R_{imp}$

<caption>Table 1</caption>

Typical parameter values for numerical models.

<table><tr><td>Planet radius</td><td>R</td><td>2000–4000 km</td></tr><tr><td>Diapir (protocol) radius</td><td>R<sub>Fe</sub></td><td>150–500 km</td></tr><tr><td>Density difference</td><td>Δρ<sub>0</sub></td><td>4500 kg m<sup>-3</sup></td></tr><tr><td>Average density</td><td>ρ<sub>0</sub></td><td>4270 kg m<sup>-3</sup></td></tr><tr><td>Iron density</td><td>ρ<sub>Fe</sub></td><td>8000 kg m<sup>-3</sup></td></tr><tr><td>Silicate density</td><td>ρ<sub>Si</sub></td><td>3500 kg m<sup>-3</sup></td></tr><tr><td>Mean thermal expansion</td><td>α</td><td>4.4 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Iron thermal expansion</td><td>α<sub>Fe</sub></td><td>1.5 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Silicate thermal expansion</td><td>α<sub>Si</sub></td><td>5 × 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Iron heat capacity</td><td>C<sub>p,Fe</sub></td><td>800 J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>Silicate heat capacity</td><td>C<sub>p,Si</sub></td><td>1000 J K<sup>-1</sup> kg<sup>-1</sup></td></tr><tr><td>Average specific heat of the impacted body</td><td>ρ<sub>0</sub>C<sub>p</sub></td><td>4 × 10<sup>6</sup> J K<sup>-1</sup> m<sup>-3</sup></td></tr><tr><td>Initial temperature</td><td>T<sub>0</sub></td><td>1250 K</td></tr><tr><td>Iron melting temperature</td><td>T<sub>melt,Fe</sub></td><td>1250 K</td></tr><tr><td>Silicate melting temperature</td><td>T<sub>melt,Si</sub></td><td>1650 K</td></tr><tr><td>Iron latent heat</td><td>L<sub>Fe</sub></td><td>4 × 10<sup>5</sup> J kg<sup>-1</sup></td></tr><tr><td>Silicate latent heat</td><td>L<sub>Si</sub></td><td>2.7 × 10<sup>5</sup> J kg<sup>-1</sup></td></tr><tr><td>Thermal diffusivity</td><td>κ</td><td>10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Mean thermal conductivity</td><td>k = k<sub>Si</sub></td><td>4 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Iron thermal conductivity</td><td>k<sub>Fe</sub></td><td>40 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Metal content</td><td>f<sub>0</sub></td><td>0.17</td></tr><tr><td>Reference viscosity</td><td>η<sub>0</sub></td><td>10<sup>22</sup> Pa s</td></tr><tr><td>Impact energy conversion coefficient</td><td>γ</td><td>0.3–0.5</td></tr><tr><td>Volume effectively heated by impact</td><td>h<sub>m</sub></td><td>2.7</td></tr><tr><td>Gravitational constant</td><td>G</td><td>6.67 × 10<sup>-11</sup> m<sup>3</sup> kg<sup>-1</sup> s<sup>-2</sup></td></tr><tr><td>Average magnetic field strength</td><td>B</td><td>2.5 mT</td></tr><tr><td>Magnetic diffusivity</td><td>ν</td><td>2 m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Magnetic permeability</td><td>μ</td><td>4π × 10<sup>-7</sup> H m<sup>-1</sup></td></tr></table>
(Monteux et al., 2009). The metal diapir sinks toward the center of the planet over time-scales from several kyr to several Myr depending on the rheology of the protoplanetary mantle. Once at the center of the undifferentiated planet, this “metallic diapir” will be referred to as a “protocore”.

## 2.2. Heat partitioning during the sinking of the diapir

The initial thermal state of the protocore is governed by its volume and the rate and location of viscous dissipation during its sinking, which depends on the viscosity ratio $\lambda$ (lower than 1) between the hot protocore and the relatively cold undifferentiated surrounding mantle (Ke and Solomatov, 2009). For a given $\lambda$ the rate of viscous dissipation depends on the shear strain rate $\dot{\epsilon}^2$, where $\dot{\epsilon} \sim V_S / R_{Fe}$ and $V_S$ is the sinking velocity. Monteux et al. (2009) showed that for small diapirs (i.e., $R_{Fe} < 70$ km), the rate of viscous dissipation is negligible and the metallic phase only cools during its sinking. In contrast, for $R_{Fe} > 70$ km and for mantle and diapir viscosities that are equal (uniform viscosity), viscous deformation within the diapir and the surrounding mantle leads to dissipation that can increase the temperature of the metal diapir over its initial value from $T_0$ to $T_{core} > T_0 + \Delta T_0$.

Fig. 1 (top left panel) shows where the viscous heating occurs when a spherical metallic diapir sinks within a Rybczynski-Hadamard analytical velocity field (Hadamard, 1911; Rybczynski, 1911) for 2 different viscosity ratios $\lambda$. For a uniform viscosity ($\lambda = 1$), viscous coupling between the metallic diapir and the surrounding material occurs and a fraction of gravitational energy is converted to heat up the diapir. However, for small viscosity ratios ($\lambda < \mathcal{O}(10)^{-2}$), viscous heating is restricted to the surrounding mantle and is concentrated at the diapir's poles where strain rates are greatest (Fig. 1, top left panel and Samuel et al., 2010). In this case, the initial temperature of the protocore is lower than in the isoviscous case and remains close to the initial temperature of the diapir (i.e. $T_{core} = T_0 + \Delta T_0$).

## 2.3. Dynamical model for the sinking of an iron diapir

### 2.3.1. Physical model

To investigate the dynamics of the metallic diapir separated after a large impact, we adapt the numerical finite volume model in spherical axisymmetric geometry of Monteux et al. (2009). In particular, neglecting the effects of adiabatic compression, conservation of energy applied to an undifferentiated planet of radius $R$ leads to

$$ \frac{DT}{Dt} = \frac{\nabla^2 T}{Ra_\chi} + D_\chi \frac{\eta}{\eta_0} \Omega. \qquad (8) $$

where $T$ and $t$ are dimensionless temperature and time ($T$ and $t$ are respectively normalized by $\Delta T_0$ and by a Stokes time $t_{Stokes} = \eta_0 / \Delta \rho g_0 R$).

$Ra_\chi = \frac{\Delta \rho_0 g_0 R^3}{\kappa \eta_0}$ is the compositional Rayleigh number and $D_\chi = \frac{\Delta \rho_0 g_0 R}{\rho_0 C_p \Delta T_0}$ is the dissipation number (with $\kappa$ the heat diffusivity and $\rho_0 C_p$ the average specific heat of the impacted body). $\Delta \rho_0$ is the density difference between metal and silicates. $g_0$ is the gravity at the surface of the impacted protoplanet. $\Delta T_0$ is the density difference between metal and silicates. $g_0$ is the gravity at the surface of the impacted protoplanet. $\Delta T_0$ is the initial impact induced temperature increase below the impact site. $\Omega$ is the dimensionless dissipation function and expresses the conversion of potential energy into heat:

$$ \Omega = \tau : \nabla \mathbf{v} \qquad (9) $$

where $\tau$ is the dimensionless deviatoric stress tensor and $\mathbf{v}$ is the non-dimensional velocity (normalized with $v_{Stokes} = R/t_{Stokes}$). Prior to impact, we assume a homogenous temperature $T_0$ of the growing planet. The viscosity is $\eta = \eta_0 \lambda^T$ where $\eta_0 = \eta(T_0)$ is the viscosity of the cold undifferentiated material at the start of the experiment (Ratcliff et al., 1997; Ziethe and Spohn, 2007).

Assuming that the protoplanet is incompressible, the other dimensionless governing equations are continuity

$$ \nabla \cdot \mathbf{v} = 0, \qquad (10) $$

and momentum conservation assuming infinite Prandtl number

$$ -\nabla P + \nabla \cdot \left( \frac{\eta}{\eta_0} \left[ \nabla \mathbf{v} + \left[ \nabla \mathbf{v} \right]^T \right] \right) + \left( \frac{T}{B} - f \right) \mathbf{r} \mathbf{e}_r = 0, \qquad (11) $$

where $P$ and $r$ are the non-dimensional pressure and radius ($P$ and $r$ are respectively normalized by $\eta_0 / t_{Stokes}$ and $R$). $\mathbf{e}_r$ is the radial unit vector and $f$ represents the volume fraction of metal. Following Monteux et al. (2009), we introduce the buoyancy ratio $B = \Delta \rho_0 / (\rho_0 \alpha \Delta T_0)$ where $\rho_0$ and $\alpha$ are the density and thermal expansion of the undifferentiated material (see Table 1 for values). Gravity depends on the radial position $r$ and $g(r) = \frac{4}{3} G \pi \rho_0 r = g_0 \frac{r}{R}$.

After an impact and the subsequent differentiation, the metallic phase and the purely silicate phase are separated within the thermal anomaly (Fig. 1, $t_3$). Hence a chemical readjustment occurs leading to the sinking of the dense iron though the lighter undifferentiated material. The buoyancy force that drives the flow of the diapir toward the center of the protoplanet increases with the volume fraction of metal $f$. In our models, we consider only two chemical phases meaning that $f$ varies between 1 (pure metal) and 0.17 (undifferentiated material with chondritic composition). During sinking, chemical diffusion is negligible. The metal volume fraction $f$ is then simply advected by the flow:

$$ \frac{\partial f}{\partial t} + \mathbf{v} \cdot \nabla f = 0. \qquad (12) $$

### 2.3.2. Simplifications of the governing equations

We make a number of approximations to simplify the analysis of the initial thermal state of the protocore. We neglect the thermal buoyancy $T/B$ in the momentum equation as the buoyancy number $B$ is very large (the density difference between metal and silicates is 140 to 560 times larger than the thermal density variations). In addition, in treating the metal-silicate separation, our model neglects multi-phase dynamics associated with core-mantle segregation (Golabek et al., 2008; Srámek et al., 2010).

Prior to the impact, the planet is undifferentiated meaning that the impact and the subsequent local differentiation lead to the first episode of core formation. Moreover, as this study focuses on the thermal state of the metallic phase once at the center of the planet, we do not consider the dynamical effects of the purely silicate phase that spreads beneath the impacted planet surface while the diapir sinks to the center of the planet (Monteux et al., 2009; Reese and Solomatov, 2006).

### 2.3.3. Numerical model

We implement a finite volume numerical model to solve Eq. 8, Eq. 10, Eq. 11 and Eq. 12 in axisymmetric spherical geometry in a 200 × 400 grid. We use a stream function formulation for the equations of motion with a direct implicit inversion method (Schubert et al., 2001). Eq. 8 and Eq. 12 are solved by an Alternating Direction Implicit (ADI) scheme (Peaceman and Rachford, 1955; Douglas, 1955). The stream function, temperature and compositional fields are described by a second-order approximation in space. To limit numerical diffusion when solving the transport equations, especially for the compositional field, we use a Total Variation Diminishing Superbee scheme (Laney, 1998; Roe, 1986) implemented in an implicit way (Srámek et al., 2010) which enables high resolution of purely advective fields. Mechanical boundary conditions are free-slip at the surface and along the symmetry axis. Thermal boundary conditions are isothermal at the surface and insulating along the symmetry axis.
# 3. Magnetic field considerations

## 3.1. Temperature of the protocore

Once at the center of the planet, core cooling depends on the temperature of the protocore and the thermal state and melting temperature of the silicate at the protocore-mantle boundary. As discussed above, the temperature of the protocore depends on the post-impact temperature difference $\Delta T_0$ (cf. Eq. 5). The maximal temperature increase after sinking is (Monteux et al., 2009):

$$ \Delta T \approx \Delta T_0 \left( 1 + a(\lambda) \left( \frac{r_0}{R} \right)^2 \frac{\Delta\Theta}{\Delta T_0} \right), \quad (13) $$

where $a(\lambda)$ is the fractional amount of viscous dissipation in the metal during the sinking of the iron toward the center of the planet. $r_0$ is the initial distance between the inertial center of the metallic diapir (noted i.c. in (Fig. 1, t$_3$)) and the center of the planet before the diapir starts to sink. $\Delta\Theta$ is the maximum temperature increase in the unrealistic case in which all the potential energy released by the diapir is used to heat up the metal phase. Assuming the impact velocity is the escape velocity of the planet (cf., Section 2) then $\Delta\Theta/\Delta T_0$ is constant, independent of the planet size and equal to 11.8 (Monteux et al., 2009) and the temperature of the metallic phase at the end of the sinking is

$$ T_{\text{core}} \approx T_0 + \Delta T = T_0 + \Delta T_0 \left( 1 + a(\lambda) \left( \frac{r_0}{R} \right)^2 \frac{\Delta\Theta}{\Delta T_0} \right). \quad (14) $$

During the sinking of the metallic diapir, its pressure increases and leads to additional heating by compressional effects. This heating process is inversely proportional to the buoyancy ratio $B$ (Monteux et al., 2009). In our models, $B$ is large enough to consider that this thermal effect is negligible compared to viscous heating (Samuel and Tackley, 2008).

## 3.2. Conditions for an early magnetic field

Three conditions are required to drive a dynamo on a growing planet:

1. The heat flow $Q$ out of the core must exceed the adiabatic value $Q_A$ such that convection can occur (Stevenson et al., 1983). This condition is

$$ Q > Q_A = \frac{k_{\text{Fe}} \alpha_{\text{Fe}} g_c T_{\text{Fe}}}{C_{p,\text{Fe}}} 4 \pi R_{\text{Fe}}^2 \quad (15) $$

where $k_{\text{Fe}}$, $\alpha_{\text{Fe}}$, $T_{\text{Fe}}$ and $C_{p,\text{Fe}}$ are respectively the thermal conductivity, the thermal expansion, the temperature and the heat capacity of the metallic phase and $g_c$ is the gravity at the surface of the protocore (see Table 1 for values). In our model, once the metallic phase has reached the center of the impacted planet, $T_{\text{Fe}}$ decreases monotonically with time due to heat conduction to the surrounding mantle.

2. The ratio of the rate at which gravitational potential energy is released by convection to the rate of ohmic dissipation, $\Phi$, must exceed a critical value (Buffett, 2002):

$$ \frac{Q}{\Phi} > \frac{1}{\epsilon_T}. \quad (16) $$

Here, $\epsilon_T$ is the Carnot-style efficiency for thermal convection. We note that a dynamo can be generated even with a subadiabatic heat flow if the compositional convection is important (Buffett, 2002; Christensen and Wicht, 2008). However, here we consider neither this mechanism or effects related to the presence of an inner core. Assuming that the characteristic length scale of the

flow leading to magnetic field generation is the radius of the protocore, $\Phi$ can be approximated as (Buffett, 2002):

$$ \Phi = \left( \frac{\nu \bar{B}^2}{\mu} \right)^{\frac{4}{3}} \pi R_{\text{Fe}}, \quad (17) $$

where $\bar{B}$ is the average strength of the magnetic field within the core, $\nu$ is the magnetic diffusivity and $\mu$ is the magnetic permeability (see Table 1 for values). In the absence of constraints on $\bar{B}$ for early planets, we assume a current Earth-like value of 2.5 mT (Kuang and Bloxham, 1997) that is independent of the protocore size or the planetary radius. The efficiency of thermal convection is then given by (Buffett, 2002)

$$ \epsilon_T = \frac{0.8\pi \alpha_{\text{Fe}} G \rho_{\text{Fe}} R_{\text{Fe}}^2}{3 C_{p,\text{Fe}}} \left(1 - \frac{Q_A}{Q}\right), \quad (18) $$

where $\rho_{\text{Fe}}$ is the density of the metallic phase and $G$ is the gravitational constant.

3. The structure of the convective motions carrying magnetic field lines must be sufficiently complicated to favor self-sustaining dynamo action. A measure of this complexity is that the magnetic Reynolds number (Christensen and Aubert, 2006)

$$ Re_m = \frac{UL}{\nu} > Re_m^{\text{crit}} = O(10-10^2). \quad (19) $$

Here, $L$ and $U$ are the characteristic length and velocity scales for the flow within the protocore and $\nu$ is the magnetic diffusivity of the metal phase. Whereas the natural length scale in the problem is the depth of the convecting iron layer, the choice of an appropriate velocity scale depends on the leading order force balance (Christensen, 2010). As the rotation rate of growing planets is potentially time-dependent and poorly constrained a convenient and reasonable choice is based on a balance between inertial and buoyancy forces and is (Stevenson, 2003a):

$$ U \sim \left( \frac{q \alpha_{\text{Fe}} g_c R_{\text{Fe}}}{\rho_{\text{Fe}} C_{p,\text{Fe}}} \right)^{1/3}, \quad (20) $$

where $q$ is the heat flux out of the core ($q = Q/(4\pi R_{\text{Fe}}^2)$). Taking $L = R_{\text{Fe}}$, $g_c = 4/3\pi G\rho_{\text{Fe}}R_{\text{Fe}}$, the combination of Eqs. 19 and Eq. 20 leads to the condition

$$ Re_m = \left( \frac{4\pi G q \alpha_{\text{Fe}} R_{\text{Fe}}^2}{3 C_{p,\text{Fe}}} \right)^{1/3} \frac{R_{\text{Fe}}}{\nu} > O(10-10^2) \quad (21) $$

Among the three criteria above, the first is typically considered a necessary condition for a thermally-driven dynamo. (Subadiabatic dynamos are possible if e.g., compositional gradients also drive convection.) As we shall see in Section 4, this condition ($Q > Q_A$) is easily met in our models. However it is not a sufficient condition for dynamo action, as indicated by the other two conditions. The second condition assumes an Earth-like dynamo: in particular it assumes an Earth-like value for the root mean square magnetic field in the core, a particular efficiency for conversion of gravitational potential energy into magnetic energy and a particular toroidal-to-poloidal field conversion efficiency. The third condition is a more general criterion for dynamo action, and whereas scalings other than Eq. 20 for the velocity field are possible, we have chosen one that does not depend upon knowing the early rotation rate of the planet or moon.
### 3.3. Heat flow across the protocore: analytic considerations

The three conditions required for an early dynamo (Eq. 15, Eq. 16 and Eq. 21) all impose conditions on the heat flow out of the protocore. At the protocore-mantle boundary there are two asymptotic heat transfer regimes. The first is a conductive transfer across a thermal boundary layer growing on the mantle side where the corresponding heat flow $Q_{cond}$ is:

$$Q_{cond} \sim \frac{4\pi k \Delta T R_{Fe}^2}{\delta} \qquad (22)$$

where $\delta \approx 2\sqrt{\kappa t}$ is the thickness of the conductive boundary layer at any time $t$. As $t$ becomes large $\delta = c_2 R_{Fe}$ (where $c_2$ is an order 1 constant (Monteux et al., 2009)), $k$ is the thermal conductivity of the surrounding material. For computational reasons, in our numerical models, we make the approximation that thermal conductivities for metallic and undifferentiated material are the same (i.e., $k = k_{Fe} = k_{Si}$, Table 1). In reality, the thermal conductivity of the core material is ≈ 10 times bigger than the conductivity for the silicates. Hence, this simplification underestimates the protocore heat flow and our criteria for dynamo action is more restrictive than in the more realistic case of unequal silicate and metal thermal conductivities.

As illustrated in Fig. 2, the temperature of the protocore might overcome the solidus temperature of the surrounding material depending on its chemical composition and on the impacted planetary radius. Hence, a melt layer may form above the protocore (see also Section 5.2). As this melt layer grows in thickness with time (and assuming it is greater than a few tens of meters thick), vigorous convection will enhance core cooling over the conductive case and lead to a heat flow that is independent of the thickness of the melt layer and of the form (Turner, 1973):

$$Q_{conv} \sim 0.8\pi R_{Fe}^2 \left( \frac{\rho_{Si}^2 g_e \alpha_{Si} k^2 C_{p.Si}}{\mu_{melt}} \right)^{1/3} (T_{core} - T_{melt.Si})^{4/3} \qquad (23)$$

where $\rho_{Si}$, $\alpha_{Si}$, $C_{p.Si}$ and $\mu_{melt}$ are respectively the density, the thermal expansion, the heat capacity and the viscosity of the molten silicate layer (see Table 1 for typical values). Eventually, the heat flux out of the core will, of course, become balanced by the heat flux across the

<figure><img src="image_157.png" /><figcaption>Fig. 2. Initial temperature of the protocore once at the center of the planet as a function of the planetary radius. Each symbol represents a numerical result for a given $R_{Fe}$ ($R_{Fe}$ in the range of 150–500 km) for a uniform viscosity. Theoretical predictions for $T_{core}$ from Eq. 14 are given for uniform viscosity (black solid line, $a(\lambda) = 15\%$) and temperature dependent viscosity (black dashed line, $a(\lambda) = 0\%$). The black solid line is without melting and the gray one is with melting. The region in between defines the range of initial protocore temperature (see Eq. 5). The dotted line shows $T = T_0$. The chondritic solidus (Agee, 1997) and the dry peridotite solidus (Zerr and Boehler, 1994) at $P = P_{CMB}$ are represented with dashed-double-dotted lines (see Section 5.2 for calculation details).</figcaption></figure>

overlying solid mantle and the melt layer will have a quasi-steady thickness.

## 4. Numerical results

### 4.1. Diapir descent: the initial protocore temperature

In our models, where we specify the viscosity of the impacted planet as $\eta_0 = 10^{22}$ Pas, metallic diapirs descend to the center of the planet relatively quickly (e.g. in less than 10 Myr for $R_{Fe} = 400$ km) (Monteux et al., 2009). During the sinking of the diapir the temperature difference between the protocore and the mantle increases rapidly in response to viscous heating depending on the fraction of gravitational energy converted into heat ($a(\lambda)$). Then this temperature difference decreases as the diapir cools once it reaches the center of the planet.

We investigate $T_{core}$ for diapirs in the range of 150–500 km and planet radii in the range of 2000–4000 km. For the diapir radii range studied here, the protocore temperature at the end of sinking is given by Eqs. 13 and 14. The results from our numerical models for uniform viscosity contrasts are shown in Fig. 2 and are in agreement with Eq. 14. Theoretical predictions of $T_{core}$ are represented in Fig. 2 with and without taking into account the effect of melting in the isobaric core after impact (gray and black solid lines respectively). The region in between defines the range of initial temperature and the difference between these two theoretical models is within the range of uncertainties of $\gamma$. For a uniform viscosity, viscous coupling between metal and undifferentiated mantle material is important and $a(\lambda) \approx 15\%$ (Monteux et al., 2009). As the viscosity contrast between the hot metallic phase and the surrounding mantle increases, the fraction of gravitational potential energy converted into heat in the metallic phase decreases (see Section 2.2). That is, as the mantle/core viscosity ratio becomes very large $a(\lambda) \to 0$ and viscous dissipation is concentrated in the mantle (Monteux et al., 2009). In this limit, the initial temperature of the protocore after the sinking is $T_{core} = T_0 + \Delta T_0$. The predictions of Eq. 14 for this asymptotic case are shown in Fig. 2 (dashed line). As expected, the reduced dissipation in the metallic diapir leads to a lower predicted $T_{core}$.

A high temperature contrast between the protocore and the surrounding undifferentiated mantle can eventually lead to the formation of a molten layer above the protocore-mantle boundary. In Fig. 2, we include the chondrite solidus temperature from Agee (1997) and the dry peridotite solidus from Zerr and Boehler (1994) at the center of the impacted planet as a function of the planetary radius. This effect will be more precisely discussed in Section 5.2.

Because of the practical difficulty in obtaining accurate solutions for very large viscosity variations, we restrict our numerical studies to the isoviscous case, which is well-explained by our theory. The agreement between our isoviscous theory and the numerical solutions gives us confidence that application of the theory for the variable viscosity case is justified. From Fig. 2 these isoviscous solutions will overestimate the core temperature by $\approx 10-40\%$, depending on the protoplanetary radius. Thus, quantitatively, the uncertainties of this enhanced heat flow for predicting dynamo action are not so large. From Eq. 13, assuming $a(\lambda) = 15\%$, the temperature difference between the protocore and the cold surrounding material in the isoviscous model is overestimated by a factor $\approx 2.8$. From Eq. 21 $R_{m} \propto Q^{1/3} \propto \Delta T^{1/3}$ and, thus, the magnetic Reynolds number will be a factor of about 1.4 too large (i.e., it is still accurate to $O(1)$).

### 4.2. Dynamo initiation in a conductive heat transfer case

To determine whether or not a dynamo is initiated in the protocore, we monitor the average heat flow $Q$ across the protocore-mantle boundary as a function of time for a wide range of $R - R_{Fe}$ conditions. We stop calculations 50 Myr after the impact because at this time the