<figure><img src="image_43.png" /><figcaption>Fig. 10. (continued)</figcaption></figure>

<caption>**Table 3**</caption>

Values of $b_1$, $b_2$, $b_3$, $b_4$ and $b_5$ obtained by fitting the results from the numerical models (Fig. 10a and b) with Eqs. (2a) and (2b). The values for models with strength are separated from those of hydrodynamic models with no strength and are shown as a function of $R_{imp}$ and $V_{imp}$.

<table><thead><tr><th>Models with strength</th><th>b1</th><th>b2</th><th>b3</th></tr></thead><tbody><tr><td>$R_{imp}$ = 50 km ($V_{imp}$ = 10 km/s)</td><td>2.16</td><td>$1 \times 10^{-2}$</td><td>$-2.4 \times 10^{-3}$</td></tr><tr><td>$R_{imp}$ = 100 km ($V_{imp}$ = 10 km/s)</td><td>2.15</td><td>$1.2 \times 10^{-3}$</td><td>$-2.4 \times 10^{-3}$</td></tr><tr><td>$R_{imp}$ = 150 km ($V_{imp}$ = 10 km/s)</td><td>2.1</td><td>$1 \times 10^{-2}$</td><td>$-2 \times 10^{-3}$</td></tr><tr><td>$R_{imp}$ = 200 km ($V_{imp}$ = 10 km/s)</td><td>2.1</td><td>$1 \times 10^{-3}$</td><td>$-2.4 \times 10^{-3}$</td></tr><tr><td>$V_{imp}$ = 5 km/s ($R_{imp}$ = 50 km)</td><td>2.17</td><td>$1.1 \times 10^{-2}$</td><td>$-2.3 \times 10^{-3}$</td></tr><tr><td>$V_{imp}$ = 10 km/s ($R_{imp}$ = 50 km)</td><td>2.16</td><td>$1 \times 10^{-2}$</td><td>$-2.4 \times 10^{-3}$</td></tr><tr><td>$V_{imp}$ = 15 km/s ($R_{imp}$ = 50 km)</td><td>2.16</td><td>$3 \times 10^{-4}$</td><td>$-1.8 \times 10^{-3}$</td></tr><tr><td>$V_{imp}$ = 20 km/s ($R_{imp}$ = 50 km)</td><td>2.16</td><td>$9 \times 10^{-3}$</td><td>$-2. \times 10^{-3}$</td></tr></tbody></table>

<table><thead><tr><th>Hydrodynamic models</th><th>b4</th><th>b5</th></tr></thead><tbody><tr><td>$R_{imp}$ = 50 km ($V_{imp}$ = 10 km/s)</td><td>1.35</td><td>$-7.5 \times 10^{-2}$</td><td>-</td></tr><tr><td>$R_{imp}$ = 100 km ($V_{imp}$ = 10 km/s)</td><td>1.3</td><td>$-4.7 \times 10^{-3}$</td><td>-</td></tr><tr><td>$R_{imp}$ = 150 km ($V_{imp}$ = 10 km/s)</td><td>1.28</td><td>$-1.2 \times 10^{-3}$</td><td>-</td></tr><tr><td>$R_{imp}$ = 200 km ($V_{imp}$ = 10 km/s)</td><td>1.33</td><td>$-6.2 \times 10^{-3}$</td><td>-</td></tr><tr><td>$V_{imp}$ = 5 km/s ($R_{imp}$ = 50 km)</td><td>1.3</td><td>$-2.2 \times 10^{-4}$</td><td>-</td></tr><tr><td>$V_{imp}$ = 10 km/s ($R_{imp}$ = 50 km)</td><td>1.35</td><td>$-7.5 \times 10^{-2}$</td><td>-</td></tr><tr><td>$V_{imp}$ = 15 km/s ($R_{imp}$ = 50 km)</td><td>1.29</td><td>$-4 \times 10^{-3}$</td><td>-</td></tr><tr><td>$V_{imp}$ = 20 km/s ($R_{imp}$ = 50 km)</td><td>1.28</td><td>$-3.8 \times 10^{-3}$</td><td>-</td></tr></tbody></table>

model with strength for $V_{imp}$ = 10 km/s and $R_{imp}$ = 50 km (top panels) and $V_{imp}$ = 10 km/s and $R_{imp}$ = 150 km (bottom panels). As mentioned before, the shock pressure is not estimated by the scaling laws where no refracted wave propagates (i.e. in the equatorial region of the core and

the lower part of the mantle beneath the shock ray that is tangent to the CMB). Hence, we cannot estimate the post-impact temperature increase in these regions (illustrated by black areas in Fig. 14). This figure shows that for a 50 km radius impactor, the region with significant impact heating is localized in a small shallow volume. For a 150 km impactor, a larger fraction of the impacted mantle is affected. Below the impact site the temperature increase can reach ≈1900 K for both the hydrodynamic model and the model with strength. The largest temperature increase is localized in Zone 1 and the temperature increase in Zone 2 is weaker. The main difference between the two rheologies is the depth of the impact temperature increase, which is deeper in the hydrodynamic model than in the model with strength, because the pressure decay is more rapid in the model with strength (Fig. 6 a, bottom left panel). Moreover, the regions with increased temperature have shapes that are more prolate than spherical, because of the influence of the $\beta$ angle on the pressure decay in Zone 1. The difference in the expression of the pressure decay between the two models as a function of $\beta$ should lead to a difference in sphericity (Fig. 6 a). However, the difference is not significant in the cases illustrated in Fig. 14.

## 5. Discussion and conclusions

We have modeled the shock pressure distributions inside the impacted hemisphere of a differentiated Mars size planet using 2D axisymmetric hydrocode simulations (iSALE-2D) for vertical impacts with impact velocities of 5–20 km/s and impactor radii ranging from 50 to 200 km. From these hydrocode simulations, we have observed a distribution of the peak shock pressure within the impacted mantle and
<figure><img src="image_44.png" /><figcaption>Fig. 11. $A_4$ and $A_5$ values (Eq. (1c): $P_{core} = A_4 (D/R_{imp})^{-A_5}$) as functions of the impactor radius (left panels) and impact velocity (right panels).</figcaption></figure>

proposed a mathematical formalism adapted to this distribution: a first zone (Zone 1) where the peak shock pressure decreases exponentially and a second zone (Zone 2) where the peak shock pressure decrease follows a power-law. With this formalism, the peak shock pressure does not go to infinity as $D/R_{imp}$ goes to 0, and the peak shock pressure decays relatively slowly from the impact site to $D/R_{imp} = 3–4$, which includes the canonical “isobaric core” with slowly decaying shock pressure with distance from the impact site and the “pressure decay regime” region proposed by Pierazzo et al. [1997]. We have shown that the ray angle $\beta$ has minor influence on the pressure decay pattern for $\beta < 65–70^\circ$.

The most important contribution to the peak shock pressure distribution is related to the rheology of the material. Adopting a model where the impacted mantle is treated using models with strength can result in a significant decrease of the peak shock pressure compared to that of the hydrodynamic models. This difference in behavior occurs principally in the shallower zone (Zone 1) while in Zone 2 the behavior is similar in both models. We have shown that the peak shock pressure is about an order of magnitude smaller at the bottom of the Zone 1 in the models with strength compared to that in the hydrodynamic models.

The characteristic of the peak shock pressure in Zone 1 has important consequences on the peak shock pressure in Zone 2, since it determines the pressure in the upper part of Zone 2. Also, the mantle rheology affects the peak shock pressure at the CMB. We have combined our shock-pressure scaling laws with the formalism proposed by Watters et al. [2009] to determine the impact heating of a differentiated Mars model by a large impactor. Our results show that the mantle rheology governs the impact heating of a planet and the impact heating may be overestimated when the material strength is not accounted for.

The Martian rheology and its evolution since Mars formation are

poorly constrained. However, the two models used in our study (hydrodynamic and with strength) enable us to cover a range of rheology models compatible with different stages of the Martian evolution. Our hydrodynamic models where the mantle is considered as an inviscid fluid may be an interesting proxy to characterize the impact heating during a magma ocean stage [Monteux et al., 2016]. On the other hand, our models with strength are more appropriate to characterize the large impact consequences within a solidified and cold Mars [Watters et al., 2009]. Our results show that the depth of the heated zone is different in the two cases (Fig. 14) with a larger thermal anomaly in the hydrodynamics models than in the models with strength as the pressure decay is more rapid in the latter case. In both cases, the shape of the impact-heated region is prolate because of the influence of $\beta$ on the pressure decay. These effects need to be accounted for in accretionary models especially when about $3 \times 10^5$ planetesimals of ~100 km in diameter are needed to build a Mars size planet.

As illustrated in Fig. 13, our parameterized approach is associated to a mean error between our scaling laws and the direct results obtained by the iSALE models that ranges from 10 to 14% within the impacted hemisphere. As a consequence, the temperature increase derived from our scaling laws is also associated to this error range. Nonetheless, in the light of the uncertainty that goes with our understanding of the early Solar System history and our knowledge of the effects of large impacts, our parameterized approach to calculate the temperature increase within a differentiated Mars-size protoplanet can be considered as a first order estimation.

As a consequence, the next step toward a better understanding of the early evolution of terrestrial planets is to implement more realistic rheological models in the hydrocode simulations. Such an effort would lead to more reliable scaling laws that could easily be implemented in accretion models to characterize the influence of large impacts on the
<figure><img src="image_45.png" /><figcaption>**Fig. 12.** Position of the peak shock pressure as a function of the time after the impact for the models with strength (top panels) and the hydrodynamic models (bottom panels).</figcaption></figure>

<caption>**Table 4**</caption>

Velocity of the peak shock pressure propagation in the mantle and in the core for the models with strength and hydrodynamic models, obtained fitting the results from Fig. 12 with a linear expression.

<table><thead><tr><th></th><th>Target mantle</th><th>Target core</th></tr></thead><tbody><tr><td colspan="3"><strong>Models with strength</strong></td></tr><tr><td>R<sub>imp</sub> = 50 km (V<sub>imp</sub> = 10 km/s)</td><td>8.84 km/s</td><td>6.38 km/s</td></tr><tr><td>R<sub>imp</sub> = 100 km (V<sub>imp</sub> = 10 km/s)</td><td>8.52 km/s</td><td>6.44 km/s</td></tr><tr><td>R<sub>imp</sub> = 150 km (V<sub>imp</sub> = 10 km/s)</td><td>8.58 km/s</td><td>6.47 km/s</td></tr><tr><td>R<sub>imp</sub> = 200 km (V<sub>imp</sub> = 10 km/s)</td><td>8.93 km/s</td><td>6.50 km/s</td></tr><tr><td>V<sub>imp</sub> = 5 km/s (R<sub>imp</sub> = 50 km)</td><td>9.04 km/s</td><td>6.33 km/s</td></tr><tr><td>V<sub>imp</sub> = 10 km/s (R<sub>imp</sub> = 50 km)</td><td>8.84 km/s</td><td>6.38 km/s</td></tr><tr><td>V<sub>imp</sub> = 15 km/s (R<sub>imp</sub> = 50 km)</td><td>8.80 km/s</td><td>6.36 km/s</td></tr><tr><td>V<sub>imp</sub> = 20 km/s (R<sub>imp</sub> = 50 km)</td><td>8.60 km/s</td><td>6.38 km/s</td></tr><tr><td colspan="3"><strong>Hydrodynamic models</strong></td></tr><tr><td>R<sub>imp</sub> = 50 km (V<sub>imp</sub> = 10 km/s)</td><td>6.96 km/s</td><td>6.36 km/s</td></tr><tr><td>R<sub>imp</sub> = 100 km (V<sub>imp</sub> = 10 km/s)</td><td>7.60 km/s</td><td>6.40 km/s</td></tr><tr><td>R<sub>imp</sub> = 150 km (V<sub>imp</sub> = 10 km/s)</td><td>8.33 km/s</td><td>6.46 km/s</td></tr><tr><td>R<sub>imp</sub> = 200 km (V<sub>imp</sub> = 10 km/s)</td><td>8.77 km/s</td><td>6.53 km/s</td></tr><tr><td>V<sub>imp</sub> = 5 km/s (R<sub>imp</sub> = 50 km)</td><td>6.76 km/s</td><td>6.41 km/s</td></tr><tr><td>V<sub>imp</sub> = 10 km/s (R<sub>imp</sub> = 50 km)</td><td>6.96 km/s</td><td>6.36 km/s</td></tr><tr><td>V<sub>imp</sub> = 15 km/s (R<sub>imp</sub> = 50 km)</td><td>8.33 km/s</td><td>6.46 km/s</td></tr><tr><td>V<sub>imp</sub> = 20 km/s (R<sub>imp</sub> = 50 km)</td><td>8.79 km/s</td><td>6.53 km/s</td></tr></tbody></table>

martian mantle dynamics [e.g. Roberts and Arkani-Hamed, 2012, 2017] or on other terrestrial planets [e.g. Roberts and Barnouin, 2012]. As large impacts are not isolated phenomena during the late regimes of planetary accretions, the scaling laws we have developed are suitable to be implemented in thermal evolution models that consider the accretion histories of growing planets [e.g. Monteux et al., 2014; Arkani-

Hamed, 2017].

The followings are some shortcomings that require improvements of our scaling laws in the future to be better suited for shock waves travelling inside an actual planet, rather than the simple two-layered planet model where a uniform mantle overlies a uniform core.

* We have considered vertical impacts, whereas majority of impacts during the accretion of a planet is oblique. The peak shock pressure produced by an oblique impact along the vertical direction beneath the impact site depends on the sine of the impact angle relative to the local horizon (see Fig. 3 of Pierazzo and Melosh, 2000a, 2000b). This implies a shock pressure reduction by a factor of ~0.71 for the most probable impact angle of 45° [Shoemaker, 1962] compared to that of a vertical impact, assuming that all other parameters are equal. Moreover, a vertical impact on a spherically symmetric planet preserves the axial symmetry. Because the shock wave velocity in the core is smaller than that in the mantle, the axial symmetry results in the convergence of the refracted wave at the CMB of the impacted hemisphere toward the axis in the core, thus artificially increasing the shock pressure near the axis. However, the near axis enhancement of the shock pressure is negligible except very close to the antipodal CMB [Ivanov et al., 2010; Arkani-Hamed and Ivanov, 2014]. We also note that this artifact does not occur in the mantle because the shock rays are diverging away from the axis of symmetry. It is desirable to adopt oblique impacts in investigating the accretion of a planet. In an oblique impact, the axial symmetry is replaced by a plane symmetry, which avoids the convergence of the waves toward the axis in the expense of introducing more complex and computationally intensive three-dimensional (3D) calculations,
<figure><img src="image_46.png" /><figcaption>**Fig. 13.** Comparison of the maximum shock pressure experienced after an impact within a differentiated Mars size planet with $R_{\text{imp}} = 50$ km and $V_{\text{imp}} = 10$ km/s. (a) is for model with strength and (b) is for Hydrodynamic model. The left panel represents the result from the iSALE model and the middle panel represents the model obtained from our scaling laws (Eqs. (1a), (1b) and (1c)). The right panel shows the error relative to the iSALE model (in %). Note that only the impacted hemisphere is shown in this figure, because the scaling laws are derived for this hemisphere.</figcaption></figure>

besides the substantial increase of the grid points in the computation domain.

* The distance travelled by a shock wave from the impact site to a given point in the planet is approximated by a straight line, D. In our simplified two-layered planet model with a uniform mantle overlapping a uniform core, the first arrival shock wave travels in quiescent medium, hence the shock rays are straight lines until they interact with the core mantle boundary where they bend according to the Snell's law. Inside a real planet the rays are actually curved due to the pre-impact non-uniform density and temperature gradients in the interior of the planet. Because of the ray curvature a wave travels longer to reach an off axis location in the interior than D. Therefore, our model slightly over estimates the shock pressure experienced in that location. Moreover, because of the ray curvature an incident angle at the CMB and the corresponding refraction angle could be slightly different than those determined by the straight line approximation. However, the core radius is about one half of the planet radius, and the angle $\beta$ is always smaller than 30°. Bearing in mind that the shock ray along the axis of symmetry is a straight line, the departure of the shock rays from the corresponding straight lines

is negligible for these $\beta$ values. As also mentioned above, the shock wave propagation weakly depends on the $\beta$ angle that is ranging between 0° (vertical) and 45°, which is already concluded by previous studies [Pierazzo et al., 1997; Monteux and Arkani-Hamed, 2016]. More importantly, providing scaling laws along straight lines, rather than along shock rays, helps to easily determine the shock pressure in the entire region.

* We have considered only direct shock waves to determine our scaling laws. Because the laws are based on shock ray theory, they can be used to determine the pressure of secondary shock waves, such as those reflected at the surface and those reflected at the antipodal CMB of the core [Arkani-Hamed and Ivanov, 2014]. In general a reflect wave arriving at a location travels longer and is weaker than a direct wave. However, because of the zero pressure boundary condition at the surface the reflected wave is rarefaction and upon interference with the direct compressional wave effectively reduced the pressure experienced at a given location close to the surface [Melosh, 1989; Arkani-Hamed, 2005; Louzada and Stewart, 2009]. The interference of the direct and secondary waves is not considered in our present study.
<figure><img src="image_47.png" /><figcaption>**Fig. 14.** Temperature increase after an impact within a differentiated Mars size planet. The first line: $R_{\text{imp}} = 50$ km and $V_{\text{imp}} = 10$ km/s, and the second line: $R_{\text{imp}} = 150$ km and $V_{\text{imp}} = 10$ km/s (Left: Model with strength, Right: Hydrodynamic model). All panels represents the model obtained from our scaling laws (Eqs. (1a), (1b) and (1c)) combined with the impact induced temperature increase model (Eqs. (4)–(6)) from Watters et al., (2009) (with $C = 7.24$ km/s, $S = 1.25$, $\rho = \rho_m$ for the mantle and $C = 4$ km/s, $S = 1.6$, $\rho = \rho_c$ for the core [Arkani-Hamed and Olson, 2010a, b]). The white dashed line represents the CMB.</figcaption></figure>

* The size of the Zone 1 depends on the impactor size. For large impactors Zone 1 directly interferes with the core, and Zone 2 does not exist. For example, there is no Zone 2 for $\beta$ values smaller than $15^\circ$ when the impactor is larger than 150 km in radius. We propose the following procedure to estimate the shock pressure distribution inside the mantle and the core of a planet using our scaling laws. The impact-induced temperature increase may then be determined using, for example Watters et al. [2009] model.

* Consider a terrestrial planet with a radius $R_p$ and a core radius $R_c$ that is impacted by an impactor of radius $R_{\text{imp}}$. We assume that the Zone 1 is a hemisphere with a radius $R_1 = 10 \times R_{\text{imp}}$ centered at the impact site. There is no interaction of Zone 1 with the core for $R_1 < (R_p - R_c)$, and the equations provided in this paper are valid. This is also the case for any shock wave specified by a ray angle $\beta$ when $\sin(\beta) > R_c/R_p$. For $R_1 > (R_p - R_c)$, Zone 1 may interact with the core. Because the shock wave velocity is greater than the particle velocity, it takes time for the CMB to deform, during which the shock wave has already entered the core, hence the incident shock wave impinges on a undeformed spherical core. Therefore, a

Zone 2 exists for the colatitudes greater than $\theta_1$, where $\cos(\theta_1) = (R_p^2 + R_c^2 - R_1^2)/2 R_p R_c$, and the equations provided in this paper are still valid. Zone 1 interacts with the core only for the colatitudes smaller than $\theta_1$.

## Acknowledgements

This research was supported by the Auvergne Fellowship program to JM, and by the Natural Sciences and Engineering Research Council (NSERC) of Canada to JAH. We gratefully acknowledge the developers of iSALE (www.isale-code.de), particularly the help we have received from Gareth S. Collins. We also would like to thank the two anonymous reviewers for their constructive comments. This paper is Laboratory of Excellence ClerVolc contribution no. 350.

## References

* Agnor, C.B., Asphaug, E., 2004. Accretion efficiency during planetary collisions. Astrophys. J. 613, L157–L160.

* Agnor, C.B., Canup, R.M., Leviston, H.F., 1999. On the character and consequences of large impacts in the late stage of terrestrial planet formation. Icarus 142, 219–237.
Ahrens, T.J., O'Keefe, J.D., 1987. Impact on the Earth, ocean and atmosphere. Int. J. Impact Eng. 5, 13–32.

Ahrens, T.J., O'Keefe, J.D., Merrill, R.B., 1977. Equations of State and Impact-Induced Shock-wave Attenuation on the Moon. Impact and Explosion Cratering. Pergamon Press, Elmsford, N.Y., pp. 639–656.

Amsden, A., Ruppel, H., Hirt, C., 1980. SALE: A Simplified ALE Computer Program for Fluid Flow at All Speeds. Los Alamos National Laboratories Report, LA-8095:101p. LANL, Los Alamos, New Mexico.

Arkani-Hamed, J., 2005. Magnetic crust of Mars. J. Geophys. Res. 110, E08005. https://doi.org/10.1029/2004JE002397.

Arkani-Hamed, J., 2017. Formation of a solid inner core during the accretion of Earth. J. Geophys. Res. Solid Earth 122. https://doi.org/10.1002/2016JB013601.

Arkani-Hamed, J., Ivanov, B., 2014. Shock wave propagation in layered planetary embryos. Phys. Earth Planet. Inter. 230, 45–59.

Arkani-Hamed, J., Olson, P., 2010a. Giant impact stratification of the Martian core. Geophys. Res. Lett. 37, L02201. https://doi.org/10.1029/2009GL04141.

Arkani-Hamed, J., Olson, P., 2010b. Giant impacts, core stratification, and failure of the Martian dynamo. J. Geophys. Res. 115, E07012.

Barr, A., Citron, R., 2011. Scaling of melt production in hypervelocity impacts from high-resolution numerical simulations. Icarus 211, 913–916.

Benz, W., Cameron, A.G.W., Melosh, H.J., 1989. The origin of the moon and the single-impact hypothesis III. Icarus 81, 113–131.

Bierhaus, M., Wünnemann, K., Elbeshausen, D., 2012. Numerical Modeling of Basin-forming Impacts: Implications for the Heat Budget of Planetary Interiors, Large Meteorite Impacts and Planetary Evolution V. LPI Contribution No. 1737.

