fraction is a strong function of the liquidus, this discontinuity happens to affect the precise depth at which the last drop of melt should solidify. However, the last part to solidify should still be located in the shallow mantle, regardless the shape of the liquidus in this region. In addition, the magma ocean lifetime $t_{MO}$ appears to be weakly affected by the choice of the model and ranges between 147 kyr (in the A-chondritic case) to 171 kyr (in the F-peridotitic case).

Still, the F-peridotitic and A-chondritic models show a significant difference on the evolution of the core temperature (Fig. 4, right). For both cases, the core temperature decreases asymptotically from its initial value to a value that ranges between 4600 K (A-chondritic case) and 4860 K (F-peridotitic) in about 0.15 Myr. Since the lowermost mantle solidifies more rapidly when considering a hotter liquidus, a F-peridotitic model for the magma ocean helps to retain some heat in the core. For both cases, the final core temperature is 100 to 150 K larger than the temperature at which the melt fraction of the lowermost mantle reaches the 40% critical value.

### 3.4. Influence of the bottom thermal boundary layer

The thickness of the bottom thermal boundary layer $e_{TBL}$ governs both the cooling rate of the core and the energy supplied to the magma ocean (Eq. (21)). In a hard-turbulent context, this thickness is difficult to determine from laboratory experiments and only theoretical models can constrain this parameter (Spiegel, 1971). Therefore, a rough estimation of $e_{TBL}$ can be obtained for $Ra \sim 10^{20}$ based on the Nusselt number calculation: $Nu \sim (RaPr)^{1/2}$ (Spiegel, 1971). Considering that $Nu \sim L/e_{TBL}$, we obtain $e_{TBL} \sim L(RaPr)^{-1/2}$. In our magma ocean context and assuming that the scaling law is still valid at Rayleigh numbers up to $10^{30}$, this corresponds to a value $e_{TBL} \sim 10^{-6}$ m. This thickness is extremely thin, however, it is compatible with previous reports (Solomatov, 2007; Lebrun et al., 2013). As a first attempt to understand the influence of bottom thermal boundary layer in the magma ocean cooling dynamics, we performed various calculations with $e_{TBL}$ fixed and ranging from 1 mm to 1 km. The lower range values can be seen as unrealistic in comparison with, for example, the core topography. However, the wide range of values considered for $e_{TBL}$ allows a better illustration of its influence on the magma ocean cooling timescales.

We find that the magma ocean lifetime does not depend largely on $e_{TBL}$ and its value remains close to 150 kyr for the whole range of $e_{TBL}$ values assumed here. Fig. 5 shows the CMB and the surface heat flows for four different values of $e_{TBL}$. For $e_{TBL} = 1000$ m (Fig. 5, b), the surface heat flow is several orders of magnitude larger than the core heat flow and the thermal coupling between these two reservoirs is inefficient. For values of thermal boundary layer thicknesses below 1 m (Fig. 5, c and d), the heat flows become comparable and the thermal coupling between the core and the molten overlying mantle becomes efficient. Ideally, $e_{TBL}$ should be an adjustable parameter in our calculation, related to the effective value of the Rayleigh number. However, this would result in extremely small $e_{TBL}$ values of the order of $10^{-6}$ m. We show in Fig. 5 that decreasing the value of $e_{TBL}$ from 100 mm to a value of 1 mm does influence neither the shape of the surface heat flow nor the magma ocean lifetime (Fig. 5, c and d). In that cases, the core heat flow rapidly reaches a maximal value and decreases within the first 20 kyr, which corresponds to the time needed by the melt fraction of the bottom of the magma ocean to reach the critical value of 40% (Fig. 6). During this brief period of time, although heat is efficiently removed from the core (Fig. 6), the CMB heat flow remains considerably smaller than the surface heat loss. This is the reason why $e_{TBL}$ does not significantly influence the magma ocean lifetime. Hence, decreasing $e_{TBL}$ to a smaller value

<figure><img src="image_285.png" /><figcaption>**Fig. 6.** Time evolution of the core temperature for different thicknesses of the bottom thermal boundary layer (assuming a A-chondritic model, $\eta_m = 100$ Pa s and $T_0^{core} = 5000$ K).</figcaption></figure>

(i.e. $e_{TBL} \ll 10^{-6}$ m) as suggested by our theoretical estimation of the TBL thickness or to relate $e_{TBL}$ to the Rayleigh number that is initially extremely high should not influence the results of our study.

Such a behavior is confirmed by the time evolution of $T_0^{core}$, which is a strong function of $e_{TBL}$ (Fig. 6). As long as $e_{TBL}$ is larger than 100 m, the initial core heat is efficiently retained and the core cooling is not influenced by the cooling of the overlying magma ocean. However, for $e_{TBL} \le 100$ m the thermal coupling between the core and the MO becomes important. For $e_{TBL} < 1$ m, the core rapidly cools down to ~4400 K, which corresponds to the temperature where the lowermost mantle reaches the critical melt fraction value of 40%. Then the core cooling efficiency strongly decreases as the lowermost mantle is becoming much more viscous. Finally, the core temperature ends up at a temperature of ~4370 K for $e_{TBL}$ values ranging between 1 mm to 10 cm.

### 3.5. Influence of the magma ocean viscosity

Measurements (Liebske et al., 2005) and ab initio calculations (Karki and Stixrude, 2010) estimate that the dynamic viscosity $\eta_m$ of peridotitic melt is in the range $10^{-2}-10^{-1}$ Pa s. At low degrees of partial melting of a peridotite, the viscosity of the generated liquid can eventually increase up to 100 Pa s (Kushiro, 1986). The viscosity of molten mafic silicate should range between $10^{-2}-10^2$ Pa s (Rubie et al., 2003). To take into account the effect of this uncertainty on the magma ocean lifetime, we perform numerical simulations considering that the fully molten magma ocean viscosity $\eta_m$ ranges between $10^{-2}-10^2$ Pa s.

Dimensional analysis of Eq. (1) indicates that $t_{MO}$ is inversely proportional to $F_{conv}$. In the hard turbulent regime relevant to a thick MO context this term scales as $\eta_m^{-3/7}$. Consequently, the lifetime of a magma ocean should scale as $\eta_m^{3/7}$. This is confirmed by our numerical results (Fig. 7), and consistent with previous work (Solomatov, 2007):

$$t_{MO}(\text{Myr}) = 0.018 \eta_m^{3/7} \qquad (23)$$

Most importantly, for realistic viscosities of the fully molten early mantle, the melt fraction drops below 50% at all mantle depths in less than 1 Myr. For the lower range of $\eta_m$, this characteristic time scale can decrease down to several kyrs, rather than 1 Myr.

## 4. Conclusion

The cooling of a thick terrestrial magma ocean is a fast process. The magma ocean lifetime is principally governed by its viscosity
<figure><img src="image_286.png" /><figcaption>**Fig. 7.** Magma ocean lifetime as a function of the magma ocean viscosity η_m (assuming a A-chondritic model, *e*<sub>TBL</sub> = 1 m and *T*<sub>0</sub><sup>core</sup> = 5000 K).</figcaption></figure>

and, for the mid range of realistic viscosities, the magma ocean reaches a 50% melt fraction at all mantle depths within 20 kyr. Depending on the thermal boundary layer at the CMB, the thermal coupling between the core and MO can either insulate the core during the MO solidification and favor a hot core (for thick TBL), or drain the heat out of the core. However, as suggested by theoretical calculations, an extremely thin thermal boundary layer in an ultraturbulent environment implies that much of the heat is removed from the core during the cooling of the overlying mantle. In this context, our F-peridotitic model (with a higher liquidus) is more willing to retain the core heat than our A-chondritic model, yielding a difference of temperature of ~170 K after the magma ocean has crystallized. In addition, the final core temperature increases by a couple hundred degrees as its early temperature increases. However, for all cases, the average core temperature at the CMB (*T*<sub>core</sub>) ends up close to the 40% melt fraction temperature of the silicate magma ocean: ±100 K depending on the initial core temperature and on the thickness of the bottom thermal boundary layer.

In contrast with previous reports (Labrosse et al., 2007), our model shows that the crystallization occurs relatively rapidly at the CMB and, after some crystallization has proceeded, the highest amount of partial melting is found at intermediate depth between the surface and the CMB (see Fig. 3). One could argue that this result is dictated by the fact that our calculation neglects the possible effects of chemical segregation during mantle cooling. On the contrary, we believe that the segregation of a melt above the core mantle boundary would not help to retain heat in the core (Labrosse, 2015; Davies et al., 2015). A melt that would accumulate just above the CMB by gravitational segregation would be depleted in refractory elements, thus with a liquidus lower than the rest of the mantle. The temperature at which this part of the mantle becomes viscous (at 40% of partial melting) would be lowered and thus the core heat would escape more easily in the presence of a basal magma ocean (Ulvrová et al., 2012; Nakagawa and Tackley, 2014), in agreement with our results displayed in Fig. 4.

These results have important consequences for the magnetic history of the Earth. Indeed, if at some point, a full magma ocean has existed on the Earth, it is likely that most of the core heat has been removed rapidly. In less than 1 Myr, a tremendous heat flow may have lead to a significant decrease of the core temperature until it reached a value that is close or slightly above the temperature at which the melt fraction of the lowermost mantle reaches the 40% critical value (i.e. ~4400 K). Assuming that the Earth-Moon system was formed by a giant impact 100 Myrs after the first solids of the Solar System (Kleine and Rudge, 2011),

and that this giant impact has completely molten the Earth's mantle (Nakajima and Stevenson, 2015), it is difficult to envision that a large amount of heat could be retained in the core to sustain the geomagnetic field by thermal convection for several Gyr (Andrault et al., 2016), in contrast with the proposition of the most recent reports (Labrosse, 2015; Davies et al., 2015).

Finally, we acknowledge that our current model neglects the effect of vertical chemical segregation. While this effect is unlikely to dominate the dynamics in a highly turbulent magma ocean (Tonks and Melosh, 1990), it could become more important when the degree of partial melting becomes close or lower than ~40% (i.e. when viscosity increases). At this point, the knowledge of the density contrast between the solid at the liquidus (the first crystal to form) and the ambient liquid becomes of major importance. Whether the melt sinks, or floats has important ramifications for understanding the first steps in the dynamic modeling of the Earth's differentiation. In the near future, a modeling effort to integrate the compositional contribution in the buoyancy calculation between liquid and solid will constitute an important step forward towards the understanding of the earliest stages of Earth's evolution.

## Acknowledgements

The authors thank J. Roberts and anonymous reviewer for thoughtful and constructive comments. The authors also thank G. Tobie and B. Cecconi for their useful help in the development of the numerical model as well as M. Le Bars and J. Lunine for fruitful discussions and help. J. Monteux and D. Andrault are funded by Agence Nationale de la Recherche (Oxydeep decision no. ANR-13-BS06-0008). J. Monteux is also funded by the Auvergne Fellowship program. H. Samuel acknowledges the support from the Deutsche Forschungsgemeinschaft (project SA 2042/3). This is Laboratory of Excellence ClerVolc contribution no. 205.

## References

* Abe, Y., 1997. Thermal and chemical evolution of the terrestrial magma ocean. Phys. Earth Planet. Inter. 100, 27-39. http://dx.doi.org/10.1016/S0031-9201(96)03229-3.

* Agnor, C.B., Canup, R.M., Levison, H.F., 1999. On the character and consequences of large impacts in the late stage of terrestrial planet formation. Icarus 142, 219-237. http://dx.doi.org/10.1006/icar.1999.6201.

* Allègre, C.J., Poirier, J.-P., Humler, E., Hofmann, A.W., 1995. The chemical composition of the Earth. Earth Planet. Sci. Lett. 134, 515-526. http://dx.doi.org/10.1016/0012-821X(95)00123-T.

* Andrault, D., Bolfan-Casanova, N., Nigro, G.L., Bouhifd, M.A., Garbarino, G., Mezouar, M., 2011. Solidus and liquidus profiles of chondritic mantle: implication for melting of the Earth across its history. Earth Planet. Sci. Lett. 304, 251-259. http://dx.doi.org/10.1016/j.epsl.2011.02.006.

* Andrault, D., Petitgirard, S., Lo Nigro, G., Devidal, J.-L., Veronesi, G., Garbarino, G., Mezouar, M., 2012. Solid-liquid iron partitioning in Earth's deep mantle. Nature 487, 354-357. http://dx.doi.org/10.1038/nature11294.

* Andrault, D., Monteux, J., Le Bars, M., Samuel, H., 2016. The deep Earth may not be cooling down. Earth Planet. Sci. Lett. 443, 195-203. http://dx.doi.org/10.1016/j.epsl.2016.03.020.

* Boyet, M., Carlson, R.W., 2005. <sup>142</sup>Nd evidence for early (>4.53 Ga) global differentiation of the silicate Earth. Science 309, 576-581. http://dx.doi.org/10.1126/science.1113634.

* Canup, R.M., 2012. Forming a Moon with an Earth-like composition via a giant impact. Science 338, 1052. http://dx.doi.org/10.1126/science.1226073.

* Carlson, R.W., Lugmair, G.W., 2000. Timescales of Planetesimal Formation and Differentiation Based on Extinct and Extant Radioisotopes. Ariz. Univ. Press, Tucson, pp. 25-44.

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

* Čuk, M., Stewart, S.T., 2012. Making the Moon from a fast-spinning Earth: a giant impact followed by resonant despinning. Science 338, 1047. http://dx.doi.org/10.1126/science.1225542.

* Davies, C., Pozzo, M., Gubbins, D., Alfè, D., 2015. Constraints from material properties on the dynamics and evolution of Earth's core. Nat. Geosci. 8, 678-685. http://dx.doi.org/10.1038/ngeo2492.
Dziewonski, A.M., Anderson, D.L., 1981. Preliminary reference Earth model. Phys. Earth Planet. Inter. 25, 297–356. http://dx.doi.org/10.1016/0031-9201(81)90046-7.

Elkins-Tanton, L.T., Parmentier, E.M., Hess, P.C., 2003. Magma ocean fractional crystallization and cumulate overturn in terrestrial planets: implications for Mars. Meteorit. Planet. Sci. 38, 1753–1771.

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.

Fiquet, G., Auzende, A.L., Siebert, J., Corgne, A., Bureau, H., Ozawa, H., Garbarino, G., 2010. Melting of peridotite to 140 gigapascals. Science 329, 1516. http://dx.doi.org/10.1126/science.1192448.

Ghosh, A., McSween, H.Y., 1998. A thermal model for the differentiation of asteroid 4 Vesta, based on radiogenic heating. Icarus 134, 187–206. http://dx.doi.org/10.1006/icar.1998.5956.

Hamano, K., Abe, Y., Genda, H., 2013. Emergence of two types of terrestrial planet on solidification of magma ocean. Nature 497, 607–610. http://dx.doi.org/10.1038/nature12163.

Herzberg, C., Zhang, J., 1996. Melting experiments on anhydrous peridotite klb-1: compositions of magmas in the upper mantle and transition zone. J. Geophys. Res., Solid Earth 101 (B4), 8271–8295. http://dx.doi.org/10.1029/96JB00170.

Javoy, M., Kaminski, E., Guyot, F., Andrault, D., Sanloup, C., Moreira, M., Labrosse, S., Jambon, A., Agrinier, P., Davaille, A., Jaupart, C., 2010. The chemical composition of the Earth: enstatite chondrite models. Earth Planet. Sci. Lett. 293, 259–268. http://dx.doi.org/10.1016/j.epsl.2010.02.033.

Karato, S.-I., Wu, P., 1993. Rheology of the upper mantle – a synthesis. Science 260, 771–778. http://dx.doi.org/10.1126/science.260.5109.771.

Karki, B.B., Stixrude, L.P., 2010. Viscosity of MgSiO₃ liquid at Earth's mantle conditions: implications for an Early magma ocean. Science 328, 740. http://dx.doi.org/10.1126/science.1188327.

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

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

Kleine, T., Rudge, J.F., 2011. Chronometry of meteorites and the formation of the Earth and Moon. Elements 7, 41–46. http://dx.doi.org/10.1016/j.gca.2006.06.004.

Kleine, T., Touboul, M., Bourdon, B., Nimmo, F., Mezger, K., Palme, H., Jacobsen, S.B., Yin, Q-Z., Halliday, A.N., 2009. HF-W chronology of the accretion and early evolution of asteroids and terrestrial planets. Geochim. Cosmochim. Acta 73, 5150–5188. http://dx.doi.org/10.1016/j.gca.2008.11.047.

Kushiro, I., 1986. Viscosity of partial melts in the upper mantle. J. Geophys. Res. 91, 9343–9350. http://dx.doi.org/10.1029/JB091iB09p09343.

Labrosse, S., 2015. Thermal evolution of the core with a high thermal conductivity. Phys. Earth Planet. Inter. 247, 36–55. http://dx.doi.org/10.1016/j.pepi.2015.02.002.

Labrosse, S., Hernlund, J.W., Coltice, N., 2007. A crystallizing dense magma ocean at the base of the Earth's mantle. Nature 450, 866–869. http://dx.doi.org/10.1038/nature06355.

Lebrun, T., Massol, H., Chassefière, E., Davaille, A., Marcq, E., Sarda, P., Leblanc, F., Brandeis, G., 2013. Thermal evolution of an early magma ocean in interaction with the atmosphere. J. Geophys. Res., Planets 118, 1155–1176. http://dx.doi.org/10.1002/jgre.20068.

Liebsche, C., Schmickler, B., Terasaki, H., Poe, B.T., Suzuki, A., Funakoshi, K.-i., Ando, R., Rubie, D.C., 2005. Viscosity of peridotite liquid up to 13 GPa: implications for magma ocean viscosities [rapid communication]. Earth Planet. Sci. Lett. 240, 589–604. http://dx.doi.org/10.1016/j.epsl.2005.10.004.

Martin, H., Moyen, J.-F., Guitreau, M., Blichert-Toft, J., Le Pennec, J.-L., 2014. Why Archaean TTG cannot be generated by MORB melting in subduction zones. Lithos 198, 1–13. http://dx.doi.org/10.1016/j.lithos.2014.02.017.

Miller, G.H., Stolper, E.M., Ahrens, T.J., 1991. The equation of state of a molten komatiite: 2. Application to komatiite petrogenesis and the Hadean Mantle. J. Geophys. Res. 96, 11,849. http://dx.doi.org/10.1029/91JB01203.

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.

Monteux, J., Jellinek, A.M., Johnson, C.L., 2011. Why might planets and moons have early dynamos? Earth Planet. Sci. Lett. 310, 349–359. http://dx.doi.org/10.1016/j.epsl.2011.08.014.

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

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

Moyen, J.-F., Martin, H., 2012. Forty years of TTG research. Lithos 148, 312–336. http://dx.doi.org/10.1016/j.lithos.2012.06.010.

Nakagawa, T., Tackley, P.J., 2014. Influence of combined primordial layering and recycled MORB on the coupled thermal evolution of Earth's mantle and core. Geochem. Geophys. Geosyst. 15, 619–633. http://dx.doi.org/10.1002/2013GC005128.

Nakajima, M., Stevenson, D.J., 2015. Melting and mixing states of the Earth's mantle after the Moon-forming impact. Earth Planet. Sci. Lett. 427, 286–295. http://dx.doi.org/10.1016/j.epsl.2015.06.023.

Neumann, W., Breuer, D., Spohn, T., 2014. Differentiation of Vesta: implications for a shallow magma ocean. Earth Planet. Sci. Lett. 395, 267–280. http://dx.doi.org/10.1016/j.epsl.2014.03.033.

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

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

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

Ringwood, A.E., 1966. Chemical evolution of the terrestrial planets. Geochim. Cosmochim. Acta 30, 41–104. http://dx.doi.org/10.1016/0016-7037(66)90090-1.

Rizo, H., Boyet, M., Blichert-Toft, J., Rosing, M.T., 2013. Early mantle dynamics inferred from $^{142}$Nd variations in Archean rocks from southwest Greenland. Earth Planet. Sci. Lett. 377, 324–335. http://dx.doi.org/10.1016/j.epsl.2013.07.012.

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

Safronov, V.S., 1978. The heating of the earth during its formation. Icarus 33, 3–12. http://dx.doi.org/10.1016/0019-1035(78)90019-2.

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

Shuvalov, V., 2009. Atmospheric erosion induced by oblique impacts. Meteorit. Planet. Sci. 44, 1095–1105. http://dx.doi.org/10.1111/j.1945-5100.2009.tb01209.x.

Simon, F., Glatzel, 1929. Fusion-pressure curve. Z. Anorg. Allg. Chem. 178, 309.

Solomatov, V., 2007. Magma Oceans and Primordial Mantle Differentiation. Schubert, G. (Ed.), Treatise of Geophysics, vol. 9. Elsevier.

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

Spiegel, E.A., 1971. Convection in stars: I. Basic Boussinesq convection. Annu. Rev. Astron. Astrophys. 9, 323. http://dx.doi.org/10.1146/annurev.aa.09.090171.001543.

Stixrude, L., de Koker, N., Sun, N., Mookherjee, M., Karki, B.B., 2009. Thermodynamics of silicate liquids in the deep Earth. Earth Planet. Sci. Lett. 278, 226–232. http://dx.doi.org/10.1016/j.epsl.2008.12.006.

Svetsov, V.V., 2005. Numerical simulations of very large impacts on the Earth. Space Sci. 53, 1205–1220. http://dx.doi.org/10.1016/j.pss.2005.04.011.

Tackley, P.J., 2012. Dynamics and evolution of the deep mantle resulting from thermal, chemical, phase and melting effects. Earth-Sci. Rev. 110, 1–25. http://dx.doi.org/10.1016/j.earscirev.2011.10.001.

Thomas, C.W., Asimow, P.D., 2013. Direct shock compression experiments on pre-molten forsterite and progress toward a consistent high-pressure equation of state for CaO-MgO-Al₂O₃-SiO₂-FeO liquids. J. Geophys. Res., Solid Earth 118, 5738–5752. http://dx.doi.org/10.1002/jgrb.50374.

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

Tonks, W.B., Melosh, H.J., 1990. The Physics of Crystal Settling and Suspension in a Turbulent Magma Ocean. Oxford University Press, pp. 151–174.

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

Tosi, N., Plesa, A.-C., Breuer, D., 2013. Overturn and evolution of a crystallized magma ocean: a numerical parameter study for Mars. J. Geophys. Res., Planets 118, 1512–1528. http://dx.doi.org/10.1002/jgre.20109.

Touboul, M., Puchtel, I.S., Walker, R.J., 2012. W-182 evidence for long-term preservation of early mantle differentiation products. Science 335 (6072), 1065–1069. http://dx.doi.org/10.1126/science.1216351.

Ulvrová, M., Labrosse, S., Coltice, N., Ráback, P., Tackley, P.J., 2012. Numerical modelling of convection interacting with a melting and solidification front: application to the thermal evolution of the basal magma ocean. Phys. Earth Planet. Inter. 206, 51–66. http://dx.doi.org/10.1016/j.pepi.2012.06.008.

Čižková, H., van den Berg, A.P., Spakman, W., Matyska, C., 2012. The viscosity of Earth's lower mantle inferred from sinking speed of subducted lithosphere. Phys. Earth Planet. Inter. 200, 56–62. http://dx.doi.org/10.1016/j.pepi.2012.02.010.

Yoshino, T., Walter, M.J., Katsura, T., 2003. Core formation in planetesimals triggered by permeable flow. Nature 422, 154–157.
# A mushy Earth's mantle for more than 500 Myr after the magma ocean solidification