Bierhaus, M., Wünnemann, K., Ivanov, B.A., 2013. Affect of Core Rheology on Shock Wave Propagation in Planetary Scale Impacts, 43rd Lunar and Planetary Science Conference. (Abstract # 2174.pdf).

Bjorkman, M.D., Holsapple, K.A., 1987. Velocity scaling impact melt volume. Int. J. Eng. Lett. 5, 155–163.

Bottke, W.F., Nolan, M.C., Greenberg, R., Kolvoord, R.A., 1994. Collisional lifetimes and impact statistics of near-Earth asteroids. In: Tucson, Gehrels T. (Ed.), Hazards Due to Comets and Asteroids. The University of Arizona Press, Arizona, pp. 337–357.

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

Chambers, J.E., 2013. Late-stage planetary accretion including hit-and-run collisions and fragmentation. Icarus 224, 43–56.

Chambers, J.E., Wetherill, G.W., 1998. Making the terrestrial planets: N-body integrations of planetary embryos in three dimensions. Icarus 136, 304–327.

Collins, G.S., Melosh, H.J., 2014. Improvements to ANEOS for multiple phase transitions. In: 45th Lunar and Planetary Science Conference, pp. 2664.

Curran, D.R., Shockey, D.A., Seaman, L., Austin, M., 1977. Mechanics and models of cratering in earth media. In: Roddy, D.J., Pepin, R.O., Merrill, R.B. (Eds.), Impact and Explosion Cratering. Pergamon, New York, pp. 1057–1087.

Davison, T.M., Collins, G.S., Ciesla, F.J., 2010. Numerical modelling of heating in porous planetesimal collisions. Icarus 208, 468–481. https://doi.org/10.1016/j.icarus.2010.01.034.

Frey, H., 2008. Ages of very large impact basins on Mars: implications for the late heavy bombardment in the inner Solar System. Geophys. Res. Lett. 35, L13203.

Gault, D.E., Heitowit, E., 1963. The partition of energy for hypervelocity impact craters formed in rock. In: Proceedings of the 6th Hypervelocity Impact Symposium. vol. 2. pp. 419–456.

Ghods, A., Arkani-Hamed, J., 2011. Effects of the Borealis impact on the mantle dynamics of Mars. Phys. Earth Planet. Inter. 188, 37–46.

Ivanov, B.A., Deniem, D., Neukum, G., 1997. Implementation of dynamic strength models into 2D hydrocodes: applications for atmospheric breakup and impact cratering. Int. J. Impact Eng. 20, 411–430.

Ivanov, B.A., Melosh, H.J., Pierazzo, E., 2010. Basin-forming impacts: reconnaissance modeling. Geol. Soc. Am. Spec. Pap. 465, 29–49.

Kaula, W.M., 1979. Thermal evolution of Earth and moon growing by planetesimal impacts. J. Geophys. Res. 84, 999–1008.

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.

Kokubo, E., Genda, H., 2010. Formation of the terrestrial planets form protoplanets under a realistic accretion condition. Astrophys. J. Lett. 714, L21–L25.

Kokubo, E., Ida, S., 2000. Formation of protoplanets from planetesimals in the solar nebula. Icarus 143, 15–27.

Kraus, R.G., Senft, L.E., Stewart, S.T., 2011. Impacts onto H2O ice: scaling laws for melting, vaporization, excavation, and final crater size. Icarus 214, 724–738.

Leinhardt, Z.M., Stewart, S.T., 2009. Full numerical simulations of catastrophic small body collisions. Icarus 199, 542–559.

Louzada, K.L., Stewart, S.T., 2009. Effects of planet curvature and crust on the shock pressure field around impact basins. Geophys. Res. Lett. 36, L15203.

Melosh, H.J., 1989. Impact Cratering: A Geologic Process. Oxford Univ. Press, New York.

Mitani, N.K., 2003. Numerical simulations of shock attenuation in solids and reevaluation of scaling law. J. Geophys. Res. 108 (E1), 5003. https://doi.org/10.1029/2000JE001472.

Monteux, J., Arkani-Hamed, J., 2014. Consequences of giant impacts in early Mars: Core merging and Martian dynamo evolution. J. Geophys. Res. 119, 480–505.

Monteux, J., Arkani-Hamed, J., 2016. Scaling laws of impact induced shock pressure and particle velocity in planetary mantle. Icarus 264, 246–256.

Monteux, J., Tobie, G., Choblet, G., Le Feuvre, M., 2014. Can large icy moons accrete undifferentiated? Icarus 237, 377–387.

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

Ohnaka, M., 1995. A shear failure strength law of rock in the brittle-plastic transition regime. Geophys. Res. Lett. 22, 25–28.

Osinski, G.R., Pierazzo, E., 2012. Impact Cratering: Processes and Products. Blackwell Publishing Ltd (330 pp.).

Pierazzo, E., Melosh, H.J., 2000a. Melt production in oblique impacts. Icarus 145, 252–261.

Pierazzo, E., Melosh, H.J., 2000b. Understanding oblique impacts from experiments, observations, and modeling. Annu. Rev. Earth Planet. Sci. 28, 141–167.

Pierazzo, E., Vickery, A.M., Melosh, H.J., 1997. A reevaluation of impact melt production. Icarus 127, 408–423.

Pierazzo, E., et al., 2008. Validation of numerical codes for impact and explosion cratering: impacts on strengthless and metal targets. Meteorit. Planet. Sci. 43 (12), 1917–1938.

Potter, R.W.K., et al., 2012. Constraining the size of the South Pole-Aitken basin impact. Icarus 220, 730–743.

Rafikov, R.R., 2003. The growth of planetary embryos: orderly, runaway, or oligarchic. Astron. J. 125, 942–961.

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

Ricard, Y., Srámek, O., Dubuffet, F., 2009. A multi-phase model of runaway core-mantle segregation in planetary embryos. Earth Planet. Sci. Lett. 284, 144–150.

Ricard, Y., Bercovici, D., Albarède, F., 2017. Thermal evolution of planetesimals during accretion. Icarus 285, 103–117.

Rivoldini, A., Van Hoolst, T., Verhoeven, O., Mocquet, A., Dehant, V., 2011. Geodesy constraints on the interior structure and composition of Mars. Icarus 213, 451–472.

Robbins, S.J., Hynek, B.M., Lillis, R.J., Bottke, W.F., 2013. Large impact crater histories of Mars: the effect of different model crater age techniques. Icarus 225, 173–184.

Roberts, J.H., Arkani-Hamed, J., 2012. Impact-induced mantle dynamics on Mars. Icarus 218, 278–289.

Roberts, J.H., Arkani-Hamed, J., 2017. Effects of basin-forming impacts on the thermal evolution and magnetic field of Mars. Earth Planet. Sci. Lett. 478, 192–202. https://doi.org/10.1016/j.epsl.2017.08.031.

Roberts, J.H., Barnouin, O.S., 2012. The effect of the Caloris impact on the mantle dynamics and volcanism of Mercury. J. Geophys. Res. Planets 117, E02007.

Roberts, J.H., Lillis, R.J., Manga, M., 2009. Giant impacts on early Mars and the cessation of the martian dynamo. J. Geophys. Res. 114, E04009. https://doi.org/10.1029/2008JE003287.

Ruedas, T., 2017. Globally smooth approximations for shock pressure decay in impacts. Icarus 289, 22–33.

Senshu, H., Kuramoto, K., Matsui, T., 2002. Thermal evolution of a growing Mars. J. Geophys. Res. 107 (E12), 5118. https://doi.org/10.1029/2001JE001819.

Shoemaker, E.M., 1962. Interpretation of Lunar Craters. Academic Press, San Diego, pp. 283–359.

Solomatov, V.S., Canup, R.M., Righter, K., 69 collaborating authors, 2000. Fluid dynamics of a terrestrial Magma Ocean. In: Origin of the Earth and Moon. University of Arizona Press, Tucson, pp. 323–338.

Thompson, S.L., 1990. ANEOS analytic equations of state for shock physics codes input manual. In: Sandia National Laboratories Report SAND89-2951.

Thompson, S., Lauson, H., 1972. Improvements in the CHART D Radiation-Hydrodynamic Code III: Revised Analytic Equations of State. Sandia National Laboratory Report SC-RR-71 0714 (113 p).

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

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

Watters, W.A., Zuber, M.T., Hager, B.H., 2009. Thermal perturbations caused by large impacts and consequences for mantle convection. J. Geophys. Res. 114, E02001.

Wünnemann, K., Collins, G.S., Melosh, H.J., 2006. A strain-based porosity model for use in hydrocode simulations of impacts and implications for transient-crater growth in porous targets. Icarus 180, 514–552.

Wünnemann, K., Collins, G.S., Osinski, G.R., 2008. Numerical modelling of impact melt production in porous rocks. Earth Planet. Sci. Lett. 269 (3-4), 530–539.
## 3 Séparation manteau/noyau à l'intérieur de la Terre Primitive

La formation du noyau terrestre est un processus hautement énergétique qui s'est déroulé pendant ou juste après les principales étapes d'accrétion de la Terre (Stevenson, 1981). Ainsi cette phase primitive de la Terre est probablement intimement liée aux processus d'impacts. Lors de la collision avec un impacteur de rayon supérieur à 100 km, l'énergie apportée par la collision peut s'ajouter à l'énergie fossile issue des radioactivités des éléments à courte période et peut générer des zones de fusion à plus ou moins grandes échelles à l'intérieur desquelles le fer de l'impacteur est dispersé sous formes de gouttelettes (Figure 8). Lors des stades primitifs de formation des planètes, de tels évènements ont pu provoquer une ségrégation locale entre le fer dense et les silicates plus légers dans la zone chauffée (Tonks and Melosh, 1992; Deguen et al., 2011). Une fois rassemblés au fond du reservoir fondu par impact, de grands volumes de fer ont ensuite pu migrer vers le centre de la Terre et contribuer à la formation du noyau terrestre par percolation (Shannon and Agee, 1998) ou par diapirisme grâce à des instabilités de type Rayleigh-Taylor (Stevenson, 1981) (Figure 8).

<figure><img src="image_48.png" /><figcaption>FIGURE 8 – Modèle de formation du noyau terrestre à partir d'un océan magmatique profond (À partir de Wood (2011)).</figcaption></figure>

L'analyse géochimique des roches du manteau terrestre montre un excès en éléments sidé-
rophiles (ex. Ni, Co, Cu) par rapport à une différenciation à l'équilibre de matériel chondritique (Becker et al., 2006; Lorand et al., 2008, 2012). Les modalités de différenciation et notamment les volumes caractéristiques de fer (gouttelettes centimétriques ou diapir pluri-kilométriques) mis en jeu après impact entre deux corps différenciés conditionnent le degré d'équilibre chimique entre le manteau et le noyau terrestre (Rubie et al., 2003; Dahl and Stevenson, 2010). En effet si les volumes caractéristiques de fer en jeu lors de la formation du noyau sont importants (tailles kilométriques), les vitesses de séparation seront très élevées et l'équilibre thermo-chimique entre la phase métallique et la phase silicatée sera partiel. En revanche, si la séparation fer/silicates se fait via des gouttelettes de fer centimétriques, les échanges entre les deux phases seront facilités et l'équilibre thermo-chimique sera atteint. Une différenciation hors-équilibre est une hypothèse émise pour expliquer l'anomalie de concentrations en éléments sidérophiles dans les roches terrestres (Kleine and Rudge, 2011). Ainsi, les processus d'accrétion et de différenciation ont pu être enregistrés dans la composition chimique et isotpique du manteau des planètes telluriques. Cette composition chimique peut refléter comment le fer a plongé vers le noyau lors de la différenciation. Les processus de séparation fer/silicates à l'intérieur de la planète gouvernent la taille caractéristique de l'écoulement du fer dans le manteau. Ainsi, dans un océan magmatique turbulent, le fer plonge sous forme de gouttelettes centimétriques avant de s'accumuler à sa base. Ensuite, de larges volumes de fer peuvent être mobilisés par instabilités kilométriques, leur taille dépendant de la rhéologie du matériel sousjacent. L'efficacité des échanges thermo-chimiques entre la phase métallique de l'impacteur et le manteau de l'impacté dépend fortement de la taille caractéristique du fer et des son temps de résidence dans le manteau. Ainsi, contraindre les processus de fragmentation siégeant dans le manteau primitif peut permettre de mieux comprendre les signatures géochimiques observées aujourd'hui.

**3.1 Fragmentation du fer de l'impacteur dans un océan magmatique**

Après un large impact, le noyau de métal de l'impacteur se disperse dans la zone fondue et plonge ensuite dans cet océan de magma moins dense sous la forme d'une pluie de fer ("Iron rain scenario"). Au cours de ce processus, deux fragmentations ont lieu : une fragmentation par impact (Kendall and Melosh, 2016) puis une fragmentation dans le réservoir fondu (Wacheul et al., 2014).

**3.1.1 Fragmentation par impact**

Lors d'un impact impliquant deux objets différenciés, l'impacteur va subir une déformation intense entraînant non seulement la fragmentation de son manteau silicaté mais aussi de son noyau métallique. Genda et al. (2017) ont modélisé les conséquences d'un impact avec un boîlide différencié (R = 1500 km) après la formation du système Terre-Lune. Pour une collision oblique statistiquement probable (avec un angle d'impact de 45°), le noyau de l'impacteur
s'allonge et se désintègre ensuite en fragments de tailles décamétriques. Certains de ces fragments frappent la Lune en orbite, mais la plupart retombent sur Terre sous forme d'impacteurs secondaires entrainant une fragmentation supplémentaire. Genda et al. (2017) montrent qu'un seul impacteur oblique peut fournir une quantité suffisante d'éléments hautement sidérophiles (HSE) aux réservoirs de silicates terrestres primordiaux au début de l'Hadéen. Le processus d'impact permet aussi d'étirer la phase métallique favorisant ainsi la dispersion du fer liquide dans un volume beaucoup plus important dans la zone fondue par impact. Kendall and Melosh (2016) ont développé des modèles d'impact SPH (voir section 2) intégrant des traceurs lagrangiens qui suivent le noyau de fer initialement intact avec un rayon de 100 km au fur et à mesure que l'impact l'étire et le disperse. La distance de déplacement finale des paires de traceurs initialement les plus proches donne une mesure de l'étirement du noyau. Les statistiques d'étirement impliquent un mélange qui sépare le noyau de fer en feuillets, en ligaments et en fragments plus petits, sur une échelle de 10 km ou moins (Kendall and Melosh, 2016).

### 3.1.2 Fragmentation dans l'océan magmatique

Les fragments de noyau dispersés par l'impact dans les réservoirs terrestres subissent en-suite un nouveau mélange par entraînement turbulent lorsque les fragments de fer en fusion traversent l'océan de magma et se déposent plus profondément dans la planète. Une fois l'impact terminé soit environ une heure après la collision, cette deuxième fragmentation va former des gouttes de fer encore plus petites (avec des rayons inférieurs au mètre). Cette fragmentation a lieu dans le réservoir fondu par impact ce qui implique un milieu avec une viscosité très faible (0.01-100 Pa.s) sur des profondeurs pouvant aller jusqu'à la profondeur de la CMB dépendamment de l'angle d'impact et de la taille de l'impacteur (Nakajima et al., 2021).

Ainsi le réservoir fondu par impact est un milieu potentiellement très turbulent. La phase métallique a pu s'y déformer puis se fragmenter en gouttes de différentes tailles entraînant des échanges thermo-chimiques plus ou moins efficaces. Les fragmentations ont pu avoir lieu en cascade jusqu'à ce que les gouttelettes atteignent une taille caractéristique qui va dépendre essentiellement de la viscosité du milieu environnant. Pour les viscosités probables des liquides silicatés, le diamètre des gouttelettes stables est estimé à environ 1 cm (Rubie et al., 2003). C'est ce qui a amené à la formulation du "iron rain scenario", scénario largement utilisé pour quantifier l'efficacité des échanges thermo-chimiques entre le fer et les silicates dans un océan magmatique (Rubie et al., 2003).

Nous avons exploré ce scenario grâce à une série d'expériences en laboratoire dans le cadre de la thèse de J.B. Wacheul (Wacheul et al., 2014) (Figure 9). Nous avons utilisé le gallium liquide et des mélanges d'eau et de glycérol respectivement comme analogues du fer et des silicate fondus. Cela nous a notamment permis d'étudier l'influence du rapport de viscosité entre le fer et les silicates sur la dynamique de fragmentation. Alors que le scénario classique de la pluie de fer considère une population de gouttes purement sphériques avec un rayon caractéristique unique qui tombent vers le fond de l'océan magmatique à une vitesse unique, nos expériences ont montré une variété de formes stables pour les gouttes de métal liquide, une
large distribution de tailles et de vitesses, et une dynamique interne intense au sein du nuage avec la superposition de nouvelles fragmentations et d'assemblage de gouttes.

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

FIGURE 9 – Haut : Modélisation analogique de la fragmentation d'un diapir de 31 mm de rayon initial. Nous représentons l'évolution temporelle (de la gauche vers la droite) pour un cas où le fluide ambiant est 500 fois plus visqueux que le diapir en fragmentation. Bas, gauche : gouttes de gallium liquide de différentes tailles tombant à travers du glycérol pur. Bas, droite : Histogramme du rayon équivalent des gouttelettes formées par la rupture d'un diapir de 30 mm de rayon initial tombant à travers un mélange de glycérol et d'eau (en bleu) et distribution Gamma la mieux ajustée (en rouge) (À partir de Wacheul et al. (2014)).

Ces résultats ont notamment servi de précurseurs aux travaux de thèse de B. Qaddah dont M. Le Bars (I.R.P.H.E. Marseille) et moi étions les encadrants. Les premiers travaux de cette thèse portaient sur la modélisation de la dynamique de fragmentation d'une goutte de fer dans un océan magmatique (Qaddah et al., 2019). Nous avons pour cela développé des simulations numériques de la dynamique de descente d'une goutte de fer liquide initialement sphérique dans une phase de silicate en fusion, jusqu'à son éventuelle fragmentation. Nous avons considéré une géométrie 2D cylindrique axisymétrique. Nous avons fait varier la viscosité des silicates fondus entre 0.05 Pa.s et 100 Pa.s et le rayon initial de la goutte de fer entre 1 mm et 350 mm. Par conséquent, nous avons considéré des nombres de Reynolds compris entre 0.027 et 85600 et des nombre de Weber compris entre 0.073 et 7480. Notre modèle numérique a permis de contraindre la morphologie, la dynamique et la stabilité de la goutte de fer en fonction des nombres de Weber et de Reynolds sans dimension ainsi que du rapport de viscosité entre les silicates fondus et la goutte de fer liquide. En particulier, nous avons montré que le
rayon maximal de la goutte stable et le nombre de Weber critique sont des fonctions monotones croissantes de la viscosité de l'océan magmatique. L'épaisseur de la couche limite de quantité de mouvement dépend principalement du rayon de la goutte et légèrement de la viscosité de l'océan magmatique. L'augmentation de la viscosité de la phase silicatée empêche les oscillations de la phase métallique et limite la surface d'échange. A l'inverse, l'augmentation du rayon initial de la goutte de fer accroît sa déformation et augmente sa surface d'échange relative. Au-dessus du nombre de Weber critique, nous avons confirmé que la fragmentation du fer liquide se produit à une distance de chute comprise entre 3.5 et 8 fois le rayon initial de la goutte dans la gamme explorée du nombre de Weber modéré.

### 3.1.3 Évolution thermo-chimique des gouttelettes de fer

Nous avons ensuite développé des modèles numériques afin de suivre les échanges de chaleur entre une goutte de fer et l'océan magmatique pour différentes tailles initiales et viscosités ambiantes (Qaddah et al., 2020) (Figure 10). Nous avons exploré des nombres de Reynolds compris entre 0.05 et 48, des rapports de viscosité compris entre 50 et 4000, des nombres de Weber compris entre 0.04 et 5, et des nombres de Peclet compris entre 70 et 850. Nous avons montré que l'échange de chaleur entre les deux phases se produit principalement dans la partie avant de la goutte. Notre étude systématique et paramétrique a montré que l'épaisseur de la couche limite thermique, la profondeur et le temps d'équilibrage, le nombre de Nusselt et le volume de l'océan magmatique affecté par les échanges thermiques suivent tous des lois de puissance en fonction du nombre de Peclet. En raison des distorsions de la goutte, ces lois d'échelle s'écartent des bilans classiques qui ne considèrent que la diffusion de la chaleur à travers une couche limite thermique laminaire. Enfin, en considérant une viscosité du fluide ambiant dépendant de la température, nous avons montré qu'une couche de faible viscosité entoure la goutte, ce qui influence l'évolution thermique des gouttes non déformables à faible nombre de Reynolds uniquement et diminue la distance de rupture pour certains modes de rupture.

Enfin dans le cadre du post-doctorat de Vincent Clési que je co-supervisais avec A. Bouhifd, nous avons caractérisé les échanges chimiques lors la descente de gouttelettes métalliques dans un océan magmatique à l'aide de modèles numériques COMSOL Multiphysics (Clesi et al., 2020) (Figure 11). En faisant varier les viscosités du silicate et du métal (entre 0.1 et 1000 Pa.s pour chaque phase) ainsi que le coefficient de partage entre le métal et le silicate ($D_{met/sil}$ variant entre 1 et 1000), nous avons dérivé des relations permettant de quantifier le degré d'équilibre entre le métal et le silicate fondus, dans un régime caractérisé par de faibles nombres de Weber ($We < 10$) et de faibles nombres de Reynolds ($10^{-3} < Re < 10^2$). Nous avons montré que les principaux paramètres contrôlant l'équilibre pour un élément sidérophile sont la viscosité de la phase silicatée et le coefficient de partage de cet élément. Nous avons appliqué notre paramétrisation pour le Ni et le Co dans le contexte de l'accrétion tardive sur Terre afin de quantifier la variation du rapport Ni/Co après un large impact en fonction de la viscosité de l'océan magmatique, pour un scénario de pluie de gouttes de fer dans un océan
<figure><img src="image_50.png" /><figcaption>FIGURE 10 – Évolution thermique d'une goutte de métal fortement déformable en fonction du temps. L'échelle de couleur représente l'anomalie de température normalisée par le contraste de température initial. La ligne noire continue représente la surface de la goutte de métal. (À partir de Qaddah et al. (2020)).</figcaption></figure>

magmatique. En considérant le modèle de Canup (2004), nous avons montré que la formation de la Lune avait un effet sur le rapport Ni/Co actuel. En fonction du rayon du noyau de Théia et de la viscosité de l'océan magmatique produit après l'impact entre la proto-Terre et Théia, la formation de la Lune pourrait être responsable de 0.45% à 3% du rapport Ni/Co actuel pour des viscosités d'océan magmatique de 0.1 à 100 Pa-s, respectivement.

## 3.2 Formation du noyau par diapirisme géant

Un impact géant peut entrainer la fusion d'importants volumes du manteau de la planète impacté formant ainsi des océan magmatiques localisés en dessous du point d'impact. Les réajustements mécaniques peuvent ensuite permettre à ces océans magmatiques de s'étaler et de
<figure><img src="image_51.png" /></figure>

FIGURE 11 – Fraction du manteau chimiquement contaminée par un élément modérément sidérophile (avec différentes valeurs de coefficient de partage $D^{met/sil}$) en fonction du rayon du noyau de l'impacteur. Dans cette figure, nous considérons que le noyau se fragmente en gouttelettes de fer avec un rayon $R_{Fe} = 1$ cm et que les viscosités du fer et des silicates fondus sont de 0.1 Pa.s. (À partir de Clesi et al. (2020)).

recouvrir l'ensemble de la surface de la Terre primitive. Ce phénomène est amplifié lorsque plusieurs impacts géants se succèdent sur des intervalles de temps suffisamment courts pour que le matériel fondu n'ait pas le temps de recristalliser. Le fer va se rassembler à la base de ces océans magmatiques (voir section précédente) car il est plus dense que les silicates fondus desquels il se différencie. Ce fer peut soit provenir du noyau de l'impacteur, soit provenir du corps impacté lui même si la séparation fer/silicates n'avait pas encore eu lieu dans la zone fondue par impact.

Ces grands volumes de fer sont instables gravitairement car ils sont eux-mêmes plus denses que le matériel sous-jacent solide et composé de silicates (si la planète est différenciée) ou d'un mélange Fer/Silicates (si la planète est encore indifférenciée). Ainsi des instabilités de type Rayleigh Taylor peuvent se former provocant la chute de diapirs de fer de tailles kilométriques vers le centre de la Terre. Honda et al. (1993) a notamment modélisé la formation du noyau terrestre en considérant une couche de fer de plusieurs centaines de kilomètres d'épaisseur recouvrant l'intégralité d'un noyau indifférencié solide. D'après leurs modèles, Honda et al. (1993) ont montré que pour former le noyau en moins d'un milliard d'années, comme l'exigent diverses contraintes géophysiques et géochimiques, la viscosité du proto-noyau silicaté et indifférencié devrait être inférieure à $10^{26}$ Pa.s, ce qui est plus grand que la viscosité actuelle du
manteau, mais beaucoup plus petit que la viscosité estimée pour le proto-noyau silicaté à basse température et à haute pression. Cette étude a notamment mis en évidence que le chauffage du proto-noyau par l'énergie gravitationnelle libérée pendant la formation du noyau pouvait contribuer à l'accélération du processus de formation du noyau. Cependant il semble peu probable qu'une couche de fer aussi épaisse puisse se former et se stabiliser avant de plonger vers le centre de la Terre.

Senshu et al. (2002) a développé des modèles numériques 1D permettant de simuler l'histoire thermique au cours de l'accrétion de Mars à la manière de Kaula (1979). Il a cependant rajouté la contribution de la séparation métal-silicate dans cette planète en croissance. Dans leurs modèle, la perturbation thermique causée par les impacts de planétésimaux est calculée pour chaque événement d'impact et un diapir métallique est considéré à la base d'un bassin de magma si le site d'impact est chauffé au-dessus de la température de fusion. Le plongeon des diapirs métalliques est retracé en supposant une vitesse de Stokes. Leurs résultats numériques suggèrent que sur Mars, (1) aucun océan magmatique global n'est formé pendant l'accrétion, (2) la séparation métal-silicate a lieu sans fusion à l'échelle globale, et (3) au lieu d'un noyau métallique, une couche riche en métal se forme à la fin du processus d'accrétion. Ce modèle de formation de noyau planétaire a ensuite été développé en 3D par Golabek et al. (2009) afin de mettre en évidence l'influence des variations de viscosité au cours de ce processus.

Ces travaux m'ont amené à m'intéresser pendant ma thèse au plongeon d'un diapir de fer collecté dans la zone chauffée après un large impact météoritique (Figure 12). Pour étudier l'évolution d'un tel diapir métallique, nous avons proposé un modèle analytique de réajustement thermochimique que nous avons comparé à des simulations numériques en géométrie sphérique axisymétrique et à viscosité variable (Monteux et al., 2009). Pour cela, j'ai développé un code numérique permettant de résoudre les équations de Navier-Stokes en volumes finis et en formulation en fonctions de courant. Par rapport au modèle précédent (Monteux et al., 2007), ce modèle se place en géométrie axisymétrique sphérique dans l'approximation de Boussinesq étendu et permet de calculer la dissipation visqueuse ainsi que le chauffage par compression/décompression adiabatique. Ce nouveau modèle numérique permet non seulement de caractériser la dynamique de rééquilibrage thermique mais aussi l'évolution de la composition chimique en fonction du temps. Afin de pouvoir traiter le transport d'un champ de composition chimique, il est nécessaire de pouvoir résoudre convenablement (i.e. sans oscillation numérique) un champ purement advectif. Pour cela, mon code numérique inclut des schémas T.V.D. (Total Variation Diminishing) (Roe, 1986) parmi lesquels nous utilisons le schéma Superbee qui est non seulement non-oscillant mais permet aussi de limiter les flux diffusifs numériques.

Grâce à ces modèles, nous avons montré que la phase métallique s'enfonce avec une vitesse de l'ordre de la vitesse de Stokes comme considéré par Senshu et al. (2002). Autre point important que nous avons contraint : l'énergie thermique libérée par la ségrégation du métal
<figure><img src="image_52.png" /></figure>

FIGURE 12 – Température adimensionnée (demi-cercle gauche) et composition chimique (demi-cercle droite) (pour le fer $f = 1$, pour le matériel indifférencié $f = 0.17$, pour les silicates $f = 0$) au temps $t = 0$ (a), $t = 1.4$ Ma (b), $t = 3.8$ Ma (c) et $t = 546$ Ma (d) (dans un cas où la viscosité est uniforme, $R = 4000$ km, $R_{imp} = 600$ km et avec une résolution de $200 \times 200$ points de grille) (À partir de Monteux et al. (2009)).

est comparable à l'énergie thermique dissipée lors de l'impact. Comme nous l'avons vu précédemment l'énergie de l'impact est distribuée potentiellement dans un grand volume de matériel indifférencié (Figure 4). L'énergie potentielle dissipée lors de la chute d'un important diapir de fer est quant à elle potentiellement libérée dans un volume beaucoup plus petit (le diapir de fer et ses environs proches). Un échauffement important du diapir métallique peut se produire lors de la chute, augmentant son excès de température d'un facteur 2 ou 3 au maximum. Lorsque l'on considère une viscosité de la phase métallique qui diminue avec la température, la proportion d'énergie thermique transférée à la matière indifférenciée environnante augmente au détriment de l'énergie dissipée dans la phase métallique. Ceci a pour effet de permettre d'une
part de former un proto-noyau à une température proche de celle de la zone d'impact et d'autre part de provoquer de la fusion et donc de la séparation Fer/silicates supplémentaire dans le manteau indifférencié environnant.

Ce processus pourrait ainsi être à l'origine d'une séparation Fer/silicates en cascade dans une planète initialement indifférenciée. En effet l'énergie potentielle dissipée en chaleur lors de la chute d'un importants volume de fer peut entraîner des phénomènes de fusion supplémentaires qui contribuent à séparer plus de fer et dissipent ainsi encore de l'énergie potentielle etc... Šrámek et al. (2010) ont développé des modèles bi-phasiques de simulations de la différenciation noyau-manteau montrant que le premier impact qui fait fondre la phase métallique près de la surface est potentiellement capable de déclencher toute la ségrégation noyau-manteau dans un phénomène d'emballement. Ils ont contraint le seuil de cette instabilité en fonction de la taille de l'impacteur et de la planète, ainsi que de la température initiale de la planète. La ségrégation du métal se produit par un mécanisme qui est intermédiaire entre l'instabilité diapir habituelle et une onde de porosité. Šrámek et al. (2010) ont contraint le rôle des paramètres les plus importants, tels que les viscosités du solide et du métal ou la dépendance en profondeur de la gravité.

Références personnelles

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

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., Volume 403, Pages 236-245. LIEN HAL

Bouhifd 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. Chemical Geology vol.461, p.128 - 139. LIEN HAL

Qaddah B., Monteux, J., Clesi V., Bouhifd A., Le Bars M., 2019, Dynamics and stability of an iron drop falling in a magma ocean. Physics of the Earth and Planetary Interiors. vol.289, p.75-89. LIEN HAL

Clesi V., Monteux, J., Qaddah B., Le Bars M., Wacheul J.B., Bouhifd A., 2020, Dynamics of core-mantle separation : Influence of viscosity contrast and metal/silicate partition coefficients on the chemical equilibrium. Physics of the Earth and Planetary Interiors, vol.306, p.106547. LIEN HAL

Qaddah B., Monteux, J., Le Bars M., 2020, Thermal evolution of a metal drop falling in a less dense, more viscous fluid. Physical Review Fluids vol.5, p.053801. LIEN HAL
# A model of metal-silicate separation on growing planets

J. Monteux*, Y. Ricard, N. Coltice, F. Dubuffet, M. Ulvrova

Université de Lyon, Lyon, F-69003, France

Université Lyon 1, Lyon, F-69003, France

Ecole Normale Supérieure de Lyon, Lyon, F-69364, France

CNRS, UMR5570, Laboratoire de Sciences de la Terre, Villeurbanne, F-69622, France

## ARTICLE INFO

Article history:

Received 26 March 2009

Received in revised form 7 August 2009

Accepted 11 August 2009

Available online xxxx

Editor: T. Spohn

Keywords:

core formation

meteorological impacts

early earth

numerical modeling

differentiation

## ABSTRACT

The thermal evolution of planets during their accretionary growth is strongly influenced by impact heating. The temperature increase following a collision takes place mostly below the impact location in a volume a few times larger than that of the impactor. Impact heating depends essentially on the radius of the impacted planet. When this radius exceeds ~1000 km, the metal phase melts and forms a shallow and dense pool that penetrates the deep mantle as a diapir. To study the evolution of a metal diapir we propose a model of thermo-chemical readjustment that we compare to numerical simulations in axisymmetric spherical geometry and with variable viscosity. We show that the metallic phase sinks with a velocity of order of a Stokes velocity. The thermal energy released by the segregation of metal is smaller but comparable to the thermal energy buried during the impact. However as the latter is distributed in a large undifferentiated volume and the former potentially liberated into a much smaller volume (the diapir and its close surroundings) a significant heating of the metal can occur raising its temperature excess by at most a factor of 2 or 3. When the viscosity of the hot differentiated material decreases, the proportion of thermal energy transferred to the undifferentiated material increases and a protocore is formed at a temperature close to that of the impact zone.

© 2009 Elsevier B.V. All rights reserved.

## 1. Introduction

Core formation is the most important differentiation event that occurred during Earth's history. Metal/silicates separation is a rapid event (<60 My) (Kleine et al., 2002; Touboul et al., 2007; Yin et al., 2002) contemporaneous with Earth accretion and involving gravitational mechanisms such as percolation, negative diapirism and Rayleigh-Taylor instabilities (Honda et al., 1993; Stevenson, 1990). In the homogeneous accretion hypothesis, metal segregation and thereby core formation need significant heating to exceed the melting temperature of iron alloys or of silicates. During the early stages of planetesimals formation, heating by decay of short lived radionuclides is a potential energy source to enhance early differentiation (Yoshino et al., 2003). As a planetesimal grows, its gravity increases and it will increasingly attract the other surrounding planetesimals. The dissipation of the kinetic energy of the impacts provides a later shallow source of heat.

Impacts of large planetesimals have strongly influenced the late accretionary and thermal state of nearly fully-formed planetary bodies (Senshu et al., 2002; Tonks and Melosh, 1992). During an impact, when the relative velocity between a planet and an impactor

overcomes the seismic velocities of the medium, a shock wave develops. The shock pressure is nearly uniform in a spherical region next to the impact (the isobaric core), and strongly decays away from it (Croft, 1982; Pierazzo et al., 1997). In this isobaric core, the kinetic energy of the impact is dissipated and leaves a temperature anomaly of several hundred degrees on Moon to Mars size bodies (Monteux et al., 2007; Senshu et al., 2002). The temperatures reached are mostly related to the properties (density and radius) of the impacted body, and only weakly to those of the impactor (Monteux et al., 2007). The melting temperature of iron alloys is lower than the silicates solidus (Agee, 1997; Fei et al., 1997; Ghosh and McSween, 1998). On large impacted planets, a local differentiation may occur between heavy metal and light silicates in the heated anomaly (Tonks and Melosh, 1992). Hence, a thermo-chemical readjustment follows, associated with the sinking of the metallic component toward the center of the impacted protoplanet (Fig. 1).

For large planets, gravitational energy release due to core formation can induce melting of the whole planet (Stevenson, 1989; Ricard et al., 2009). This subsequent melting depends on the mechanisms of the metal descent (Golabek et al., 2008; Samuel and Tackley, 2008). The aim of this study is to determine the thermal evolution of metal during descent and the thermal state of the core.

First, we propose analytical and numerical isoviscous models of segregation of a purely spherical iron diapir. As the viscosity contrast between molten metal and undifferentiated cold material can reach

<footnote>* Corresponding author. Université de Lyon, Lyon, F-69003, France.
E-mail address: jmonteux@eos.ubc.ca (J. Monteux).</footnote>
<figure><img src="image_53.png" /><figcaption>**Fig. 1.** Schematic view of the chemical equilibration following a large impact on an undifferentiated protoplanet. In the isobaric core resulting from the dissipation of the shock wave (a,b), the temperature increase melts the metal (c) that segregates rapidly (d), then sinks toward the planetary embryo center by a diapiric instability (e).</figcaption></figure>

several orders of magnitude, we then focus on more realistic models of segregation of metal after a large impact with temperature dependent rheologies. We show that the size of impactors and viscosities involved largely determine the inner thermal state of a young planet.

## 2. Thermo-chemical state after large impact

### 2.1. Thermal state

After a meteorological impact, heating is localized in a spherical region called the isobaric core just beneath the impact site. The radius of the isobaric core $R_{ic}$ is comparable to the radius of the impactor $R_{imp}$ and depends on energy conversion during the shock. With a minimal set of assumptions, we get $R_{ic} = 3^{1/3}R_{imp}$ following (Senshu et al., 2002) and (Pierazzo et al., 1997). Just after the adiabatic pressure release, the isobaric core is isothermal and we call $\Delta T_0$ the shock induced temperature increase. The lower script 0 indicates that we consider this instant as the origin of our time variable. Outside the isobaric core, the temperature anomaly decays as $\Delta T_0(r) = \Delta T_0(R_{ic}/r)^m$ with $m \sim 4.4$ as proposed by (Senshu et al., 2002). Assuming that the kinetic energy of the impactor is controlled by the escape velocity of the impacted body and that impactor and impacted body have the same densities (i.e., $\rho_{ic} = \rho_{imp} \equiv \rho_0$), a simple energy balance (Monteux et al., 2007), indicates that

$$ \Delta T_0 = \frac{4\pi \gamma \rho_0^2 GR^2}{9 h(m) \rho C_p}, \qquad (1) $$