J. Monteux,<fnref n="1" /> D. Andrault,<fnref n="1" /> M. Guitreau,<fnref n="1" /> H. Samuel and S. Demouchy<fnref n="1" />,<fnref n="3" />

Accepted 2020 January 24. Received 2020 January 22; in original form 2019 July 17

## SUMMARY

In its early evolution, the Earth mantle likely experienced several episodes of complete melting enhanced by giant impact heating, short-lived radionuclides heating and viscous dissipation during the metal/silicate separation. After a first stage of rapid and significant crystallization (Magma Ocean stage), the mantle cooling is slowed down due to the rheological transition, which occurs at a critical melt fraction of 40–50%. This transition first occurs in the lowermost mantle, before the mushy zone migrates toward the Earth’s surface with further mantle cooling. Thick thermal boundary layers form above and below this reservoir. We have developed numerical models to monitor the thermal evolution of a cooling and crystallizing deep mushy mantle. For this purpose, we use a 1-D approach in spherical geometry accounting for turbulent convective heat transfer and integrating recent and solid experimental constraints from mineral physics. Our results show that the last stages of the mushy mantle solidification occur in two separate mantle layers. The lifetime and depth of each layer are strongly dependent on the considered viscosity model and in particular on the viscosity contrast between the solid upper and lower mantle. In any case, the full solidification should occur at the Hadean–Eoarchean boundary 500–800 Myr after Earth’s formation. The persistence of molten reservoirs during the Hadean may favor the absence of early reliefs at that time and maintain isolation of the early crust from the underlying mantle dynamics.

**Key words:** Numerical modelling; Heat flow; Rheology: Mantle; Dynamics: convection currents, and mantle plumes; Heat generation and transport.

## 1. INTRODUCTION

After the giant impact which led to the formation of the Earth–Moon system, the Earth’s mantle was likely completely molten (e.g. Nakajima & Stevenson 2015). During the subsequent Magma Ocean (MO) stage, the Earth’s early mantle undergone a rapid global cooling until its melt fraction decrease to a critical value $\varphi_c$ (≈40%) associated with a major increase of its viscosity (Solomatov 2007; Monteux et al. 2016). This step was followed by a slow cooling stage that triggered a complete crystallization of the silicate mantle (Solomatov 2007). This second cooling stage could have lasted hundreds of Ma, or even a couple of Ga as its dynamics was governed by the rheology of the slowly deforming solid-like mantle, in contrast to the first one which was driven by the magma viscosity (Solomatov 2007; Ulvrová et al. 2012; Monteux et al. 2016).

Recent experimental results have shown that the upper mantle solidus is at lower temperature than previously expected for a chondritic composition (Andrault et al. 2018). According to this study, such a solidus associated with a hotter earlier mantle would enable the presence of a deep and persistent molten layer in the

Archean mantle. The progressive solidification of this melt layer could have enhanced the mechanical coupling between the lithosphere and the asthenosphere. Such a change might explain the transition from surface dynamics dominated by a stagnant lid to modern plate tectonics with deep-slab subductions. Assuming that the intersect between the mantle solidus and an adiabat temperature profile with a potential surface temperature corresponding to a surface melt fraction of $\varphi_c \approx 40\%$ could constrain the depth of the bottom of the remaining partially molten layer, Solomatov (2007) obtained a depth $\approx 300$ km. In the case of bottom-up solidification of a mushy mantle, the depth of the last remaining partially molten layer should, hence, be smaller than 300 km and the full mantle crystallization should occur within $\approx 1$ Ga (Sleep et al. 2014).

The depth at which the full crystallization is reached is likely governed by (1) the solidus temperature which controls the depth and temperature of solidification, (2) the thickness of the top thermal boundary layer (TBL), where heat convecting from the deep mantle is transferred to the surface by conduction and (3) the temperature contrast on both sides of the TBL which controls the efficiency of heat evacuation. In the mushy regime, the mantle viscosity is a

<footnote n="1">¹Université Clermont Auvergne, CNRS, IRD, OPGC, Laboratoire Magmas et Volcans, F-63000 Clermont-Ferrand, France. E-mail: julien.monteux@uca.fr</footnote>

<footnote n="2">²Institut de Physique du Globe de Paris, Université Sorbonne Paris Cité, 75005 Paris, France</footnote>

<footnote n="3">³Géosciences Montpellier – Université Montpellier & CNRS, F-34095 Montpellier, France</footnote>
key parameter, which governs mantle dynamics and the formation of TBLs. This day viscosity of the Earth's mantle is difficult to constrain, particularly in the lower mantle (Čižková et al. 2012), and the value of the viscosity of the Earth's early mantle prior to the occurrence of major differentiation events is even more hypothetical. Constraining the value of this physical quantity is of first importance, since deep mantle viscosity likely governs the initiation of major geological features, such as plate tectonics and mantle plumes (Sleep 2014; Foley et al. 2014).

We have developed a numerical model to monitor the cooling, and crystallization of an isochemical mantle, starting from a partially molten stage. We aim to characterize the influence of the viscosity of the solid early mantle on its cooling dynamics. In particular, we tested the impact of variation in activation energy and the viscosity prefactor, as well as monitored the cooling and crystallization processes to determine the solidification timescales and the depths at which the last melt fraction solidifies.

## 2. MODEL

We considered the cooling of a partially molten magma ocean with an initial depth of 2900 km, which we modelled using a thermal evolution described below.

### 2.1 Thermal evolution model

#### 2.1.1 Heat transfer model

We used a 1-D spherical approach (e.g. Abe 1997; Laneuville et al. 2018) accounting for turbulent convective heat transfer (e.g. Abe 1997; Monteux et al. 2016; Bower et al. 2018). This approach is relevant for the ranges of low viscosities and high Rayleigh numbers expected within a partially molten planetary mantle which are difficult to numerically resolve in 2-D and 3-D spatial domains. Indeed, even if Ra numbers are lower during the mushy stage than during the magma ocean stage, Ra numbers are still too high (up to $10^{30}$) to correctly solve the TBLs in such a dynamic reservoir. Moreover, molten reservoirs may survive during the early thermal evolution of the mushy mantle. In these regions, the melt fraction can reach values larger than 40% with very low associated viscosity of the mushy material, making the local Rayleigh number too high for 2-D or 3-D computational domains.

Our numerical model solves the following heat equation:

$$ \rho C_p \frac{\partial T}{\partial t} = \nabla \cdot (k \nabla T) + \rho H, \qquad (1) $$

where $\rho$ is the density, $C_p$ is the heat capacity, $T$ is the temperature and $H$ is the heat production from radiogenic sources. $k$ is an 'effective' conductivity defined as:

$$ k = k_c + k_v, \qquad (2) $$

with $k_v$ the effective conductivity relative to thermal convection of the mushy material is:

$$ k_v = \frac{F_{\text{conv}} L}{\Delta T}, \qquad (3) $$

$F_{\text{conv}}$ is the convective heat flux accounting for thermal buoyancy, $L$ is the thickness of the Earth's mantle and $k_c$ the intrinsic thermal conductivity of the material ($k_c = 5 \text{ W.m}^{-1}.\text{K}^{-1}$).

For mantle convection to occur, the temperature gradient must be larger than the adiabatic gradient:

$$ \left(\frac{dT}{dr}\right)_S = \frac{-\alpha g T}{C_p}. \qquad (4) $$

When the temperature gradient is subadiabatic, the mantle heat is transported only by conduction, and $k = k_c$. When the temperature gradient is superadiabatic, the mantle is convecting. The convective heat flux $F_{\text{conv}}$ depends on the local Rayleigh number $Ra$:

$$ Ra = \frac{\alpha g C_p \rho^2 \Delta T L^3}{k_c \eta}, \qquad (5) $$

where $\alpha$ is the thermal expansion coefficient of the mushy material, $g$ is the gravitational acceleration assumed to be constant through the whole mantle and $\eta$ is the local dynamic viscosity. In eq. (5), $k_c$ is constant, $C_p$ is a function of the melt fraction, $\alpha$ and $\rho$ vary with depth and melt fraction. $\eta$ varies with depth, temperature and melt fraction. $\Delta T$ is the thermal driving force for the convection, therefore the temperature difference between the surface and the core-mantle boundary (CMB) minored by the increase of temperature along the mantle adiabat. Mantle dynamics and cooling is also governed by the Prandtl number, $Pr$:

$$ Pr = \frac{\eta C_p}{k_c}, \qquad (6) $$

which is the ratio of the momentum diffusivity over the thermal diffusivity. $Pr$ is calculated at each depth using the local viscosity value $\eta$. Depending on the values of $Pr$ and $Ra$, two flow regimes arise and as a consequence two convective heat fluxes:

(1) a soft turbulent regime where the corresponding convective heat flux is (Solomatov 2007; Monteux et al. 2016):

$$ F_{\text{conv}} = \frac{0.089 k_c \Delta T R a^{1/3}}{L} \quad \text{if } Ra < 10^8 Pr^{5/3}; \qquad (7a) $$

(2) and a hard turbulent regime (following Solomatov 2007; Monteux et al. 2016) where:

$$ F_{\text{conv}} = \frac{0.22 k_c \Delta T R a^{2/7} Pr^{-1/7}}{L} \quad \text{if } Ra > 10^8 Pr^{5/3}. \qquad (7b) $$

In our numerical model, we compare the temperature gradient to the adiabatic gradient. If the temperature gradient is larger, $k_v$ is calculated according to eq. (3). If the temperature gradient is lower, $k = k_c$.

#### 2.1.2 Boundary and initial conditions

The large impacts, radiogenic heating and the energy dissipated during metal-silicate separation control the early thermal state of the core. The core temperature at the end of the magma ocean stage is governed by the heat accumulated in this reservoir during its formation but also by the efficiency of the magma ocean to retain heat within the core by forming thick TBLs (Monteux et al. 2016). The core heat flow at the CMB can be expressed as:

$$ F_{\text{core}} = \frac{k_c (T^{\text{core}} - T^{\text{mantle}})}{\delta_{\text{TBL,bot}}}, \qquad (8) $$

where $T^{\text{core}}$ is the average core temperature just below the CMB (i.e. here, a TBL within the core is not considered) and $T^{\text{mantle}}$ is the mantle temperature above the CMB. $\delta_{\text{TBL,bot}}$ is the thickness of the TBL at the bottom of the mantle where the heat is extracted from the core by conduction. $T^{\text{mantle}}$ is calculated from eq. (1) whereas $T^{\text{core}}$
is obtained by the integration of the following differential equation:

$$V_{\text{core}} \rho_{\text{Fe}} C_{p,\text{Fe}} \frac{dT^{\text{core}}}{dt} = S_{\text{core}} F_{\text{core}}, \quad (9)$$

where $V_{\text{core}}$ is the core volume, $S_{\text{core}}$ is the core surface, $\rho_{\text{Fe}}$ is the core density, $C_{p,\text{Fe}}$ is the core heat capacity. This formulation does not consider the increase in adiabatic temperature within the core but allows following the evolution of core temperature right below the CMB as a function of time, based on the CMB heat flux. The error associated with this simplification (which in the case of an Earth-like body would amount to less than 10%) is small compared to uncertainties in other model parameters.

The efficiency of mantle cooling also depends on the heat transfer at the surface. During the magma ocean cooling, heat is efficiently radiated toward space, but the formation of a primitive atmosphere may significantly slow down the cooling. With a primitive atmosphere composed of 300 bars H₂O and 100 bars CO₂ overlaying the magma ocean, Lebrun et al. (2013) estimated that the surface temperature remains constant throughout the entire duration of the mushy stage at $T_{\text{surf}} \approx 500$ K. Sleep et al. (2014) also showed that the surface temperature during the mushy stage was maintained at $\approx 500$ K for 1000 bars and 100 bars atmospheres in equilibrium with bulk silicated magmas.

### 2.1.3 TBLs parametrization

Monteux et al. (2016) showed that the initial core heat can only be efficiently retained within the core when the bottom TBL thickness ($\delta_{\text{TBL,bot}}$) is larger than $\approx 100$ m. For $\delta_{\text{TBL,bot}} < 100$ m the thermal coupling between the core and the MO is important and the core's heat is efficiently transferred to the mantle during the short timescales of the MO cooling. For $\delta_{\text{TBL,bot}} < 1$ m, the core rapidly cools down to $\approx 4400$ K, which corresponds to the core-mantle boundary temperature at a critical melt fraction $\varphi = \varphi_{\text{crit}}$. Actually, the thicknesses of both the bottom and top TBL are governed by the cooling dynamics of the mantle (Solomatov 1995). The formation of a TBL is induced by the velocity field decrease close to the boundary. In 1-D models, the TBL cannot form numerically by themselves because the velocity field is not calculated and, hence, it has to be parametrized. Moreover, as we consider a compressible fluid with properties changing with depth, the top and bottom TBL do not have the same thicknesses. In our models, we consider that the thermal boundary thickness scale as (Grott & Breuer 2008):

$$\delta_{\text{TBL}} = L \left( \frac{Ra_{\text{crit}}}{\overline{Ra}} \right)^{1/3}, \quad (10)$$

with $Ra_{\text{crit}} = 450$ (Choblet & Sotin 2000) and $\overline{Ra}$ is the Rayleigh number value calculated using eq. (5) and the characteristic parameters corresponding to either the top or the bottom of the magma ocean, for upper and lower TBL, respectively. It results in a TBL thicker at the top of the mantle, than at the bottom of the mantle. Strictly, the value of $Ra_{\text{crit}}$ should be different when considering the upper or lower TBL (Thiriet et al. 2019). The value of 450 is adapted for the upper mantle. For the lower mantle, it should be expressed as a function of the internal Rayleigh number (Deschamps & Sotin 2000). We have implemented such a parametrization in our models. Our numerical tests (not shown here) show that the bottom TBL thickness derived from Deschamps & Sotin (2000) leads to a decrease of the bottom TBL thickness, but does not change significantly the solidification depth and time scales. Hence, we used $Ra_{\text{crit}} = 450$ for both top and bottom TBL calculations. Within the TBLs, the heat is transferred by conduction and, again, $k = k_c$. In

all the models presented here, the TBL thickness has a thickness larger than 10 km, therefore given our 1-km grid spacing, at least 10 gridpoints are used to identify and characterize the heat transfer in the TBL.

## 2.2 Geochemical model and derived parameters

Despite the difficulty to characterize the chemical composition of the Earth's early mantle, a consensus has emerged that Earth's mantle should be of chondritic composition (e.g. Mc Donough & Sun 1995; Javoy et al. 2010; Palme & O'Neill 2014). Still, chondrites present a large diversity in major, minor and trace element compositions. According to several isotopic tracers, Earth has accreted from a large majority of building blocks typical of high-enstatite chondrites (EH, Javoy et al. 2010). Then, late-accretion processes and core-mantle differentiation (Rubie et al. 2011) have induced a drift of the bulk mantle composition to an MgO-enriched composition, compared to EH (e.g. Mc Donough & Sun 1995; Palme & O'Neill 2014). In the following section, we detail our chemical model for the mushy mantle following the magma ocean stage.

### 2.2.1 Radiogenic heating

Short timescales inferred for the duration of the magma ocean stage appear comparable to the timescales for the decay of short-lived radionuclides such as $^{26}$Al. However, large bodies such as the proto-Earth, or the Theia Earth's impactor, appeared only much later, from the accretion of pre-differentiated planetesimals. Therefore, only the long-lived radiogenic elements, such as $^{238}$U, $^{235}$U, $^{232}$Th and $^{40}$K can provide heat at the long time scale corresponding to the cooling of the mushy mantle. In our models, we consider the radiogenic heating from these radionuclides in eq. (1) assuming that the abundance of these elements in the primitive mantle is similar to the concentration in EH-chondrites (Javoy 1999). The radiogenic heat production rate is computed as:

$$H = \sum A_i [i] H_i \exp(-\lambda_i (t - t_0)), \quad (11)$$

where the meaning of $A_i$, $E_i$, $H_i$, $\lambda_i$ is detailed in Table 1.

### 2.2.2 Solidus and liquidus temperatures

The chemical composition of the mantle governs its melting properties, based on its solidus and liquidus profiles. The latter plays a major role in the early thermal evolution of the magma ocean, because it defines the temperature and the depth at which crystallization starts. Laboratory experiments have constrained the liquidus and solidus of mantle-like material up to pressures compatible with the CMB conditions (Fiquet et al. 2010; Andrault et al. 2011). We performed calculations using the melting curves derived from chondritic-type mantle composition from Andrault et al. (2011, Fig. 1). The experimental solidus and liquidus profiles are fitted with a modified Simon and Glatzel equation (Simon & Glatzel 1929). For pressures $P$ below 24 GPa, we use solidus and liquidus temperatures of chondritic mantle reported from Andrault et al. (2018):

$$T_{\text{sol}} = 1373 \left( \frac{P}{0.82 \times 10^9} + 1 \right)^{(1/6.94)}, \quad (12a)$$

$$T_{\text{liq}} = 1983.4 \left( \frac{P}{6.48 \times 10^9} + 1 \right)^{(1/5.35)}. \quad (12b)$$
<table><caption>**Table 1.** Radiogenic heat sources and characteristics (from Javoy 1999).</caption><thead><tr><th>Isotope/element</th><th>Half-life (λ<sub>i</sub>) (yr)</th><th>Heat production per unit mass of isotope (H<sub>i</sub>) (W kg<sup>-1</sup>)</th><th>Natural Abundance (A<sub>i</sub>) (%)</th><th>Present-day concentrations ([i]) (ppm wt.)</th></tr></thead><tbody><tr><td><sup>238</sup>U</td><td>4.46 × 10<sup>9</sup></td><td>9.17 × 10<sup>-5</sup></td><td>99.28</td><td></td></tr><tr><td><sup>235</sup>U</td><td>7.04 × 10<sup>8</sup></td><td>5.75 × 10<sup>-4</sup></td><td>0.72</td><td></td></tr><tr><td>U</td><td></td><td></td><td></td><td>0.20</td></tr><tr><td><sup>232</sup>Th</td><td>1.4 × 10<sup>11</sup></td><td>2.56 × 10<sup>-5</sup></td><td>100</td><td></td></tr><tr><td>Th</td><td></td><td></td><td></td><td>0.069</td></tr><tr><td><sup>40</sup>K</td><td>1.26 × 10<sup>9</sup></td><td>2.97 × 10<sup>-5</sup></td><td>0.0117</td><td></td></tr><tr><td>K</td><td></td><td></td><td></td><td>270</td></tr></tbody></table>

<figure><img src="image_287.png" /><figcaption>**Figure 1.** Solidus (green) and liquidus (red) (computed from eq. 12) as a function of the Earth's radius. Black solid lines show adiabatic temperature profiles (computed from eqs 4 and 18) with three different potential temperatures (1600, 1750 and 2000 K). Note that these three temperature profiles are arbitrary and do not result from our numerical model. The blue dashed line separates the upper and lower mantle. The dashed lines correspond to the hypothetic conductive temperature profiles in the top thermal boundary layer (i.e. recalculated at each step of our modelling procedure). The two ellipsoids illustrate critical points where deep and shallow last molten layers should solidify (SML, shallow mantle layer and DML, deep mantle layer).</figcaption></figure>

For pressures larger than P = 24 GPa, we use results from Andrault *et al.* (2011):

$$T_{\text{sol}} = 1334.5 \left( \frac{P}{9.63 \times 10^9} + 1 \right)^{(1/2.41)}, \qquad (12c)$$

$$T_{\text{liq}} = 1862 \left( \frac{P}{21.15 \times 10^9} + 1 \right)^{(1/2.15)}. \qquad (12d)$$

Mantle solidification may induce some chemical fractionation. In such case, the melting curves may evolve (Andrault *et al.* 2017). Major changes concern the liquidus temperature, which increases with the MgSiO₃-content in the mantle. On the other hand, the mantle's solidus is almost independent of composition. These effects are not accounted in our study, because there is an insufficient knowledge on melting properties as a function of pressure, temperature and mantle composition.

### 2.2.3 Thermodynamic parameters

Thermodynamic parameters of the magma ocean depend of its chemical composition. Volumetric and elastic parameters of silicate liquids have been characterized up to a pressure of 140 GPa using shock compression experiments (Mosenfelder *et al.* 2007, 2009;

Thomas *et al.* 2012; Thomas & Asimow 2013). Here we assume a multicomponent system with a chondritic-type composition (62% enstatite + 24% forsterite + 8% fayalite + 4% anorthite + 2% diopside). Using fourth-order Birch-Murnaghan/Mie-Grüneisen equation of state fits for molten silicate liquids from Thomas & Asimow (2013), we obtain the melt density ρ<sub>m</sub>, the volumetric thermal expansion α as a function of pressure as well as the specific heat C<sub>p</sub> of the molten material for a chondritic multicomponent assemblage. The density of the solid phase is then calculated as:

$$\rho_s = \rho_m + \Delta\rho, \qquad (13)$$

with Δρ being the density difference between solid and liquid phases which is fixed to 64 kg m⁻³ (Monteux *et al.* 2016).

For a partially molten material, the density ρ', the coefficient of volumetric thermal expansion α' and the specific heat C'<sub>p</sub> are given as follows (Solomatov 2007):

$$\frac{1}{\rho'} = \frac{1 - \varphi}{\rho_s} + \frac{\varphi}{\rho_m}, \qquad (14)$$

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

$$C'_{p} = C_{p} + \frac{\Delta H}{T_{\text{liq}} - T_{\text{sol}}}, \qquad (16)$$

where ΔH is the latent heat released during solidification, and ϕ is the melt fraction:

$$\varphi = \frac{T - T_{\text{sol}}}{T_{\text{liq}} - T_{\text{sol}}}. \qquad (17)$$

### 2.2.4 Adiabats

In vigorously convecting systems such as magma oceans, the temperature distribution is nearly adiabatic (Solomatov 2007). For one-phase systems, such as a completely molten or a completely solid layer, eq. (4) gives the equation for an adiabat. In two-phase systems (liquid + solid), the effects of phase changes need to be considered (Solomatov 2007). The equation for such adiabat follows:

$$\left(\frac{dT}{dr}\right)_s = \frac{-\alpha'gT}{C'_{p}}. \qquad (18)$$

This results in an increase of the adiabat gradient at depth where the two phases coexist, compared to the purely liquid or solid one-phase adiabats (Solomatov 2007). Fig. 1 compares three adiabatic temperature profiles and the melting curves used in our study. The adiabatic temperature profiles are calculated by numerical integration of eqs (4) and (18) using a fourth-order Runge-Kutta method
(Press et al. 1993). These adiabatic temperature profiles are used to calculate at each depth, and when it is super-adiabatic, the temperature difference ΔT from eq. (5).

### 2.2.5 Assumptions

During the solidification of the early mantle (i.e. the magma ocean and mushy mantle stages) chemical fractionation may occur between compatible and incompatible elements that partition preferentially into solid and liquid phases, respectively. Initially, the bridgmanite grains are denser than the liquid and they could fall toward the core-mantle boundary. Then, after a significant fraction of the MO is crystallized, the liquid could become denser as iron is a relatively incompatible element in mantle minerals. This could produce late mantle overturns (e.g. Boukaré et al. 2015). Such chemical differentiations could also induce heterogeneous distribution of radiogenic elements due to their incompatible behaviour. Still, the early chemical fractionation of the Earth's mantle history remains debated, based on contradictory geodynamical (Solomatov 2000) and geochemical (Mc Donough & Sun 1995; Boyet & Carlson 2005; Palme & O'Neill 2014) arguments. Therefore, we do not consider chemical differentiation in the solidifying mushy mantle in this study. Hence, melting curves, density and concentration of radiogenic elements are considered unchanged along the cooling process.