where $\overline{\rho C_p}$ is the average heat capacity of the impacted body that is plausibly a mixture of silicate and metal, G is the gravitational constant, $\rho_0$ is the density of the undifferentiated material, R is the radius of the impacted planet and where the function $h(m)$ represents the volume effectively heated normalized by the volume of the isobaric core (typically $h(m) \sim 2-3$ (Monteux et al., 2007)). The empirical coefficient $\gamma$ is the fraction of the kinetic energy of the impactor dissipated as heat. From shock experiments, $\gamma$ ranges between 0.2 and 0.4 depending on material properties and shock velocities (O'Keefe and Ahrens, 1977) (i.e., 20 to 40% of the kinetic energy is buried at depth, the rest rapidly radiated away during or shortly after the impact). The shock-induced temperature excess, $\Delta T_0$, strongly increases with the radius of the impacted body. According to the set of parameters of Table 1, $\Delta T_0(K) = 4.7 \times 10^{-5} R^2(\text{km})$; for a Moon size body $\Delta T_0$ is 140 K while it is 1925 K for an Earth size body.

The thermal state of a protoplanet before an impact depends on its growth history and on its initial heating caused by short lived

<table><caption>**Table 1**

Typical parameter values for numerical models.</caption><tr><td>Planet radius</td><td>R</td><td>1000–4000 km</td></tr><tr><td>Impactor radius</td><td>R<sub>imp</sub></td><td>100–400 km</td></tr><tr><td>Silicate density</td><td>ρ<sub>Si</sub></td><td>3500 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>Density difference</td><td>Δρ<sub>0</sub> = ρ<sub>Fe</sub> - ρ<sub>Si</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>Heat capacity</td><td>ρ<sub>Cp</sub></td><td>4 × 10<sup>3</sup> kJ K<sup>-1</sup> m<sup>-3</sup></td></tr><tr><td>Heat diffusivity</td><td>κ</td><td>10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Thermal conductivity</td><td>k</td><td>4 W m<sup>-1</sup> K<sup>-1</sup></td></tr><tr><td>Initial temperature</td><td>T<sub>0</sub></td><td>K</td></tr><tr><td>Metal content</td><td>f<sub>0</sub></td><td>0.17</td></tr><tr><td>Viscosity</td><td>η<sub>0</sub></td><td>10<sup>22</sup> Pa s</td></tr><tr><td>Viscosity factor</td><td>λ</td><td>2.5 × 10<sup>-2</sup> – 1</td></tr><tr><td>Gravity</td><td>g<sub>0</sub> = 4πGρ<sub>0</sub>R<sup>3</sup></td><td>m s<sup>-2</sup></td></tr><tr><td>Stokes velocity scale</td><td>Δρ<sub>0</sub>g<sub>0</sub>R<sup>2</sup>/η<sub>0</sub></td><td>~ 100 m yr<sup>-1</sup></td></tr><tr><td>Time scale</td><td>η<sub>0</sub>/Δρ<sub>0</sub>g<sub>0</sub>R</td><td>~ 20 kyr</td></tr><tr><td>Rayleigh number Ra<sub>χ</sub></td><td>ρ<sub>Cp</sub>Δρ<sub>0</sub>g<sub>0</sub>R<sup>3</sup>/η<sub>0</sub>k</td><td>~ 10<sup>8</sup></td></tr><tr><td>Buoyancy B</td><td>Δρ<sub>0</sub>/ρ<sub>0</sub>g<sub>0</sub>T<sub>0</sub></td><td>25–250</td></tr><tr><td>Dissipation number D<sub>χ</sub></td><td>Δρ<sub>0</sub>g<sub>0</sub>R / ρ<sub>Cp</sub>ΔT<sub>0</sub></td><td>36.6</td></tr><tr><td>Impact energy conversion coefficient</td><td>γ</td><td>0.3</td></tr><tr><td>Volume effectively heated by impact</td><td>h(m)</td><td>2.7</td></tr><tr><td>Stokes velocity coefficient</td><td>c<sub>1</sub></td><td>0.1–0.2</td></tr><tr><td>Heat diffusion coefficient</td><td>c<sub>2</sub></td><td>0.3–1.05</td></tr></table>

radionuclides like <sup>26</sup>Al and <sup>60</sup>Fe. This early radioactive heating can eventually cause melting and differentiation of planetesimals that have quickly grown (Yoshino et al., 2003). The impact heating superimposed to a sufficiently hot protoplanetary interior can trigger melting of the Fe-FeS system (the eutectic temperature is close to 1250 K at 1 bar) (Fei et al., 1997) and potentially of silicates (solidus temperature is around 1500 K at 1 bar) (Agee, 1997). In these cases, a fraction of the thermal energy is converted to latent heat during the phase transformations.

### 2.2. Compositional state

An impact on a large enough undifferentiated protoplanet composed of a mixture of metal alloys and silicates can trigger phase transformations and initiate differentiation. The first component that melts is the metal phase. In the region where metal melting occurs, the liquid metal can percolate through the solid silicate matrix. Percolation is only possible for small dihedral angles (<60°) or for large melt volume fraction above a percolation threshold. The dihedral angle of liquid iron alloy within silicates is large (~100°) in the upper mantle but decreases with increasing pressure (Shannon and Agee, 1996). However, the volume fraction of liquid alloy is typically larger than 10% if melting is complete, which overcomes the percolation threshold (Von Bargen and Waff, 1986). On Earth the core represents 17% of the volume of the planet, Mars has likely a slightly smaller core but Mercury's core is 43% of the planet. The metal is collected at the bottom boundary of the melted zone forming a diapir that ultimately sinks within the interior of the impacted protoplanet (Ricard et al., 2009).

If the temperature exceeds the silicate solidus and eventually the liquidus, the separation of metal and silicates can occur as a metal rainfall through a turbulent magma (Höink et al., 2005; Stevenson, 1990). Small droplets of heavy metal sediment at the bottom of the melted region. This scenario may not be the generic one, as it would imply that a planet embryo maintains a melted metal component without differentiating until the silicates start melting. It has been suggested that the metal may segregate per percolation, as soon as it melts, while the silicates are still mostly solid (Ricard et al., 2009). Locally, however, the impact of an undifferentiated planetesimal on an already differentiated large planetary embryo, may of course, be energetic enough to melt (or even vaporize) the silicate and metal contents of the impactor and the silicates of the impacted body inside the isobaric core.
The two processes (percolation or metal “rain”) lead to a local differentiation within the melted region between light silicates and heavy metals on a short timescale compared to that of the slow viscous deformation (Tonks and Melosh, 1992). The melted region is as large as or a few times larger than the isobaric core (Pierazzo et al., 1997). Here, we identify the initially differentiated zone to the isobaric core, metal being overlaid by pure silicates shortly after the impact (see Fig. 1).

### 3. Dynamic model of differentiation

The setting described in the previous section is gravitationally unstable and the metal phase sinks toward the center of the impacted planetesimal while the silicates (lighter than undifferentiated material) spread underneath the surface. To study the global dynamics of this differentiation event, we develop a thermo-mechanical model in spherical axisymmetric geometry, of viscous flow with three chemical components. Using a viscous and linear rheology during the segregation of the core is clearly a large approximation. The large deviatoric stress generated by the metallic diapirs should lead to a non-linear rheology (Samuel and Tackley, 2008), elasto-plastic deformations (Gerya and Yuen, 2007) or even to hydrofracturation if they exceed the ultimate strength of rocks which is ~1–2 GPa (Davies, 1982). Pressure dependence of the rheology can also influence the metal sinking time but is not considered here since we focus on small growing planets. During the early stages of accretion, the interior of the growing planets may have been colder or hotter than the outer layers depending on the ratio of radioactive and impact heating and on the history of accretion. For simplicity, we assume in our models a homogenous temperature on the growing planet before the impact.

#### 3.1. Physical model

Sinking occurs under the action of gravity in a spherical homogeneous protoplanetary body. We neglect for simplicity the changes of gravity during the differentiation. Hence gravitational acceleration $g(r)$ increases linearly with radius $r$:

$$g(r) = \frac{4}{3} G \rho_0 r = g_0 \frac{r}{R}, \qquad (2)$$

where $g_0$ the surface gravity. The density of undifferentiated material is $\rho_0=f_0\rho_{Fe}+(1-f_0)\rho_{Si}$ where $f_0$ is the volume fraction of metal and $\rho_{Fe}$, $\rho_{Si}$, the densities of the metallic phase and the pure silicates, respectively (see typical numerical values in Table 1.)

The dynamics of segregation potentially involves a series of multiscale physical processes, especially to take the effects of melting into account and a realistic multiphase dynamics (Golabek et al., 2008; Ricard et al., 2009). No numerical models can handle simultaneously all these complexities and as a consequence, we follow the approach of (Samuel and Tackley, 2008) and consider a thermo-chemical system with infinite Prandtl limit, with no possible subsequent phase separation within the undifferentiated material except that caused by the impact (e.g., the volumes of pure metal and pure silicates remain constant during the simulations and equal to 17% and 83% of the initial isobaric core).

The necessary approximations are somewhat different from the classic treatment of thermal convection (Ricard, 2007). We non-dimensionalize the lengths by the planetary radius $R$, the velocities by a Stokes velocity $\Delta\rho_0 g_0 R^2/\eta_0$ (where $\Delta\rho_0=\rho_{Fe}-\rho_{Si}$ and $\eta_0$ is the reference viscosity of cold material far from the impact site), the temperature by $\Delta T_0$ (see Eq. (1)). The governing mechanical non-dimensional equations are the conservation of mass

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

and the conservation of momentum

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

where $\mathbf{v}$, $P$, $T$ and $r$ are the non-dimensional velocity, pressure, temperature and radius, $\eta$ the viscosity and $\mathbf{e}_r$ the radial unit vector. The buoyancy ratio $B$ (Christensen and Yuen, 1985) is:

$$B = \frac{\Delta\rho_0}{\rho_0 \alpha \Delta T_0}. \qquad (5)$$

The downward buoyancy force that drives the flow increases with the volume fraction of metal $f$ that varies between 0 (pure silicates) and 1 (pure metal), 0.17 being that of undifferentiated material. A depth dependent and constant in time gravity has been used in Eq. (4) although, in principle, gravity should have been computed self-consistently from the time-dependent density distribution. We assume a temperature dependent viscosity such as $\eta=\eta_0\lambda^T$ with $\lambda$ being the viscosity factor (lower than 1) which is equivalent to the viscosity ratio between the hottest and coldest material at the start of the experiment. Such a viscosity decreases sharply with temperature and is simpler to implement than the usual Arrhenius law (Ratcliff et al., 1997; Ziethe and Spohn, 2007).

The conservation of energy writes

$$\frac{\mathrm{D}T}{\mathrm{D}t} = \frac{\nabla^2 T}{Ra_\chi} + D_\chi \frac{\eta}{\eta_0} \Phi + \frac{1}{B} \frac{\Delta\rho}{\rho} \frac{T}{\Delta T_0} D_\chi \frac{\mathrm{D}P}{\mathrm{D}t}. \qquad (6)$$

The importance of diffusion is controlled by the compositional Rayleigh number $Ra_\chi$,

$$Ra_\chi = \frac{\Delta\rho_0 g_0 R^3}{\kappa \eta_0}, \qquad (7)$$

the chemical dissipation number is

$$D_\chi = \frac{\Delta\rho_0 g_0 R}{\overline{\rho C_p} \Delta T_0}, \qquad (8)$$

considering for simplicity that $\overline{\rho C_p} = \rho_{Fe}C_p^{Fe} = \rho_{Si}C_p^{Si}$ (truly, see Table 1, $\rho_{Fe}C_p^{Fe}=4 \times 10^3 \mathrm{kJ} \mathrm{K}^{-1} \mathrm{m}^{-3}$, $\rho_{Si}C_p^{Si}=3.85 \times 10^3 \mathrm{kJ} \mathrm{K}^{-1} \mathrm{m}^{-3}$, and we use $\overline{\rho C_p} = 4 \times 10^3 \mathrm{kJ} \mathrm{K}^{-1} \mathrm{m}^{-3}$). As $g_0$ is proportional to $R$ and $\Delta T_0$ to $R^2$, see Eq. (1), the chemical dissipation is independent of the planet radius and amounts to 36.6 (see Table 1).

An important energy source is provided by the dimensionless dissipation function $\Phi$ that expresses the conversion of potential energy into heat

$$\Phi = 2 \varepsilon : \varepsilon. \qquad (9)$$

where $\varepsilon$ is the dimensionless strain rate tensor. For simplicity, we make the approximation that the thermal conductivities of the metal, silicates and undifferentiated materials are the same (truly $k_{Fe}=10 \mathrm{W} \mathrm{m}^{-1} \mathrm{K}^{-1}>k_{Si}=3 \mathrm{W} \mathrm{m}^{-1} \mathrm{K}^{-1}$).

The metal volume fraction is then simply advected by the flow,

$$\frac{Df}{Dt} = 0. \qquad (10)$$

#### 3.2. Model approximations

The equations of momentum and energy conservations, Eqs. (4) and (6), are similar to those classically used for mantle convection simulations but a number of differences should be discussed. As the buoyancy number $B$ is very large (the density difference between
metal and silicates is 40 to 620 times larger than the thermal density variations), the thermal buoyancy $T/B$ can be safely neglected in the momentum equation.

Neglecting the terms in $1/B$ implies to omit the adiabatic heat transfer (the term in $(DP/Dt)$) in Eq. (6) but to keep the dissipation term $D_\chi(\eta/\eta_0)\Phi$. The differentiation of the planet liberates a large amount of potential energy converted into heat by the dissipation term but the adiabatic heating remains small. This is very different from the typical convection situation in which there is no time variation of the potential energy, and where the dissipation is on average balanced by the work due to compression and expansion over the convective cycle (Hewitt et al., 1975).

### 3.3. Numerical model

We implement a finite volume numerical model to solve Eqs. (3), (4), (6) and (10) in axi-symmetric spherical geometry. We use a stream function formulation for the equations of motion with a direct implicit inversion method (Schubert et al., 2001). Eqs. (6) and (10) are solved by an Alternating Direction Implicit (ADI) scheme (Douglas, 1955; Peaceman and Rachford, 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 (Sramek, 2007) which enables a high resolution of pure advection fields. We use at least 200 × 200 grid points. Velocity 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. We benchmark the viscous flow solver with variable viscosity and the transport scheme against several analytical solutions (Monteux, 2009).

## 4. Thermal evolution of sinking metallic diapir: analytical considerations

Before showing the results of complex numerical simulations with temperature dependent rheologies, we develop a simple model describing the thermal evolution of the sinking metal diapir, by approximating the metal diapir by a spherical drop falling into undifferentiated medium of uniform viscosity with a Stokes-like velocity. The radius of the metal drop $R_{Fe}$ can be related to the radius $R_{ic}$ of the volume initially differentiated after impact heating by $R_{Fe}^3 = f_0 R_{ic}^3$ and to the radius of the impactor by $R_{Fe}^3 = 3f_0 R_{imp}^3$.

### 4.1. Sinking velocity

The velocity $V$ of the metallic diapir in an undifferentiated medium is comparable to the Stokes velocity of a sphere of similar volume. The density difference between the metal and the undifferentiated material is a function of temperature and composition but the temperature contribution is minor. Hence, we consider $\Delta\rho = (1 - f_0)\Delta\rho_0$. Because gravity is a linear function of depth, the velocity of the sphere decreases during sinking as

$$V = \frac{dr}{dt} = -c_1(1 - f_0) \frac{\Delta\rho_0 g_0 R_{Fe}^2 r}{\eta_S R}. \quad (11)$$

In Eq. (11), the dimensionless constant $c_1$ depends on the geometry of the system and on the viscosity contrast between the falling sphere and the surrounding medium.

The viscosity of the surrounding undifferentiated material $\eta_S$ controls the sinking velocity. In the case of a sphere sinking in an infinite medium, the coefficient $c_1$ is given by the Hadamard-Rybczynski equation and varies from $4/15 = 0.27$ (isoviscous) to $1/3 = 0.33$ for an inviscid sphere

(Hadamard, 1911; Rybczynski, 1911). In the situation described in this paper, the boundary conditions are applied at a finite distance (the planetary surface) and the Hadamard and Rybczynski equation is thus only an approximation (Honda et al., 1993; Samuel and Tackley, 2008). The exact value of the constant $c_1$ will be obtained later through numerical experiments.

The position of the metallic drop obtained by solving Eq. (11) varies from an initial position $r_0$ ($r_0 = R - R_{Fe} \sim R$) as

$$r(t) = r_0 \exp\left(-\frac{t}{\tau_S}\right), \quad (12)$$

with a characteristic time equal to

$$\tau_S = \frac{\eta_S R}{c_1 \Delta\rho_0 g_0 \frac{1}{(1-f_0)R_{Fe}^2}}. \quad (13)$$

As $g_0$ is proportional to the planetary radius $R$ (Eq. (2)), the time $\tau_S$ is independent of the planetary radius but depends only on the diapir size $R_{Fe}$. Of course, no segregation occurs, i.e., $\tau_S \to +\infty$, for a planet of pure silicates ($f_0 = 0$ which means $R_{Fe} = 0$) or of pure metal ($f_0 = 1$). This characteristic sinking time is strongly dependent of the viscosity of the surrounding undifferentiated material which is poorly constrained. With the typical values of Table 1, this time can be computed from the size $R_{imp}$ of the impactor and we find $\tau_S(\text{kyr}) = 2.7 \times 10^9 (\eta_S/\eta_0) R_{imp}^{-2} (\text{km})$.

### 4.2. Global energy conversion

As we assume that gravity remains constant with time (albeit non-uniform), the energy equation Eq. (6) integrated over the whole planet with the use of the momentum equation Eq. (4) and neglecting the adiabatic decompression of the planet during the core segregation is simply

$$\frac{d}{dt}(\Delta E_p + \Delta E_T) = F, \quad (14)$$

where the total potential and thermal energies changes are

$$\Delta E_p = \int_{\Omega} \frac{1}{2} [\rho(\mathbf{r},t) - \rho(\mathbf{r},0)] g_0 \frac{r^2}{R} dV, \quad (15)$$

($\Omega$ is the planetary volume),

$$\Delta E_T = \int_{\Omega} \overline{\rho c_p [T(\mathbf{r},t) - T(\mathbf{r},0)]} dV, \quad (16)$$

and the heat flux $F$ is

$$F = \int_{\Sigma} k \frac{\partial T}{\partial r} dS, \quad (17)$$

($\Sigma$ is the planetary surface).

As we neglect the term in $1/B$ in the energy equation Eq. (4), the Eq. (14) misses the energy variation $\Delta E_a$ due to the changes in pressure (the subscript $a$ means that this term is related to changes in adiabatic compression)

$$\frac{d\Delta E_a}{dt} = \int_{\Omega} \alpha T \frac{\partial P}{\partial t} dV \sim \alpha T_0 \frac{d}{dt} \int_{\Omega} [P(\mathbf{r},t) - P(\mathbf{r},0)] dV \quad (18)$$

where the last approximation assumes that the temperature remains close to $T_0$. The difference of pressure between a homogeneous and a differentiated planet is easy to compute analytically and is of order $\alpha T_0 \Delta E_p$, i.e., a few percent of the changes in potential energy. This confirms that the energy change due to pressure changes is a minor effect.
### 4.3. Maximum temperature

The maximum temperature that the sinking metal can reach can be estimated by assuming that the whole variation of potential energy is only used to heat up the metal, without any heat transfer to the surrounding material.

Let us consider a melted zone of radius $R_{ic}$ underneath and tangent to the planetary surface that differentiates ultimately forming a metallic core of volume $V_{Fe}$ and radius $R_{Fe}$ (with $R_{Fe}^3 = f_0 R_{ic}^3$) and a silicate layer of volume $V_{Si}$ within a shell surrounding the whole planet with inner shell radius $R_S$ and outer shell radius $R$ i.e., $R_S^3 = R^3 - (1 - f_0) R_{ic}^3$. The change of potential energy is according to Eq. (15) (Flasar and Birch, 1973):

$$ \Delta E_p = \frac{2\pi}{5R} g_0 ((\rho_{Fe} - \rho_0) R_{Fe}^5 + (\rho_{Si} - \rho_0) (R^5 - R_S^5)). \quad (19) $$

Assuming $R_{ic} \ll R$, a Taylor expansion of Eq. (19) leads to

$$ \begin{aligned} \Delta E_p &\sim -\frac{1}{2} \Delta \rho_0 g_0 R f_0 (1 - f_0) V_{ic} = -\frac{1}{2} (\rho_{Fe} - \rho_0) g_0 R V_{Fe} \\ &= -\frac{1}{2} (\rho_0 - \rho_{Si}) g_0 R V_{Si}, \end{aligned} \quad (20) $$

where $V_{ic}$ is the volume of the isobaric core. The change of potential energy is thus equivalent to that released by the sinking of the isobaric volume $V_{ic}$ and excess density $f_0(1 - f_0)\Delta\rho_0$. Alternatively it corresponds to the energy released by a metal sphere of volume $V_{Fe}$ sinking, or of a silicate sphere rising, through undifferentiated material. If only the metal heats up, the change of thermal energy according to Eq. (16) is $\Delta E_T = \overline{\rho C_p f_0 \Delta\Theta} V_{ic}$ where $\Delta\Theta$ is the temperature increase (just after the impact, the metal temperature is $T_0 + \Delta T_0$, then it reaches at most $T_0 + \Delta T_0 + \Delta\Theta$). A scaling value for the temperature increase during segregation is thus

$$ \Delta\Theta = \frac{1 - f_0}{2} \frac{1}{\rho C_p} \Delta\rho_0 g_0 R. \quad (21) $$

As $g_0$ is proportional to $R$ (Eq. (2)), the core segregation can increase the temperature by a quantity proportional to $R^2$ (Flasar and Birch, 1973; Ricard et al., 2009). The ratio of $\Delta\Theta$ to the post impact temperature $\Delta T_0$ is, according to Table 1 and Eq. (1),

$$ \frac{\Delta\Theta}{\Delta T_0} = \frac{3h(m)}{2\gamma} (1 - f_0) \frac{\Delta\rho_0}{\rho_0} \sim 11.8 \quad (22) $$

or, $\Delta\Theta(K) = 11.8\Delta T_0(K) = 5.6 \times 10^{-4} R^2(\text{km})$ which rapidly becomes a large quantity as $R$ increases. Of course, in a real situation not all energy will remain within the metal, and we will see that, when the metal diapir is too small, the metal can even cool off rather than warm up during its motion.

### 4.4. Thermal regime of the metallic sphere

While the hot metallic sphere is sinking, it warms up by shear heating but it also cools down by diffusion. In the reference frame of the sinking drop, the conservation of energy integrated over the volume $V_{Fe}$ of the metallic drop (or through its surface $S_{Fe}$) indicates that

$$ \overline{\rho C_p} V_{Fe} \frac{d\Delta T}{dt} = -k \frac{\Delta T}{\delta} S_{Fe} + \tau_{\perp} : \nabla \mathbf{v} V_{Fe}, \quad (23) $$

where we assume that the temperature and the dissipation are at first order uniform in the metal. The difference $\Delta T$ is the difference between the diapir and the undifferentiated material. We assume that $\Delta T = T - T_0$, i.e., that the hot diapir sinks into a medium that keeps its initial temperature outside the boundary thickness $\delta$. Even when the

diapir viscosity is low and when the dissipation occurs significantly outside it, our numerical simulations show that the maximum temperature is reached inside the diapir.

The thickness $\delta$ over which the temperature diffuses should be written as $R_{Fe}$ times a dimensionless function $c_2$ of the various parameters of the problem. The thickness of the diffusive boundary layer, $c_2$, should decrease with the sinking velocity of the diapir (i.e., with the Peclet number $VR_{Fe}/\kappa$) as a power law with exponent $-1/2$ or $-1/3$, depending on the viscosity ratio between the metal and the undifferentiated material (see e.g., Ribe, 2007). We can also write the dissipation $\tau$: $\nabla\mathbf{v} = \eta_e V^2 / R_{Fe}^3$ where $\eta_e$ is the effective viscosity of the region where dissipation occurs. In this case, Eq. (23) using the expressions of the time dependent position, Eq. (12), and of the maximum temperature increase, Eq. (21), can be recast as

$$ \frac{d\Delta T}{dt} = -\frac{\Delta T}{\tau_D} + 2a \left(\frac{r_0}{R}\right)^2 \frac{\Delta\Theta}{\tau_S} \exp\left(-2\frac{t}{\tau_S}\right), \quad (24) $$

where the dimensionless constant

$$ a = c_1 \frac{\eta_e}{\eta_S} \quad (25) $$

characterizes the proportion of heat effectively dissipated in the metal and $\tau_D$ the characteristic time of diffusion

$$ \tau_D = \frac{c_2 R_{Fe}^2}{3\kappa}, \quad (26) $$

where $c_2$, measuring in terms of $R_{Fe}$ the thickness of the thermal boundary layer around the metal, $\delta = c_2 R_{Fe}$, is a dimensionless number.

Eq. (24) cannot be used predictively in a complex situation as it requires the knowledge of various parameters $c_1$, $c_2$ and $a$. The dependences of these parameters with more fundamental quantities (mostly with the temperature dependence of the viscosity) have to be determined empirically. We will see however, that for a given choice of the rheology, Eq. (24) captures the evolution of the metallic diapir temperature as a function of time and the dependence of this temperature with the diapir size. For example, Eq. (24) suggests that the diffusion term decreases with $R_{Fe}$ (as $R_{Fe}^{-2}$ if one considers $c_2$ as a constant) while the dissipation term increases with $R_{Fe}^2$. We can also use Eq. (24) qualitatively by assuming $a \sim c_1 \sim 4/15$ (using Stokes law) and $c_2 \sim 1$.

Eq. (24) shows that the temperature is not necessarily an increasing function of time. More precisely, according to Eq. (24) the metal temperature increases just after the impact ($t \sim 0$), if

$$ -\frac{\Delta T_0}{\tau_D} + 2a \left(\frac{r_0}{R}\right)^2 \frac{\Delta\Theta}{\tau_S} = 0 \quad (27) $$

Using the expressions for the temperature increase upon impact $\Delta T_0$ (see Eq. (1)), the maximum temperature increase during segregation $\Delta\Theta$ (see Eq. (21)) and for the two time constants $\tau_S$ and $\tau_D$ (see Eqs. (13) and (26)), this condition implies that dissipative heating overcomes the conductive diffusion when

$$ R_{Fe} > R_{Fe,\min} \quad (28) $$

where $R_{Fe,\min}$ involves the properties of the planet, but not its radius since $\Delta T_0$ is proportional to $R^2$

$$ R_{Fe,\min}^4 = \frac{9}{8\pi} \left(\frac{r_0}{R}\right)^2 \frac{1}{c_1 c_2 a} \frac{\Delta T_0}{\Delta\Theta} \frac{\eta_S \kappa}{G \rho_0 (1 - f_0) \Delta\rho_0}. \quad (29) $$

According to the set of parameters shown in Table 1, $R_{Fe,\min} \sim 45$ km (using $c_1 \sim a \sim 4/15$, $r_0 \sim R$ and $c_2 \sim 1$). However, using values fitted from experiments does not change this radius very much for the moderate
level of viscosity variations used in our simulations hereafter. Such a diapir corresponds to an impactor of radius $R_{imp} \sim 60$ km ($R_{imp} = R_{ic}/3^{1/3} = R_{Fe}/(3f_0)^{1/3}$). Therefore, only impactors larger than $R_{imp} = 60$ km generate metallic diapirs that heat up during sinking, although their initial temperature set by the impact is not dependent on the size of the impactor.

Integration of Eq. (24) leads to

$$ \Delta T = \Delta T_0 \exp\left(-\frac{t}{\tau_D}\right) + a \left(\frac{r_0}{R}\right)^2 \Delta\Theta \frac{2\tau_D}{2\tau_D - \tau_S} \left(\exp\left(-\frac{t}{\tau_D}\right) - \exp\left(-2\frac{t}{\tau_S}\right)\right). \quad (30) $$

The initial temperature anomaly $\Delta T_0$ decreases exponentially with time while the interplay between diffusion and dissipation controls the general temperature evolution. For the diapir to heat up, the heating time $\tau_S/2$ must be shorter than the diffusive time $\tau_D$. Typically $r_0 \sim R$ and in the regime where the diapir heats up, the dissipation occurs before the diffusion, $\tau_S/2 \ll t \ll \tau_D$; the temperature rapidly increases to $\Delta T = \Delta T_0 + a(r_0/R)^2 \Delta\Theta$, and the physical interpretation of $a$ is therefore the percentage of heat dissipated inside the metal. According to Eq. (25), $a$ should be lower than the coefficient $c_1$ of the Rybczinski–Hadamard velocity as the effective viscosity of the hot diapir $\eta_e$ is likely lower than the average viscosity $\eta_S$. For a numerical application we take however $a \sim c_1 = 4/15 \sim 0.27$ as obtained for the isoviscous Rybczinski–Hadamard velocity. As $\Delta\Theta$ and $\Delta T_0$ are simultaneously proportional to $R^2$, the maximum temperature of the diapir is at most $\Delta T = 4.2\Delta T_0$ and is independent of the planet size.

Dissipation decreases as $\exp(-2t/\tau_s) = (r/r_0)^2$ according to Eq. (12). Hence, the dissipation term in Eq. (24) decreases with depth. When a diapir heats up, its temperature increases therefore to the maximum $\Delta T_{max}$ reached at the radius $r$ that satisfies $d\Delta T/dt = 0$ or

$$ 0 = -\frac{\Delta T_{\max}}{\tau_D} + 2a \frac{\Delta\Theta}{\tau_S} \left(\frac{r}{R}\right)^2, \quad (31) $$

which implies

$$ \left(\frac{r}{r_0}\right)^2 = \frac{\Delta T_{\max}}{\Delta T_0} \left(\frac{R_{\text{Fe,min}}}{R_{\text{Fe}}}\right)^4. \quad 32 $$

The factor $\Delta T_{max}/\Delta T_0$ varies between 1 (no heating) and 4.2 (maximum estimated temperature). As an example, an impactor of radius 120 km, generates a metallic diapir of 96 km (two times $R_{Fe,min}$) that heats up until it reaches half the radius of the impacted planet. The expression Eq. (32) is only valid when $R_{Fe} > R_{Fe,min}$, otherwise the diapir temperature simply decreases.

## 5. Numerical simulations

We compare the predictions of the analytical model to spherical axisymmetric calculations of a sinking metallic drop, especially to extract the diffusive and sinking times $\tau_D$ and $\tau_S$ and the fraction of heat trapped in the metallic phase (e.g., the constants $c_1, c_2$ and $a$, that we expect to be close to 4/15, 1 and 4/15). We then compare these results to more complex numerical experiments where a compositional anomaly is generated in the isobaric core after a large impact. The effect of variable viscosity is also studied in these models.

### 5.1. Numerical models of sinking metallic drops

#### 5.1.1. Sinking velocity

We solve numerically a set of problems in which we introduce metallic spheres ($f_0 = 1$) of different sizes, tangent to the surface, in undifferentiated planets ($f_0 = 0.17$) of various radii. From this set of experiments, we compare the temporal evolution of the sphere position to what is predicted by Eq. (12). The calculations presented

here are isoviscous for simplicity but variable viscosity will be introduced in more complex cases. Fig. 2 shows that the values of $\tau_S$ obtained by fitting the center of the diapir position to an exponential in the numerical models, vary as $1/R_{Fe}^2$ as expected from the analytical model, with $c_1 = 0.187$ (almost 70% of the Hadamard–Rybczynski velocity for a homogenous viscosity $4/15 = 0.27$). For large sphere radii, boundary effects are stronger and the sinking times are slightly larger.

#### 5.1.2. Temperature evolution

Large sinking diapirs heat up before cooling down by diffusion when the velocity of the metal decreases sufficiently towards the center. Our theoretical predictions given by Eq. (30) are in good agreement with the computed evolutions using the value $c_1$ obtained previously. Fig. 3 shows the consistency between the numerical results and the theory when the parameters $c_2$ and $a$ are fitted ($c_2 = 0.72, a = 0.2$ which is reasonably close to $c_1 = 0.187$). The value of $a$, indicates that 20% of the released heat is trapped in the metal. The maximum temperature value, $2.2\Delta T_0$, is in rough agreement with the estimate $\Delta T = \Delta T_0 + a(r_0/R)^2 \Delta\Theta = 2.88\Delta T_0$. This value is obtained for sufficiently large impactors (> 200 km) since smaller ones can cool off very early upon sinking as seen from Eq. (32).

We monitor the temperature evolution for various diapir radii. Fitting the temperature evolution with Eq. (30) leads to values of $\tau_D$ and $a$ for each diapir radius. The corresponding characteristic diffusive times are plotted in Fig. 4. These times are consistent with analytical

<figure><img src="image_54.png" /><figcaption>Fig. 2. Characteristic sinking time $\tau_S$ as a function of $1/R_{Fe}^2$ normalized with the impacted planet radius. Results from numerical experiments (with uniform viscosity $\eta_S = 10^{22}$ and $R = 1000$ km) are represented with black circles. Theoretical fit from Eq. (13) is shown by the dashed line with $c_1 = 0.187$.</figcaption></figure>

<figure><img src="image_55.png" /><figcaption>Fig. 3. Temperature evolution (black line) of a metallic sphere ($R_{Fe} = 130$ km) falling in an undifferentiated planet with $R = 1000$ km. Theoretical evolution from Eq. (30) is shown with a dashed line ($c_1 = 0.187, c_2 = 0.72$ and $a = 20\%$).</figcaption></figure>
predictions from Eq. (26) and increase with the square of the diapir size. For all the experiments, the fraction of heat a trapped in the metal is therefore reasonably constant (~22 ± 5%) and close to $c_1$.

To verify Eq. (32) that predicts the radius for which dissipation overcomes diffusion, we computed the rate of heating or cooling of metallic spheres as a function of their radius and depths. Various planetary radii have been used and, as predicted, the heating always occur in the external part of the planet (filled symbols). Near the center of the impacted planet, when the gravity decreases, diffusion dominates (open symbols) and the temperature of the sinking metallic phase decreases. As shown in Fig. 5, the transition between heating and cooling occurs consistently within the shaded area predicted by the analytical expressions Eq. (32). For small diapirs (i.e. for $R_{Fe} \le 45$ km), diffusion dominates and prevents heating. Large diapirs reach their maximum temperature and start cooling near the high temperature estimate of the analytical model.

## 5.2. Application to global evolution after an impact

The thermo-chemical initial conditions after an impact differ from a simple hot metallic sphere sinking within an undifferentiated material. Indeed, the denser metallic pond collected at the bottom of

<figure><img src="image_56.png" /><figcaption>**Fig. 4.** Non-dimensional characteristic time of diffusion $\tau_D$ as a function of the non-dimensionalized metallic sphere radius. Results from numerical experiments (with a uniform viscosity and $R=1000$ km) are represented with black circles. Theoretical fit from Eq. (26) is shown in dashed line with $c_2 = 1.01$.</figcaption></figure>

<figure><img src="image_57.png" /><figcaption>**Fig. 5.** Thermal behaviour of a sinking metal sphere in an undifferentiated media as function of position and sphere radius. Each symbol represents the instantaneous thermal behaviour of a hot metallic sphere with radius $R_{Fe}$ for a given initial position. Filled symbols represent numerical experiments with viscous heating and open symbols represent numerical experiments with only cooling. Different symbols characterize different planets radii. The analytical transition between heating and cooling is predicted within the shaded area and the boarders of this area are defined with $\Delta T_{max}/\Delta T_0$ between 1 (no heating) and 4.2 (maximum heating) (see, Eq. (32)).</figcaption></figure>

the isobaric core is not spherical and above it, a volume of light silicates rises and spreads underneath the surface until it covers the entire surface of the planet. These deviations from our analytical model potentially modify the results obtained from the sinking metallic drop model. Here we show numerical simulations of segregation after an impact and compare them to the analytical model previously developed.

Fig. 6 depicts the thermal and compositional evolution after an impact of a large impactor ($R=4000$ km, $R_{imp}=600$ km and $R_{Fe}=480$ km). The four rows correspond to real time snapshots at 0, 1.4, 3.8 and 546 Myrs. The temperature field is depicted in the left column, and the composition in the right column (undifferentiated material in light blue, metal in red, silicates in green). The metallic pond sinks towards the center of the planet while heating. This heating is in agreement with our previous findings that dissipation is larger than diffusion for large impacts. However, the metal develops a tail through sinking and is significantly deformed. In the meantime, the light silicates rise upward and heat up as well, while stretching laterally to cover the whole surface of the planet. Of course, the diffusion of heat out of the silicate layer near the surface, is much faster than that out of the deep protocore and this shallow hot silicate layer cools rapidly. On a much longer time scale (assuming unrealistically that no other impact occurs, hot thermal plumes should start from the proto core–mantle boundary and deliver the protocore heat to the surface (Behounkova and Choblet, 2009).

Fig. 7 illustrates the evolution of the conversion from potential to thermal energy with time. During the thermo-chemical reequilibration, the potential energy (thick line) decreases as the metal approaches the center and as the silicates spread beneath the surface. Viscous heating induces an increase of thermal energy (grey line). Once the metal has reached the center of the impacted protoplanet, the thermal energy can only decrease. During this whole process, heat

<figure><img src="image_58.png" /><figcaption>**Fig. 6.** Non-dimensional temperature (left) and composition (right) at times $t=0$ (first line), $t=1.4$ Myr (second line), $t=3.8$ Myr (third line) and $t=546$ Myr (fourth line) (computed for a uniform viscosity with $R=4000$ km, $R_{imp}=600$ km and $200 \times 200$ grid points).</figcaption></figure>
<figure><img src="image_59.png" /><figcaption>**Fig. 7.** Non-dimensionalized potential (solid black line) and thermal (solid grey line) energies and time integrated surface heat flow (dotted black line) as functions of time. The sum of these three quantities times 100 is shown in dashed black line. Its difference to zero is indicative of the accuracy of the energy conservation of the numerical code (for R = 2000 km, R<sub>imp</sub> = 300 km and R<sub>Fe</sub> = 240 km and uniform viscosity).</figcaption></figure>

is slowly removed by diffusion through the surface of the planet and the cumulative heat flux (dotted line) balances the total energy budget. This global balance (sum of potential energy, thermal energy and cumulative heat flux (see Eq. (14))) is closely satisfied which illustrates the good accuracy of the numerical code.

We now introduce a temperature-dependence of the viscosity in the calculations. Experimental results suggest that the viscosity contrast between melt iron and solid silicates can reach 20 orders of magnitude (Vocadlo et al., 2000). Such a viscosity contrast is difficult to handle numerically and we use much smaller values.

In our models, the viscosity varies as $\eta = \eta_0 \lambda^T$ and as the temperature of metal may increase while sinking by a factor up to 2, it implies maximum viscosity contrasts up to $1/\lambda^2$ orders of magnitude between cold and hot materials. Using a composition dependent viscosity would have been more realistic but viscous fronts are too difficult to handle numerically. We compare the thermo-

<figure><img src="image_60.png" /><figcaption>**Fig. 8.** The four rows depict the temperature (left) and the composition (right) at t = 3.2 Myr (with R = 2000 km and R<sub>imp</sub> = 300 km), for a uniform viscosity (top) and for variable viscosities (contrast of ~16 (second row), ~100 (third row) and ~1600 (bottom row)). As expected, the sinking velocity of the metallic diapir and the rising velocity of the silicates, both increase when their viscosity is decreased.</figcaption></figure>

<figure><img src="image_61.png" /><figcaption>**Fig. 9.** Position of the inertia center of the metal phase as a function of time for a uniform viscosity (black line) and for temperature-dependent viscosities with $\lambda = 0.25$ (dashed dotted line), $\lambda = 0.1$ (grey line) and $\lambda = 2.5 \times 10^{-2}$ (dotted line) (R = 2000 km and R<sub>imp</sub> = 300 km). Thin dashed lines correspond to simple exponential fittings from which the sinking times are extracted (see Table 2).</figcaption></figure>

chemical states at the same time, t = 3.2 Myr for different viscosity factors in Fig. 8. We use $\lambda = 0.25$ (Fig. 8 second row), $\lambda = 0.1$ (Fig. 8 third row) and $\lambda = 2.5 \times 10^{-2}$ (Fig. 8 bottom row), the top row being the reference isoviscous case.

Increasing the temperature-dependence of the viscosity softens the surrounding material around the metallic drop and the metallic diapir, at a given time, is closer to the center when its viscosity is decreased, as shown in Fig. 8. However, this effect remains small. Because the metallic pond becomes less viscous, its shape becomes more spherical and the tail developed in isoviscous experiments becomes thinner. Increasing the sinking velocity increases the rate of shear heating but not the total release of thermal energy which is only related to the change in gravitational energy. Lowering the viscosity in the surrounding material and within the metallic pond has also the effect of diminishing $\eta_e$. The dissipation is therefore increased in the undifferentiated material and decreased in the hot and less viscous metallic diapir. This effect combined with the faster spreading of the hot silicate that removes the heat more rapidly lead to lower maximum temperatures (see Fig. 8).

We monitor the position of the inertia center of the metallic diapir as a function of time and compute the sinking times $\tau_S$ (see Fig. 9). The position of the diapir obeys reasonably the exponential law predicted by Eq. (12). In the isoviscous case, the observed normalized time is $\tau_S = 563$ which is twice longer than what is predicted by Eq. (13). This is due to the fact that the initial diapir shape is not spherical and to the presence of the rising volume of silicates. When the viscosity decreases with temperature the sinking is faster, $\tau_S = 249, 170$ and 114, for $\lambda = 0.25, 0.1$ and $2.5 \times 10^{-2}$ (see Fig. 9 and Table 2). This is due to two effects: the reduction of viscosity inside the metal (the Rybczinski–Hadamard formula predicts an increase of the velocity factor $c_1$ from 0.27 to 0.33 when the interior viscosity of the diapir decreases) and the decrease of viscosity of the heated surrounding material.

In the experiments depicted in Figs. 6 and 8, the metal temperature increases and reaches a value close to twice the initial temperature of the isobaric core (Fig. 10). However, heating within the metal is less

<caption>**Table 2**</caption>

Values obtained by fitting numerical experiments with theoretical predictions (Eqs. (12) and (30)) for different values of $\lambda$ (with R = 2000 km and R<sub>imp</sub> = 300 km).

<table><thead><tr><th></th><th>λ = 1</th><th>λ = 0.25</th><th>λ = 0.1</th><th>λ = 2.5 × 10<sup>-2</sup></th></tr></thead><tbody><tr><td>τ<sub>S</sub></td><td>563</td><td>249</td><td>170</td><td>114</td></tr><tr><td>τ<sub>D</sub></td><td>20054</td><td>16520</td><td>13316</td><td>8974</td></tr><tr><td>a</td><td>19%</td><td>14.7%</td><td>11%</td><td>7%</td></tr></tbody></table>
<figure><img src="image_62.png" /><figcaption>**Fig. 10.** Temperature evolution of the metal phase as a function of time for a uniform viscosity (solid black line) and for temperature-dependent viscosities with $\lambda = 0.25$ (dashed dotted line), $\lambda = 0.1$ (grey line) and $\lambda = 2.5 \times 10^{-2}$ (dotted line) ($R = 2000$ km and $R_{\text{imp}} = 300$ km). Thin dashed lines correspond to theoretical results from Eq. (30) from which the diffusive times and the proportion of energy heating the metal diapir are extracted (see Table 2).</figcaption></figure>

pronounced with variable viscosity and decreases with the viscosity contrast. Fitting the computed temperature evolutions in the metallic diapirs with our theoretical model gives values of $c_2$ in the isoviscous case and the variable viscosity cases (see Table 2 and Fig. 10). The thickness of the thermal boundary measured by $c_2$ decreases with the sinking velocity (the Peclet number). The values of $c_2$ and of $\tau_D$ are therefore related to $Pe^{-n} \propto \tau_s^n$ with an exponent $\sim 1/3$ in the range of values, $n = 1/2 - 1/3$ predicted in Ribe (2007).

When the temperature dependence of the viscosity increases, the proportion of energy heating the metal diapir, $a$, decreases (see Table 2). As a consequence, the heat release of the gravitational energy becomes increasingly efficient in the surrounding undifferentiated material. This suggests that a diapir of very small viscosity does not heat much during its motion while most of the release of gravitational energy occurs in the undifferentiated materials. A low viscosity diapir keeps basically its initial temperature because its characteristic diffusive time is larger than its sinking time and also because of the buffering effect of the temperature dependent viscosity (i.e., a too large cooling would increase the viscosity and would bring back the dissipation within the diapir itself).

**6. Discussion and conclusion**

Core formation events induced by meteorological impacts play a major role in determining the early thermo-chemical state of growing planets. Large meteorological impacts can trigger a local differentiation between metal and silicates in a spherical zone above the surface called the isobaric core. The segregation of dense and light phases through the undifferentiated material of the impacted protoplanet induces a large viscous heating.

We followed the dynamics of the metal phase after a large impact with numerical experiments in axisymmetrical spherical geometry. The sinking velocity of the metal phase is Stokes-like and is function of the viscosity contrast between the metal phase and the undifferentiated crossed media. The velocity increases when viscous heating decreases the viscosity of the surrounding material. A stress dependent viscosity (not considered here) would also increase this velocity (Samuel and Tackley, 2008). The sinking process in a planet with a cold interior compared to its surface would eventually imply higher viscosity contrasts between the metal and the surrounding material and would lead to longer sinking times.

The gravitational energy release during the segregation is converted into viscous heating in the metal and in the silicates. Our results show that a net viscous heating of the metallic phase only occurs for large metallic diapirs ($R_{\text{Fe}} > 45$ km). This metallic volume at

the bottom of the isobaric core would be produced by an impactor of order $R_{\text{imp}} > 60$ km. This result underlines the importance of accretion conditions on the inner thermal state of planetary bodies. Small metallic diapirs cool while sinking and may ultimately bring the metal in a solid state to the core of the impacted planet.

The heat repartition between the metal phase, the silicates and the undifferentiated material is not only a function of the size of the metallic diapir but also of the rheology of the various phases. For low viscosity of the metal and of the sheared zone around the metallic diapir, the metal phase is weakly heated. Hence, gravitational energy release will mainly lead to the heating of the surrounding undifferentiated material and ultimately to its differentiation.

The viscosity variations that we explore in our simulations are of order $\lambda^{\Delta T_{\text{max}}}$ which in the most extreme cases reach about four orders of magnitude over very short distances. This is certainly modest relative to the viscosity contrasts of 20 orders of magnitude that exist between liquid metal and solid silicates (Vocadlo et al., 2000). Viscosity contrasts based on composition rather than temperature would be more realistic but would have occurred on even shorter distances (the computation grid itself) that could not be resolved with classical numerical methods. Our model is therefore an end-member of possible models on heating modes during core formation. However, the description of the physics of the processes would still be valid for larger viscosity contrasts.

As soon as a growing planet reaches a few thousand kilometers in radius $R$, the heating by impacts becomes significant (the temperature increase varies as $R^2$ and reaches 400 K for $R = 3000$ km, (Monteux et al., 2007)). This temperature increase superimposed on the fossil temperature $T_0$ from short half-life radionucleides ($^{26}$Al and $^{60}$Fe) and previous impacts can lead to a temperature larger than the melting temperature of the metallic phase. Our analytical models confirmed by numerical experiments show that the metallic drop reaches the planet center in a time depending on the size of the metallic drop and the background viscosity of the planet but not of its radius (see Eq. (13)). Even in the case where the impacted planet is relatively cold and with a high viscosity of $10^{22}$ Pa s, this time is smaller than a few million years for an impactor of 300 km. The sinking timescales obtained in our models are comparable to those obtained with an Arhenius rheology (Ziethe and Spohn, 2007) and within the timeframe required for an early core formation (<60 My). The temperature increase in the undifferentiated material localized along the sinking path of the metallic diapir could provide a preferential low viscosity channel for the following differentiation events.

Proposing predictive models for the thermal consequences of differentiation after an impact is fundamental in order to understand the thermal state of the interior of growing planets. As shown in Ricard et al. (2009), core formation of terrestrial protoplanets could be the consequence of a runaway segregation induced by a large enough impact on undifferentiated material. These results also underline the importance of accretionary conditions (size and temporal repartition of impacts) on the thermal energy repartition and, hence, on the magnetic history of growing planets (Elkins-Tanton et al., 2005).

**Acknowledgements**

The authors thank H. Samuel, S. Labrosse and M. Jellinek for helpful discussions. The authors also thank H. Schmeling, T. Spohn and anonymous reviewers for constructive comments. This project was funded by the A.N.R. ETHER.

**References**

Agee, C.B., 1997. Melting temperatures of the Allende meteorite: implications for a Hadean magma ocean. Phys. Earth Planet. Int. 100, 41–47.

Behounkova, M., Choblet, G., 2009. Onset of convection in a basally heated spherical shell, application to planets. Earth Planet. Sci. Lett., submitted.

Christensen, U.R., Yuen, D.A., 1985. Layered convection induced by phase transitions. J. Geophys. Res. 90, 10291–10300.
Croft, S.K., 1982. A first-order estimate of shock heating and vaporization in oceanic impacts. Geological Implications of Impacts of Large Asteroids and Comets on Earth: In: Silver, T.L., Schultz, P.H. (Eds.), Spec. Pap. Geol. Soc. Am., vol. 190.

Davies, G.F., 1982. Ultimate strength of solids and formation of planetary cores. Geophys. Res. Lett. 9, 1267–1270.

Douglas, J., 1955. On the numerical integration of $\frac{\partial^2 u}{\partial x^2} + \frac{\partial^2 u}{\partial y^2} = \frac{\partial u}{\partial t}$ by implicit methods. J. Soc. Ind. Appl. Math. 3, 42–65.

Elkins-Tanton, L.T., Zaranek, S.E., Parmentier, E.M., Hess, P.C., 2005. Early magnetic field and magmatic activity on Mars from magma ocean cumulate overturn. Earth Planet. Sci. Lett. 236, 1–12.

Fei, Y., Bertka, C.M., Finger, L.W., 1997. High-pressure iron–sulfur compound, $Fe_3S_2$, and melting relations in the FeFeS system. Science 275, 1621–1623.

Flasar, F.M., Birch, F., 1973. Energetics of core formation: a correction. J. Geophys. Res. 78, 6101–6103.

Gerya, T.V., Yuen, D.A., 2007. Robust characteristics method for modelling multiphase visco-elastic thermo-mechanical problems. Phys. Earth Planet. Int. 163, 83–105.

Ghosh, A., McSween, H.Y., 1998. A thermal model for the differentiation of asteroid 4 vesta, based on radiogenic heating. Icarus 134, 187–206.

Golabek, G.J., Schmeling, H., Tackley, P.J., 2008. Earth's core formation aided by flow channelling instabilities induced by iron diapirs. Earth Planet. Sci. Lett. 271, 24–33.

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

Hewitt, J.M., McKenzie, D.P., Weiss, N.O., 1975. Dissipative heating in convective flows. J. Fluid Mech. 68, 721–738.

Höink, T., Schmalzl, J., Hansen, U., 2005. Formation of compositional structures by sedimentation in vigorous convection. Phys. Earth Planet. Int. 153, 11–20.

Honda, R., Mizutani, H., Yamamoto, T., 1993. Numerical simulation of Earth's core formation. J. Geophys. Res. 98, 2075–2090.

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

Laney, C.B., 1998. Computational gasdynamics. Cambridge University Press, Cambridge.

Monteux, J., 2009. Modélisation numérique de la formation du noyau terrestre : contribution des impacts météoritiques. Ph.D. thesis, Université Claude Bernard Lyon 1.

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

O'Keefe, J.D., Ahrens, T.J., 1977. Impact-induced energy partitioning, melting, and vaporization on terrestrial planets. In: Merril, R.B. (Ed.), Lun. Planet. Sci. Conf., vol. 8, pp. 3357–3374.

Peaceman, D.W., Rachford, H.H., 1955. The numerical solution of parabolic and elliptic differential equations. J. Soc. Ind. Appl. Math. 3, 28–41.

Pierazzo, E., Vickery, A.M., Melosh, H.J., 1997. A reevaluation of impact melt production. Icarus 127, 408–423.

Ratcliff, J.T., Tackley, P.J., Schubert, G., Zebib, A., 1997. Transitions in thermal convection with strongly variable viscosity. Phys. Earth Planet. Int. 102, 201–212.

Ribe, N.M., 2007. In: Schubert, G. (Ed.), Analytical Methods in Mantle Dynamics, Treatise of Geophysics, vol. 7. Elsevier.

Ricard, Y., 2007. In: Schubert, G. (Ed.), Physics of Mantle Convection, Treatise of Geophysics, vol. 7. Elsevier.

Ricard, Y., Sramek, O., Dubuffet, F., 2009. A multi-phase model of runaway core–mantle segregation in planetary embryos. Earth Planet. Sci. Lett., in press.

Roe, P.L., 1986. Characteristic-based schemes for the Euler equations. Annu Rev Fluid Mech 18, 337–365.

Rybczynski, W., 1911. über die fortschreitende bewegung einer flüssigen kugel in einen medium. Bull. Acad. Sci. Cracovie 1, 40–46.

Samuel, H., Tackley, P.J., 2008. Dynamics of core formation and equilibration by negative diapirism. Geochem. Geophys. Geosyst. 9, 6011–6026.

Schubert, G., Turcotte, D.L., Olson, P., 2001. Mantle convection in the Earth and planets. Cambridge University Press. Sep.

Senshu, H., Kuramoto, K., Matsui, T., 2002. Thermal evolution of a growing Mars. J. Geophys. Res. 107, 1–13.

Shannon, M.C., Agee, C.B., 1996. High pressure constraints on percolative core formation. Geophys. Res. Lett. 23, 2717–2720.

Sramek, O., 2007. Modèle d'écoulement biphasé en sciences de la terre: Fusion partielle, compaction et différenciation. Ph.D. thesis, Ecole Normale Supérieure de Lyon.

Stevenson, D.J., 1989. Formation and early evolution of the Earth. In: Peltier, W.R. (Ed.), Mantle convection and plate tectonics.

Stevenson, D.J., 1990. Fluid dynamics of core formation. In: Newsom, H.E., Jones, J.H. (Eds.), Origin of the Earth. Oxford Univ, New York.

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

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

Vocadlo, L., Alfè, D., Price, G.D., Gillan, M.J., 2000. First principles calculations on the diffusivity and viscosity of liquid Fe–S at experimentally accessible conditions. Phys. Earth Planet. Int. 120, 145–152.

Von Bargen, N., Waff, H.S., 1986. Permeabilities, interfacial areas and curvatures of partially molten systems: results of numerical computations of equilibrium microstructures. J. Geophys. Res. 91, 9261–9276.

Yin, Q., Jacobsen, S.B., Yamashita, K., Blichert-Toft, J., Télouk, P., Albarède, F., 2002. A short timescale for terrestrial planet formation from Hf-W chronometry of meteorites. Nature 418, 949–952.

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

Ziethe, R., Spohn, T., 2007. Two-dimensional stokes flow around a heated cylinder: a possible application for diapirs in the mantle. J. Geophys. Res. 112, 1–13.
<figure><img src="image_63.png" /></figure>

Laboratory experiments on the breakup of liquid metal diapirs

Jean-Baptiste Wacheul<sup>a, b</sup>, Michael Le Bars<sup>a, c, *</sup>, Julien Monteux<sup>d</sup>, Jonathan M. Aurnou<sup>a</sup>

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

CrossMark

ARTICLE INFO

Article history:

Received 21 January 2014

Received in revised form 14 May 2014

Accepted 26 June 2014

Available online 24 July 2014

Editor: C. Sotin

Keywords:

metal diapir

breakup

core formation

magma ocean

metal-silicate equilibration

laboratory experiments

ABSTRACT

The validity of the iron rain scenario, i.e. the widely accepted model for the dynamics of iron sedimentation through a magma ocean during the latest stage of the Earth's accretion, is explored via a suite of laboratory experiments. Liquid gallium and mixtures of water and glycerol are used as analogs of the iron and the molten silicate respectively. This allows us to investigate the effects of the viscosity ratio between iron and silicate and to reproduce the relevant effects of surface tension on the fragmentation dynamics. While the classical iron rain scenario considers a population of purely spherical drops with a single characteristic radius that fall towards the bottom of the magma ocean at a unique velocity without any further change, our experiments exhibit a variety of stable shapes for liquid metal drops, a large distribution of sizes and velocities, and an intense internal dynamics within the cloud with the superimposition of further fragmentations and merging events. Our results demonstrate that rich and complex dynamics occur in models of molten metal diapir physics. Further, we hypothesize that the inclusion of such flows into state of the art thermochemical equilibration models will generate a similarly broad array of complex, and likely novel, behaviors.

© 2014 Elsevier B.V. All rights reserved.

# 1. Introduction

Differentiation of Earth into a core and a mantle was likely completed within the first tens million years after its accretion (e.g. Stevenson, 1990; Boyet et al., 2003; Boyet and Carlson, 2005). Numerical simulations (Neumann et al., 2012) and geochemical data on meteorites (Yoshino et al., 2003) also show that small planetesimals could have differentiated even earlier when accounting for heating by decay of short-lived radionuclides. There is also strong evidence that the Earth's late accretion is due to collisions with large planetesimals (a tenth to a third of Earth mass), when both the impactor and the proto-Earth were already differentiated (Morbidelli et al., 2012). During accretion, the Earth and other planets in formation underwent several mechanisms of heating: 1) the decay of relatively abundant radioactive elements with short half-life (Merk et al., 2002; Walter and Tronnes, 2004), 2) the conversion of gravitational potential energy by viscous forces during differentiation (Rubie et al., 2007; Monteux et

al., 2009; Samuel, 2012), 3) the collisions themselves with the conversion of huge amount of kinetic energy (Safronov, 1978; Kaula, 1979; Reese and Solomatov, 2006; Monteux et al., 2007), these impacts alone being able to generate a local melting resulting in a shallow magma ocean (Tonks and Melosh, 1992). In addition, the primitive atmosphere was certainly much more opaque to IR radiation, so the effect of thermal blanketing was highly efficient (Abe and Matsui, 1985). Thus, according to the simulations, Earth has probably had one or several episodes of global magma ocean, with depths possibly reaching thousands of kilometers (Tonks and Melosh, 1993). In this context, further impacts of differentiated planetesimals would require, in order for the cores of the Earth and the meteorites to merge, that the latter flows through the magma ocean (Fig. 1). This process can be seen as a secondary step of mixing between core and mantle, since it could lead to partial or complete thermo-chemical equilibration between the sinking metal and the molten surrounding silicates, depending on the characteristics of the flow of the core material through the mantle.

The importance of this exchange is an issue for the interpretation of numerous geochemical proxies, such as the tungsten 182 signal. Hafnium (¹⁸²Hf) disintegrates in tungsten (¹⁸²W) with a relatively short half-life of 9 My, comparable with the time scale of core differentiation. They are both refractory but tungsten is siderophile whereas hafnium is lithophile. This is why the

<footnote><sup>a</sup> Department of Earth, Planetary and Space Sciences, University of California, Los Angeles, CA 90095-1567, USA</footnote>

<footnote><sup>b</sup> École Normale Supérieure, Paris, France</footnote>

<footnote><sup>c</sup> CNRS, Aix-Marseille Université, École Centrale Marseille, IRPHE UMR 7342, Marseille, France</footnote>

<footnote><sup>d</sup> Laboratoire de Planétologie et de Géodynamique, LPG Nantes, CNRS UMR 6112, Université de Nantes, France</footnote>

<footnote>* Corresponding author at: CNRS, Aix-Marseille Université, École Centrale Marseille, IRPHE UMR 7342, Marseille, France. Tel.: +33 4 13552020; fax: +33 4 13552001.</footnote>

<footnote>E-mail address: lebars@irphe.univ-mrs.fr (M. Le Bars).</footnote>
<figure><img src="image_65.png" /><figcaption>**Fig. 1.** Left: schematic of the metal/silicate separation after the impact of a differentiated planetesimal with the early Earth in reference to the work of Stevenson (1990). The equilibration by diffusion between the impactor's core and the magma ocean during the metal rainfall and later when sinking through the solid mantle as a large diapir is poorly constrained, and strongly depends on the fluid dynamics of the iron sedimentation. Right: schematic of our experiment.</figcaption></figure>

¹⁸²Hf–¹⁸²W system has been used in multiple studies to approximate the age of terrestrial bodies' core (Lee and Halliday, 1996; Harper and Jacobsen, 1996; Kleine et al., 2004). However due to the late impacts and possible mixing that could occur, there are very poor constraints on which event – late impacts or core differentiation – is relevant for the interpretation of the ratio ¹⁸²W/¹⁸⁴W. Then the age given by this proxy could be any intermediary between the initial differentiation of the proto-Earth and the most recent giant impact that it endured, depending on how much ¹⁸²W has been absorbed by the asteroid's core during its passing through the mantle (Kleine et al., 2004). The same kind of interrogation can be held against interpretations of the U/Pb proxy, and for the coefficients of partition between metal and silicate, which strongly depend on the details of the small-scale processes at the iron-silicate interface during sedimentation. Actually, these interfacial dynamics influence every mechanism of equilibration by diffusion, such as diffusion of heat and diffusion of momentum by viscosity, both leading to indetermination on the initial thermal state of the mantle and the core, and on the reparation of the energy between these two (Monteux et al., 2009; Samuel et al., 2010). Thus, in order to model the evolution of both Earth's core and mantle, it is important to understand the fluid dynamics at the drop scale during the iron sedimentation (Solomatov, 2000). Towards this goal, we present here novel laboratory experiments investigating the fluid dynamics of sedimenting liquid metal droplets.

## 2. Parameters controlling the fluid dynamics of the iron sedimentation

The equilibration between the iron and molten silicate strongly depends on the typical size of the metal entities. Indeed, for a given volume of metal, a single large diapir would fall rapidly through the magma ocean with a relatively small surface of exchanges, while the fragmentation of the same volume of iron through a large number of small structures broadens the surface area of exchanges and slows down the falling velocity, hence extending the time during which iron and silicate equilibrate. Note that in the present paper, we generically use the term diapir to designate any large blob of fluid moving through an ambient fluid via the action of buoyancy forces.

Several approaches have been developed in order to give a physically coherent description of what happens when a liquid iron diapir falls through a magma ocean, and ultimately to provide a time scale for the equilibration. At first order, the shape of the falling diapir is dominated by two forces. The surface tension tends to stabilize a spherical shape, while the dynamic pressure

deforms the diapir and tears it apart. Let us assume for instance, a typical diapir with a radius $R_0 = 10$ km falling at the inviscid, free fall Newtonian velocity $V$ valid for a rigid sphere

$$V \approx \sqrt{\frac{\Delta \rho}{\rho_a}} g R_0, \quad (1)$$

where $\rho_a$ is the silicate density ("a" standing for "ambient"), $\Delta\rho$ is the density difference between iron and silicate, and $g$ is the gravitational acceleration. Assuming an Earth with more than half its final mass, $V$ is close to 1 km/s. The Re number for the flow in the mantle, which estimates the ratio of the inertial and viscous forces, is

$$Re_a = \frac{\rho_a V R_0}{\mu_a} \approx 10^{10}, \quad (2)$$

where $\mu_a$ is the silicate dynamic viscosity. This large Re value provides an *a posteriori* validation of the inviscid velocity estimate provided by Eq. (1). It also implies that dynamic pressure scales as the inertia. On the other hand, the characteristic strength of surface tension is directly linked to the radius of curvature of the surface, so its order of magnitude roughly depends on the radius of the spheroid diapir $R_0$. A good estimation of the stability of a diapir is given by the Weber number, which is the ratio of the inertial and surface tension forces:

$$We_a = \frac{\rho_a V^2 R_0}{\sigma}, \quad (3)$$

where $\sigma$ is the coefficient of surface tension. For $We \gg 1$, diapirs are unstable and break-up. Below some threshold of order 1 (e.g. $We_c = 6$ for rain drops, see Villermaux and Bossa, 2009), surface tension and inertia compensate, and the diapir is stable. This widely used breakup criterium (e.g. Tonks and Melosh, 1992; Rubie et al., 2003; Dahl and Stevenson, 2010; Deguen et al., 2011) allows a calculation of the maximal radius for stable diapirs, given some hypothesis regarding its falling speed. For the simple Newtonian velocity scaling given by (1), the maximal radius corresponds to

$$R_{cap} \approx \sqrt{We_c \frac{\sigma}{\Delta\rho g}}, \quad (4)$$

which is about 1 to 2 cm for the Earth's iron-silicate system. Such a criterium is well known in the case of water drops in the air, for which it has been confirmed by experiments (Villermaux and Bossa, 2009). It has also been supported by a recent numerical study designed for the case of an iron diapir in molten
silicate (Samuel, 2012). In this study, the initial spherical diapir flattens, then breaks up within a distance of no more than 10 di- apir radii. The diapir’s sons keep breaking up, cascading downward to the scale $R_{cap}$ where surface tension compensates the pressure forces. However, in these axisymmetric simulations, the breakup is not actually computed, but is supposed to occur shortly af- ter the topological change from pancake to torus. Thus, the size of the resulting drops remains unknown. According to this sce- nario, there is no way for a 10 km radius diapir corresponding to a Weber number above $10^{14}$ to remain entirely intact during its fall. It should rapidly break up into a cloud of drops of radius $R_{cap}$, i.e. the so-called “iron rain”. Then, most models of metal- silicate equilibration (e.g. Rubie et al., 2003; Ulvrová et al., 2011; Samuel, 2012) assume that all iron drops have an identical and fixed spherical shape with a radius equal to $R_{cap}$ and a unique sedimentation velocity based on free fall models.

On the other hand, for diapirs with a radius of 10 km and above, the Weber number is so large (above $10^{14}$) that its dynam- ics are controlled by the inertia of the flow only, allowing one to completely ignore the effects of surface tension. The resulting mis- cible models lead to interesting findings. For instance, Deguen et al. (2011) supposed that the diapir rapidly becomes a cloud of drops that then evolves in a coherent manner as a buoyant thermal, whose radius grows linearly with depth because of turbulent mix- ing with the outside. They found that for a shallow magma pond created by an impact as for a deep magma ocean (Deguen et al., 2014), there should be enough dilution between metal and silicate for them to reach chemical equilibration. On the contrary, Dahl and Stevenson (2010) considered a diapir with an almost constant shape eroded by Raleigh–Taylor and Kelvin Helmholtz instabilities. Their conclusions were very different: erosion was found to be in- sufficient to allow for full equilibration for diapirs with a radius above 10 km.

Open questions thus remain on all stages of the iron sedimenta- tion, from the large-scale dynamics to the behavior at the smallest scales, where surface tension modifies the modalities of diffusive exchanges. In this context, the numerical study by Ichikawa et al. (2010), in agreement with the results for water drops in the air (Villermaux and Bossa, 2009), pointed out that the cloud of drops contains a whole distribution of sizes and not just a single one (Karato and Murthy, 1997). Furthermore, the distribution of speeds and the relation between size and speed of drops are not triv- ial, supporting the idea that the iron rain scenario with a single size and speed may be inaccurate. However, the conclusions of Ichikawa et al. (2010) are supported mostly by 2D numerical simu- lations, where surface tension is not properly rendered. Hence they still need confirmation.

In addition, it is worth noting that all previous studies of the stability of diapirs have neglected the influence of the viscosity ratio between the ambient fluid and the metal

$$r_{\mu} = \frac{\mu_{a}}{\mu_{m}}. \qquad (5)$$

When the magma ocean cools down, or as the diapir goes deeper into the mantle, the magma becomes closer to the solidus and contains more crystals. Therefore, the viscosity of the magma is expected to increase from $10^{-3}$ to $10^2$ Pa s, corresponding to $r_{\mu}$ ranging from 1 to $10^5$ (Deguen et al., 2011; Samuel, 2012). Such a large viscosity ratio is known to have a significant effect on the shape of falling drops (e.g. Bonometti and Magnaudet, 2006; Ohta et al., 2010). Indeed, from the continuity of the velocity and the stress tensor at the drop interface, a large $r_{\mu}$ in a large Reynolds number flow implies intense internal recirculations in- side the drop, which then deforms and may be dynamically stabi- lized. This has not been investigated in a geophysical context.

<table><caption>**Table 1**
List of relevant parameters and typical Earth’s and experimental values.</caption><thead><tr><th>Symbol</th><th>Parameter</th><th>Value for the Earth</th><th>Value in our experiment</th></tr></thead><tbody><tr><td>ρ<sub>a</sub></td><td>ambient fluid density</td><td>3000 kg m<sup>-3</sup></td><td>997–1260 kg m<sup>-3</sup></td></tr><tr><td>ρ<sub>m</sub></td><td>liquid metal density</td><td>7000 kg m<sup>-3</sup></td><td>6095 kg m<sup>-3</sup></td></tr><tr><td>μ<sub>a</sub></td><td>ambient fluid viscosity</td><td>10<sup>-3</sup>–10<sup>2</sup> Pa s</td><td>10<sup>-3</sup>–1 Pa s</td></tr><tr><td>μ<sub>m</sub></td><td>liquid metal viscosity</td><td>10<sup>-3</sup> Pa s</td><td>1.9 × 10<sup>-3</sup> Pa s</td></tr><tr><td>r<sub>μ</sub></td><td>viscosity ratio</td><td>1 to 10<sup>5</sup></td><td>0.5 to 500</td></tr><tr><td>σ</td><td>surface tension</td><td>1 J m<sup>-2</sup></td><td>0.7 J m<sup>-2</sup></td></tr></tbody></table>

Here we report the first experiments designed to simulate the last stages of the fragmentation process with realistic values of the iron–silicate viscosity ratio and relevant behaviors regarding the ef- fects of surface tension on the fragmentation dynamics. The set-up is presented in Section 3. In Section 4, we investigate the variety of stable shapes for iron drops, as well as their mutual interactions and dynamics after the initial breakup. The measured distributions of size and velocity are presented in Sections 5 and 6, highlighting the influence of the viscosity ratio. The relevance of the complex fluid dynamics shown by our experiments to the Earth is illus- trated in Section 7 by a simple equilibration model based on our results. Conclusions and open questions are given in Section 8.

### 3. Set-up and methods

Our experimental set-up is sketched in Fig. 1. As an analogue for the magma ocean, we use a 160 cm high cylindrical tank, with a diameter of 19 cm, filled with a mixture of water and glyc- erol. The glycerol (resp. water) has a dynamic viscosity of 1.08 Pa s (resp. 0.00093 Pa s) at room temperature (23 °C): the mixture of the two allows us to explore a range of 3 orders of magnitude for the viscosity of the ambient fluid $\mu_a$, with a density rang- ing from $\rho_a = 1260$ kg m$^{-3}$ for pure glycerol to $\rho_a = 997$ kg m$^{-3}$ for pure water. As an analogue for the liquid iron diapir, we use liquid gallium. It has a viscosity $\mu_m = 1.9 \times 10^{-3}$ Pa s and a density $\rho_m = 6095$ kg m$^{-3}$ (King and Aurnou, 2013). The gal- lium is initially contained in a latex balloon at the top of the set-up; the balloon is then popped by a syringe needle at the beginning of the experiment. This method has two advantages: the amount of gallium is precisely known by weighing the bal- loon, and since the retraction of the balloon occurs within about 1/50 s, the diapir has no initial speed and its initial shape is the one imposed by the balloon. The fall of the diapir is then recorded by a high speed USB camera at 136 frames per sec- ond, with a resolution of 196 × 1280 pixels (examples of these videos can be viewed at https://www.youtube.com/watch?v=g- AmGmWWK1o&list=UU7u8FUioww0nKsdXPO6emhA). In addition to this camera, higher resolution videos of the lower part of the experiment are taken at 60 frames per second with a resolution of 1280 × 720 pixels. Beyond direct visualization, the videos are used to recover the shape and velocity of the droplets, after removing the background and after binarization. Each droplet in a selected frame is detected using the Matlab Image Processing toolbox, and an equivalent radius is retrieved by measuring its apparent area A and applying the formula $r = \sqrt{A/\pi}$. This is a lower estimate since the drops are oblate at different degrees. We also construct space–time diagrams by extracting the same horizontal line from all frames of a chosen video. The resulting image then gives us the horizontal radius of each droplet reaching the selected depth as well as its arrival time, hence its mean fall velocity.

Relevant parameters are listed in Table 1 in comparison with Earth’s values. Experiments were run with 6 different sizes of the initial diapir ranging from an equivalent radius of 14 mm to 30 mm, and with 4 viscosities of the ambient fluid 0.001, 0.01, 0.1 and 1 Pa s. In the inertial regime, the characteristic velocity of the
<figure><img src="image_66.png" /><figcaption>**Fig. 2.** In blue, area of the parameter space associated with the fall of iron diapirs in the context of Earth's accretion. Variations are due to changes in the viscosity of the magma ocean and in the diapir radius, ranging from its initial value to the radius corresponding to a Weber number of 1. For comparison, the area explored by previous studies (Ichikawa et al., 2010; Samuel, 2012; Deguen et al., 2014) is shown in green, whereas the parameter space explored by the present study is shown in red. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

flow scales as the Newtonian velocity introduced in Eq. (1). Then, the dynamics is entirely characterized by four dimensionless pa- rameters: the density ratio $\Delta\rho/\rho_a$, the viscosity ratio $r_\mu = \mu_a/\mu_m$, and the Newtonian Reynolds and Weber numbers $Re_a$ and $We_a$ (based on the Newtonian velocity). In our experiment, $\Delta\rho/\rho_a$ remains almost constant at about 5 (the ambient density only marginally varies with the proportion of water), which means that as in the geophysical setting, density changes are of order 1, such that the Boussinesq approximation cannot be applied. With the ac- cessible range of initial diapir radii, we are able to produce initial Reynolds numbers from 10 to $4 \times 10^4$ and initial Newtonian We- ber numbers from 14 to 64. As shown in Fig. 2, this is obviously limited compared to the possible geophysical values, especially for diapirs with large initial radius. But even if the dynamic similitude between the experiment and reality is not exact, the experiments are capable of reaching the relevant dynamical regime, with a fully turbulent flow and a Weber number above the critical value for breakup. Furthermore, the similitude is exact for diapir radii around the decimeter scale, hence for the final stages of the fully- developed iron rain. Note finally that our experiment is the first

one to take into account the effects of the viscosity ratio, which spans the range 0.5–500 in our set-up.

### 4. Zoology of the breakup

The results of 3 experiments in pure glycerol with different ini- tial radius of the diapir are shown in Figs. 3A, 3B and 4. Fig. 5 shows a closer look of the various possible shapes of gallium drops obtained in the experiment presented in Fig. 4. According to the classical iron rain scenario applied to our experimental set- up, any diapir with a Weber number larger than 6 should rapidly give rise to a cloud of spherical drops with a single characteris- tic radius $R_{cap} = 1$ cm that should fall with a Newtonian velocity $V = 74$ cm s$^{-1}$, without any further dynamical change. This is not the case.

For the lowest initial radius (Fig. 3A), the diapir is stable even if its Weber number is above the known threshold 6 (note that the effective Weber number is $\sim 9$, calculated using the measured falling velocity instead of the scaling given by Eq. (1)). One can also notice that its equilibrium shape is not spherical: the drop takes the form of a cap, where the intense internal recirculation driven by viscous coupling with the ambient fluid stabilizes the non- spherical shape of the molten metal droplet (see also the sketch in Fig. 5B). The viscosity ratio strongly influences the shape and the condition for stability of a single structure, as found in the stud- ies of Bonometti and Magnaudet (2006), Ohta et al. (2010) using axisymmetric numerical simulations. There is a clear tendency for large viscosity ratios to stabilize the drops, and the diapir shown in Fig. 3A is indeed unstable for $r_\mu$ below 50 (not shown here). According to results in Bonometti and Magnaudet (2006), Ohta et al. (2010), which are coherent with our experimental observations, the critical Weber number can actually be more than one order of magnitude larger than the classically used value $We_c = 6$, depend- ing on the viscosity ratio and the initial shape of the drop.

The slightly larger diapir shown in Fig. 3B rapidly breaks up into three large caps plus some smaller drops. A striking point here is that after this initial break-up, the dynamics is not frozen: the three caps interact and two of them finally merge to re-build a larger diapir. For comparison, the same experiment but with an ambient viscosity 100 times less viscous is shown in Fig. 3C. The initial behavior of the diapir is similar but the breakup dynam- ics is clearly different, even if the Weber number is the same in both experiments. The Reynolds number is greater in Fig. 3C be- cause of the smaller ambient viscosity; but comparing both series of pictures, one can notice that the falling velocities are close in the two experiments: both cases are in the Newtonian regime, and changes in $Re_a$ cannot explain changes in the fragmentation be- havior. Thus, we argue that the differences between Figs. 3B and 3C are due to viscosity ratio effects, which allow for the dynami- cal stabilization of larger drops by strong internal recirculations in case 3B.

The dynamics of initially large diapirs systematically follows the series of stages illustrated in Fig. 4: after a short acceleration, we observe waves forming on the surface of the diapir, qualitatively similar to the description of Dahl and Stevenson (2010). But al- most simultaneously, the diapir flattens as the wave amplitudes rise: it evolves towards a thin wavy sheet where the axisymme- try is fully broken. It breaks up shortly after this stage: holes appear in the sheet, the transient ligaments retract and break-up. The drops resulting from this burst have various sizes and shapes, and the biggest ones continue to break up in the same way un- til a steady-state is reached. The whole process can be seen as a downward cascade toward small scales where surface tension is important. However, we also observe multiple coalescence of droplets several times in a row, corresponding to a non-linear in- verse cascading process. From the experimental videos, we see that
<figure><img src="image_67.png" /><figcaption>**Fig. 3.** Successive snapshots with a fixed time interval for 3 experiments. A: Fall of a 14 mm initial radius diapir in pure glycerine, corresponding to $Re_a = 12$, $We_a = 14$ and a viscosity ratio of 500. The diapir has a stable cap shape. B: same as A but for a 23 mm initial radius diapir, corresponding to $Re_a = 24$ and $We_a = 37$. The diapir breaks up in 3 main cap-shaped stable drops of close sizes. C: same as B but for an ambient fluid 100 times less viscous, corresponding to $Re_a = 2400$, $We_a = 37$ and a viscosity ratio of 5.</figcaption></figure>

small drops accelerate and merge when they are near bigger drops whose rear recirculation engulfs them. One can argue that the narrowness of our experiment limits the spreading of the metallic droplet cloud, thereby enhancing the typical frequency of droplet collisions. But as noted above, coalescence happens even if there are very few drops. We thus believe that this inverse cascade process is also relevant in larger domains, as for instance in magma oceans.

Fig. 5A shows a closer instantaneous view of the droplet cloud once the steady-state is reached. A whole distribution of shapes and sizes is observed. Large drops have a cap shape stabilized by internal recirculation (see Fig. 5B), and the smallest drops adopt a spherical shape; a continuity of flattened ellipsoidal shapes is observed in between.

We have also observed a novel metastable structure (Fig. 5C) formed during the breakup of the biggest diapirs. These structures

look like a hot air balloon and fall slowly relative to other droplets of comparable size. From slow motion videos of their formation and disappearance, we have inferred that they are bubbles of galium enclosing glycerol (see Fig. 5D), similar to bubbles of soapy water in air. This bubble structure explains their relative stability as well as their anomalously low settling velocity for their rather large size.

**5. Distribution of sizes and influence of the viscosity ratio**

Fig. 6 shows the cumulative distribution of droplets sized obtained from the breakup of the largest class of diapir for a viscosity ratio of $r_μ = 50$. This corresponds to $Re_a = 368$ and $We_a = 64$. Taking into account the resolution of the video (9.0 pixels cm$^{-1}$), we are able to detect the size of drops larger than 0.6 mm in radius, the smaller ones being detected as drops of 0.6 mm. The
<figure><img src="image_68.png" /><figcaption>**Fig. 4.** Successive snapshots with a fixed time interval $\Delta t = 0.18$ s for a 30 mm initial radius diapir falling through pure glycerol, corresponding to $Re_a = 37$, $We_a = 64$ and a viscosity ratio of 500. One can notice that the mean speed of the front of the diapir is rather constant.</figcaption></figure>

breakup of metal diapirs does not create droplets of one single size but a whole distribution of equivalent radii, in agreement with the simulations of Ichikawa et al. (2010). Because the formation of drops results from the generic process of ligament rupture, their size distribution is well fitted by a Gamma distribution, whose probability density function $p$ writes:

$$p(x) = \frac{x^{k-1}e^{-x/\theta}}{\theta^k \Gamma(k)}, \qquad (6)$$

where $k$ is the shape of the Gamma distribution, $\theta$ its scale, and $\Gamma(k)$ is the Gamma function evaluated at $k$. Such a Gamma distribution of droplet sizes is similar to the one obtained in the case of water drops in the air (Villermaux and Bossa, 2009). Note however that in this last case, the viscosity ratio is reverse (i.e. the more viscous fluid is inside) and the shapes of the obtained drops are very different. Our distribution is tightened around a mean radius of 4 mm. This value can be related to a breakup criterion, now understood in a statistical sense: surface tension sets the characteristic length scale of the distribution, the mean radius. Using our experimental results, the critical Weber number corresponding to this radius is $We_c \simeq 1$. The distributions obtained for diapirs with different initial sizes are similar to the one shown in Fig. 6, and so is the measured mean radius, provided that these diapirs are large enough to create a distribution of sizes that converges statistically. This condition is verified for the 4 biggest classes of diapirs that we have produced.

Fig. 7 shows series of snapshots from 4 experiments with the same initial diapir but different viscosity ratios. With our present set-up, because a large number of drops are superimposed on the video images, it was not possible to detect their individual contours for a viscosity ratio smaller than 50. Hence we could not perform a systematic quantitative study of the sizes distribution as a function of $r_\mu$. But relying on direct observation, we see that there is a clear tendency for large viscosity ratios to stabilize bigger drops, as already noticed in Section 4 for single structures. In

all cases, we expect to systematically recover a Gamma-type distribution for the equivalent radii. This means that the distributions always have the same shape, with a peak at a small scale corresponding to $We_c \simeq 1$, and an exponential tail. But we expect the slope of this tail to be significantly more gentle when the viscosity ratio increases. This corresponds to a decreasing value of the shape of the Gamma distribution: for instance, Villermaux and Bossa (2009) found a shape value of 4 for the breakup of water in the air (viscosity ratio $2 \times 10^{-2}$), while we find a shape of 2.2 for a viscosity ratio of 50, as shown in Fig. 6.

## 6. Simultaneous distributions of sizes and velocities

The interactions between the droplets lead to a wide range of sizes and velocities in our experiments. Fig. 8 shows the distribution of sizes and speeds for a viscosity ratio $r_\mu = 50$ and a diapir with an initial radius of 23 mm, using the values obtained from a space-time diagram at a distance of 140 cm from the initial position of the center of mass of the diapir. It is plotted in the same way as in Ichikawa et al. (2010). The fact that velocities are calculated from the travel time through the tank averages out a large part of the variability due to raw turbulence and allows us to then consider the mean structure of the flow. Interestingly, the drops' velocities do not follow a Newtonian scaling based on their individual radius, even when adjusting the pre-factor (see Fig. 8). This result seems to validate the entrainment hypothesis described by Deguen et al. (2011, 2014): after the breakup, the drops fall as an interacting cloud whose velocity is determined by the inertia of the whole flow, related to the initial mass of the diapir. Additional fluctuations are related to the turbulent mixing and interactions between drops.

## 7. Typical equilibration length

We can estimate a rough length scale of equilibration following the same reasoning as in Samuel (2012), but using the distribution
<figure><img src="image_69.png" /><figcaption>**Fig. 5.** A: drops of liquid gallium of various sizes falling through pure glycerol. Several intermediates between cap shape and spherical shape are present. B: expected streamlines of the relative flow within and around a falling diapir with a cap shape (shown in gray). C: gallium bubble enclosing glycerol formed during the breakup of a 30 mm radius diapir in pure glycerol. Due to its reduced buoyancy compared to full gallium drops, it falls slowly. These structures last from one tenth of a second to several minutes. D: schematic of a gallium bubble.</figcaption></figure>

of sizes and velocities found in our experiments. Note that the relatively simple equilibration models shown below are meant only for illustration of potential consequences of the complex dynamics exhibited in our experiments. The question of geochemical equilibration clearly deserves more detailed studies.

In our model, we neglect the breakup distance and consider that the distribution of sizes and speeds has reached a steady state. From our observations, from simulations (Samuel, 2012) as well as from previous experiments (Deguen et al., 2014), breakup occurs within a typical length scale of a few initial radii, so the previous hypothesis is valid for initial entities at least 10 times smaller than the mantle depth: we choose here an initial diapir of radius 10 km sedimenting in a magma ocean with a minimum depth of 100 km. We suppose that the metal concentration $C_m$ in a chemical element of interest (e.g., Ni, Co, W, Hf) is uniform inside each drop and that it evolves following a simple Fick's law

$$ \frac{4\pi}{3} r^3 \frac{dC_m}{dt} = -4\pi r^2 \kappa \frac{C_m - C_{m,eq}}{\delta_{BL}}, \quad (7) $$

where r is the drop radius, κ the chemical diffusivity coefficient of typical order $10^{-8}$ m² s⁻¹, and $C_{m,eq}$ the equilibrium concentration. We further assume that the thickness of the chemical boundary layer $\delta_{BL}$ scales as $\delta_{BL} = \sqrt{2\kappa r/v}$, where v is the local speed

<figure><img src="image_70.png" /><figcaption>**Fig. 6.** Histogram of the equivalent radius of the droplets created by the breakup of a 30 mm initial radius diapir falling through a mixture of glycerol and water (in blue) and best fitting Gamma distribution (in red, with shape 2.2 and scale 1.9). Dimensionless parameters of the experiment are $Re_a = 368$, $We_a = 64$ and $r_\mu = 50$. The low cutoff radius is set by the resolution of the video to 0.6 mm. A large number of measurements is necessary for obtaining converged statistics. This cumulative distribution was thus constructed from 6 runs of the same experiment and using the last 10 frames of each video, i.e. once a statistically steady state is reached and before the fastest droplet touches the bottom of the tank. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>
<figure><img src="image_71.png" /><figcaption>**Fig. 7.** Series of snapshots of the diapirs evolution for 4 runs of the experiment with a 31 mm initial radius diapir. From top to bottom, only the viscosity of the ambient fluid increases; the corresponding viscosity ratios are 0.5, 5, 50, and 500. Rigorously, changing the ambient viscosity also changes the Reynolds number. But it may be noticed from the snapshots that the mean sedimentation velocities in the four cases are comparable, and, hence, largely independent of the ambient fluid viscosity, as expected in the Newtonian limit. We argue then that the various observed dynamics are primarily related to changes in the viscosity ratio.</figcaption></figure>
<figure><img src="image_72.png" /><figcaption>**Fig. 8.** Distribution of the drops equivalent radii (in mm, along the x-axis) and speeds (non-dimensionalized by the Newtonian speed of the initial diapir, along the y-axis). The initial radius of the diapir is 23 mm and the viscosity ratio is 50, which gives $Re_a = 244$ and $We_a = 37$. Three curves of the corresponding Newtonian velocity scaling (1) with a prefactor respectively equal to 1, 1/2 and 1/3 from top to bottom, are shown in green. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

of the flow around the droplet, assumed to be constant (following Samuel, 2012). We define the degree of equilibration as

$$C^*(t) = \frac{C_m(t) - C_m(t=0)}{C_{m,eq} - C_m(t=0)}. \qquad (8)$$

Eq. (7) then leads to an exponential solution

$$C^*(t) = 1 - e^{-t/\tau} \quad \text{where } \tau = \frac{r^{3/2}}{3} \sqrt{\frac{2}{\kappa \nu}}. \qquad (9)$$

Following our experimental results, we evaluate the degree of equilibration using three different values for the local speed of the drops and the speed at which the cloud of drops is entrained. The reference case corresponds to the standard iron rain scenario with a cloud of drops with a uniform radius $R_{cap}$, and where the local and global speeds are the corresponding Newtonian speed. Two other scenarios based on our experimental results are investigated. The local speed is given in one case by the Newtonian speed for each drop whose radius follows a Gamma distribution, and in the other case by the size-speed distribution results presented in Fig. 8, where we are assuming that the distribution does not depend on the large scale parameters such as the mean diameter and velocity of the drop cloud. In both cases, the global speed is the Newtonian velocity of the initial 10 km radius diapir. This estimate agrees with our present experimental results, as shown, for example, in Fig. 7. The global sedimentation speed, however, probably decreases over time because of the progressive entrainment of ambient fluid, as shown by Deguen et al. (2011, 2014). This effect is not seen here, possibly because of the limited size of our container. Corresponding equilibration results for the three models are shown in Fig. 9.

In all cases, the entire distribution of drops fully equilibrates before reaching the bottom of the magma ocean. However, the predicted depth of equilibration is 1 to 3 orders of magnitude larger when considering the scenarios derived from our fluid mechanics experiments compared to the idealized iron rain. This is the result of two combined effects highlighted by our laboratory experiments: droplet velocities are significantly larger than in the classical iron rain, and the drop size distribution puts a significant fraction of the metal phase in drops larger than $R_{cap}$, which implies a smaller surface of exchange between iron and silicate. As a result, equilibration will integrate a broader range of conditions regarding pressure and temperature, and the energy and

<figure><img src="image_73.png" /><figcaption>**Fig. 9.** Degree of equilibration of a cloud of droplets resulting from the fragmentation of an initially 10 km radius diapir and falling through a magma ocean. The classical iron rain scenario is shown in black for reference. Two scenarios derived from the statistics of our experiments with an initial radius $R = 30$ mm and $r_\mu = 50$ are also shown, where the global speed of the drop cloud is taken as the Newtonian velocity of the initial 10 km radius diapir. In red, the local speed of the drops is derived from the distributions shown in Fig. 8. In blue, the local speed is given by the Newtonian scaling law for each drop radius. Note that the present model neglects the breakup distance of the initial diapir. For an initial diapir of radius 10 km, this implies a minimum depth of the magma ocean of 100 km, shown by the vertical green line. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

mass exchanged by diffusion during the sinking of a particular metal droplet can occur substantively deeper within the magma ocean. Besides, for several elements, the partition coefficient depends on the pressure of equilibration. Expected changes in the depth of equilibration will also change their final repartition between core and mantle.

## 8. Conclusions and open questions

Our laboratory experiments on the fragmentation of gallium drops falling through glycerol shed new light on the dynamics of the iron sedimentation process that likely occurs during a terrestrial body's final accretionary stages. The classical iron rain scenario considers a population of spherical drops with a single characteristic radius that fall towards the bottom of the magma ocean at a unique velocity without any further change. In contrast, here we find that the fragmentation of an initially large diapir actually gives rise to a variety of stable shapes, a large distribution of sizes and velocities, and an intense internal dynamics within the cloud with the superimposition of further fragmentations and merging. Thus, we argue that previous models of chemical and thermal equilibration processes can only be accurate in a statistical mean sense (e.g. Rubie et al., 2003; Samuel et al., 2010). On the basis of our simple equilibration model, we still predict a complete equilibration before reaching the core, but at a significantly deeper depth. One should also keep in mind that the probability of a “strange” event, i.e. an anomalously large diapir or an anomalously slow falling velocity, is statistically possible, especially for large viscosity ratios.

Additional experiments are necessary to complement the conclusions drawn here, notably with a larger tank to avoid confinement effects and to allow for initially larger diapirs. Also it now appears necessary to take into account, in more evolved models of equilibration, the complex internal dynamics between drops inside the clouds, including the observed inverse cascade and the global sedimentation dynamics. It is also necessary, in addition to the first study for fixed-shape spherical drops by Ulvrová et al. (2011), to account for the very intense internal dynamics inside and outside large drops, which both stabilizes and deforms them, and should increase the equilibration process by advection.
All these effects clearly deserve further studies that are beyond the scope of this paper. Finally, our experiments have highlighted the importance of the viscosity ratio on the fluid dynamics of the fragmentation and sedimentation processes. As in Earth, the heat brought by the conversion of gravitational and kinetic energy during accretion is not negligible (see e.g. Tonks and Melosh, 1992; Samuel et al., 2010), it would now be interesting to study the strong coupling between the heating by viscous damping of the intense flows caused by the fall of iron diapirs, the changes in the ambient viscosity induced by this thermal evolution, and the corresponding evolution of the drop size distribution.

## Acknowledgements

We thank Professor E. Villermaux for helpful discussions. M.L.B. acknowledges financial support from the European Commission, Research Executive Agency, Marie Curie Actions (project FP7-PEOPLE-2011-IOF-298238). J.M. is funded by Agence Nationale de la Recherche (Accretis decision No. ANR-10-PDOC-001-01). J.M.A. and the UCLA SpinLab are supported by the NSF Geophysics Program.

## References

Abe, Y., Matsui, T., 1985. The formation of an impact-generated H₂O atmosphere and its implications for the thermal history of the Earth. J. Geophys. Res. 90, 545–559.

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

Boyet, M., Carlson, R.W., 2005. ¹⁴²Nd evidence for early (>4.53 Ga) global differentiation of the silicate Earth. Science 309, 576–581.

Boyet, M., Blichert-Toft, J., Rosing, M., Storey, M., Télouk, P., Albarède, F., 2003. ¹⁴²Nd evidence for early Earth differentiation. Earth Planet. Sci. Lett. 214, 427–442.

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

Deguen, R., Landeau, M., Olson, P., 2014. Turbulent metal-silicate mixing, fragmentation, and equilibration in magma oceans. Earth Planet. Sci. Lett. 391, 274–287. http://dx.doi.org/10.1016/j.epsl.2014.02.007.

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

Harper, C., Jacobsen, S., 1996. Evidence for ¹⁸²Hf in the early solar system and constraints on the timescale of terrestrial accretion and core formation. Geochim. Cosmochim. Acta 60, 1131–1153.

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

Karato, S.-I., Murthy, V., 1997. Core formation and chemical equilibrium in the Earth—I. Physical considerations. Phys. Earth Planet. Inter. 100, 61–79.

Kaula, W.M., 1979. Thermal evolution of Earth and moon growing by planetesimal impacts. J. Geophys. Res. 84, 999–1008.

King, E.M., Aurnou, J.M., 2013. Turbulent convection in liquid metal with and without rotation. Proc. Natl. Acad. Sci. USA 110 (17), 6688–6693.

Kleine, T., Mezger, K., Münker, C., Palme, H., Bischoff, A., 2004. ¹⁸²Hf–¹⁸²W isotope systematics of chondrites, eucrites, and martian meteorites: chronology of core formation and early mantle differentiation in Vesta and Mars. Geochim. Cosmochim. Acta 68, 2935–2946.

Lee, D.-C., Halliday, A.N., 1996. Hf-W isotopic evidence for rapid accretion and differentiation in the early solar system. Science 274, 1876–1879.

Merk, R., Breuer, D., Spohn, T., 2002. Numerical modeling of ²⁶Al-induced radioactive melting of asteroids considering accretion. Icarus 159, 183–191.

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

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

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

Neumann, W., Breuer, D., Spohn, T., 2012. Differentiation and core formation in accreting planetesimals. Astron. Astrophys. 543, A141.

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

Reese, C., Solomatov, V., 2006. Fluid dynamics of local martian magma oceans. Icarus 184, 102–120.

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

Rubie, D.C., Nimmo, F., Melosh, H.J., 2007. Formation of Earth's core. In: Treatise on Geophysics. Elsevier, pp. 51–90.

Safronov, V.S., 1978. The heating of the Earth during its formation. Icarus 33, 3–12.

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

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

Solomatov, V.S., 2000. Fluid dynamics of a terrestrial magma ocean. In: Origin of the Earth and Moon. University of Arizona Press, pp. 323–338.

Stevenson, D.J., 1990. Fluid dynamics of core formation. In: Origin of the Earth. Oxford University Press, pp. 231–249.

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

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

Ulvrová, M., Coltice, N., Ricard, Y., Labrosse, S., Dubuffet, F., Velimský, J., Srámek, O., 2011. Compositional and thermal equilibration of particles, drops, and diapirs in geophysical flows. Geochem. Geophys. Geosyst. 12, 1–11.

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

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

Yoshino, T., Walter, M., Katsura, T., 2003. Core formation in planetesimals triggered by permeable flow. Nature 422, 154–157.
# Silicate melts during Earth's core formation

M.A. Bouhifd a,*, V. Clesi a, A. Boujibar b, N. Bolfan-Casanova a, C. Cartier c, T. Hammouda a, M. Boyet a,
G. Manthilake a, J. Monteux a, D. Andrault a

## ARTICLE INFO

**Article history:**

Received 3 March 2016
Received in revised form 17 December 2016
Accepted 21 December 2016
Available online xxxx

**Keywords:**

Magma ocean
Earth's core
Magma ocean crystallization
Bridgmanite
Liquidus
Helium
Nb-Ta paradox
Redox state

## ABSTRACT

Accretion from primordial material and its subsequent differentiation into a planet with core and mantle are fundamental problems in terrestrial and solar system. Many of the questions about the processes, although well developed as model scenarios over the last few decades, are still open and much debated. In the early Earth, during its formation and differentiation into rocky mantle and iron-rich core, it is likely that silicate melts played an important part in shaping the Earth's main reservoirs as we know them today.

Here, we review several recent results in a deep magma ocean scenario that give tight constraints on the early evolution of our planet. These results include the behaviour of some siderophile elements (Ni and Fe), lithophile elements (Nb and Ta) and one volatile element (Helium) during Earth's core formation. We will also discuss the melting and crystallization of an early magma ocean, and the implications on the general feature of core-mantle separation and the depth of the magma ocean. The incorporation of Fe²⁺ and Fe³⁺ in bridgmanite during magma ocean crystallization is also discussed. All the examples presented here highlight the importance of the prevailing conditions during the earliest time of Earth's history in determining the composition and dynamic history of our planet.

© 2016 Elsevier B.V. All rights reserved.

## 1. Introduction

Silicate melts have played an important role in Earth's evolution from the core formation 4.5 billion years ago to present-day volcanic activity. Modern planetary formation models show that terrestrial planets likely formed through three modes of growth that were ordered in time in the protoplanetary disk. In step I, km-sized planetesimals formed by the coagulation of gas and dust in less than ~10⁴ years. In step II, the collisional evolution of the planetesimals population led to the growth of Moon to Mars-sized planetary embryos formed by runaway growth through collisions among the planetesimals in <10⁶ years. In step III, giant collisions among the planetary embryos resulted in the formation of a few terrestrial planets that swept up all the other bodies. This stage may have taken up to ~10⁸ years (e.g. Chambers, 2004; Morbidelli et al., 2012). At least step II and step III were energetic enough to have likely melted a large fraction of the silicate Earth (e.g. Canup, 2012). Theory and observations point in fact to the occurrence of magma oceans in the early evolution of terrestrial planets and in many planetesimals (e.g. Greenwood et al., 2005; Elkins-Tanton, 2012). In addition a magma ocean scenario provides a favorable environment in which metal-silicate segregation can occur rapidly and efficiently (e.g. Stevenson, 1990).

It is now well established that for metal and silicate to segregate rapidly on a planetary scale requires that at least the metal is molten (Stevenson, 1990). On the other hand, when the silicate is also largely molten, liquid metal can segregate extremely efficiently (e.g. Rubie et al., 2003). The short time scale of Earth's core formation, ranging from 30 to 100 Ma after the beginning of solar system accretion (e.g. Kleine et al., 2009) can be also used as an indirect evidence for large scale melting during the first 100 Ma of Earth's history. There are three main sources of energy that can produce the melting that is required for core formation: (1) the decay of short-lived radioactive nuclides (²⁶Al and ⁶⁰Fe) during the first 1–3 Ma of the solar system evolution (Yoshino et al., 2003); (2) the energy delivered by impacts that can generate local or global melting especially during the late stages of Earth accretion (Tonks and Melosh, 1993); (3) the conversion of potential energy into heat via viscous dissipation during the segregation of the metallic phase that can increase the silicate temperature by 1000–3000° (Ke and Solomatov, 2009; Monteux et al., 2009; Samuel et al., 2010; more details can be found in the review paper by Rubie et al., 2015b).

Melting events are common in the first stages of Solar System history and start very rapidly after the beginning of solar nebula condensation, dated at 4.567–4.568 Ga, age of refractory inclusions, known as CAIs (see Amelin et al., 2002; Bouvier and Wadhwa, 2010; Connelly et al., 2012). Oldest ages measured in iron meteorites show that metal-silicate separation began on planetesimals within 1 to 1.5 Ma of the start of

<footnote>a Laboratoire Magmas et Volcanes, Université Blaise Pascal, CNRS UMR 6524, Campus Universitaire des Cézeaux, 6 Avenue Blaise Pascal, 63178 Aubière Cedex, France</footnote>

<footnote>b NASA Johnson Space Center, 2101 NASA Parkway, Houston, TX 77058, USA</footnote>

<footnote>c Department of Geology (B20), Université de Liège, B-4000 Sart Tilman, Belgium</footnote>

<footnote>* Corresponding author.
E-mail address: A.Bouhifd@opgc.univ-bpclermont.fr (M.A. Bouhifd).</footnote>
Solar System formation (Kruijer et al., 2012, 2014). The study of basaltic
achondrites (angrites, eucrites) coming from different parent bodies
shows that small planetesimals have experienced magma ocean stages
during the first million years of their history (e.g. Brennecka and
Wadhwa, 2012; Touboul et al., 2015). This chronology shows that the
building blocks of forming planets implied mainly differentiated objects
instead of primitive and undifferentiated materials. In planets, the last
stages of accretion end with very large collisions. The most spectacular
event for Earth corresponds to the Moon-forming giant impact event.
The study of short-lived systematics in terrestrial samples gives evi
dence for successive magma-ocean stages on Earth. Noble gas systemat-
ics trace the mantle's degassing and the Xenon isotope signature is
consistent with extensive loss of volatile elements within the first
100 Ma of Earth's existence (Allègre et al., 1983; Yokochi and Marty,
2005; Mukhopadhyay, 2012). The high ³He/²²Ne fractionation ob-
served in the Earth's mantle suggests at least two separate magma
ocean steps during the Earth's accretion (Tucker and Mukhopadhyay,
2014). Excesses in ¹⁸²W (produced by the decay of ¹⁸²Hf) measured in
2.7 Ga-old komatiites are explained by a large-scale magmatic differen-
tiation during the first 30 Ma of the Solar System's history (Touboul et
al., 2012). Since the giant Moon-forming impact is dated between 40
and 200 Ma (Boyet et al., 2015 and references therein), this collision
may not have induced a complete homogenization of the Earth's man-
tle. The combined measurements of ¹⁴⁶,¹⁴⁷Sm-¹⁴²,¹⁴³Nd systematics of
undisturbed Eoarchean metabasalts (from the southwestern part of
Greenland, Isua Supracrustal Belt) suggest differentiation of the mantle
source at ~4.47 Ga (Rizo et al., 2011). Chemical fractionation would
have occurred during the crystallization of a large magma ocean
which extended down into the lower mantle. This chronology is consis-
tent with the presence of zircon formed in a terrestrial crust at ~4.4 Ga
(Wilde et al., 2001; Valley et al., 2014).