## 2.3 Viscosity model

Viscosity governs mantle cooling dynamics, which is strongly dependent on the melt fraction $\varphi$. In our study we consider that $\varphi$ is a linear function of the temperature difference between the liquidus and the solidus (eq. 17). In the following section, we detail the parametrization used to compute the viscosity in our numerical models.

### 2.3.1 Liquid fraction viscosity ($\varphi = 1$)

During the cooling of a mushy mantle, the melt fraction globally decreases; however, locally, mantle layers may remain largely molten for a long period of time. Therefore, a mushy mantle may locally be extremely turbulent because of the low viscosity of the molten mantle material (Cochain et al. 2017). For the fully molten mantle (i.e. when $\varphi = 1$), we consider that its viscosity is equal to the viscosity reported for liquid MgSiO$_3$ (Karki & Stixrude 2010):

$$ \eta = \eta_l = \exp \left( -7.75 + 0.005 P(\text{GPa}) - 0.00015 P(\text{GPa})^2 + \frac{5000 + 135 P(\text{GPa}) + 0.23 P(\text{GPa})^2}{T - 1000} \right). \quad (19a) $$

### 2.3.2 Solid fraction viscosity ($\varphi = 0$)

The viscosity of the solid fraction within the early Earth's deep mantle is a key parameter, which governs its cooling efficiency during the mushy stage. However, such a quantity is poorly constrained for a chondritic mantle. Instead, the viscosity of a deep 'bridgmanite-bearing' mantle has become increasingly documented (e.g. Boioli et al. 2017; Reali et al. 2019). Due to the absence of seismic anisotropy in the current lower mantle, diffusion creep was generally considered to be the dominant deformation mechanism at these depths (e.g. Karato & Li 1992), however, recent results advocate for diffusion-driven pure dislocation climb creep as a main

deformation mechanism for bridgmanite (e.g. Boioli et al. 2017; Reali et al. 2019). During the early Earth's history, the mantle temperatures were hotter, and ionic diffusion (thus diffusion creep) is expected to have been even more important in the solid mantle fraction (e.g. Frost & Ashby 1982) for both upper and lower mantle. Therefore, we assume here that the deformation of the interconnected solid phase occurs via diffusion creep only. We neglect the possible effect of polymineralic aggregates as one phase is expected to be volumetrically abundant (olivine or bridgmanite in the upper or lower mantle, respectively, see Ji et al. 2001; Huet et al. 2014). Hence for $\varphi = 0$:

$$ \eta = \eta_s = \frac{1}{A_{\text{diff}}} \exp \left( \frac{E_{\text{diff}} + P V_{\text{diff}}}{RT} \right), \quad (19b) $$

with P the pressure, R the gas constant (= 8.314 J K$^{-1}$ mol$^{-1}$) and T the absolute temperature. $A_{\text{diff}}$ is the viscosity pre-factor, which includes grain-size sensitivity. Here the grain size is kept constant, as well as the grain size exponent (i.e. equals to 3). $E_{\text{diff}}$ is the activation energy, and $V_{\text{diff}}$ is the activation volume for diffusion creep.

In addition, we considered $E_{\text{diff}}$, $V_{\text{diff}}$ and $A_{\text{diff}}$ values based on two requirements: (1) the values of the rheological parameters must be compatible with those derived from experiments (e.g. Hirth & Kohlstedt 2003 for dry olivine/upper mantle, and Xu et al. 2011 for the bridgmanite/lower mantle), (2) the calculated viscosity profile corresponding to a realistic present-day Earth mantle geotherm [for an adiabatic temperature profile with $T_p = 1600$ K (Tackley 2012) and references therein] and a PREM pressure profile in eq. (19b), must be compatible with viscosity profiles constrained by geoid and postglacial rebound [see Čižková et al. (2012) and references therein, Fig. 2 and Table 2].

In our calculations, we also investigate the potential role of the upper mantle, which presents a different mineralogy and, therefore, distinct rheological properties. For the sake of simplicity, we did not implement a transition zone (composed of wadsleyite (410–520 km) and ringwoodite (520–660 km). The mineralogical transition from olivine to bridgmanite significantly increases the viscosity of the mantle, which in turn could affect the ability of the mantle to lose its primordial heat. The rheological parameters for the upper mantle (whenever considered) are listed in Table 2.

### 2.3.2 Viscosity of the partially molten mantle ($0 \le \varphi \le 1$)

During the solidification of the mushy mantle, the fraction of solid material increases until reaching a threshold ($\varphi = \varphi_{\text{crit}}$), which separates the turbulent regime from viscous regime (Solomatov 2015). For $\varphi_{\text{crit}} < \varphi < 1$, the viscosity of the highly molten material scales with the viscosity of the molten mantle $\eta_l$ (Roscoe 1952):

$$ \eta = \frac{\eta_l}{\left(1 - \left(\frac{1-\varphi}{1-\varphi_{\text{crit}}}\right)\right)^{2.5}} \quad (19c) $$

As soon as the melt fraction threshold is reached at any mantle depth, the cooling efficiency of the primitive mantle significantly reduces, even if the mantle remains partially molten at other depths (Monteux et al. 2016). In a mushy mantle context where most of the material is solid, the viscosity is still strongly influenced by the fraction of the molten material $\varphi$. For $0 < \varphi < \varphi_{\text{crit}}$ the partially molten viscosity scales with the solid mantle viscosity $\eta_s$:

$$ \eta = \eta_s \exp(-\alpha_n \varphi), \quad (19d) $$
<figure><img src="image_288.png" /><figcaption>**Figure 2.** The different mantle viscosity profiles considered in our models for the solid phase, computed using eq. (19b). Upper panels, $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ (with reference case highlighted in yellow) and lower panels $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. Grey domain viscosity models from geoid inversion and post-glacial rebound [Čížková et al. (2012) and references therein]. Left-hand frame: no viscous dichotomy between the upper and lower mantle is considered. Right-hand frame: when a different upper mantle is considered. We performed calculations for $A_{\text{diff}}$ and $V_{\text{diff}}$ ranging from $2 \times 10^{-17}$ to $2 \times 10^{-13} \text{ Pa}^{-1}\text{s}^{-1}$ and $10^{-6}$ to $2 \times 10^{-6} \text{ Pa}^{-1}\text{s}^{-1}$, respectively (Table 2). In the B-frame, we have also studied the influence of the upper mantle viscosity with $A_{\text{diff}}$ ranging from $6 \times 10^{-10}$ to $6 \times 10^{-9} \text{ Pa}^{-1}\text{s}^{-1}$. The viscosity profiles (red, green and blue lines) are calculated considering an adiabatic temperature profile with $T_p = 1600 \text{ K}$ from the surface of the Earth to the CMB, thus neglecting the presence of top and bottom boundary layers for the figure readability.</figcaption></figure>

<table><caption>**Table 2.** Values used in eq. (19b) to calculate the solid mantle viscosity.</caption><thead><tr><th colspan="3">Lower 'bridgmanite-like' mantle rheology:</th></tr></thead><tbody><tr><td>A<sub>diff</sub></td><td>Pre-exponential parameter for diffusion creep</td><td>2 &times; 10<sup>-17</sup> to 2 &times; 10<sup>-13</sup></td><td>Pa<sup>-1</sup>s<sup>-1</sup></td></tr><tr><td>E<sub>diff</sub></td><td>Activation energy for diffusion creep</td><td>200 (Čížková et al. 2012)<br>300 (Xu et al. 2011)</td><td>kJ mol<sup>-1</sup></td></tr><tr><td>V<sub>diff</sub></td><td>Activation volume for diffusion creep</td><td>1 &times; 10<sup>-6</sup> to 2 &times; 10<sup>-6</sup></td><td>m<sup>3</sup> mol<sup>-1</sup></td></tr><tr><td colspan="4">Upper mantle rheology (dry olivine, grain size = 15 &micro;m):</td></tr><tr><td>A<sub>diff</sub></td><td>Pre-exponential parameter for diffusion creep</td><td>6 &times; 10<sup>-9</sup> to 6 &times; 10<sup>-10</sup></td><td>Pa<sup>-1</sup>s<sup>-1</sup></td></tr><tr><td>E<sub>diff</sub></td><td>Activation energy for diffusion creep</td><td>375 (Hirth &amp; Kohlstedt 2003)</td><td>kJ mol<sup>-1</sup></td></tr><tr><td>V<sub>diff</sub></td><td>Activation volume for diffusion creep</td><td>2.5 &times; 10<sup>-6</sup> (Hirth &amp; Kohlstedt 2003)</td><td>m<sup>3</sup> mol<sup>-1</sup></td></tr></tbody></table>

with $\alpha_n$ the coefficient in melt fraction-dependent viscosity. The latter equals 26 for deformation via olivine diffusion creep mechanism under anhydrous conditions (Mei et al. 2002).

## 2.4 Numerical model

We model the thermal evolution of a 2900-km-thick isochemical silicate mantle overlying an iron core by solving for the conservation of energy (eq. 1) in a 1-D, spherically symmetric domain (with a radius

ranging from 3500 to 6400 km). To this end, we used a modified version of the numerical model developed in Monteux et al. (2016). Eq. (1) is discretized using a semi-implicit predictor-corrector Finite Difference scheme, of second-order in both space and time (Press et al. 1993). Our numerical scheme was successfully benchmarked against steady and unsteady analytical solutions for diffusion problems (Crank 1975). We have also successfully benchmarked our physical model with 3-D spherical calculations at both steady and transient states from the models developed by Wagner et al. (2019)
<table><caption>**Table 3.** Parameter values used to calculate the viscosity of the lower mantle. Viscosity layering means that we consider an upper mantle in our model.</caption><thead><tr><th>Series #</th><th>Viscosity layering</th><th>E<sub>diff</sub> (kJ mol<sup>-1</sup>)</th><th>V<sub>diff</sub> (m<sup>3</sup> mol<sup>-1</sup>)</th><th>A<sub>diff</sub> (Pa<sup>-1</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>1</td><td>No</td><td>200</td><td>1 &times; 10<sup>-6</sup> &minus; 2 &times; 10<sup>-6</sup></td><td>2 &times; 10<sup>-17</sup> &minus; 1 &times; 10<sup>-15</sup></td></tr><tr><td>2</td><td>Yes</td><td>200</td><td>1 &times; 10<sup>-6</sup> &minus; 2 &times; 10<sup>-6</sup></td><td>2 &times; 10<sup>-17</sup> &minus; 1 &times; 10<sup>-15</sup></td></tr><tr><td>3</td><td>No</td><td>300</td><td>1 &times; 10<sup>-6</sup> &minus; 2 &times; 10<sup>-6</sup></td><td>5 &times; 10<sup>-15</sup> &minus; 2 &times; 10<sup>-13</sup></td></tr><tr><td>4</td><td>Yes</td><td>300</td><td>1 &times; 10<sup>-6</sup> &minus; 2 &times; 10<sup>-6</sup></td><td>5 &times; 10<sup>-15</sup> &minus; 2 &times; 10<sup>-13</sup></td></tr></tbody></table>

considering a relatively lower Ra number and a smaller viscosity contrast (see the Appendix). The mantle is discretized using 2900 equally spaced gridpoints resulting in a constant spatial resolution $\delta r = 1$ km. The variable time step is set as $\delta t = \min(\delta r^2/\kappa(r))$, where $\kappa(r) = k/(\rho C_p)$ is the effective diffusivity. The boundary conditions in our models are those described in Section 2.1.2: isothermal at the surface with $T_{surf} = 500$ K and variable heat flux accounting for heat transfer between the core and the mantle at the CMB. In all our models, the core temperature below and just above the CMB are initialized to the same value ($T_{0,core} = 4370$ K).

# 3. RESULTS: COOLING AND SOLIDIFICATION DYNAMICS

## 3.1 A reference case

We followed the thermal evolution of a deep mushy ocean with an initial temperature profile corresponding to a melt fraction of 40%

throughout the whole mantle. As a reference case, we considered the chondritic-type mantle from Series 1 with $A_{diff} = 10^{-15}$ Pa$^{-1}$s$^{-1}$ and $V_{diff} = 10^{-6}$ m$^3$ mol$^{-1}$ through the whole mantle (see Table 3). This reference case represents a lower bound in terms of viscosity (Fig. 2). The temperature initially decreases rapidly from the surface where heat is efficiently extracted by conductive cooling, and where a thin TBL initially forms (Fig. 3). In the deepest part of the mantle, the temperature profile bends towards an adiabatic temperature profile, which is more vertical than the solidus profile. As a consequence, the solidification front starts from the lowermost mantle. After 100 Myr, most of the lower mantle temperatures lie below the solidus, nevertheless two molten reservoirs remain (named hereafter SML and DML). Fig. 3 (bottom panels) shows that 270 Myr after the beginning of our simulation, the deeper one (DML) is located at a depth centred at 650 km, and the depth of shallower one (SML) ranges between 20 and 60 km. Full solidification of DML occurs prior to that of SML. Finally, after 900 Myr of cooling, the entire temperature profile is below the solidus, but remains super adiabatic,

<figure><img src="image_289.png" /><figcaption>**Figure 3.** Upper panels: temperature time evolution as a function of depth. Lower panels: melt fraction time evolution as a function of depth. In this reference case (Series 1, $A_{diff} = 10^{-15}$ Pa$^{-1}$s$^{-1}$ and $V_{diff} = 10^{-6}$ m$^3$ mol$^{-1}$), no dichotomy in the viscosity model is considered between the upper and lower mantle. The blue dashed line separates the upper and lower mantle. The right-hand panels represent close-up views of the left-hand panels. When $t > 90$ Myr, the partially molten layer is separated in 2 layers: SML and DML.</figcaption></figure>
<figure><img src="image_290.png" /><figcaption>**Figure 4.** Thickness of the top (circles) and bottom (squares) thermal boundary layer when the mushy mantle is fully solidified as a function of $A_{\text{diff}}$ and $E_{\text{diff}}$ (i.e. $E_{\text{diff}} = 200$ kJ mol$^{-1}$ for top figures or $E_{\text{diff}} = 300$ kJ mol$^{-1}$ for bottom figures). The A, B, C, D panels correspond to Series 1, 2, 3, 4, respectively (see Tables 2 and 3 for complete set of parameter values). The cases with $V_{\text{diff}} = 10^{-6}$ m$^3$ mol$^{-1}$ are illustrated with red (bottom TBL) and black (top TBL) symbols and cases with $V_{\text{diff}} = 2 \times 10^{-6}$ m$^3$ mol$^{-1}$ are illustrated with green symbols. Dashed lines represent power law fits to the thermal evolution data points.</figcaption></figure>

especially in the mid-top mantle. This solidification timescale is in good agreement with the timescale proposed by Solomatov (2000), where the complete crystallization of the shallow early mantle could last more than $10^8$ yr.

## 3.2 Influence of the viscosity parameters on the mushy mantle solidification

### 3.2.1 TBL thicknesses

From the initial thermal state, two TBLs rapidly form above and below the convecting portion of the mantle. Upon cooling, the Rayleigh number within the convecting mantle decreases and the two boundary layers thicken following the scaling used in eq. (10). Therefore, the thickness of the boundary layers scales with $Ra^{-1/3}$ and as a consequence scales with $\eta^{1/3}$ and $A_{\text{diff}}^{-1/3}$. In Fig. 4, we plotted both the bottom (red) and top (black) boundary layer thicknesses at the end of the mushy stage (i.e. as soon as the mantle reaches complete solidification) as a function of the viscosity exponential pre-factor ($A_{\text{diff}}$) and for two different values of the activation energy ($E_{\text{diff}}$) and activation volume ($V_{\text{diff}}$). Our results show that the

evolution of the TBL thickness is strongly dependent on the value of the activation energy ($E_{\text{diff}} = 200$ or $300$ kJ mol$^{-1}$) and on $A_{\text{diff}}$ as illustrated in Fig. 4. For Series 1 (Fig. 4a), both the top and bottom boundary layer thicknesses scale with $A_{\text{diff}}^{-0.27}$, which is close to the theoretical scaling of $A_{\text{diff}}^{-1/3}$ for an entirely solid mantle. This indicates that the viscosity of the solid mantle governs the thickness of the two TBL. At the end of the mushy stage, the top boundary layer thickness ranges between 80 and 250 km, whereas the bottom boundary layer thickness ranges between 45 and 140 km. For Series 3 (Fig. 4c), the behaviour of the bottom TBL thickness is similar to the Series 1 cases. The bottom TBL thickness decreases as $A_{\text{diff}}^{-0.26}$ with corresponding values ranging from ≈40 to ≈100 km, and the top TBL thickness decreases as $A_{\text{diff}}^{-0.27}$, with corresponding values ranging from ≈60 to ≈160 km.

When an upper mantle is considered (Figs 4b and d), the influence of the lower mantle viscosity on both the top and bottom TBL thicknesses vanishes. The bottom TBL thickness decreases as $A_{\text{diff}}^{-0.084}$ for Series 2 and as $A_{\text{diff}}^{-0.11}$ for Series 4. For both values of $E_{\text{diff}}$ we used for the lower mantle viscosity, the top TBL thickness decreases to a value of ≈60 km at the end of the mushy stage, independently of the value of $E_{\text{diff}}$.
<figure><img src="image_291.png" /><figcaption>**Figure 5.** Depth at which the last layers of partially molten material solidify as a function of $A_{\text{diff}}$. Open symbols represent the solidification depth of the upper molten layer (SML) whereas solid symbols illustrate the solidification depth of the lower molten layer (DML) (See also Figs 3c and d). The a, b, c, d panels correspond to Series 1, 2, 3, 4, respectively (see Tables 2 and 3 for parameter values). The cases with $V_{\text{diff}} = 10^{-6}$ m$^3$ mol$^{-1}$ are in black symbols and cases with $V_{\text{diff}} = 2 \times 10^{-6}$ m$^3$ mol$^{-1}$ are in green symbols. Dashed and dotted lines represent power law fits of the numerical data for SML and DML, respectively.</figcaption></figure>

Also, the top TBL thickness is weakly dependent on the value of $V_{\text{diff}}$, as shown in Fig. 4. However, this parameter strongly influences the bottom TBL thickness. Indeed, increasing $V_{\text{diff}}$ from $1 \times 10^{-6}$ to $2 \times 10^{-6}$ m$^3$ mol$^{-1}$ results into a viscosity increase by at least a factor 2 to 3 for both $E_{\text{diff}} = 200$ kJ mol$^{-1}$ and $E_{\text{diff}} = 300$ kJ mol$^{-1}$ cases, and for cases considering an upper mantle and its influence in rheology or not. This result illustrates the influence of $V_{\text{diff}}$ on the viscosity, and can be understood when comparing the red and green viscosity profiles from Fig. 2. As the value of $V_{\text{diff}}$ governs how the viscosity increases with depth from a reference value, increasing $V_{\text{diff}}$ does not change significantly the viscosity close to the Earth's surface. However, increasing $V_{\text{diff}}$ increases significantly the viscosity in the lowermost mantle, leading to a significant thickening of the TBL above the core mantle-boundary. In our models, right after the solidification of the molten layers (SML and DML), the viscosity above the bottom TBL for cases with $V_{\text{diff}} = 2 \times 10^{-6}$ m$^3$ mol$^{-1}$ is larger than the viscosity above the bottom TBL for $V_{\text{diff}} = 1 \times 10^{-6}$ m$^3$ mol$^{-1}$ by a factor 10–30. This important increase in the lower mantle viscosity explains the increase in TBL thickness illustrated in Fig. 4 for our range of $V_{\text{diff}}$ values as the TBL scales with $\eta^{1/3}$.

### 3.2.2 Depth of final melt layer

During the cooling and the solidification of the mushy mantle, the melt fraction decreases from a global value of 0.4 to 0 (Figs 3c and d). Depending on the solid viscosity parameters used for the early mushy mantle, two layers can remain molten before full solidification: a deep one (DML) and a shallower one (SML, see also Fig. 3). The depths at which the two last layers of melt solidify as a function of $A_{\text{diff}}$ for two different values of $E_{\text{diff}}$ and $V_{\text{diff}}$ is reported in Fig. 5. This figure shows that the crystallization mechanism strongly depends on the presence of a viscosity dichotomy between the upper and lower mantle. When no dichotomy is considered (Figs 5a and c), the behaviour is similar for $E_{\text{diff}} = 200$ kJ mol$^{-1}$ and $E_{\text{diff}} = 300$ kJ mol$^{-1}$. For both $E_{\text{diff}}$ values, the solidification depth of the upper molten layer (SML) decreases as the viscosity decreases (scaling with $A_{\text{diff}}^{-0.25}$ and $A_{\text{diff}}^{-0.35}$, respectively) whereas the solidification depth of the deep molten layer (DML) is constant and equals 660 km (i.e. the depth of the transition between the upper and lower mantle). In the later case, the transition is not the consequence of rheological properties but is related to the change of the solidus slope (eqs (12a) and (12c) and Fig. 3), which is steeper in the lower mantle than in the upper mantle. These
<figure><img src="image_292.png" /><figcaption>**Figure 6.** Duration for full solidification of the mushy mantle (initially with $\varphi = 0.4$). Open symbols represent the solidification time of the upper molten layer (SML) whereas solid symbols illustrate the solidification time of the lower molten layer (DML). The a, b, c, d panels correspond to Series 1, 2, 3, 4, respectively (see Tables 2 and 3 for parameter values). The cases with $V_{\text{diff}} = 10^{-6} \text{ m}^3 \text{ mol}^{-1}$ are illustrated with black symbols and cases with $V_{\text{diff}} = 2 \times 10^{-6} \text{ m}^3 \text{ mol}^{-1}$ are illustrated with green symbols. Dashed and dotted lines represent power law fits of the numerical data for SML and DML, respectively.</figcaption></figure>

changes in the melting properties with depth coupled with the slope of the temperature profile computed from our models explain this particular behaviour. In these cases, the depth of the deep molten layer is insensitive to the value of the viscosity parameters.

However, a mineralogical dichotomy between the upper and lower mantle is likely to appear rapidly during the solidification of the mushy mantle, due to the high-pressure polymorphism. When considering rheologically distinct upper and lower mantles (Figs 5b and d), the solidification depth of the deep molten layer is no longer tied to a depth of 660 km, but now depends on the viscosity of the lower mantle. The depth at which DML solidifies decreases when $A_{\text{diff}}$ increases (i.e. when the lower mantle viscosity decreases) and scales with $A_{\text{diff}}^{-0.06-0.1}$ (Figs 5b and d). This results in a deep solidification stage occurring at depth decreasing from 1250 to 800 km when $A_{\text{diff}}$ increases within the range envisioned in our study (i.e. when the deep mantle viscosity decreases). In contrast, the depth of final upper molten layer SML remains nearly constant (≈30–50 km) for the whole range of lower mantle viscosities considered in Fig. 5. This illustrates the fact that the depth at which the last upper layer of melt solidifies is not governed by the viscosity of the lower mantle but rather by the rheological properties of the upper mantle (we tested this hypothesis in Section 3.2.4).

### 3.2.3 Mushy stage timescale

The influence of $A_{\text{diff}}$ on the time required to fully solidify a partially molten mantle and for two different values of $E_{\text{diff}}$ and $V_{\text{diff}}$ is reported in Fig. 6. A quick inspection of eq. (1) indicates that this solidification timescale should be inversely proportional to the convective heat flux $F_{\text{conv}}$. In the hard-turbulent regime, this term scales as $\eta^{-3/7}$ whereas in the soft-turbulent regime, this term scales as $\eta^{-1/3}$. The eq. (19b) implies that an increase of either $E_{\text{diff}}$ or $V_{\text{diff}}$ yields an increase of the viscosity. On the contrary, an increase of $A_{\text{diff}}$ yields a decrease of the solid viscosity at given P and T conditions scaling with $A_{\text{diff}}^{-1}$. Consequently, if the viscosity of its solid fraction governs the characteristic solidification timescale, this time should scale as $A_{\text{diff}}^{-n}$ with n ranging between 1/3 and 3/7. This is confirmed by our numerical results (Figs 6a and c). The time required for the upper molten layer (SML) to fully solidify scales with $A_{\text{diff}}^{-0.39}$ for Series 1 (Fig. 6a) and with $A_{\text{diff}}^{-0.29}$ for Series 3 (Fig. 6c). For the lower molten layer (DML), the influence of mantle viscosity is even stronger and the time required for DML to fully solidify scales with $A_{\text{diff}}^{-0.55}$ for Series 1 (Fig. 6A) and with $A_{\text{diff}}^{-0.49}$ for Series 3 (Fig. 6c).

When no upper/lower mantle dichotomy is considered, the duration of the complete mushy mantle solidification ranges between
<figure><img src="image_293.png" /><figcaption>**Figure 7.** TBL thickness (a), depth at which the last layers (upper and lower) of partially molten material solidify (b) and time to fully solidify the last layers (upper and lower) of partially molten (c) as a function of the value of $A_{\text{diff}}$ in the upper mantle. In these figures, $A_{\text{diff}} = 10^{-15} \text{ Pa}^{-1}\text{s}^{-1}$, $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ and $V_{\text{diff}} = 10^{-6} \text{ m}^3 \text{ mol}^{-1}$ in the lower mantle whereas $E_{\text{diff}} = 375 \text{ kJ mol}^{-1}$ and $V_{\text{diff}} = 2.5 \times 10^{-6} \text{ m}^3 \text{ mol}^{-1}$ in the upper mantle. Open symbols represent the solidification time of the upper molten layer (SML) whereas solid symbols illustrate the solidification time of the lower molten layer (DML). Dashed and dotted lines represent power law fits of the numerical data for SML and DML, respectively.</figcaption></figure>

900 Myr and 4.4 Gyr for Series 1 (Fig. 6a). For Series 3 (Fig. 6c), the solidification duration ranges between 1.6 and 4.6 Gyr. Our results also show that the solidification of the deeper DML occurs prior to that of the shallower SML, with times ranging between 375 Myr and 3.2 Gyr for $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ or between 590 Myr and 3.7 Gyr for $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. Fig. 6 also shows that our solidification timescale is nearly insensitive to the value of $V_{\text{diff}}$.

When considering an upper mantle rheologically different from the lower mantle (Figs 6b and d), several changes occur, compared to models without viscous dichotomy. The solidification timescale for the upper molten layer (SML) is less affected by changes in the value of $A_{\text{diff}}$ than DML. Thus, the time required to fully solidify the whole mantle (SML and DML) exhibits a narrower range of values (between 460 and 770 Myr for Series 2 and between 570 and 700 Myr for Series 4). On the contrary, the solidification time of the deeper molten layer DML strongly depends on the viscosity of the lower mantle and scales with $A_{\text{diff}}^{-0.57}$ for Series 2 and with $A_{\text{diff}}^{-0.4}$ for Series 4. Again, the solidification time of DML is faster than SML (between 37 and 365 Myr for Series 2 or between 76 and 322 Myr for Series 4). Our results suggest that the solidification of the SML (i.e. the final episode of solidification of the early mantle in our models) is weakly sensitive to the viscosity of the lower mantle but is mostly governed by the viscosity of the upper mantle. On the other hand, the solidification time of the DML is in comparison faster, and the viscosity of the lower mantle governs the time delay. This is certainly related to the relatively low upper mantle viscosity used in this calculation (see Fig. 2).

### 3.2.4 Influence of the upper mantle viscosity

We then investigated the influence of the upper mantle viscosity on the characteristic time and length scales of mushy terrestrial mantle crystallization. Thus we considered the reference case detailed in Fig. 3 with a dichotomy in the viscosity between the upper and lower mantle. In this section, we consider constant values for $A_{\text{diff}}$, $E_{\text{diff}}$ and $V_{\text{diff}}$ for the lower mantle, and we used three different values for $A_{\text{diff}}$ for the upper mantle ranging between $6 \times 10^{-10}$ and $6 \times 10^{-9} \text{ Pa}^{-1}\text{s}^{-1}$.

The results given in Fig. 7 show that the viscosity of the upper mantle influences both the shallow characteristic time and length scales. We recall here that increasing the value of $A_{\text{diff}}$ results in a viscosity decrease. Fig. 7(a) shows that both the top and the bottom TBL thicknesses now decrease when $A_{\text{diff}}$ in the upper mantle increases. The top TBL is more sensitive to variations in the $A_{\text{diff}}$

values in the upper mantle than the bottom TBL (power exponent -0.19 compared to -0.9 in top and bottom TBL, respectively). In Fig. 7(b), the results show that the depth at which the shallow molten layer SML solidifies decreases with $A_{\text{diff}}^{-0.22}$, whereas the depth at which the deep molten layer DML solidifies increases with $A_{\text{diff}}^{0.14}$. This means that decreasing the upper mantle viscosity (i.e. increasing $A_{\text{diff}}$ in the upper mantle) favors a deeper solidification of DML, while favoring the solidification of SML closer to the surface. Finally, results from Fig. 7(c) shows that both the time at which SML and DML solidify decrease when increasing $A_{\text{diff}}$ in the upper mantle. Interestingly, the viscosity of the upper mantle has a stronger influence on the solidification time for the DML than on the shallow SML (power exponent -0.4 compared to -0.22). Hence, by controlling the heat loss in the shallower part of the early Earth, the upper mantle viscosity strongly influences the characteristic solidification time and length scales.

Figs 5(b), (d) and 6(b), (d) show that, when a dichotomy in viscosity is considered, the depth and solidification time of the shallow molten layer are weakly dependent on the viscosity of the deep mantle. However, Fig. 7 illustrates that the viscosity of the upper mantle plays a key role on the time and depth of solidification of the shallow molten layer and is more important than the influence of the lower mantle on this time. Concerning the DML, a decrease of either upper or lower mantle viscosity leads to a decrease of the solidification time of this layer. Nevertheless, this time is more influenced by the viscosity of the deep mantle than by the viscosity of the upper mantle (power exponent -0.57 compared to -0.4). An interesting behaviour arises from the depth at which the deep molten layer DML solidifies. Indeed, Figs 5(b)–(d) shows that this depth decreases when $A_{\text{diff}}$ in the lower mantle increases (i.e. when the viscosity decreases) for a fixed value of $A_{\text{diff}}$ in the upper mantle. On the contrary, the depth at which DML solidifies increases when $A_{\text{diff}}$ in the upper mantle increases for a fixed value of $A_{\text{diff}}$ in the lower mantle. Our results show that the rheological parameters of both the upper and lower mantle govern the deep processes of solidification in the lower mantle, whereas the shallower solidification processes are governed only by the properties of the upper mantle.

### 3.2.5 Summary

We have developed a numerical approach to constrain the characteristic depth and time of solidification of a mushy mantle. We have identified two persistent molten layers (SML) and (DML). We show that the cooling and solidification dynamics are very sensitive to the
Rayleigh number that increases with decreasing viscosity. Hence, an increase of the pre-exponential factor $A_{\text{diff}}$ (i.e. a decrease of the viscosity) systematically leads to a decrease of the TBL thicknesses (in agreement with eq. 10) and as a consequence to the depth of solidification of the last layer of molten material. As higher Ra values lead to a more efficient cooling of the early mantle, the timescale of complete solidification of the mantle also decreases with decreasing $A_{\text{diff}}$. Within a moderately convecting viscous mantle, the TBL thickness and the cooling timescale are expected to scale with $Ra^{-1/3}$ while within a turbulent reservoir they are expected to scale with $Ra^{-2/7}$. In our models were important changes in the parameters occur with temperature, pressure and melt fraction, the value of the exponent in the power law is slightly different (from -0.25 to -0.55 when no viscous dichotomy is considered) but the behaviour is similar.

We have characterized the influence of the solid mantle viscosity with or without a rheological contrast between the upper and lower mantle. Our parametrical study shows that the viscosity of the deep mantle influences the solidification of the DML. This result is not surprising since the bottom TBL thickness is related to the viscosity of the deep material. Hence, one can expect that the solidification characteristics (depth and time scales) of the DML to be strongly influenced by the values of $A_{\text{diff}}$ in the lower mantle. Our results show that the same reasoning can be applied to the solidification characteristics of the SML that is governed by the values of $A_{\text{diff}}$ in the upper mantle.

Our results show that the viscosity of the upper mantle affects the DML solidification characteristics. This feature illustrates that the upper mantle governs the global mantle dynamics by acting as a thermal blanket that reduces the efficiency of heat loss. Hence, a decrease in the upper mantle viscosity leads to an increase of the surface heat flux, to a more vigorous internal convection associated with a thinner bottom TBL (Fig. 7a), and to a more rapid solidification (Fig. 7c). Conversely, the viscosity of the lower mantle does not influence significantly the SML solidification characteristics (Figs 4–6, right-hand panels). Again, this illustrates that the viscosity of the upper mantle mostly controls the cooling and solidification dynamics. A low viscosity lower mantle enhances the heat transfer from the core toward the mantle but the mantle heat loss is limited by the viscous properties of the upper mantle.

# 4. DISCUSSION

## 4.1 Geological constraints

A mineralogical dichotomy and a subsequent transition of the viscous behaviour between the upper and lower mantle are likely to appear rapidly during the solidification of the mushy mantle. In the following discussion we only consider the results from the models that account for a viscous dichotomy between the upper and lower mantle (i.e. left-hand column, b and d graphics in Figs 4–6). Our model results show that the top melt layer (SML) crystallizes at the Hadean-Eoarchean boundary (500–800 Myr after Earth’s formation; Fig. 6), regardless of the model, and the crystallization proceeds at relatively shallow depths of 35–45 km (Fig. 5). On the contrary, the bottom melt layer (DML) crystallizes at deeper levels (800–1000 km; Fig. 5) and earlier (40–400 Myr after Earth’s formation; Fig. 6). If the upper mantle viscosity is considered separately, these time and depth estimates are only slightly decreased or increased (Fig. 7).

The presence of molten material and the resulting rheological weakening may have profound effect on the evolution of the early crust, on its ability to deform and on how orogens develop (Sawyer *et al.* 2011). The persistence of a melt layer at shallow depth during the Hadean and its final crystallization around the Hadean–Eoarchean boundary could prevent the formation of elevated orogenic formations due to fast isostatic compensation of any created reliefs and development of large-scale tectonic fault and shear-zones. The absence of reliefs would, in turn, result in a water-world with most of the Earth being covered by shallow water. Major faults and shear-zones could represent pathways for liquid water to penetrate to lower crustal levels and, in turn, induce intense hydrothermal activity. In addition, this weak layer at the depth of the lower crust could possibly isolate the crust from the underlying mantle. Doglioni *et al.* (2011) proposed that a stable partial melt layer between the asthenosphere and the lithosphere could induce an effective viscous decoupling between the two layers and explain the lifetime of cratonic roots. At the Hadean–Eoarchean, the viscous coupling between the mantle and the crust could have induced the beginning of large-scale Hadean crust reworking and the formation of stable Archean crustal blocks. The persistence of SML could, hence, account for the absence of Hadean crustal fragments in geological record and at the beginning of the Archean geological record. They are solid outputs from our geodynamic models and, therefore, they should have affected the dynamics of our planet early in its history. Hence, linking the timing of major differentiation events in the geological record with SML and DML crystallization could help understanding early shallow processes.

Little is known about the Hadean period since we do not have the rock record at the Hadean–Eoarchean boundary (e.g. Goodwin 1996; Guitreau *et al.* 2012). Yet, some detrital zircon crystals, formed during the Hadean, survived until today. They offer a window into the Earth’s infancy (e.g. Froude *et al.* 1983; Cavosie *et al.* 2019). In addition, relics of global chemical fractionation that occurred during the Hadean are recorded by extinct radionuclides, such as <sup>142</sup>Nd and <sup>182</sup>W (e.g. Boyet *et al.* 2003; Touboul *et al.* 2012). The <sup>182</sup>Hf-<sup>182</sup>W system operated during the first 50 Myr after Solar System formation, and it is, hence, unlikely to have recorded processes depicted in our model. In contrast, the lifetime of <sup>146</sup>Sm-<sup>142</sup>Nd system matches very well the timescale for DML crystallization and is, hence, very appropriate to help constrain the physical parameters of our models. Interestingly, most <sup>142</sup>Nd signatures point to major differentiation event(s) of the Earth around 4.3–4.4 Ga (i.e. 150–250 Ma after the Earth’s formation, e.g. Saji *et al.* 2018, and references therein, Guitreau *et al.* 2019), also consistent with detrital zircon ages (e.g. Cavosie *et al.* 2019). On the other hand, the disappearance of SML would correspond to the start of the rock record (i.e. preservation of stable crustal blocks) between 4.0 and 3.8 Ga.

Considering that the SML crystallization is correlated with the end of major resurfacing on Earth, the comparison with Venus is quite appealing. Based on the crater population, it was suggested that the surface of Venus seems uniformly young. With absence of plate tectonics, this observation suggested that catastrophic resurfacing occurs episodically on Venus (Phillips & Hansen 1998; Harris & Bédard 2014; Smrekar *et al.* 2018). The available geodynamic models point out the importance of radioactive heating in the Venusian mantle which, correlated to the presence of a rigid stagnant lid, could have resulted in an increase of the mantle potential temperature with geological time, especially in the first 1–2 Ga (O’Rourke & Korenaga 2012; Tosi *et al.* 2017). The mantle potential temperature could still be today above 1800 K on Venus, thus at a similar
<figure><img src="image_294.png" /><figcaption>**Figure 8.** For the lower mantle, range of $A_{\text{diff}}$ values considered in our study for (left-hand panel) $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ and (right-hand panel) $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. The grey domain marks dislocation creep regime, below is diffusion creep regime. The white domains represent $A_{\text{diff}}$ values not considered in our study. The red dashed domains on the colour bar represent the values derived from the geological constraints.</figcaption></figure>

level than it was early in the Earth's history (e.g. Herzberg et al. 2010). A logical conclusion is that partial melting still takes place today at shallow depths in the Venusian mantle.

## 4.2 Refinement of mantle's rheological parameters

Following the idea that SML and DML final crystallization correspond to identified Hadean geological events on Earth, the exact timing of these events can help refine most realistic values of $A_{\text{diff}}$, $E_{\text{diff}}$ and $V_{\text{diff}}$. In our models, the upper mantle viscosity does not significantly influence the solidification time of the last global molten layers, and should not strongly affect the timing of the geological events discussed above. We cannot estimate the best pair of $E_{\text{diff}}$ and $V_{\text{diff}}$, since $V_{\text{diff}}$ has very little influence on the timing of crystallization of SML and DML. Nevertheless, we can propose a couple of solutions for fixed values of $E_{\text{diff}}$. SML crystallization is essentially insensitive to $A_{\text{diff}}$ values and we, hence, cannot use it to estimate $A_{\text{diff}}$ values. On the other hand, the crystallization of DML is sensitive to $A_{\text{diff}}$ values. In order to explain the ages of 4.3–4.4 Ga inferred from $^{142}$Nd signatures, the lower mantle $A_{\text{diff}}$ values should range between $3 \times 10^{-17}$ and $7 \times 10^{-17} \text{ Pa}^{-1}\text{s}^{-1}$ for Series 2, and $7 \times 10^{-15} \text{ Pa}^{-1}\text{s}^{-1}$ to $2 \times 10^{-14}$ for Series 4. Assuming that SML accounts for the start of the geological record (i.e. preservation of stable crustal blocks) between 4.0 and 3.8 Ga, the upper-mantle $A_{\text{diff}}$ values should range between $1 \times 10^{-9}$ and $4 \times 10^{-9} \text{ Pa}^{-1}\text{s}^{-1}$. These ranges of values obtained for $A_{\text{diff}}$ in the lower and upper mantle are pretty narrow given that reference viscosities are generally unknown to multiple orders of magnitude. Moreover, the values inferred for both the lower and the upper mantle are within the range of those proposed for the mantle (Cížková et al. 2012).

## 4.3 Chemical weakening and grain size

Among the parameters used to compute the solid-state mantle viscosity the pre-exponential factor exhibits a large range of plausible values (typically two orders of magnitude). While $A_{\text{diff}}$ is called the material constant, its variability expresses the grain size sensitivity and the potential influence of chemical weakening. For olivine, this influence has been experimentally characterized for hydrogen at crustal and upper mantle pressures (e.g. Mackwell et al. 1985; Mei & Kohlstedt 2000a, b, Demouchy et al. 2012; Girard et al. 2013; Tielke et al. 2017), for iron (Zhao et al. 2009; Hansen et al. 2012),

and for titanium, (Faul et al. 2016). It can be expressed as:

$$ \frac{1}{A_{\text{diff}}} = 2 \left( \frac{A A'_{\text{CW}}}{\mu} \right)^{-1} \left( \frac{b}{d} \right)^{-3}, \qquad (20) $$

where A is thus a material constant ($A = 8.17 \times 10^{15} \text{ s}^{-1}$), $\mu$ is the shear modulus ($\mu = 80 \text{ GPa}$), b is the magnitude of the Burgers vector ($b = 0.55 \times 10^{-9} \text{ m}$) and d is the grain size (Karato & Wu 1993). $A'_{\text{CW}}$ is a dimensionless parameter characterizing the influence of the potential chemical weakening due to the incorporation of, for example Al$^{3+}$, Fe$^{2+}/^{3+}$, Ti$^{4+}$ in mantle minerals. $A'_{\text{CW}}$ can be envisioned as a stress factor in the sense that an increase of its value leads to an increase of $A_{\text{diff}}$ and, as a consequence, to a viscosity decrease. Note that hydrogen is expected to have only very minor to negligible effect on lower mantle properties, since hydrogen can barely be embedded in bridgmanite as a point defects (Bolfan-Cavanova, Keppler & Rubie 2003) and since the hydrogen solubility in periclase remains very limited (Bolfan-Casanova et al. 2002, see Bolfan-Casanova 2005, for a review). Therefore, ‘water’ weakening in the lower mantle is discarded in this study.

In our models, we have considered different values for $A_{\text{diff}}$ ranging between $2 \times 10^{-17} \text{ Pa}^{-1}\text{s}^{-1}$ and $10^{-15} \text{ Pa}^{-1}\text{s}^{-1}$ for $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ and between $5 \times 10^{-15} \text{ Pa}^{-1}\text{s}^{-1}$ and $2 \times 10^{-13} \text{ Pa}^{-1}\text{s}^{-1}$ for $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. According to eq. (20), each value of $A_{\text{diff}}$ corresponds to a set of values for the pair d and $A'_{\text{CW}}$. In Fig. 8, we plotted $A_{\text{diff}}$ as a function of d and $A'_{\text{CW}}$. From Fig. 8, we can estimate the range of plausible values for d and $A'_{\text{CW}}$ corresponding to the values of $A_{\text{diff}}$ considered in our models. The transition from diffusion to dislocation creep is expected to occur for d larger than 100 µm in the lower mantle (Boioli et al. 2017). Considering only a domain where the diffusion creep scaling applies, Fig. 8 illustrates that d ranges between 10 and 100 µm for $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ and between 1 and 100 µm for $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. In the mean time, $A'_{\text{CW}}$ ranges between $10^{-7}$ and $10^{-4}$ for $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ and between $10^{-7}$ and $10^{-2}$ for $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. Fig. 8 shows that when increasing the grain size by a factor 10, the stress factor associated to chemical weakening has to be increased by a factor 1000 to maintain a constant pre-exponential factor $A_{\text{diff}}$.

Our arguments developed in previous sections suggest the following range for $A_{\text{diff}}$ values within the lower mantle: between $3 \times 10^{-17}$ and $7 \times 10^{-17} \text{ Pa}^{-1}\text{s}^{-1}$ for $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$, and between $7 \times 10^{-15} \text{ Pa}^{-1}\text{s}^{-1}$ and $2 \times 10^{-14} \text{ Pa}^{-1}\text{s}^{-1}$ for $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$,
based on the comparison between our models and geological features. When reporting these two ranges of values in Fig. 8, we illustrate that these values correspond to the material colored within the red dashed box. Using eq. (20), we can relate our preferred values for $A_{\text{diff}}$ and the grain size that prevailed in the deep mantle during this period. Hence we can constrain $d$ between 30 and 100 µm when $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ and between 3 and 100 µm when $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. In the meantime, from our models, the dimensionless parameter characterizing the influence of the chemical weakening $A'_{\text{CW}}$ would range between $10^{-7}$ and $10^{-5}$ when $E_{\text{diff}} = 200 \text{ kJ mol}^{-1}$ and between $10^{-7}$ and $10^{-3}$ when $E_{\text{diff}} = 300 \text{ kJ mol}^{-1}$. Therefore, constraining the $E_{\text{diff}}$ values in the deep mantle from geological features could help to reduce the range of potential stress factor related to chemical weakening and grain size within the deep mantle.

## 4.4 Model limitations

1-D models are appropriate to characterize the first order cooling and solidification dynamics of a compressible mushy mantle where viscosity contrasts between solid mantle and magmas can reach up to 20 orders of magnitude. However, we acknowledge that the timescales for SML and DML crystallization may be affected by 3-D lateral variations. The early geotherms obtained from our models can be implemented in geometrically more realistic models developed to characterize the global geodynamic regime that operated prior to the onset of modern plate tectonics (Rozel et al. 2017; Agrusta et al. 2018).

The surface temperature and the ability of the top boundary layer to facilitate the heat loss are key parameters that also control the cooling dynamics of the mushy mantle. In our study we have assumed that the surface temperature was constant (=500 K) according to Lebrun et al. (2013) and Sleep et al. (2014). However, this surface temperature may be overestimated, as the dissolution of CO₂ in the condensed water is not accounted. In addition, depending on the buoyancy of the residual melt, partial melting in the shallow melt layer is likely to rise towards the surface through magma conduits and cool the interior efficiently by volcanic heat transport (Ricard et al. 2014; Kankanamge & Moore 2019). Both these shallow processes may enhance the heat evacuation and affect the timescales and depths obtained in our study.

The chemical differentiation that is not accounted in our models could affect the results obtained from our models by different ways.

(1) The vertical chemical segregation due to compatible/incompatible elements separation could be an important feature during this early cooling event (Ballmer et al. 2017) and could delay the solidification. Such a vertical segregation depends on the possible occurrence of gravitational fractionation of the Fe-enriched melt and the solid in a turbulent magma ocean.

(2) Considering a chemical differentiation between compatible/incompatible elements would also affect the partitioning of radiogenic heat producing elements. Indeed, K, Th and U are incompatible elements which will accumulate preferentially in the molten reservoirs. Hence the chemical segregation is likely to affect the heat partitioning within the mantle and as a consequence its cooling dynamics.

(3) During the chemical segregation, the composition of the molten phase will evolve towards an enrichment in incompatible elements. This chemical evolution will affect the melting curves (Andrault et al. 2017)

## 5. CONCLUSION

We have performed 1-D numerical simulations to monitor the temperature and melt fraction evolutions of an initially 40% molten early mantle. In our models, we have implemented recent and solid experimental constraints from mineral physics. We have considered a range of solid fraction viscosity compatible with the knowledge of the current lowermost mantle viscosity. Our models illustrate the influence of the solid mantle fraction viscosity on the cooling of a deep mushy mantle, and in particular on the characteristic time and depth at which complete solidification is achieved. Considering that deformation occurs via diffusion creep, the cooling dynamics is mainly governed by the pre-exponential factor $A_{\text{diff}}$. Our models highlight two molten layers (SML and DML) whose complete crystallization are separated both in time and space: DML solidifying earlier (between 40 and 400 Ma) and deeper (between 800 and 1200 km) than SML, whose solidification occurs during the first 400–800 Ma, and at depth ranging between 30 and 50 km. The viscosity of the upper mantle plays a key role on the time and depth of solidification of the shallow SML, whereas the viscosity of the deep mantle governs the duration of DML solidification.

The solidification timescales derived from our models suggest a full crystallization of the early mantle at the Hadean-Eoarchean boundary. A shallow molten layer stable during 150–250 Ma after the Earth’s formation could favor the absence of early reliefs and isolate the early crust from the underlying mantle dynamics. We associate the crystallization of SML and DML to major events in the geological record. Then, we use the timing of these records to refine a preferred set of parameters defining the mantle viscosity. Our models favor the highest values of the range of deep mantle viscosities derived from geoid inversion and post-glacial rebound (Čížková et al. 2012, and references therein). From our preferred viscosity models, our study suggests a 3–100 µm range for grain size and a chemical weakening parameter ranging between $10^{-7}$ and $10^{-3}$. Our 1-D approach is relevant when characterizing a fully mushy early mantle where large viscosity contrasts associated to large Rayleigh numbers can persist during the whole solidification of the reservoir. The temperature profiles and the TBLs characteristics inferred from our study can easily be incorporated in thermochemical evolution models of a solid terrestrial mantle.

## ACKNOWLEDGEMENTS

The authors thank Maud Boyet, Craig Storey and Hugo Moreira for fruitful discussions. The authors also thank J. Roberts and anonymous reviewer for thoughtful and constructive comments. The authors also thank G. Tobie and B. Cecconi for their useful help in the development of the numerical model. This is ClerVolc contribution number 400. H. S. acknowledges the support from the “Deutsche Forschungsgemeinschaft” (grant number SA 2042/3-1).

## References

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

* Agrusta, R., van Hunen, J. & Goes, S., 2018. Strong plates enhance mantle mixing in early Earth, *Nature Communications*, **9**, 2708.

* Andrault, D., Bolfan-Casanova, N., Lo Nigro, G., Bouhifd, M.A., Garbarino, G. & Mezouar, M., 2011. Melting curve of the deep mantle applied to properties of early magma ocean and actual core-mantle boundary, *Earth planet. Sci. Lett.*, **304**, 251–259.

* Andrault, D., et al., 2018. Deep and persistent melt layer in the Archaean mantle, *Nat. Geosci.*, **11**, 139–143, doi:10.1038/s41561-017-0053-9.
Andrault, D. et al., 2017. Toward a coherent model for the melting behaviour of the deep Earth's mantle, *Phys. Earth planet. Inter.*, **265**, 67–81.

Ballmer, M.D., Lourenço, D.L., Hirose, K., Caracas, R. & Nomura, R., 2017. Reconciling magma-ocean crystallization models with the present-day structure of the Earth's mantle, *Geochem. Geophys. Geosyst.*, **18**, 2785–2806.

Boioli, F., Carrez, P., Cordier, P., Devincre, B., Gouriet, K., Hirel, P., Kraych, A. & Ritterbex, S., 2017. Pure climb creep mechanism drives flow in the Earth's lower mantle, *Sci. Adv.*, **3**, e1601958, doi: 10.1126/sci-adv.1601958.

Bolfan-Casanova, N., 2005. Water in the Earth's mantle, *Mineral. Mag.*, **69**, 3, 229–257.

Bolfan-Casanova, N., Keppler, H. & Rubie, D., 2003. Water partitioning at 660 km depth and evidence for very low water solubility in magnesium silicate perovskite. *Geophys. Res. Lett.*, **30**, doi: 10.1029/2003GL017182.

Bolfan-Casanova, N., Mackwell, S., Keppler, H., McCammon, C. & Rubie, D., 2002. Pressure dependence of h solubility in magnesiowustite up to 25 GPa: Implications for the storage of water in the Earth's lower mantle, *Geophys. Res. Lett.*, **29**, 1449, 1–4.

Boukare, C.E., Ricard, Y. & Fiquet, G., 2015. Thermodynamics of the MgO-FeO-SiO₂ system up to 140 GPa: application to the crystallization of Earth's magma ocean, *J. geophys. Res.*, **120**, 6085–6101.

Bower, D.J., Sanan, P. & Wolf, A.S., 2018. Numerical solution of a non-linear conservation law applicable to the interior dynamics of partially molten planets, *Phys. Earth planet. Inter.*, **274**, 49–62.

Boyet, M., Blichert-Toft, J., Rosing, M., Storey, M., Télouk, P. & Albarède, F., 2003. 142Nd evidence for early Earth differentiation, *Earth planet. Sci. Lett.*, **214**, 427–442.

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

Cavosie, A.J., Valley, J.W. & Wilde, S.A., 2019. The oldest terrestrial mineral record: Thirty years of research on Hadean zircon from Jack Hills, Western Australia, in: *Earth's Oldest Rocks*, 2nd edn, Chapter 12, pp. 255–278, eds Van Kranendonk, M.J. Bennett, V.C. & Hoffmann, J.E., Elsevier.

Choblet, G. & Sotin, C., 2000. 3D thermal convection with variable viscosity: can transient cooling be described by a quasi-static scaling law?, *Phys. Earth planet. Inter.* **119**, 321–336.

Cochain, B., Sanloup, C., Leroy, C. & Kono, Y., 2017. Viscosity of mafic magmas at high pressures, *Geophys. Res. Lett.*, **44**(2), 818–826.

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

Čížková, H., van den Berg, A.P., Spakman, W. & Matyska, C., 2012. The viscosity of Earth's lower mantle inferred from sinking speed of subducted lithosphere, *Phys. Earth Planet. Inter.*, **200–201**, 56–62.

Demouchy, S., Tommasi, A., Barou, F., Mainprice, D. & Cordier, P., 2012. Deformation of olivine in torsion under hydrous conditions, *Phys. Earth planet. Inter.*, **202–203**, 56–70.

Deschamps, F. & Sotin, C., 2000. Inversion of two-dimensional numerical convection experiments for a fluid with a strongly temperature-dependent viscosity, *Geophys. J. Int.*, **143**, 204–218.

Doglioni, C., Ismail-Zadeh, A., Panza, G. & Riguzzi, F., 2011. Lithosphere-asthenosphere viscosity contrast and decoupling, *Phys. Earth planet. Inter.*, **189**, 1–8.

Faul, U.H., Cline, C.J., David, E.C. & Jackson, I., 2016. Titanium-hydroxyl defect-controlled rheology of the Earth's upper mantle, *Earth planet. Sci. Lett.*, **452**, 227–237.

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

Foley, B.J., Bercovici, D. & Elkins-Tanton, L.T., 2014. Initiation of plate tectonics from post-magma ocean thermochemical convection, *J. geophys. Res.*, **119**, 8538–8561.

Frost, H.J. & Ashby, M.F., 1982, *Deformation-Mechanism Maps: The Plasticity and Creep of Metals and Ceramics*, Pergamon, Oxford.

Froude, D.O., Ireland, T.R., Kinny, P.D., Williams, I.S. & Compston, W., 1983. Ion microprobe identification of 4,100-4,200 Myr-old terrestrial zircons, *Nature*, **304**, 616–618.

Girard, J., Chen, J., Raterron, P. & Holyoke, C.W., 2013. Hydrolytic weakening of olivine at mantle pressure: Evidence of [100](010) slip system softening from single-crystal deformation experiments. 1–9.

Goodwin, A.M., 1996. *Principles of Precambrian Geology*, pp. 327, Academic Press.

Grott, M. & Breuer, D., 2008, The evolution of the martian elastic lithosphere and implications for crustal and mantle rheology, *Icarus*, **193**, 503–515.

Guitreau, M., Blichert-Toft, J., Martin, H., Mojzsis, S.J. & Albarède, F., 2012. Hafnium isotope evidence from Archean granitic rocks for deep mantle origin of continental crust, *Earth planet. Sci. Lett.*, **337–338**, 211–223.

Guitreau, M., Boyet, M., Paquette, J.-L., Gannoun, A., Konc, Z., Benbakkar, M., Suchorski, K. & Hénot, J.-M., 2019, Hadean protocrust reworking at the origin of the Archean Napier Complex (Antarctica), *Geochem. Perspect. Lett.*, **12**, 7–11.

Hansen, L.N., Zimmerman, M.E. & Kohlstedt, D.L., 2012. The influence of microstructure on deformation of olivine in the grain-boundary sliding regime, *J. geophys. Res.*, **117**, B09201, doi: 10.1029/2012JB009305.

Harris, L. & Bédard, J., 2014. Modern approaches in solid earth sciences, in *Evolution of Archean Crust and Early Life*, Vol. 7, Chapter 9, pp. 215–288, eds Dilek, Y. & Furnes, H., Springer.

Herzberg, C., Condie, K. & Korenaga, J., 2010. Thermal history of the Earth and its petrological expression, *Earth planet. Sci. Lett.*, **292**, 79–88.

Hirth, G. & Kohlstedt, D.L., 2003. Rheology of the upper mantle and the mantle wedge: a review from the experimentalists, in *Inside the Subduction Factory*, Geophys. Monogr. Ser., Vol. **138**, pp. 83–105, ed. Eiler, J., AGU.

Huet, B., Yamato, P. & Grasemann, B., 2014. The minimized power geometric model: an analytical mixing model for calculating polyphase rock viscosities consistent with experimental data, *J. geophys. Res.*, **119**, 3897–3924.

Javoy, M., 1999. Chemical earth models, *C. R. Acad. Sci. Paris*, **329**, 537–555.

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

Ji, S., Wang, Z. & Wirth, R., 2001. Bulk flow strength of Forsterite-Enstatite composites as a function of forsterite content, *Tectonophysics*, **341**, 69–93.

Kankanamge, D.G.J. & Moore, W.B., 2019. A parameterization for volcanic heat flux in heat pipe planets, *J. geophys. Res.*, **124**, 114–127.

Karato, S., & Li, P., 1992. Diffusion creep in perovskite: Implications for the rheology of the lower mantle, *Science*, **255**, 1238–1240, doi:10.1126/science.255.5049.1238.

Karato, S.I. & Wu, P., 1993. Rheology of the upper mantle: a synthesis, *Science*, **260**, 771–778.

Karki, B.B. & Stixrude, L.P., 2010. Viscosity of MgSiO₃ liquid at Earth's mantle conditions: implications for an early magma ocean, *Science*, **328**, 740–742.

Laneuville, M., Hernlund, J., Labrosse, S. & Guttenberg, N., 2018. Crystallization of a compositionally stratified basal magma ocean, *Phys. Earth planet. Inter.*, **276**, 86–92.

Lebrun, T., Massol, H., Chassefiere, E., Davaille, A., Marcq, E., Sarda, P., Leblanc, F. & Brandeis, G., 2013. Thermal evolution of an early magma ocean in interaction with the atmosphere, *J. geophys. Res.*, **118**, 1155–1176.

Mackwell, S.J., Kohlstedt, D.L. & Paterson, M.S., 1985. The role of water in the deformation of olivine single crystals, *J. Geophys. Res.*, **90**, 11319–11333.

McDonough, W. & Sun, S.-s., 1995. The composition of the Earth, *Chem. Geol.*, **120**(3-4), 223–253.

Mei, S., Bai, W., Hiraga, T. & Kohlstedt, D.L., 2002. Influence of melt on the creep behavior of olivine-basalt aggregates under hydrous conditions, *Earth planet. Sci. Lett.*, **201**, 491–507.

Mei, S. & Kohlstedt, D.L., 2000a. Influence of water on plastic deformation of olivine aggregates 2. Dislocation creep regime, *J. geophys. Res.*, **105**, 21471–21481.

Mei, S. & Kohlstedt, D.L., 2000b. Influence of water on the plastic deformation of olivine aggregates: 1. Diffusion creep regime, *J. geophys. Res.*, **105**, 21457–21469.
Monteux, J., Andrault, D. & Samuel, H., 2016. On the cooling of a deep terrestrial magma ocean, *Earth planet. Sci. Lett.*, **448**, 140–149.

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

Mosenfelder, J.L., Asimow, P.D., Frost, D.J., Rubie, D.C. & Ahrens, T.J., 2009. The MgSiO₃ system at high pressure: thermodynamic properties of perovskite, post-perovskite, and melt from global inversion of shock and static compression data, *J. geophys. Res.*, **114**, B01203.

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

O'Rourke, J.G. & Korenaga, J., 2012. Terrestrial planet evolution in the stagnant-lid regime: size effects and the formation of self-destabilizing crust, *Icarus*, **221**, 1043–1060.

Palme, H. & O'Neill, H., 2014. Cosmochemical estimates of mantle composition, in ed. *Treatise on Geochemistry*, 2nd edn, Vol. 3. pp. 1–39, ed. Carlson, R.W., Elsevier.

Phillips, R.J. & Hansen, V.L., 1998, Geological evolution of Venus: rises, plains, plumes and plateaus, *Science*, **279**, 1492–1497.

Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P., 1993. *Numerical Recipes in FORTRAN the Art of Scientific Computing*, 2nd edn. Cambridge Univ. Press.

Reali, R., van Orman, J., Pigott, J., Jackson, J.M., Boioli, F., Carrez, P. & Cordier, P., 2019. The role of diffusion-driven pure climb creep on the rheology of bridgmanite under lower mantle conditions, *Sci. Rep.*, **9**, 2053.

Ricard, Y., Labrosse, S. & Dubuffet, F., 2014. Lifting the cover of the cauldron: convection in hot planets, *Geochem. Geophys. Geosyst.*, **15**, 4617–4630.

Roscoe, R., 1952. The viscosity of suspensions of rigid spheres, *Br. J. Appl. Phys.*, **3**, 267–269.

Rozel, A.B., Golabek, G.J., Jain, C., Tackley, P.J. & Gerya, T., 2017. Continental crust formation on early Earth controlled by intrusive magmatism, *Nature*, **545**, 332–335.

Rubie, D.C. et al., 2011. Heterogeneous accretion, composition and core-mantle differentiation of the Earth, *Earth planet. Sci. Lett.*, **301**, 31–42.

Saji, N.S., Larsen, K., Wielandt, D., Schiller, M., Costa, M.M., Whitehouse, M.J., Rosing, M.T. & Bizzarro, M., 2018. Hadean geodynamics inferred from time-varying 142Nd/144Nd in the early Earth rock record, *Geochem. Perspect. Lett.*, **7**, 43–48, doi:10.7185/geochemlet.1818.

Sawyer, E.W., Cesare, B. & Brown, M., 2011. When the continental crust melts, *Elements*, **7**, 229–234.

Simon, F. & Glatzel, G., 1929. Fusion-pressure curve, *Z. Anorg. Allg. Chem.*, **178**, 309.

Sleep, N.H., Zahnle, K.J. & Lupu, R.E., 2014. Terrestrial aftermath of the Moon-forming impact, *Philos. Trans. R. Soc. A*, **372**, doi:10.1098/rsta.2013.0172.

Smrekar, S., Davaille, A. & Sotin, C., 2018. Venus interior structure and dynamics, *Space Sci. Rev.*, **214**(5), 88, 34 pp.

Solomatov, V., 1995. Scaling of temperature-and stress-dependent viscosity convection, *Phys. Fluids*, **7**, 266.

Solomatov, V., 2007. Magma Oceans and Primordial Mantle Differentiation, in *Treatise of Geophysics*, Vol. 9, ed. Schubert, G., Elsevier.

Solomatov, V.S., 2000. Fluid dynamics of terrestrial magma ocean, in *Origin of the Earth and Moon*, pp. 323–338, eds Canup, R.M. & Righter, K., The University of Arizona Press.

Solomatov, V.S., 2015, Magma oceans and primordial mantle differentiation, in *Treatise on Geophysics*, 2nd edn., pp. 81–104, ed. Schubert, G., Elsevier, Amsterdam.

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

Thiriet, M., Breuer, D., Michaut, C. & Plesa, A.-C., 2019. Scaling laws of convection for cooling planets in a stagnant lid regime, *Phys. Earth planet. Inter.*, **286**, 138–153.

Thomas, C.W. & Asimow, P.D., 2013. Direct shock compression experiments on premolten forsterite and progress toward a consistent high-pressure

equation of state for CaO-MgO-Al₂O₃-SiO₂-FeO liquids, *J. geophys. Res.*, **118**, 5738–5752.

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

Tielke, J.A., Zimmerman, M.E. & Kohlstedt, D.L. 2017. Hydrolytic weakening in olivine single crystals, *J. geophys. Res.*, **122**, 3465–3479.

Tosi, N. et al., 2017. The habitability of a stagnant-lid Earth, *Astron. Astro-phys.*, **605**, A71.

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

Ulvrová, M., Labrosse, S., Coltice, N., Råback, P. & Tackley, P.J., 2012. Numerical modelling of convection interacting with a melting and solidification front: application to the thermal evolution of the basal magma ocean, *Phys. Earth Planet. Inter.*, **206**, 51–66, http://dx.doi.org/10.1016/j.pepi.2012.06.008.

Wagner, F.W., Plesa, A.C. & Rozel, A.B., 2019. Calibrating mixing-length theory for thermal convection in rocky planets, *Geophys. J. Int.*, **217**, 75–89.

Xu, J., Yamazaki, D., Katsura, T., Wu, X., Remmert, P., Yurimoto, H. & Chakraborty, S., 2011. Silicon and magnesium diffusion in a single crystal of MgSiO₃ perovskite, *J. geophys. Res.*, **116**, B12205, doi: 10.1029/2011JB008444.

Zhao, Y.H, Zimmerman, M.E. & Kohlstedt, D.L. 2009. Effect of iron content on the creep behavior of olivine, 1: anhydrous conditions, *Earth planet. Sci. Lett.*, **287** (1), 229–240.

## APPENDIX: NUMERICAL BENCHMARK OF THE ENERGY BALANCE EQUATION

We validated our numerical approach by benchmarking our model with the analytical non-steady state diffusion problem (eq 6.18 in the section 6. Diffusion in a sphere from Crank 1975) and obtained a relative error below $10^{-7}$. We also benchmarked our physical approach by comparing our model with those developed for modeling convection in rocky planets from Wagner *et al.* (2019). In particular, we considered the same setting used to obtain their Fig. 8(a) where they $Ra = 10^7$ and they impose a viscosity contrast of 100 between the top and the bottom of the mantle. This comparison is illustrated in Fig. A1.

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

**Figure A1.** Comparison of the temperature profiles obtained by Wagner *et al.* 2019 (Fig. 8a from their study) (black dots, red dashed line and solid blue line) and the temperature profile obtained from our 1-D model with the same parametrization ($Ra = 10^7$ and a viscosity contrast of 100).
Fig. A1 shows that our results reasonably agree with spherical calculations displayed in the study quoted above, which confirms that our approach can reproduce the thermal evolution in non-symmetric spherical geometry. Moreover, Fig. A1 shows that our results are in agreement with the 1-D (Mixing Length Theory) approach and the evolutions computed in spherical geometry by Wagner *et al.* (2019). Differences between our approach and the MLT approach by Wagner *et al.* (2019) are due to different parametrization details. Moreover, as we impose a conductive heat flux within the two

TBLs, it leads to a change in the temperature profile between the convective mantle, and the conductive TBL that is less smooth than in the models from Wagner *et al.* (2019). In the two TBL, the error between our models and the models from Wagner *et al.* (2019) can reach ≈20% while in the central parts of the mushy mantle, the error is less than 5%. However, for larger values of Ra numbers that are more relevant to our study, the TBL are considerably thinner, therefore these differences are expected to vanish.
JGR Planets

RESEARCH ARTICLE

10.1029/2023JE007805

**Key Points:**

* We develop numerical models to monitor the crystal fraction evolution within a convecting magma ocean

* We illustrate that during the early crystallization of a magma ocean, suspension should be the dominant process compared to sedimentation

* We observe that bridgmanite segregation is more likely to occur at the bottom of a MO enriched in SiO₂ compared to pyrolite

**Correspondence to:**

J. Monteux,
julien.monteux@uca.fr

**Citation:**

Monteux, J., Qaddah, B., & Andrault, D. (2023). Conditions for segregation of a crystal-rich layer within a convective magma ocean. *Journal of Geophysical Research: Planets*, 128, e2023JE007805. https://doi.org/10.1029/2023JE007805

Received 28 FEB 2023
Accepted 11 APR 2023

# Conditions for Segregation of a Crystal-Rich Layer Within a Convective Magma Ocean

J. Monteux¹, B. Qaddah¹,², and D. Andrault¹

¹Laboratoire Magmas et Volcans, Université Clermont-Auvergne, CNRS, IRD, Clermont-Ferrand, France, ²Institut Jean Lamour, Campus Artem, Université de Lorraine, Nancy, France

**Abstract** The vigor of the thermal convection in a terrestrial magma ocean tends to prevent the sedimentation of the solid grains. Understanding of the overall dynamical behavior of this solid phase segregation is required to anticipate the solidification mechanisms in the early Earth mantle. We develop numerical models to monitor the crystal fraction evolution of a convecting magmatic reservoir. Our models show that the ability of the crystal fraction to disperse or sediment within the domain strongly depends on the crystal size, the density difference and the magma viscosity. Our models show that the critical value of the convection/buoyancy stress ratio separating sedimentation/suspension regimes can be smaller than 0.1. Hence, during the early crystallization of a magma ocean, suspension should be the dominant process. We then investigate the possibility of bridgmanite segregation by comparing the density difference between the MO and the solid crystals for different compositions and the critical density contrast above which crystal segregation is likely to occur. We define the relevant set of parameters; including the P-V-T equations of state of coexisting melt and bridgmanite in the mushy MO. We observe that bridgmanite grains are unlikely to segregate in a mantle of pyrolite composition. However, bridgmanite segregation is more likely to occur at the bottom of a MO enriched in SiO₂, compared to pyrolite. When a solidifying layer contains 60% of bridgmanite and 40% of melt, we observe a significant SiO₂ enrichment with increasing mantle depth in a primitive mantle compatible with seismic and geochemical observations.

**Plain Language Summary** Following the giant impact that formed the Earth-Moon system, the Earth's mantle probably underwent a major melting episode forming a magma ocean several hundred kilometers thick. During the solidification of this magma ocean, the first crystals either settled at depth or were efficiently mixed with the magma. This process played a determining role in the chemical composition of the early mantle. We have developed numerical models to monitor the evolution of a thin layer of crystals in a magma ocean. Our models characterize the influence of crystal size, magma viscosity and density difference between crystals and magma. Our models support efficient mixing of crystals in the magma ocean. However, in some settings such as a SiO₂-rich magma ocean, bridgmanite grains may separate from the magma and form a solid layer at the base of the Earth's mantle early during cooling. This process may generate SiO₂ enrichment with depth in agreement with seismic and geochemical observations.

## 1. Introduction

During its early history, the Earth experienced several episodes during which a significant volume fraction of its mantle was molten. These melting events were enhanced by radiogenic heating (Yoshino et al., 2003), viscous dissipation during core formation (Monteux et al., 2009; Samuel et al., 2010) and kinetic energy converted into heat during large impacts (Tonks & Melosh, 1993). The giant impact leading to the Earth/Moon system could even have resulted in complete melting of the Earth mantle (Nakajima & Stevenson, 2015) and the subsequent formation of a nearly 3,000 km-thick magma ocean.

Extremely vigorous convection is likely to occur within such a molten reservoir (Monteux et al., 2016; V. Solomatov, 2015). The variations in temperature are strong enough to generate significant natural convective flows. Moreover the potentially large thickness of a molten-silicate reservoir coupled with low viscosity material (Karki & Stixrude, 2010; Xie et al., 2021) leads to Rayleigh numbers (comparing advective and diffusive timescales) ranging from $10^{20}$ to $10^{30}$ (Patočka et al., 2020) compared to the current Rayleigh number of the Earth's mantle of $10^6-10^8$ (Ricard, 2007). The cooling and solidification dynamics are difficult to monitor within such

© 2023. American Geophysical Union.
All Rights Reserved.
an environment. However, recent numerical models agree on the characteristic cooling timescales that can be constrained to several thousands of years (Lebrun et al., 2013; Monteux et al., 2016).

The study of particle sedimentation in vigorously convecting fluids and, by extension, the study of crystal behavior in convecting magmas is a highly debated subject in both fluid dynamics (Lavorel & Le Bars, 2009) and Earth Sciences (Jaupart & Tait, 1995). A large number of studies have investigated crystal settling in magma chambers using laboratory experiments (Martin & Nokes, 1988; Sturtz et al., 2021), numerical modeling (Höink et al., 2005; Verhoeven & Schmalzl, 2009) or both (Weinstein et al., 1988). All these studies point out the importance of particle deposition and the subsequent differentiation mechanisms that might occur within large plutonic reservoirs. In magma oceans, the physical context is quite different, with potentially more vigorous convection involving much thicker reservoirs. Recent numerical models from Patočka et al. (2020) constrained the settling of inertial particles in turbulent Rayleigh-Bénard convection without taking into consideration the re-entrainment process. In particular, Patočka et al. (2020) identified different settling regimes and derived characteristic residence timescales. For their range of particle sizes and densities, their results show that the settling should be rapid, much shorter than the typical timescales for the solidification of a whole-mantle terrestrial magma ocean, thus supporting the idea of fully fractional crystallization. However, using analog models, V. S. Solomatov et al. (1993) introduced the conditions for particle re-entrainment after sedimentation at the bottom of a convecting reservoir. At sufficiently high Rayleigh numbers, their models illustrate that the particles might be re-entrained by the viscous stress produced by thermal plumes, emphasizing that entrainment of crystals cumulated in dunes was likely at the bottom of magma oceans. Such a re-entrainment process could therefore favor equilibrium crystallization. More recent laboratory experiments from Lavorel & Le Bars (2009) confirmed the findings of V. S. Solomatov et al. (1993) and characterized the influence of the density ratio between fluid and crystals and the temperature difference driving thermal convection.

Several parameters likely influence the settling dynamics: (a) the density difference between the settling crystal and the surrounding convecting magma ocean, (b) the crystal size, which is strongly dependent on the magma ocean cooling rate, and (c) the viscosity of the magma ocean whose value may span several orders of magnitude depending on the pressure/temperature conditions within the early magma ocean. Recently, Caracas et al. (2019) proposed that the accumulation of crystals might occur near the depth of neutral buoyancy between crystals and the coexisting melt. They determined this neutral buoyancy depth for different values of iron partition coefficients between melt and solid fractions. They showed that the crystals can cumulate and form a mushy layer within the magma ocean, thus separating a shallow MO from a basal MO. Here, we monitor the stability of a thin mushy layer cumulated within a magma ocean. We evaluate the effects of grain size, density difference between crystals and the surrounding magma ocean, and viscosity values for the magma ocean derived from petrological experiments. The main purpose of our models is to determine if the crystal layer is sedimented at the bottom of the reservoir or is efficiently mixed into the convective reservoir.

The hypothesis of a fully molten mantle is difficult to reconcile with geochemical data and the need for hidden reservoirs (Boyet & Carlson, 2005). However, the cooling and solidification processes occurring within large partially molten silicate reservoirs may lead to significant chemical differentiation events and thus to chemical segregation. The depth at which solidification initiates is governed by the intersection between the early geotherm and the melting curves (Caracas et al., 2019) proposed that crystal settling would occur toward neutral buoyancy depths where crystals are paradoxically more easily entrained by the convective flow. Both chemical composition and temperature contribute to the buoyancy of a crystal layer within a convecting reservoir. The chemical differentiation that happens following a large melting event is strongly related to the ability of the crystal fraction to sink or float within the convective reservoir. This ability is itself governed by the ability of the crystallized minerals to integrate heavy elements such as Fe into their structure (Andrault et al., 2012; Nomura et al., 2011). Depending on the crystal buoyancy relative to the molten magma, different scenarios may involve either a bottom up solidification (Monteux et al., 2016) or a mid-mantle solidification, in which case a basal magma ocean would be possible (Labrosse et al., 2007). In our study, we discuss the implications of our model on the mechanism of magma ocean solidification after a major Moon-forming impact on an Earth-like planet.

In this study, we have developed a numerical model to monitor the cooling, segregation and chemical evolution of the early Earth mantle. Our study is laid out as follows. In Section 2, we present our physical and numerical models. In Section 3, we present the main numerical results from our systematic parametric study emphasizing the influence of melt viscosity, crystal size and density difference between crystals and melt. Section 4 details a
<table><caption>**Table 1**

Symbol Definitions and Values of the Physical Parameters Used in This Study</caption><thead><tr><th></th><th>Symbol</th><th>Value or range</th></tr></thead><tbody><tr><td>Magma ocean density</td><td>ρ<sub>c,0</sub></td><td>3,000 kg/m<sup>3</sup></td></tr><tr><td>Crystals density</td><td>ρ<sub>d,0</sub></td><td>ρ<sub>c,0</sub> + Δρ</td></tr><tr><td>Density difference</td><td>Δρ</td><td>5–100 kg/m<sup>3</sup></td></tr><tr><td>Magma ocean viscosity</td><td>η<sub>c</sub></td><td>0.1–10 Pa.s</td></tr><tr><td>Crystal diameter</td><td>D<sub>d</sub></td><td>10<sup>-4</sup>–10<sup>-3</sup> m</td></tr><tr><td>Heat capacity</td><td>C<sub>p</sub></td><td>1,000 J.kg<sup>-1</sup>.K<sup>-1</sup></td></tr><tr><td>Thermal conductivity</td><td>k</td><td>5 W.m<sup>-1</sup>.K<sup>-1</sup></td></tr><tr><td>Thermal diffusivity</td><td>κ = k/(ρ<sub>c,0</sub>C<sub>p</sub>)</td><td>1.7 × 10<sup>-6</sup> m<sup>2</sup>.s<sup>-1</sup></td></tr><tr><td>Thermal expansion coefficient</td><td>α</td><td>10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Reference temperature</td><td>T<sub>0</sub></td><td>4,000 K</td></tr><tr><td>Temperature difference</td><td>ΔT</td><td>100 K</td></tr></tbody></table>

criterion for the stability of a crystallizing layer. Finally, Section 5 describes the evolution of mushy material in the context of progressive magma ocean crystallization. The results and conclusion are presented in Section 6.

## 2. Physical and Numerical Models

### 2.1. Governing Equations

We monitored the stability of a horizontal mushy layer composed of solid crystals within an initially convecting fully molten magma reservoir. For this purpose, we used the Euler-Euler model which is an accurate dispersed multiphase flow model describing both the crystal and liquid phases. We defined two non-miscible phases in our models: the continuous phase (subscript c) to represent the liquid magma ocean and the dispersed phase (subscript d) that represents the crystal phase. Both phases behaved as Newtonian and incompressible fluids. The local volume fraction of the fluid continuous phase $\phi_c$ and the dispersed solid phase $\phi_d$ were the local average volume of melt and crystals respectively ($\phi_c = 1 - \phi_d$).

The dynamic and thermal evolution of the mushy layer and ambient liquid are governed by the Navier-Stokes and heat transfer equations. We initially solved the equations for the continuous phase heat in 2D rectangular geometry using the Rayleigh-Bénard system. Once it reached a steady state, we solved the dynamic evolution of the mushy layer within the convective reservoir. The Euler-Euler flow model defines one set of Navier-Stokes equations for each phase:

* Mass conservation is represented by the continuity equations (Crowe et al., 1998):

$$ \frac{\partial \phi_d}{\partial t} + \nabla \cdot (\phi_d \mathbf{u}_d) = 0 \qquad (1) $$

$$ \nabla \cdot ((1 - \phi_d) \mathbf{u}_c + \phi_d \mathbf{u}_d) = 0, \qquad (2) $$

with $\mathbf{u}_d$ and $\mathbf{u}_c$ representing the dispersed phase and the continuous phase fluid velocity vector respectively (m.s<sup>-1</sup>). We assumed that the mass transfer between the two phases was zero. In our study, we did not consider either solidification or melting processes but choose to develop a parametric study (see Section 3) in which the parameters such as particle volume fraction, density contrast and particle radii vary within a range of values compatible with particle sedimentation for magma ocean conditions.

* Conservation of momentum is represented by the Navier-Stokes equations. These equations represent a balance between the inertial, pressure, viscous, gravitational and drag forces. In the case of an incompressible Newtonian fluid, this yields

$$ \rho_{c,0} \left( \frac{\partial \mathbf{u}_c}{\partial t} + \mathbf{u}_c \cdot \nabla \mathbf{u}_c \right) = \nabla \cdot \left[ -P \mathbf{I} + \eta_c \left( \nabla \mathbf{u}_c + (\nabla \mathbf{u}_c)^T \right) \right] + \rho_c g + \frac{F_{m,c}}{\phi_c} \qquad (3) $$

$$ \rho_{d,0} \left( \frac{\partial \mathbf{u}_d}{\partial t} + \mathbf{u}_d \cdot \nabla \mathbf{u}_d \right) = \nabla \cdot \left[ -P \mathbf{I} + \eta_d \left( \nabla \mathbf{u}_d + (\nabla \mathbf{u}_d)^T \right) \right] + \rho_d g + \frac{F_{m,d}}{\phi_d} \qquad (4) $$

with $\rho_{c,0}$ the reference melt density, $\rho_{d,0}$ the reference solid phase density (see Table 1 for values), $\rho_c$ the melt density (kg.m<sup>-3</sup>), $\rho_d$ the solid phase density, $\eta_c$ the melt viscosity (Pa.s), $\eta_d$ the solid phase viscosity (Pa.s), $t$ the time (s), $P$ the fluid pressure (Pa), $g$ the gravitational acceleration (m.s<sup>-2</sup>) and $\mathbf{I}$ the identity matrix. $F_{m,c}$ and $F_{m,d}$ are the drag forces (N/m<sup>3</sup>).

* Heat conservation: Assuming an incompressible fluid with no viscous dissipation and no internal heat sources, the internal energy balance can be written as

$$ \rho_{c,0} C_p \left[ \frac{\partial T}{\partial t} + \mathbf{u} \cdot \nabla T \right] + \nabla \cdot \left[ -k \nabla T \right] = 0 \qquad (5) $$

with $C_p$ the heat capacity at constant pressure (J.kg<sup>-1</sup>.K<sup>-1</sup>), $T$ the fluid temperature (K), $\mathbf{u}$ the fluid velocity vector (m.s<sup>-1</sup>), $k$ the thermal conductivity (W.m<sup>-1</sup>.K<sup>-1</sup>). No heat source was considered in our models (e.g.,
radiogenic, tidal or latent heat). For simplicity, we assigned the same heat capacities and thermal conductivities for both the continuous and dispersed phases.

For each fluid parcel, the viscosity ranges between the viscosity of the fully molten phase $\eta_c$, which is imposed for each simulation (see Table 1 for values), and the mixture viscosity $\eta$, which depends on the dispersed solid phase volume fraction $\phi_d$. In our magma ocean context, it is the liquid phase viscosity that governs the flow (Monteux et al., 2016), and the solid phase viscosity is not used in our calculations. For the mixture, we applied a Krieger-type dynamic viscosity which is appropriate for monitoring dispersed phases:

$$ \eta = \eta_c \left( 1 - MIN \left( \frac{\phi_d}{\phi_{d,max}}, 0.999 \right) \right)^{-2.5\phi_{d,max}} \quad (6) $$

where $\phi_{d,max}$ is the maximum packing concentration. $\phi_{d,max}$ depends on the geometry of the crystals within the system and on the shear rate. For mono-dispersed systems with spheres and low shear rates, $\phi_{d,max}$ ranges between 0.5 and 0.7 (Dörr et al., 2013). In our models $\phi_{d,max}$ was set to 0.6 which is the critical packing crystal fraction at which it becomes rheologically impossible for a magma to erupt (Marsh, 1981). In our calculations, there is no requirement for $\phi_d$ to remain between 0 and $\phi_{d,max}$. The maximum value $\phi_{d,max}$ was only used in the expression of the viscosity (Equation 6), which controls the behavior of the flow. When $\phi_d \ge \phi_{d,max}$, the solid viscosity should govern the flow (e.g., Monteux et al., 2020) instead of the liquid viscosity. Hence the viscosity calculation in Equation 6 should account for this discrepancy. Since the cases of our study fall mainly around the suspension/ sedimentation criteria, the crystal fraction never overcome maximum packing concentration.

The drag force acts in opposition to the relative motion of the crystals in the surrounding magma. It is governed by the particle shape and the relative Reynolds number (comparing inertial and viscous forces). In reality the crystal shapes are various and complex. To compute the drag force $F_m$, we considered the Schiller-Naumann drag model, a widely used model in the multiphase flow literature that is particularly relevant for dispersed rigid spheres and for the range of Reynolds number values in our study (lower than $10^4$) (Ibrahim & Meguid, 2020; Schiller & Naumann, 1935). $F_m$ is defined as:

$$ F_{m,c} = -F_{m,d} = \beta(\mathbf{u}_d - \mathbf{u}_c) \quad (7) $$

with $\beta$ the momentum transfer coefficient (in kg m$^{-3}$ s$^{-1}$) between fluid and solid particles (Ibrahim & Meguid, 2020)

$$ \beta = \frac{3}{4} \frac{\rho_{c,0} \phi_d C_D (\mathbf{u}_d - \mathbf{u}_c)}{D_d} \quad (8) $$

$D_d$ is the particle diameter (m) and $C_D$ is the drag coefficient for a single particle. For $Re < 1,000$ (Schiller & Naumann, 1935):

$$ C_D = \frac{24(1 + 0.15Re^{0.687})}{Re} \quad (9) $$

For $Re > 1,000$, $C_D = 0.44$ (Wen & Yu, 1966).

The density of the magma ocean varies linearly with temperature in the buoyancy term of the momentum equation following the Boussinesq approximation. The density calculation accounts for the contrast relating to the thermal effects and the phase influence. It was calculated as follows:

$$ \rho_c = \rho_{c,0}(1 - \alpha(T - T_0)) \quad (10) $$

where $T_0$ is the reference temperature (in K) and $\alpha$ is the thermal expansion coefficient (K$^{-1}$). The initial reference temperature for the reservoir is highly influenced by the full thermal history of the early Earth as well as its depth within the Earth mantle. Here we fixed $T_0 = 4000$ K. The particle density was calculated as:

$$ \rho_d = \rho_c + \Delta\rho \quad (11) $$

with $\Delta\rho$ the intrinsic density difference between the molten silicates and the solid crystals.

For the dynamics, we used free-slip boundary conditions at all the system boundaries for both the dispersed and continuous phases. Considering free-slip boundary conditions at the bottom of the reservoir favors re-suspension compared to cases with no-slip boundary conditions. We fixed the temperature to 4100 K at the bottom boundary
<figure><img src="image_296.png" /><figcaption>**Figure 1.** Nusselt number as a function of the Rayleigh number during the first phase where only thermal convection is considered (no settling of particles). Results from our models are shown by black circles. For comparison, we plot the scaling laws proposed by Wolstencroft et al. (2009) and Turcotte and Schubert (1982) with blue and green dashed lines respectively. The red dashed line represents the result from a power-law fit of our numerical results. Corresponding scaling laws are detailed in the legend of the figure.</figcaption></figure>

and to 4000 K at the top boundary (i.e., $\Delta T = 100$ K), with no flux conditions at lateral boundaries.

## 2.2. Non Dimensional Parameters

Before considering the stability of a crystal layer within the molten reservoir, we first solved only the equations relating to thermal convection within the reservoir. In our geophysical investigation, the least constrained parameter was the magma ocean viscosity, which ranges widely in value. Here we used viscosities for the continuous phase ($\eta_{c,0}$) ranging between 0.1 and 10 Pa.s (Karki & Stixrude, 2010; Xie et al., 2021). The convection flow induced within the magma ocean was characterized by the thermal Rayleigh number (Ra) describing the vigor of convection as a consequence of basal heating:

$$Ra = \frac{\alpha \rho_{c,0} g \Delta T H^3}{\eta_c \kappa} \qquad (12)$$

with $\kappa$ the heat diffusivity of the convective fluid (see Table 1 for value) and $H$ the thickness of the reservoir in our models. We used a reservoir thickness $H = 1$ m which led to Ra values ranging between $1.8 \times 10^6$ and $1.8 \times 10^8$. In reality, the magma ocean is likely to involve the full thickness of the Earth's early mantle (i.e., a thickness of 2,900 km). In these conditions the Rayleigh number could reach values as high as $10^{30}$ and the associated dynamics would be extremely turbulent, playing a role on all scales down that of the crystals (V. S. Solomatov, 2000). However, crystal settling is likely to occur

close to mechanical boundaries (CMB or solidified wall) where the velocity is lower than in the middle of the magma ocean. Hence our choice of parameters is relevant for monitoring the crystal settling at the bottom of a cooling magma ocean.

From our purely thermally convective models, we monitored the surface heat flux $F$ which is proportional to the average value of the temperature gradient at the surface of the reservoir. We obtained $F$ as a function of time for the range of Ra values considered in our study. We stopped our thermally convective simulations once the surface heat flux had converged to a constant value and used this value to calculate the corresponding Nusselt number Nu:

$$Nu = \frac{HF}{k\Delta T} \qquad (13)$$

Figure 1 shows the Nu-Ra relation obtained from our models and compares our numerical results with the scalings obtained from Turcotte and Schubert (1982) and Wolstencroft et al. (2009). Our results show an agreement with these two scalings and we obtain $Nu = 0.14Ra^{0.34}$. For each Ra value used in our models, we obtained a velocity field for the continuous phase and a temperature field that were used as the initial temperature and velocity fields for the Euler-Euler model when we investigated the stability of a thin mushy layer with different characteristics ($D_d$ and $\Delta \rho$). The initial velocity field for the dispersed phase used in the Euler-Euler model is null.

## 2.3. Numerical Model

We used the commercial software COMSOL Multiphysics (version 5.4) that has been previously validated for two phase flow applications (Qaddah et al. (2019, 2020)). Equations 1–5 were solved using the multiphase-flow tool-suite of COMSOL Multiphysics 5.4 (Computational Fluid Dynamics and Heat transfer modules) (For more informations see the Code Availability section). We used a 2D domain with a thickness $H$ and a width of $4H$ (see Figure 2). Applying the model to a 2D geometry involves making assumptions, which may lead to discrepancies when compared with 3D models (e.g., V. S. Solomatov et al., 1993). The difference between 2D and 3D Rayleigh-Bénard simulations varies with the Prandtl number (Pr) which is the ratio between the momentum and the thermal diffusivities. For the range of Ra values used in our study, the Nu – Ra and the Nu – Pr scalings are similar for Pr values larger than 1 but differ when Pr values are close to 1 (van der Poel et al., 2013). In our simulations, the Prandtl numbers were much greater than 1 which makes our 2D models a reliable means of characterizing the
<figure><img src="image_297.png" /><figcaption>**Figure 2.** Schematic representation of the computational domain of our study. The blue material represents the continuous phase (molten silicates) and the green material represents the initially mixed continuous and dispersed phases. A zoomed view is illustrated in the top-right panel.</figcaption></figure>

$h = H/10$, $\phi_{d,0} = 0.5$, $\eta_c = 10$ Pa.s (i.e., $Ra = 1.8 \times 10^6$) and $\Delta\rho = 30$ kg.m$^{-3}$. Figure 3 shows the average value $\phi_{d,average}$ along a horizontal profile located at the bottom of the reservoir at $t = 600$ s when the mushy layer is fully sedimented. $\phi_{d,average}$ increases with decrease in the mesh size from 0.05 to 0.01 m. However, below a mesh size value of 0.01 m, $\phi_{d,average}$ converges to a constant value of 0.45. In order to save computational time, we used a mesh size of 0.008 which gave an error of 7.7% relative to the converged value obtained for the smallest mesh size. In addition, for the most vigorously convective contexts, we decreased the size of the mesh by a factor of two either in the lowermost part of the reservoir where crystal-rich patches are likely to settle or within the whole reservoir. Therefore, the calculation domain was composed of 100,000–180,000 mesh elements. We list the details of our models in Table 2. Each run presented in this study represents (for the Euler-Euler flow model only) 12–72 hr computation time on a bi-processor, 8-cores, 3.2-GHz workstation.

particle sinking dynamics in a thermally convecting reservoir. The computational cost of 2D models being significantly lower than for 3D models, we considered a 2D Cartesian geometry in our models as a first step toward more realistic models (Patočka et al., 2020). 3D models can be developed with COMSOL Multiphysics and will be the subject of a future study.

As a first attempt to characterize the stability of a thin newly formed mushy layer within a convective magma ocean, we arbitrarily considered that the thickness of this mushy layer was $h = H/10$. The mush thickness was set as a constant in all our models. Within the mushy layer, the crystals were considered as spheres with a fixed diameter $D_d$. In our models, this diameter ranges between $10^{-4}$ and $10^{-3}$ m. The initial volume fraction of crystal $\phi_{d,0}$ was also set to a constant value ($\phi_{d,0} = 0.1$).

The COMSOL Multiphysics software applies a finite-element method to a grid mesh to solve the set of equations (Equations 1–5). In order to study the sensitivity analysis of mesh size on the main results of this paper, we ran several models varying the size of the mesh (between 1.6 and 50 mm). In these resolution tests we looked at the overturn of a mushy layer with

## 2.4. Reference Case

In this subsection, we detail the case corresponding to simulation no. 47 in Table 2. For this reference case, we used $Ra = 4.4 \times 10^6$ (corresponding to $\eta_c = 4$ Pa.s), $\Delta\rho = 50$ kg.m$^{-3}$ ($\Delta\rho/\rho_{c,0} = 1.7\%$) and $D_d = 10^{-3}$ m. Figure 4 illustrates the thermal and dynamic time evolutions of the convective magma ocean and the mushy crystal layer. Because of negative buoyancy, the crystal layer sinks rapidly. Figure 4 shows that a rapid overturn occurs concentrating the crystals in the lower part of the reservoir during the first hundreds of seconds after the formation of the layer. First, the crystal fraction is diluted in the lower reservoir. Then it starts to separate out into patches in which the concentration is much higher than the initial concentration. This early stage limits the heat transfer from the bottom surface to the reservoir, leading to a transient cooling of the magma.

To monitor the evolution of the average concentration $C$ (particles/m$^3$) of crystals in the upper mid reservoir for the reference case, we defined:

$$e = \frac{\int_{H/2}^{H} \int_{0}^{4H} C dxdz}{\int_{0}^{H} \int_{0}^{4H} C dxdz} \quad (14)$$

Höink et al. (2005) refers to this parameter as the relative entrainment. Figure 5 illustrates the time evolution of $e$ for our reference case and shows that for a period of 2,000 s, the upper reservoir is depleted in crystals. After nearly 10,000 s, large hot plumes form at the bottom of the reservoir (see Figure 4) and transport the crystals up to the upper reservoir. Hence the crystal concentration in the upper reservoir increases. This second phase lead to

<figure><img src="image_298.png" /><figcaption>**Figure 3.** Dispersed fraction averaged along an horizontal profile at the bottom of the computational domain as a function of the mesh size. In these models, $h = H/10$, $\phi_{d,0} = 0.5$, $\eta_c = 10$ Pa.s, and $\Delta\rho = 30$ kg.m$^{-3}$. The red dashed line represents the asymptotic value at which $\phi_{d,average}$ converges for high-resolution (i.e., small mesh size) models.</figcaption></figure>
<table><caption>**Table 2**

Parameters for All Simulations Used in This Study</caption><thead><tr><th></th><th>η<sub>c</sub> (Pa.s)</th><th>Ra</th><th>D<sub>a</sub> (m)</th><th>Δρ (kg/m<sup>3</sup>)</th><th>Sh</th></tr></thead><tbody><tr><td>#1</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>5</td><td>4.16</td></tr><tr><td>#2</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>10</td><td>2.08</td></tr><tr><td>#3</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>20</td><td>1.04</td></tr><tr><td>#4</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>30</td><td>0.69</td></tr><tr><td>#5</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>40</td><td>0.52</td></tr><tr><td>#6</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>50</td><td>0.42</td></tr><tr><td>#7</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>5</td><td>0.41</td></tr><tr><td>#8</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>10</td><td>0.2</td></tr><tr><td>#9</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>20</td><td>0.1</td></tr><tr><td>#10</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>30</td><td>0.07</td></tr><tr><td>#11</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>40</td><td>0.052</td></tr><tr><td>#12</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>50</td><td>0.042</td></tr><tr><td>#13</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>80</td><td>0.026</td></tr><tr><td>#14</td><td>0.1</td><td>1.8 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>100</td><td>0.02</td></tr><tr><td>#15</td><td>0.2</td><td>9. × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>50</td><td>0.052</td></tr><tr><td>#16</td><td>0.2</td><td>9. × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>60</td><td>0.044</td></tr><tr><td>#17</td><td>0.4</td><td>4.4 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>5</td><td>0.65</td></tr><tr><td>#18</td><td>0.4</td><td>4.4 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>10</td><td>0.33</td></tr><tr><td>#19</td><td>0.4</td><td>4.4 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>20</td><td>0.16</td></tr><tr><td>#20</td><td>0.4</td><td>4.4 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>30</td><td>0.11</td></tr><tr><td>#21</td><td>0.4</td><td>4.4 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>40</td><td>0.082</td></tr><tr><td>#22</td><td>0.4</td><td>4.4 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>50</td><td>0.065</td></tr><tr><td>#23</td><td>0.7</td><td>3. × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>40</td><td>0.1</td></tr><tr><td>#24</td><td>0.7</td><td>3. × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>50</td><td>0.081</td></tr><tr><td>#25</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>5</td><td>8.9</td></tr><tr><td>#26</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>10</td><td>4.45</td></tr><tr><td>#27</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>20</td><td>2.22</td></tr><tr><td>#28</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>30</td><td>1.48</td></tr><tr><td>#29</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>40</td><td>1.1</td></tr><tr><td>#30</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>50</td><td>0.089</td></tr><tr><td>#31</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>10</td><td>0.44</td></tr><tr><td>#32</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>20</td><td>0.22</td></tr><tr><td>#33</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>30</td><td>0.148</td></tr><tr><td>#34</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>35</td><td>0.127</td></tr><tr><td>#35</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>40</td><td>0.111</td></tr><tr><td>#36</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>50</td><td>0.09</td></tr><tr><td>#37</td><td>1.</td><td>1.8 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>100</td><td>0.044</td></tr><tr><td>#38</td><td>4.</td><td>4.4 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>50</td><td>1.4</td></tr><tr><td>#39</td><td>4.</td><td>4.4 × 10<sup>6</sup></td><td>3 × 10<sup>-4</sup></td><td>50</td><td>0.47</td></tr><tr><td>#40</td><td>4.</td><td>4.4 × 10<sup>6</sup></td><td>6 × 10<sup>-4</sup></td><td>50</td><td>0.23</td></tr><tr><td>#41</td><td>4.</td><td>4.4 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>50</td><td>0.14</td></tr><tr><td>#42</td><td>4.</td><td>4.4 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>70</td><td>0.1</td></tr><tr><td>#43</td><td>4.</td><td>4.4 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>80</td><td>0.0875</td></tr></tbody></table>

efficient mixing of the whole reservoir and the crystal concentration returns to its initial value, with slight oscillations. Figure 5 shows that the initial concentration in the upper reservoir is slightly higher than its theoretical value of 0.01. This is due to the mesh properties that encompass a larger volume of crystal leading to an initial value of $e = 0.0105$.

### 3. Parametric Study

In the following section, we look at the influence of the density difference between the crystals and the magma, the crystal size, and the magma viscosity on the settling dynamics within a convecting magma ocean.

#### 3.1. Influence of the Density Difference Δρ

During the solidification of the magma ocean, the density of the crystal phase varies relative to the molten phase based on its ability to integrate dense chemical elements into its structure. Within the fully molten mantle, the first crystal to form is bridgmanite. This first crystal is Fe-poor, and achieves neutral chemical buoyancy close to the core mantle boundary (Caracas et al., 2019). The relative Fe content then increases in the melt, potentially leading to the formation of a basal magma ocean overlain by a thick layer of bridgmanite. The level of neutral buoyancy moves toward the shallow mantle as solidification proceeds (Boukaré et al., 2015; Caracas et al., 2019; V. S. Solomatov, 2000). This scenario could lead to fractional crystallization of the early mantle (Patočka et al., 2020). However, if re-entrainment enhanced by vigorous convection is taken into account, crystal dispersion might be an efficient process that would favor batch (or equilibrium) crystallization of the early mantle (V. S. Solomatov et al., 1993). The melt-crystal density crossover in both these scenarios is discussed by Caracas et al. (2019). Here, we monitor the influence of the chemical density difference on the settling dynamics considering relatively small chemical density differences ($10 - 100 \text{ kg.m}^{-3}$), such as are relevant during the early stages of mantle solidification. The chemical density difference is fixed in each simulation and is superimposed onto the thermal density difference which can evolve during cooling. We note that in our models, the density difference related to temperature differences ($\alpha\rho\Delta T = 3 \text{ kg.m}^{-3}$) is smaller than the density difference between the magma and the crystals.

In Figure 6, we illustrate the reference case presented in Figure 4, but in which the density difference has been increased from $\Delta\rho = 50 \text{ kg.m}^{-3}$ to $\Delta\rho = 100 \text{ kg.m}^{-3}$. As already illustrated in Figures 4 and 6 show that because of negative buoyancy, the crystal layer sinks rapidly. After this fast settling at the bottom of the reservoir, particle entrainment occurs from the dunes (10,000 s after initiation of the settling). However, because of a larger density difference between the dispersed phase and the continuous phase, efficient mixing does not occur in this case. Particle entrainment dwindles until the layer is fully sedimented at the bottom of the reservoir (28,000 s after initiation of the settling). Again, the sedimented layer acts as an insulating layer by limiting heat transfer from the bottom of the reservoir. As a consequence, the mean temperature is significantly lower due to this insulating layer at the bottom when full sedimentation is achieved.

Figure 5 compares the time evolution of the concentration of particles in the upper part of the reservoir for the two density differences illustrated in Figures 4 and 6. Throughout the simulation, the particle concentration for $\Delta\rho = 50 \text{ kg.m}^{-3}$ remains higher than the concentration for $\Delta\rho = 100 \text{ kg.m}^{-3}$. This feature is in agreement with Lavorel and Le Bars (2009) who monitored the relative number of particles in suspension as a function of the fluid/
<caption>**Table 2**</caption>

*Continued*

<table><thead><tr><th></th><th>η<sub>c</sub>(Pa.s)</th><th>Ra</th><th>D<sub>d</sub> (m)</th><th>Δρ (kg/m<sup>3</sup>)</th><th>Sh</th></tr></thead><tbody><tr><td>#44</td><td>4.</td><td>4.4 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>100</td><td>0.07</td></tr><tr><td>#45</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>5</td><td>19</td></tr><tr><td>#46</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>10</td><td>9.51</td></tr><tr><td>#47</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>20</td><td>4.76</td></tr><tr><td>#48</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>30</td><td>3.17</td></tr><tr><td>#49</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>40</td><td>2.37</td></tr><tr><td>#50</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>50</td><td>1.9</td></tr><tr><td>#51</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-4</sup></td><td>100</td><td>0.95</td></tr><tr><td>#52</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>5</td><td>1.9</td></tr><tr><td>#53</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>10</td><td>0.95</td></tr><tr><td>#54</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>20</td><td>0.47</td></tr><tr><td>#55</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>30</td><td>0.31</td></tr><tr><td>#56</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>40</td><td>0.24</td></tr><tr><td>#57</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>50</td><td>0.19</td></tr><tr><td>#58</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>70</td><td>0.136</td></tr><tr><td>#59</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>100</td><td>0.095</td></tr><tr><td>#60</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>125</td><td>0.076</td></tr><tr><td>#61</td><td>10</td><td>1.8 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>200</td><td>0.048</td></tr></tbody></table>

*Note.* We calculate the Sh value from Equation 15 by combining Equations 12 and 13 and the expression Nu = 0.14Ra<sup>0.34</sup> obtained from Figure 1 to compute the heat flow F.

particle density difference during the first 1,000 s from a starting scenario of uniform distribution within the reservoir. Figure 5 also shows that increasing the density difference decreases both the minimum value of the particle concentration in the upper reservoir and the time taken to reach this minimum.

## 3.2. Influence of the Crystal Size

Upon solidifying, crystals form following the nucleation and growth mechanisms (V. Solomatov, 2015). The crystal size reached at equilibrium depends on the nature of the mineral and on the cooling rate within the convecting reservoir. The grain size ranges between $10^{-4}$ and $10^{-3}$ m in typical magma ocean conditions (V. Solomatov, 2015). Here, we monitor the influence of the crystal size in our models. It was not possible for us to consider grain sizes larger than our grid mesh size ($8 \times 10^{-3}$ m). Indeed, when the mesh size and crystal size are similar, the Euler-Euler description is no longer relevant because (a) the solid fraction $\phi_d$ as a smooth local average of the volume fraction of the solid phase becomes ill-defined at the scale of the crystal size, and (b) given that the crystals are of comparable size to the flow features, taking their shape into consideration becomes important for solving the flow itself. Hence, we compare the reference case (with $D_d = 10^{-3}$ m) to cases with smaller sizes (down to $D_d = 10^{-4}$ m).

Figure 7 shows the time evolution of the crystal fraction in the upper part of the reservoir for the four particle diameters $D_d$ considered. It shows that for the range of crystal sizes considered ($10^{-4} \le D_d \le 10^{-3}$ m), the time needed for the relative entrainment $e$ to reach its minimum value is very little affected the particle diameter $D_d$ and typically ranges between 305 and 315 s. Moreover the final behavior of the reservoir in terms of sedimentation/suspension is not dependent on the crystal size for the range of sizes considered here. For our range of particle diameters and in the range of controlling parameters corresponding to the suspension regime, the overturn process is rapidly followed by a suspension regime in which the particles are efficiently mixed within the convecting reservoir.

## 3.3. Influence of the Viscosity

The viscosity of the molten phase ($\eta_c$) is a key parameter governing the cooling dynamics of the magma ocean and in particular the solidification timescale (Monteux et al., 2016). According to Karki and Stixrude (2010), $\eta_c$ for anhydrous MgSiO$_3$ melt may decrease down to 0.048 Pa.s at mid-mantle conditions ($P = 70$ GPa and $T = 4000$ K along the magma ocean adiabat). For molten peridotite, Xie et al. (2021) obtained a viscosity ranging between 0.017 and 0.038 Pa.s. It should be noted that vertical velocities are low and the flow is much less turbulent adjacent to horizontal solid boundaries rather than in the center of a deep convecting magma ocean (Elkins-Tanton, 2012). In the context of this study (i.e., close to solid boundaries), the flow is laminar. Hence the characteristic Reynolds number has to be maintained below a value of ~$10^3$ and the range of viscosity values in our models must be chosen accordingly. To maintain laminar flow, we consider $0.1 \le \eta_c \le 10$ Pa.s.

The initial conditions for the continuous phase velocity and temperature fields change as a function of the viscosity. To ensure that the differences we observe between the various viscosity values are not due to differences in the initial conditions, we ran each model for a fixed viscosity value with different time values corresponding to different initial conditions. Figure 8 shows that the continuous phase viscosity governs the ability of the solid phase to mix or settle at the bottom of the reservoir. What is interesting is non linearity of this process. In Figure 8 we can see that for the lowest viscosity, the crystal layer rapidly mixes with the convecting reservoir. However, for $0.4 \le \eta_c \le 1$, the rapid overturn process is not followed by a mixing process, and the final concentration within the upper part of the reservoir does not reach the value typical of efficient mixing. Finally, for $\eta_c > 1$ Pa.s, the overturn process is followed by a vigorous mixing process similarly to the case with $\eta_c = 0.1$ Pa.s. Figure 8 shows that the molten phase viscosity also governs the timing of the early overturn (i.e., the time needed for the particle concentration in the upper part of the reservoir to reach the minimal value). This overturn time increases from nearly 40 s (for $\eta_c = 0.1$ Pa.s) up to 2,000 s (for $\eta_c = 10$ Pa.s).
<figure><img src="image_299.png" /><figcaption>**Figure 4.** Time evolution (from top to bottom) of the dispersed volume concentration (left) and the temperature (right) for the reference case (with Ra = 4.4 × 10⁶ (corresponding to η_c = 4 Pa.s), Δρ = 50 kg.m⁻³, and D_d = 10⁻³m).</figcaption></figure>

<figure><img src="image_300.png" /><figcaption>**Figure 5.** Time evolution of the particles concentration in the upper part of the reservoir with Ra = 4.4 × 10⁶ (corresponding to η_c = 4 Pa.s) and D_d = 10⁻³ m. The red line illustrates the case where Δρ = 50 kg.m⁻³ and the green line the case where Δρ = 100 kg.m⁻³. Vertical dashed lines illustrate the time when the minimum value of e is obtained.</figcaption></figure>

4. A Criterion for the Stability of a Crystallizing Layer

We have illustrated the separate influences of the density difference between crystals and melt, grain size and melt viscosity on the behavior of a crystallizing layer within a magma ocean. We now look at the ability of the flow generated by convection within the magma ocean to initiate movement and suspension of sedimented crystals and to lead to efficient mixing within the magma ocean. We introduce the Shields number Sh which is the ratio of the tangential stress from convective flow over the buoyancy stress relative to the density difference between the crystals and the melt (V. S. Solomatov et al., 1993):

$$ Sh = \left( \frac{\eta_c \alpha g F}{C_p} \right)^{1/2} \frac{1}{\Delta \rho g D_d} \qquad (15) $$

Large Sh values favor a suspension regime while small Sh values favor a sedimentation regime. The two regimes are separated by a critical Shields number value Sh_c. From their analog experiments, V. S. Solomatov et al. (1993) estimated a critical value Sh_c = 0.1 – 0.2. We have run a total of 61 numerical models (see Table 2 for values) with different sets of D_d, Δρ, and η_c whose value ranges are given in Table 1. Each model corresponds to an Sh value
<figure><img src="image_301.png" /><figcaption>**Figure 6.** Time evolution (from top to bottom) of the dispersed volume concentration (left) and the temperature (right) for the case with Ra = 4.4 × 10⁶ (corresponding to η_c = 4 Pa.s), Δρ = 100 kg.m⁻³, and D_d = 10⁻³ m.</figcaption></figure>

<figure><img src="image_302.png" /><figcaption>**Figure 7.** Time evolution of the particle concentration in the upper part of the reservoir with Ra = 4.4 × 10⁶ (corresponding to η_c = 4 Pa.s) and Δρ = 50 kg.m⁻³. The red line illustrates the case with D_d = 10⁻³ m, the green line the case with D_d = 10⁻⁴ m, the blue line the case with D_d = 3 × 10⁻⁴ m and the orange line the case with D_d = 6 × 10⁻⁴ m. Vertical dashed lines illustrate the time when the minimum value of e is obtained.</figcaption></figure>

calculated using a combination of Equations 12 and 13 and the expression Nu = 0.14Ra⁰.³⁴ obtained from Figure 1 to compute the heat flow F (see Table 2 for values). Monitoring the concentration of grains in the upper reservoir (e(t)) in each model makes it possible to determine whether the models end up in the suspension or sedimentation regime.

In Figure 9, we illustrate the Sh values as a function of the molten phase viscosity η_c. Two regimes are considered: one in which the grain concentration in the upper part of the reservoir (e) does not return to its initial value, known as the sedimentation regime; another one in which the grain concentration in the upper part of the reservoir does return to its initial value, known as the suspension regime. As shown in Figure 8, some cases may end up in an intermediate zone where the final concentration is 0.005 < e < 0.01. We consider that these intermediate cases fall into the sedimentation regime. Sedimentation and suspension regimes are separated by the critical Shields number value (Sh_c).

According to Figure 9, the Sh_c value is a non-linear function of the continuous molten phase viscosity. For η_c > 1, Sh_c seems to be constant and independent of the η_c value. However, when η_c < 1, the Sh_c value decreases with decreasing viscosity values following Sh_c = 0.1η_c⁰.⁴⁵. As a consequence Sh_c is not influenced by the increase in Rayleigh number up to Ra = 1.8 × 10⁷. For Ra values higher than 1.8 × 10⁷ the increase in kinetic energy induced by the vigorous convection in the system results in the following decrease in Sh_c:
<figure><img src="image_303.png" /><figcaption>**Figure 8.** Time evolution of the particles concentration in the upper part of the reservoir with $\Delta\rho = 50 \text{ kg.m}^{-3}$ and $D_d = 10^{-3} \text{ m}$. The different color lines illustrate cases with different viscosities for the continuous phase ranging from 0.1 to 10 Pa.s. Vertical dashed lines illustrate the time when the minimum value of $e$ is obtained. Note that we tested the case with $\eta_c = 0.1 \text{ Pa.s}$ with two different meshes (6 and 8 mm mesh) and the difference between the two curves is not significant.</figcaption></figure>

<figure><img src="image_304.png" /><figcaption>**Figure 9.** Shields number value (Sh) as a function of the molten continuous phase viscosity ($\eta_c$) (or Ra for the top x-axis). Square symbols correspond to the models run in this study (see Table 2 for values). Empty black symbols represent cases where the suspension regime is reached. Filled green symbols represent cases where the sedimentation regime is reached. The two regimes are separated by dashed lines (green for $\eta_c \le 1$, blue for $\eta_c \ge 1$) illustrating the critical Shields number (Sh$_c$). Red symbols and the red dotted line illustrate the values for the cases represented in Figures 5, 7, and 8 respectively. We calculate the Sh values by combining Equations 12 and 13 and the expression $\text{Nu} = 0.14\text{Ra}^{0.34}$ obtained from Figure 1 to compute the heat flow $F$ within Equation 15.</figcaption></figure>
$$ Sh_c = 184 \cdot Ra^{-0.45} \qquad (16) $$

The non-linear dependence of Sh_c on viscosity explains why the crystal settling behavior can vary back and forth between sedimentation and suspension regimes when the viscosity increases, as shown with a dotted red line in Figure 9, which shows the path followed by the models detailed in Figure 8. Moreover, Figure 9 shows that Sh_c is smaller than the range of values proposed by V. S. Solomatov et al. (1993) (0.1 – 0.2). This feature is particularly visible for the lowest range of viscosities used in our models. For a magma ocean viscosity value of 0.1 Pa.s, the Sh_c value may decrease down to 0.04 which makes it difficult to envision crystal settling in such a convecting environment, especially given the typical viscosity values for a magma ocean (Karki & Stixrude, 2010; Xie et al., 2021). Decreasing the viscosity down to the lowest values proposed by Xie et al. (2021) may enhance this feature. Finally, we have tested the cases with Ra = 1.8 × 10⁶ (i.e., η_c = 10 Pa.s) (cases referred as #52–#61 in Table 2) considering no-slip boundary conditions. From these models (not shown here), we obtain 0.19 < Sh_c < 0.31 while, when considering free-slip boundary conditions, we obtain 0.095 < Sh_c < 0.136. This result illustrates that free-slip boundary conditions favor suspension regimes compared to no-slip boundary conditions.

## 5. Fate of Mushy Material in the Context of a Progressive Magma Ocean Crystallization

In the previous sections, we showed the influence of the viscosity, crystal size and density difference between melt and crystals on the crystals settling dynamics within a convecting magma ocean. We proposed a criterion (Sh_c) to determine the stability of a crystallizing layer based on the convection/buoyancy stress ratio. We showed that the transition between suspension and sedimentation could occur for Sh_c values lower than 0.1. In the Sh_c expression (see Equation 15), the buoyancy relation between newly formed crystals and the magma ocean itself is a major parameter for determining how the crystal settling and the early mantle differentiation is likely to have occurred. Caracas et al. (2019) have constrained the evolution of the neutral buoyancy depth during magma ocean solidification and the progressive enrichment in Fe content of the residual melt. In particular, they showed that the neutral buoyancy depth moved toward shallower depths during magma ocean solidification, thus favoring the formation of a basal magma ocean. In the following sections, we aim to constrain the magma ocean crystallization dynamics by comparing (a) the critical density difference separating the suspension and sedimentation regimes derived from our detailed numerical approach described above and (b) the density difference between crystals and melt calculated from a plausible range of chemical compositions.

### 5.1. Parameters Controlling the Solidification of a Magma Ocean on an Earth-Like Planet

We investigated the implications of our model on the mechanism of magma ocean (MO) solidification on Earth. More precisely, we investigated the possibility of bridgmanite (Bg) segregation based on its buoyancy (positive or negative). For this, we compared the density difference between the MO and the solid crystal Δρ_MO calculated for different compositions and the critical density contrast Δρ_crit above which crystal segregation is likely to occur. Δρ_crit is derived from Equation 15 by calculating:

$$ \Delta\rho_{crit} = \left( \frac{\eta_c \alpha g F}{C_p} \right)^{1/2} \frac{1}{Sh_c g D_d} \qquad (17) $$

and considering Sh_c = 0.1 which is a conservative value for a magma ocean viscosity of 0.1 Pa.s or less (Xie et al., 2021). At a given depth (z), a Δρ_MO higher than the critical value Δρ_crit would lead to efficient Bg sedimentation at the bottom of the MO. By extension, we also considered that Bg would efficiently segregate into a solid layer above a basal magma ocean for Δρ_MO/Δρ_crit of less than -1. Calculating Δρ_MO requires a number of assumptions:

1. We used pyrolite as a reference model for the mantle composition. Pyrolite is obtained by adding a crustal component to the present-day depleted (peridotitic) upper mantle (Ringwood, 1975). To retain some flexibility in terms of possible mantle composition, we subtracted or added bridgmanite (Bg) to pyrolite, producing a more or less mafic mantle material, respectively. Below, we argue that the composition of the primitive upper mantle (PUM) should be close to that of pyrolite, as is widely accepted.

2. The composition of the bulk silicate Earth (BSE) is still not well established, because a number of deep mantle reservoirs could have a different composition to pyrolite: The large low seismic velocity provinces (LLSVP)
seem to present a different mineralogy (e.g., Vilella et al., 2021), with the deep lower mantle possibly being more like bridgmanitic than pyrolitic (Murakami et al., 2012) and there could be relics of a basal magma ocean above the core-mantle boundary (CMB) (Labrosse et al., 2007) or of early mantle overturns (e.g., Boukaré et al., 2015). In our models, we allowed a lot of variation in terms of the size of the possible Bg-rich or mafic reservoirs that might have segregated in the deep mantle during the solidification of the MO. Variance of the BSE composition was achieved by adding or subtracting up to 40% of Bg to the pyrolite mantle reference. This is a much broader range of chemical variation than the current uncertainty on the major element concentration in the BSE.

3. The Bg composition was calculated by equilibrating it with the melt in the mushy MO, using partitioning coefficients $K(Al)^{Bg/Melt} = 0.5$ and $K(Fe)^{Bg/Melt} = 0.1$ or 0.25. The value of $K(Al)^{Bg/Melt}$ does not significantly impact our calculations. $K(Fe)^{Bg/Melt}$ remains controversial as values in the literature range from ~0.07 to ~0.5 (Andrault et al., 2012; Nomura et al., 2011). The former value might be an underestimate as it was measured in a sample showing extreme chemical segregation along a temperature gradient, which artificially heightens Fe-depletion in Bg at the highest temperatures. In the latter study, the experimental conditions corresponded to solidus-like temperatures, at which Fe is more compatible in Bg. Relevant temperatures for MO solidification cover the range from liquidus-like to intermediate between liquidus and solidus. Therefore, we chose values for $K(Fe)^{Bg/Melt}$ of 0.1 or 0.25 for our calculations.

4. As a reference model for the mushy MO, we used a Bg grain diameter ($D_{Bg}$) of 0.1–10 mm, viscosity ($\eta$) of 0.1 Pa.s (Xie et al., 2020), intrinsic conductivity ($k$) of 5 W/m/K (Hofmeister, 1999) and heat capacity ($Cp$) of 1,800 J/kg/K (Asimow, 2018). We considered that the melt had the same composition at all depths because of efficient mixing in a turbulent MO.

5. We considered that the term $\Delta T$ used in $F$ calculation from Equations 13 and 17 was the temperature difference between the Earth's surface and the bottom of the MO, after correcting for the adiabatic increase in temperature with mantle depth. We initially fixed the surface temperature ($T_{surf}$) to 500 K (Lebrun et al., 2013; Monteux et al., 2020), but later investigated $T_{surf}$ of 1000 and 1500 K. We considered that the bottom of the MO ($z_{MO}$) was the depth where mantle solidification was achieved for a mixture of 60% Bg and 40% melt. During mantle solidification, $z_{MO}$ evolves from the CMB depth toward the surface. We assumed that the temperature at a depth $z_{MO}$ followed the solidus temperature profile of pyrolite ($T_{sol}(z_{MO})$) (Pierru et al., 2022). We note that using the solidus profile yields a conservative $\Delta T$ compared to the true temperature which should lie between the solidus and the liquidus. Then, we calculated $\Delta T(z_{MO})$ of a MO extending from a depth $z_{MO}$ to the Earth's surface using the following equation:

$$ \Delta T(z_{MO}) = T_{sol}(z_{MO}) - T_{surf} - \Delta T_{ad}(z_{MO}) \quad (18) $$

where $\Delta T_{ad}(z_{MO})$ is the temperature difference between $T_{sol}(z_{MO})$ and the corresponding potential surface temperature extrapolated from $T_{sol}(z_{MO})$ to the surface along the MO adiabatic temperature profile. The adiabatic gradient in the MO was calculated using the equations of state available for different silicate melt compositions (Asimow, 2018). $\Delta T(z_{MO})$ is found to range between ~1300 K and ~2000 K for MO extending from 660 km (24 GPa) to 2,800 km (135 GPa), respectively (see Figure 10).

6. We calculated the density of Bg ($\rho_{Bg}$) along the mantle solidus, using the equations of state (EoS) for MgSiO$_3$, FeSiO$_3$, and Al$_2$O$_3$ end-members (Matas et al., 2007). The breakdown into end-members assumes Bg of ($Mg_{1-x}Fe_xAl_y$)($Si_{1-y}Al_y$)O$_3$ stoichiometry with a dominant fraction of Fe$^{2+}$ over Fe$^{3+}$ at lower mantle conditions, in agreement with previous works (e.g., Andrault et al., 2018). This assumption has little impact on the Bg density; the dominant parameter is the global Fe content in Bg, itself controlled by $K(Fe)^{Bg/Melt}$ (discussed above). At the same P-T conditions, we calculated the MO density ($\rho_{Melt}$) and thermal expansion ($\alpha_{Melt}$) using the EoS of the molten end-member MgSiO$_3$ enstatite, Mg$_2$SiO$_4$ forsterite, Fe$_2$SiO$_4$ fayalite, CaMgSi$_2$O$_6$ diopside, and CaAl$_2$Si$_2$O$_8$ anorthite (Asimow, 2018). These calculations yield the Bg-to-melt density contrast ($\Delta\rho_{MO} = \rho_{Bg} - \rho_{Melt}$, Figure 11) along the mantle solidus: $\Delta\rho$ decreases with increasing pressure regardless of the composition of the mushy mantle, because the melt is more compressible than Bg.

7. We did not attempt to model the final stages of MO crystallization, when the bottom of the MO would be shallower than 660 km depth, because the phase relations rely on the melting diagrams at upper mantle P-T conditions.
<figure><img src="image_305.png" /><figcaption>**Figure 10.** Temperature profiles corresponding to (gray thick line) the solidus of the pyrolite mantle (Pierru et al., 2022) and (colored lines) adiabatic gradients of MO for potential surface temperatures from 1800 to 2500 K.</figcaption></figure>

## 5.2. Conditions for Segregation of Bg in the Mushy Magma Ocean

### 5.2.1. Effect of Melt Composition

We then calculated the ratio of density differences $\chi_{\rho} = \Delta\rho_{MO}/\Delta\rho_{crit}$ at all lower mantle depths using the parameters defined above. First, we fixed the grain size of Bg to 1 mm, $K(Fe)^{Bg/Melt}$ to 0.1 or 0.25, and varied the MO composition from pyrolite –40% Bg to pyrolite +40% Bg (i.e., 70% pyrolite and 30% Bg). $\chi_{\rho}$ was either positive or negative (Figure 12), depending on the solid-melt density difference (Figure 11). A $\Delta\rho_{MO}$ value greater than $\Delta\rho_{crit}$ (i.e., $\chi_{\rho} > 1$) indicates efficient Bg segregation to the bottom of the MO. Our calculations show that Bg

<figure><img src="image_306.png" /><figcaption>**Figure 11.** Solid-melt density difference ($\Delta\rho_{MO} = \rho_{Bg} - \rho_{melt}$) calculated using the P-V-T equation of state for Bg (Matas et al., 2007) and MO (Asimow, 2018), after decomposition into end-members. Calculations are performed for several MO compositions by adding or subtracting a fraction of Bg to a pyrolite-based mantle reference, and for Fe partition coefficients between Bg and MO ($K(Fe)^{Bg/Melt}$) of 0.1 or 0.25. Positive values indicate lighter melt.</figcaption></figure>
<figure><img src="image_307.png" /><figcaption>**Figure 12.** Density contrast ratio $\Delta\rho_{MO}/\Delta\rho_{crit}$ in the MO. Values greater than 1, or less than -1, indicate efficient Bg segregation below, or above, a mushy MO, respectively.</figcaption></figure>

segregation happens in shallow MOs for all chemical compositions investigated. For pyrolite, this is in agreement with previous work to determine the MO viscosity (Xie et al., 2020). The higher the Bg content in the MO, the earlier (and deeper) Bg segregation occurs during the MO solidification; this is because $\Delta\rho_{MO}$ is higher when the bulk Fe-content is lower. As expected, lower $K(Fe)^{Bg/Melt}$ favors a denser melt, and thus lower $\Delta\rho_{MO}$. Still, our calculations show a moderate impact on the $K(Fe)^{Bg/Melt}$ variations, from 0.1 to 0.25. Finally, efficient Bg segregation could occur above a dense MO for an MO composition highly depleted in Bg; this is suggested by $\chi_{\rho} < -1$ for pyrolite -40%Bg at pressures greater than ~120 GPa, for a Bg grain size of 1 mm (Figure 12) or more (see below).

### 5.2.2. Effect of Grain Size

We also investigated the role of the Bg grain size, for five different mantle compositions (Figure 13). First we focused on a melt of pyrolite composition. $\chi_{\rho}$ is found to be positive at most mantle depths, because $\Delta\rho_{MO}$ is also positive (Figure 11), except in the lowermost mantle for $K(Fe)^{Bg/Melt} = 0.1$. Still, $\chi_{\rho}$ can only be greater than 1 for $D_{Bg}$ of more than 0.5 mm; Bg-grains of for example, 2 mm segregate when pressure (depth) is lower than 70 GPa (1,650 km). Grain size in a crystallizing MO remains a complicated issue (V. Solomatov, 2015): (a) as long as the MO is fully liquid close to the planetary surface, grains are often remelted within the large scale convection movements, meaning that grain size is controlled by nucleation processes; (b) When the crystallization zone reaches the surface, it enables some grains to grow indefinitely. Depending on the curves of the mantle solidus and liquidus, this change in mechanism could happen when the bottom of the MO reaches a depth of 1,000 km (40 GPa) (see Xie et al., 2020). In this article, we are aiming to model the chemical segregation in the initial stages of the MO solidification, when grain size is controlled by nucleation processes. $D_{Bg}$ in these conditions depends largely on the cooling rate of the MO, with estimated values of 0.1–1 mm (V. Solomatov, 2015). For a grain size of 1 mm, $\chi_{\rho}$ is higher than unity at pressures (depths) lower than 45 GPa (1,100 km) (Figure 12). As a result, it appears unlikely that Bg grains could segregate in a crystallizing MO of pyrolite composition at depths greater than ~1,100 km.

A BSE composition enriched in Bg yields a lighter melt during the early stages of MO solidification, and thus a higher $\Delta\rho_{MO}$. It facilitates the segregation of Bg at the bottom of the mushy MO. For a melt composition enriched by 20% or 40% of Bg, Bg grains of 1 mm could efficiently segregate at pressures (depths) lower than 55 GPa (1,350 km) or 65 GPa (1,550 km), respectively (Figures 13a and 13b). From a geochemical point of view, early segregation of Bg to the CMB requires a BSE composition enriched in Bg if the late shallow MO is of pyrolite composition; we will discuss this matter later in the article.
<figure><img src="image_308.png" /><figcaption>**Figure 13.** Density contrast ratio $\Delta\rho_{MO}/\Delta\rho_{crit}$ in the MO calculated as a function of the Bg grain size, from 0.1 to 5 mm, for 5 different melt compositions, from pyrolite –40%Bg to pyrolite +40%Bg. Values greater than 1 or less than –1 indicate efficient Bg segregation below, or above, the mushy MO, respectively. $\chi_p = 0$ is found for $\Delta\rho_{MO} = 0$ (Figure 11). It indicates an infinite value for Sh and perfect crystal suspension (Figure 9). Continuous and dashed lines indicate $K(Fe)^{Bg/Melt}$ of 0.1 and 0.25, respectively.</figcaption></figure>

Alternatively, a BSE composition depleted in Bg yields a denser melt in the mushy MO, as discussed in several works (Ballmer, Lourenço, et al., 2017; Boukaré et al., 2015; Caracas et al., 2019). Our calculations of $\chi_p$ for melt compositions of pyrolite –20% and –40% Bg show that Bg is indeed less dense than the melts over a large fraction of the lower mantle (Figure 11). Bg remains denser at low pressures, due to its very high bulk modulus. For a Bg grain size of 1 mm, values of $\chi_p$ are found to be lower than –1 at pressures (depths) greater than 120 GPa (2,600 km) for a melt composition of pyrolite –40% of Bg (Figure 13e). The pressure (depth) at which Bg could segregate upwards decreases with increasing Bg grain size. From a geochemical point of view, early segregation close to the CMB of a Bg-depleted melt requires a BSE composition depleted in Bg, if the shallow MO were to be pyrolitic.

### 5.2.3. Effect of Temperature

Another important parameter in $\Delta\rho_{crit}$ is $\Delta T$ (Equation 18), which is intrinsically related to the rate of MO cooling. During the mushy stage of MO solidification, the surface temperature is $\approx$500 K, which leads to final MO solidification 100–200 thousand years after the Moon-forming impact (Lebrun et al., 2013; Monteux et al., 2020). Presence of a blanketing atmosphere above the MO would induce a higher surface temperature and a lower rate of cooling (Salvador et al., 2017; Sleep et al., 2014). We note that the formation of a blanketing atmosphere is unlikely to affect the early stages of MO solidification for a number of reasons: (a) H$_2$O (a highly efficient blanketing compound) is only efficiently vaporised into the atmosphere after a high degree of MO solidification (Elkins-Tanton, 2008; Salvador et al., 2017). (b) A surface temperature higher than 1500 K yields an atmosphere dominated by rock-vapor, which is not efficient at stopping the surface cooling by thermal radiation (Sleep et al., 2014). For these reasons, we also performed calculations of $\Delta\rho_{MO}$ for surface temperatures of 500, 1000, and 1500 K (Figure 14).
<figure><img src="image_309.png" /><figcaption>**Figure 14.** Density contrast ratio $\Delta\rho_{MO}/\Delta\rho_{crit}$ calculated for various surface temperatures (Equation 18). Thick, intermediate and thin lines correspond to 500, 1000, and 1500 K, respectively. The calculations are also performed for various mantle compositions and $K(Fe)^{Bg/Melt}$ values.</figcaption></figure>

Calculations show that higher surface temperatures yield more extreme positive and negative values of $\chi_p$ and therefore easier Bg segregation. This is because higher $T_{surf}$ corresponds to a smaller temperature difference $\Delta T$ (Equation 18) and thus a less turbulent MO. Still, the impact remains modest; for a pyrolite mantle, a variation in $T_{surf}$ from 500 to 1500 K yields Bg segregation at pressures (depths) lower than 45 GPa (1,100 km) and 67 GPa (1,600 km), respectively.

## 5.3. Possible Scenarios for MO Solidification

### 5.3.1. Conditions for Bg Segregation at the Bottom of the MO

Geochemical scenarios of the Earth's formation involve building blocks of chondritic composition. Uncertainties remain about the composition of the primitive "chondritic-type" mantle after the core-mantle segregation. Still, it is very likely that it would have been enriched in $SiO_2$, compared to pyrolite (Allègre et al., 1995; Javoy, 1995; McDonough & Sun, 1995). Such enrichment has raised long-running discussions about how much Si the core might contain and how much the deep mantle might be enriched in $SiO_2$, compared to the present day upper mantle. A high $MgSiO_3$/MgO ratio in the lower mantle is compatible with results from inversions of seismic profiles based on elastic parameters of minerals (Matas et al., 2007; Murakami et al., 2012; Samuel et al., 2005).

Our model shows that increasing the Bg content in BSE favors the segregation of Bg in a deep MO (Figures 13a and 13b). Bg grains of 1 mm diameter segregate efficiently (i.e., $\chi_p > 1$) at pressures (depths) lower than 45 GPa (1,100 km), 55 GPa (1,350 km), and 65 GPa (1,550 km) for pyrolite, pyrolite with addition of +20%Bg and +40%Bg, respectively. Still, Bg grains would segregate during the early stages of the solidification if the grain size were 5 mm or more. Such large grains are not expected in a very thick MO (V. Solomatov, 2015). Bg segregation becomes possible when the MO is shallower, and the greater the Bg enrichment, the earlier (and deeper) the Bg-segregation occurs. For a grain size of 1 mm, a Bg layer could form at MO depths of 1,100–1,550 km, depending on the BSE composition. Such a thin MO could also favor grain growth, because the MO would be partially crystallized up to the Earth's surface (V. Solomatov, 2015; Xie et al., 2020). This Bg-layer could yield silica-enriched domains (Ballmer, Houser, et al., 2017). Unfortunately, our models are not able to assess the effects of changes in the thickness of this layer, because there are too many adjustable parameters. For a BSE enriched in Bg, segregation of Bg in the lower mantle would drive the MO toward a pyrolite mantle composition. Decreasing the Bg-content would then make the segregation of Bg grains less efficient (Figure 12).
### 5.3.2. Solidification of a Pyrolite-Type Bulk Silicate Earth

A pyrolite-like BSE would solidify with no Bg segregation from the melt, as long as the Bg grain size remained lower than 0.5 mm (Figure 13c). Bg segregation could still happen in the uppermost lower mantle if the Bg grain size were to increase. We note, however, that the grain growth of olivine and other minerals of the upper mantle may have a moderate impact on the grain size of Bg in the lower mantle. Modeling the mechanism of MO solidification in the upper mantle is beyond the scope of this article.

Interestingly, the absence of Bg segregation does not imply a chemically homogeneous mantle once it is entirely solidified. The origin of this effect is that the solidified mantle contains ~60% Bg and ~40% melt, while the overlying mushy MO contains a different fraction of the two phases. Here we assume a rate of mushy-mantle compaction (with extraction of the melt) that is significantly slower than the rate of MO solidification, following previous work (V. Solomatov, 2015). We modeled the chemical stratification resulting from homogeneous solidification as follows: (a) we chose a BSE composition close to pyrolite, albeit with an ajustable Bg content, to allow the final shallow MO (thickness of ~660 km) composition to be as close as possible to that of pyrolite. (b) At each time step of the progressive MO solidification, we estimated the fraction of Bg grains in the MO above the bottom of the MO ($\phi_{MO}$) based on our previous model of progressive MO solidification (Figure 3 of Monteux et al. (2016)). (c) We calculated the composition of the coexisting melt and Bg grains, fixing $K(Fe)^{Bg/Melt} = 0.25$. (The value of $K(Fe)^{Bg/Melt}$ has no significant impact on this calculation in any case). (d) Starting from the CMB, we cooled the mantle and calculated the composition of a mixture of 60% Bg grains and 40% melt that solidified at the bottom of the MO. The chemical content of this basal layer was then subtracted from the composition of the overlying mushy MO, before a next step of the solidification was performed within a thinner MO.

This calculation shows significant changes in the SiO₂, FeO, and CaO contents with depth in the solidified mantle (Figure 15). They occur due to Bg having higher SiO₂ and lower FeO and CaO contents compared to the melt. The shallow MO has a chemical composition close to that of pyrolite, when the BSE composition is adjusted to 80% pyrolite and 20% Bg. The SiO₂ content in the deep mantle is found to decrease from ~46% to 39%, which translates into a decrease from ~90% to 70% in Bg-content, with decreasing depth from 2,900 to 660 km, respectively. The Fe-content also increases with decreasing depth, because of its incompatible character; this effect is found even greater when using $K(Fe)^{Bg/Melt} = 0.1$, instead of $K(Fe)^{Bg/Melt} = 0.25$. However, the overlying mushy MO does not become denser than the solidified layers, because Bg is denser than the melt at all mantle depths for Bg-enriched BSE compositions (Figure 11).

Such chemical layering would only survive in a mantle that became stagnant after solidification. Instead, one could expect significant vertical mixing in a mushy mantle with a relatively low viscosity and with a major heat flux from the core to the surface (Monteux et al., 2016). Another effect that might moderate the chemical stratification is simultaneous mantle solidification over a large range of mantle depths. Upon cooling, the temperature profile becomes almost parallel to the solidus (Miller et al., 1991; Monteux et al., 2016; Stixrude et al., 2009), which could result in bulk solidification. However, these effects may not entirely suppress the chemical gradients, in particular between a shallow MO and the lower mantle. The progressive enrichment of Bg with mantle depth is in good agreement with the mineralogical models of, for example, Samuel et al. (2005); Matas et al. (2007); Murakami et al. (2012), as well as with the difference between the shallow mantle composition today (pyrolite) and a "chondritic-type" BSE (see above).

We note that the scenario of Bg-segregation in the lowermost mantle does not help explain the LLSVPs above the CMB. Several origins have been proposed for the LLSVP, in particular primitive mantle that did not melt following the Moon forming impact. Primitive mantle compositions imply a relatively low proportion of ferropericlase, which is compatible with the comparison between the LLSVP seismic signature and the elastic parameters of the minerals (Vilella et al., 2021). Still, LLSVPs represent a quite small fraction of the mantle volume, typically 1.5%-2.5% (Burke et al., 2008). Therefore, our models would not be much different if some primitive mantle were to have survived melting during the Moon-Forming impact.

### 5.3.3. Conditions for the Formation of a Primary Basal Magma Ocean

The segregation of large Bg grains at the bottom of a pyrolite mantle would yield a mafic Fe-enriched MO on top of this Bg solid layer. Following this scenario, it has been proposed that the melt could become denser than the underlying mantle and induce early mantle overturns (Ballmer, Lourenço, et al., 2017; Boukaré et al., 2015; Caracas et al., 2019). Such a reversal could result in a secondary basal MO in the lowermost mantle with a
<figure><img src="image_310.png" /><figcaption>**Figure 15.** Chemical profiles obtained for the case in which there is no Bg segregation from the melt. The chemical gradient comes from the solidification of ~60% Bg and 40% of the melt at the bottom of the MO. These profiles assume no mantle mixing once the mantle layers are eventually solidified. “BSE” and “s-MO” are bulk silicate Earth and a shallow MO extending up to 660 km depth, respectively. Here, the BSE composition is adjusted to 80% pyrolite and 20% Bg.</figcaption></figure>

composition typical of the shallow MO. Our results do not strongly support this scenario, because the lower the Bg-content in the mushy MO, the less efficient the Bg segregation (Figure 13). Therefore, bulk solidification is more likely, after a certain amount of Bg-segregation in a thick MO. Mantle overturns may still occur due to the crystallization processes in the upper mantle; this should be investigated in future work based on analysis of the melting diagrams in the upper mantle.

Formation of a basal MO during the solidification process (i.e., a primary basal MO) would require efficient segregation of Bg into an overlying solid layer (Labrosse et al., 2007). These conditions translate into $\chi_p < -1$ in our model. Such conditions are reached for a mafic BSE composition, which has a lower Bg content compared to pyrolite (Figures 13c and 13d). However, our calculations show that upward Bg segregation is only expected (a) in the lowermost part of the MO, (b) for a BSE composition depleted in more than –20% Bg compared to pyrolite, and (c) for Bg grain sizes larger than 0.5 mm. We note that low $K(Fe)^{Bg/Melt}$ values at these conditions help make the melt denser than Bg (Figure 12, Nomura et al. (2011)). For an MO composed of pyrolite –40% Bg, Bg grains of 1 mm diameter would aggregate efficiently into a Bg-rich layer at a pressure (depth) of 120 GPa (2,600 km). For an MO composed of pyrolite –20% Bg, only Bg grains of 5 mm diameter could efficiently segregate upwards (Figure 13d). The simultaneous formation of a basal MO and a solid Bg layer may ensure that the