Concerning the Earth one of the most important events during the first 100 Ma is the segregation of the metallic core. Among the several uncertainties about the formation of the Earth and other planets, the physical conditions upon which Earth differentiated are important questions that remain much debated. Evidence on the earliest times on Earth comes mainly from two sources: (1) The concentrations of core-forming elements left behind in the silicate part provide a fingerprint on the conditions under which this segregation occurs. (2) The Moon is the second important source of evidence, it is believed to have formed after the Earth's core formation as a result of a giant impact between Earth and a Mars-size body (named Theia) (Hartmann and Davis, 1975; Tonks and Melosh, 1993). The magma oceans were reservoirs where partial to complete equilibrium and segregation of the metallic phase from the impactor and the silicate phase from the target were likely to occur (Rubie et al., 2003). The thermo-chemical exchnages that occurred in the magma oceans have strongly influenced the geochemical signature within the current deepest parts of the Earth. Hence they play a key role on how we constrain the differentiation processes in terms of characteristic timescales, P - T conditions or redox conditions. The metal silicate separation processes and especially the characteristic volumes that are involved (cm-droplets to km-diamirs) during the core formation govern the degree of chemical equilibra-tion between the mantle and the core of the Earth (Samuel, 2012; Wacheul et al., 2014). During the latest stages of planetary accretions the large impactors were probably differentiated. After a collision, the impactor's core is assimilated into the partially to completely molten mantle leading to a first fragmentation due to the impact itself. This first fragmentation is poorly constrained but it is the starting point of the second fragmentation that is likely to occur within the turbulent magma ocean. During these two stages of fragmentation, the impactor's core was potentially fragmented down to cm-size droplets (Rubie et al., 2003). If the impactor's core was completely emulsified in droplets ("iron rain scenario"), the thermo-chemical exchanges between the metallic phase from the impactor and the silicates from the magma ocean were extremely efficient (Rubie et al., 2003). Many geochemical

models are based on this “iron rain scenario”. However, after a large im-
pact, the emulsification of the impactor's core might be incomplete and
large volumes of iron may sink without experiencing a complete equil-
ibration with the mantle they are passing through (Dahl and Stevenson,
2010). In this case, a limited fraction of the mantle could equilibrate
with the metallic phase. This could explain the relatively high abun-
dance of highly siderophile elements of the terrestrial mantle that
should have efficiently segregated into the core. Recently, laboratory ex-
periments have constrained the mixing (Deguen et al., 2011, 2014) and
the fragmentation dynamics (Wacheul et al., 2014) during the sinking
of a metallic volume within a magma ocean. Wacheul et al. (2014)
showed that there was a wide range of metallic droplet sizes and shapes
as well as important interactions within the sinking metal cloud in
which fragmentation and coalescence events are superimposed. For
the first time, Wacheul et al. (2014) have estimated the influence of
the viscosity contrast between liquid iron and the liquid silicates and
showed that a large viscosity contrast favors large-scale diapirs. The
characteristic size (that governs the sinking velocity) and the shape of
the metallic droplets (that governs the exchange surface) are both key
parameters to anticipate the degree of equilibration of any chemical el-
ement within a magma ocean (Rubie et al., 2003). This degree of equil-
ibration is a strong function of the partition coefficient of each element
considered (Deguen et al., 2014).

In summary, the Earth accreted over a period of at least 100 Ma with
several magma ocean episodes. During this period the main part of core
formation process occurred. In this paper we will review results involv-
ing exchanges between molten silicate and metal liquid during the first
100 Ma that give rise to strong constraints on the evolution of our plan-
et. This review includes the conditions of the Earth's core formation de-
termined from the behaviour of some siderophile elements (Ni and Fe)
and lithophile elements (Nb and Ta). Recent results also show that the
early core could have incorporated enough helium to supply deep-root-
ed plumes with high ³He/⁴He ratios throughout the age of the Earth. We
will discuss the melting and crystallization of an early magma ocean,
and the implications on the general feature of core-mantle separation
and the depth of the magma ocean. We will also discuss the incorpora-
tion of Fe²⁺ and Fe³⁺ in bridgmanite during magma ocean
crystallization.

**2. Conditions during Earth's core formation**

2.1. *Constraints from siderophile elements*

The knowledge of the pressure, temperature and oxygen fugacity (f<sub>O2</sub>) conditions prevailing during Earth's core formation is derived mainly from the behaviour of the siderophile elements (iron loving elements) that are depleted in the mantle. In this context, the behaviour of Ni and Co has been considered to provide an important clue, both being refractory elements and present in the Earth's mantle in a near-chondritic Ni to Co ratio (Newsom, 1990). In addition, Ni and Co are almost identical to Fe in cosmochemical volatility (Grossman, 1972; O'Neill and Palme, 1998). By using all existing data of metal – silicate partitioning of Ni and Co the consensus at the present time is that metal-silicate equilibration at high pressures, in the range of 40–60 GPa (corresponding to depths of 1000–1500 km), was required to produce the observed Ni and Co depletions in the mantle (see: Bouhifd and Jephcoat, 2011; Righter, 2011; Fischer et al., 2015; Siebert et al., 2012; for the most recent studies). Similar conclusions were reached based on the metal-silicate partitioning of lithophile and weakly-siderophile elements (Mann et al., 2009; Siebert et al., 2011). One should note here that the conditions of Earth's core formation derived from metal-silicate partitioning of siderophile elements cannot be used as arguments for single-stage core formation as this is highly unlikely considering that Earth's core formation occurred during a series of large impact events (e.g. Canup, 2008; O'Brien et al., 2006; Morbidelli et al., 2012; Wetherill, 1985). Such impacts added Fe-rich
metal which segregated to the Earth's proto-core. The experimental results of metal-silicate partitioning of Ni and Co suggest only that some impactor must re-equilibrate in the deep magma ocean above 40 GPa, presumably as a result of emulsification (Rubie et al., 2003) before reaching the center of our planet to form the core, and leaving the observed mantle abundances of siderophile elements.

### 2.1.1. Silicate melts structure control on metal-silicate partitioning of Ni

As noted above, Ni is the most studied element to constrain the conditions of Earth's core formation. The available experimental data of molten metal-silicate liquid partitioning of Ni that range from about 1 GPa up to 100 GPa (e.g. Li and Agee, 1996, 2001; Thibault and Walter, 1995; Ohtani and Yurimoto, 1996; O'Neill et al., 1998; Gessmann and Rubie, 2000; Bouhifd and Jephcoat, 2003, 2011; Chabot et al., 2005; Siebert et al., 2012; Fischer et al., 2015; and references therein) show a decrease of metal-silicate partition coefficients of Ni ($D_{Ni}$) with increasing pressure (at roughly constant oxygen fugacity and temperature following the liquidus of chondritic mantle). Another important feature of $D_{Ni}$ versus pressure is its marked slope change for pressures higher than 35 GPa. As noted by Sanloup et al. (2013), this slope change coincides with a change in the coordination of silicon (in a basaltic melt composition) from 4-fold coordination under ambient conditions to 6-fold coordination at about 35 GPa. One can note that similar results were obtained by studying the structure of basaltic melt at high pressures using first-principles simulations (Bajgain et al., 2015). The similarity of the transition pressure from 4- to 6-fold coordination of Si and the pressure of the marked change in the $D_{Ni}$ coefficients indicates that melt compressibility controls siderophile-element partitioning—at least for Ni (Sanloup et al., 2013). For other siderophile elements, like Cr, V or Mn, metal-silicate partition coefficients versus pressure do not show any slope change at around 35 GPa (see Fischer et al., 2015 for experiments up to 100 GPa). We reported in Fig. 1 a compilation of $D_{Fe}$ versus pressure up to about 55 GPa that shows no slope change at around 35 GPa (the reported experiments in Fig. 1 are at $f_{O_2}$ varying between IW-1 and IW-2). These examples show that melt compressibility is not affecting all siderophile elements partitioning in the same way.

<figure><img src="image_74.png" /><figcaption>**Fig. 1.** Molten metal-silicate liquid partition coefficients of Fe as a function of pressure (GPa). The temperatures of the experiments were above the liquidus of a chondritic composition, and the oxygen fugacity ranges between IW-2 and IW-1. The compiled data are taken from the authors reported in the figure.</figcaption></figure>

### 2.2. Constraints from refractory lithophile elements

How elements are distributed between the main reservoirs (core, magma ocean, atmosphere) during planetary differentiation depends on their chemical affinity. Geochemists and cosmochemists have divided the chemical elements into groups, such as the iron-loving siderophile elements, the silicate-loving lithophile elements, the gas phase-loving atmophile elements, and the sulfur-loving chalcophile elements (Goldschmidt, 1937). This classification has been constructed based on the form in which the elements are encountered at the Earth's surface, i.e. under ambient pressure and temperature, and oxidizing conditions. However, each element behaviour is likely to be modified by intensive parameters variations. This is illustrated in the following for the case of Nb and Ta distribution between molten metal and silicate, this pair of elements giving important clues on oxygen fugacity prevailing during early Earth history.

#### 2.2.1. Oxygen fugacity effect on metal-Silicate partitioning of Nb and Ta

The variation of the partitioning with the oxygen fugacity ($f_{O_2}$) of any cation in silicate liquid informs on its electronic charge in the silicate melt. In fact, metal-silicate partitioning of an element M between silicate and metal involves exchange between oxidized species ($M O_{n/2}$) in the silicate and reduced species in the metal (M). This exchange reaction must depend on oxygen fugacity and the oxidation state n of the element M in the silicate melt as reported in the following equation:

$$M O_{n/2} = M + n/4 O_2 \quad (2)$$

Working with trace elements, some simplifying assumptions can be made to express their activities and the equilibrium constant can be approximated to the Nernst partition coefficient, yielding the following equation:

$$D_M = \frac{X_M^{\text{alloy}}}{X_M^{\text{silicate}}} \approx \frac{a_M^{\text{alloy}}}{a_M^{\text{silicate}}} = \exp\left(-\frac{\Delta G^0}{RT}\right) - \frac{n}{4}(f_{O_2}) \quad (3)$$

with X representing the weight fractions, a the activities, $\Delta G^0$ the Gibbs free energy change at the standard state and $f_{O_2}$ the oxygen fugacity. According to the Eq. (3) the metal-silicate partition coefficient for any cation is inversely related to the oxygen fugacity. Additionally, it is also a function of the element valence state in the silicate melt. It can also be predicted that, in theory, any element should become siderophile, provided that the conditions are sufficiently reducing. The most reducing conditions recorded in the solar system are those of calcium aluminium refractory inclusions (CAIs) of carbonaceous chondrites, about 7 log units below Fe-FeO equilibrium (IW-7) (Grossman et al., 2008). (Note that enstatite chondrite conditions may have been even lower, of the order of IW-10 (Larimer and Buseck, 1974).

#### 2.2.2. The "Nb paradox"

As the Earth accreted and differentiated, metallic materials segregated from silicates to form the iron-rich core. Because they were excluded from the core, the proportions of the refractory lithophile elements left over in the silicate part of the Earth are thought to have remained identical to those of chondritic meteorites throughout accretion (Kargel and Lewis, 1993). Niobium and tantalum, two refractory elements, have been classified as high field strength elements (HFSE), because their ionic charge over ionic radius ratio is high. This electronic feature mostly defines Nb and Ta behaviour in magmatic systems: as their sizes differ significantly from that of the available lattices in most of the common minerals, both elements tend to partition into the liquid phase during melting (Hofmann, 1988). Moreover Nb and Ta are described as "geochemical twins": they are known to occur only in pentavalent state and octahedrally coordinated (VI-fold) in geological systems, and share the same ionic radius (0.64 Å, Shannon, 1976). For those reasons,
Nb and Ta should behave similarly during the extraction of the core. However, the Nb/Ta ratio of the bulk silicate Earth (BSE) is significantly lower (Nb/Ta<sub>(BSE)</sub> = 14 ± 0.3) than the chondritic value (Nb/Ta<sub>(chondrites)</sub> = 19.9 ± 0.6), as well as lunar rocks (Nb/Ta<sub>(Moon)</sub>) = 17.0 ± 0.8 (Münker et al., 2003)). This observation is quoted as the “niobium paradox”. To resolve this paradox, three main explanations have been proposed: the Earth has a hidden silicate reservoir characterized by a superchondritic Nb/Ta that originated from the subduction of refractory eclogite through Earth’s history (Rudnick et al., 2000); early storage of the late stage residual melt from the magma ocean (Nebel et al., 2010); the “missing Nb” has been sequestrated into the core (Cartier et al., 2014; Wade and Wood, 2001).

### 2.2.3. Nb and Ta partitioning between metal and silicate melts

Various experimental studies support the idea that the core may be the “missing Nb” reservoir. The results are based on Nb and Ta metal-silicate melts partition coefficients, which appear to be mainly dependent on oxygen fugacity, compared to pressure, temperature or silicate composition (Cartier et al., 2014). When plotted as a function of oxygen fugacity (Fig. 2), partition coefficients describe lines which slopes that depend directly on the cation valence state in the silicate melt.

Under moderately reducing conditions (IW to IW-3) data are consistent with Nb<sub>2</sub>O<sub>5</sub> and Ta<sub>2</sub>O<sub>5</sub> in the silicate melt, i.e. Nb and Ta are mainly 5+ in the silicate liquid. More surprisingly, Nb and Ta describe different slopes under very reducing conditions (for f<sub>O2</sub> > IW-3), consistent with NbO and Ta<sub>2</sub>O<sub>3</sub> in the silicate melts. This contrasting valence change of Nb and Ta under reducing conditions is responsible for the two elements decoupling in metal-silicate systems. One can note here that the observed behaviour above for Nb and Ta is not a general feature for all lithophile elements as for Sm and Nd for example we did not observe any slope break of D<sub>Sm</sub> and D<sub>Nd</sub> as a function of oxygen fugacity in the range of IW-5 and IW-1. The reported results for Sm and Nd for S-free systems are consistent with Sm and Nd being dissolved in the silicate melt as Sm<sup>3+</sup> and Nd<sup>3+</sup> in reducing conditions between IW-5 and IW-1 (Bouhifd et al., 2015). Thus, partitioning data is a good tool to reveal element valence states. However the uncertainty on the element valence state is quite large, because of experimental and analytical error, coupled with significant uncertainties on experimental samples intrinsic oxygen fugacity estimation. The slope is also a function of the interaction coefficient of the element which can vary as Si enters the metal at low f<sub>O2</sub> - this is another good reason to independently constrain the valence through spectroscopy with direct valence probing in the samples.

<figure><img src="image_75.png" /><figcaption>**Fig. 2.** Experimentally determined metal-silicate partition coefficients for Nb (a) and Ta (b) as a function of oxygen fugacity, expressed in log units relative to Fe-FeO (IW) equilibrium (open, gray and closed symbols; compilation after Cartier et al., 2014). The data have been obtained over a range of pressure between 0 and 20 GPa. The red lines delineate the domain where lithophile behaviour is observed from that of siderophile behaviour). The black lines illustrate the slope of data alignment with suggested valence state in the silicate melt (Nb: 5+ and 2+; Ta: 5+ and 3+). Yellow diamonds represent experiments, for which Nb and Ta valence states have been directly measured by x-ray spectroscopy (Cartier et al., 2015). In this case, estimated errors on valence states are ±0.6 and ±1.3 for Nb and Ta, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

### 2.2.4. Nb and Ta speciation in silicate melts

In a recent study, XANES spectra have been acquired in silicate glasses equilibrated under moderately and highly reducing conditions (Fig. 2) (Cartier et al., 2015). Under moderately reducing conditions, the studied silicate glasses mainly contain Nb and Ta in the 5+ state. Under highly reducing conditions, silicate glasses contain Nb<sup>2+</sup> and Ta<sup>3+</sup>. These results confirm the information about Nb and Ta valence obtained with partition coefficients. Thus Nb and Ta behaviour in metal-silicate melts systems is mainly controlled by redox conditions, which simplifies the modeling of their partition coefficients as a function of oxygen fugacity.

### 2.2.5. Nb/Ta ratio as an fO₂ proxy during core formation

The results presented here show that Nb/Ta ratio in the primitive magma ocean depends mainly on the prevailing oxygen fugacity. Additionally, the experiments show that pressure has no effect (Cartier et al., 2014), meaning that Nb-Ta systematics cannot be used to estimate the depth of the magma ocean, in contrast to Ni and Co for example. Nb-Ta systematics, however, can be used to shed light on the oxidation state of the planet primitive building blocks.

For this, planet growth from planetesimals has been modeled in light of Nb/Ta evolution in the magma ocean during core formation. As noted above, when observed, Nb depletion is likely due to its trapping in the metallic core. Therefore, present day Earth's mantle bears the evidence of a reduced stage at some early stage of our planet history. By applying the metal-silicate partitioning data to a heterogeneous accretion scenario (such as proposed by Rubie et al., 2011) it is possible to reach the observed subchondritic Nb/Ta ratio (14 ± 0.3) of the bulk silicate Earth (Fig. 3). According to this scenario, the first ~70 wt% of the Earth accreted from reduced material, similar to enstatite chondrites, except for its lower volatile content, yielding extremely low Nb/Ta in the silicate mantle. The last 30% consist of more oxidized material, allowing for Nb/Ta to increase, until it reaches its present day value.

In the case of the Moon, the Nb/Ta estimate is of 17.0 ± 0.8 (Münker et al., 2003). It is higher than the Earth's value, but still subchondritic. Noting that the Moon's core is very small and has a negligible effect on Nb budget, the Nb/Ta of the Moon mantle is necessarily the result of the mixing of the proto-Earth mantle with the impactor Theia. Considering that terrestrial and lunar samples are isotopically similar, a successful scenario is that at the time of impact, the Moon may have inherited the chemical signature of the Earth's upper mantle, whose Nb/Ta was higher than the modern value due to the crystallization of bridgmanite. This signature would have been subsequently erased due to the homogenization of the whole Earth's mantle. Nb/Ta ratio is