# Mixing time of heterogeneities in a buoyancy-dominated magma ocean

B Thomas, H Samuel, C G Farnetani, J Aubert, C Chauvel

► To cite this version:

B Thomas, H Samuel, C G Farnetani, J Aubert, C Chauvel. Mixing time of heterogeneities in a buoyancy-dominated magma ocean. Geophysical Journal International, 2024, 236 (2), pp.764-777. 10.1093/gji/ggad452. hal-04737219

HAL Id: hal-04737219

https://hal.science/hal-04737219v1

Submitted on 15 Oct 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

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

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

# Mixing time of heterogeneities in a buoyancy-dominated magma ocean

B. Thomas<sup>●</sup>, H. Samuel<sup>●</sup>, C. G. Farnetani, J. Aubert<sup>●</sup> and C. Chauvel

Université Paris Cité, Institut de physique du globe de Paris, CNRS, F-75005 Paris, France. E-mail: bthomas@ipgp.fr

Accepted 2023 November 14. Received 2023 November 14; in original form 2023 June 5

## SUMMARY

During the accretion stage, large impacts provided sufficient energy to melt the entire mantle into a terrestrial magma ocean. Processes occurring in the magma ocean may have led to the formation of heterogeneities still found in modern ocean island basalts. So far, no definitive mechanism exists to explain the survival of early heterogeneities for approximately 4.5 Ga. Addressing this question requires understanding the efficiency of convective mixing during both the early molten and the solid-state stages experienced by the Earth's mantle. While mixing in the solid mantle and in an essentially crystallized magma ocean has been relatively well documented, the efficiency of convective mixing in a liquid magma ocean has received less attention. In this paper we characterized the mixing efficiency of a convecting fluid in a rotating spherical shell, accounting for inertial effects, by computing finite-time Lyapunov exponents (i.e. the Lagrangian strain rate). We conducted a series of numerical experiments for a regime where the influence of the buoyancy force dominates that of rotation and we derived scaling laws to predict the mixing efficiency. We found that for a terrestrial magma ocean, in its fully liquid state, mixing time is of the order of a few minutes or less, even for initially large (~1000 km) heterogeneities. Therefore, passive early mantle heterogeneities cannot survive a fully molten magma ocean stage. This suggests that short-lived heterogeneities (e.g. $^{182}$Hf-$^{182}$W) were either created at the end of the accretional stage, or were stored in deeper regions of the Earth.

**Key words:** Numerical modelling; Planetary interiors; Magma ocean; Mixing.

## 1 INTRODUCTION

The compositionally heterogeneous nature of the Earth's mantle is revealed by geochemical differences recorded between mid-ocean ridge basalts (MORBs) and ocean island basalts (OIB). Isotopic systems with long-lived radioactive parent elements (i.e. half-life of billions of years) inform us on mantle processes that occurred over long geological timescales. For example, the long-term history of melt extraction, the subduction of ancient oceanic crust and continental sediments, along with their storage prior to their entrainment by mantle plumes (Hofmann & White 1982).

Recently, heterogeneous mantle domains have been defined based on isotopic ratios from short-lived radioactive parent elements (Touboul et al. 2012; Horan et al. 2018; Rizo et al. 2019). For these parent-daughter pairs, the parent element has such a short half-life that it became extinct soon after the formation of the solar system. Consequently, isotopic anomalies observed for these systems can only result from fractionation processes that occurred very early in the Earth's history. The Hafnium-Tungsten (Hf-W) systematics offers a clear example: $^{182}$Hf has a half-life of ~8.9 Myr, thus after 50 Myr no more daughter isotope $^{182}$W could be produced. It was

therefore a major discovery to find anomalous $\mu^{182}$W (i.e. ppm deviations relative to a terrestrial standard) in present-day lavas from Hawaii, Samoa, Iceland, Galápagos and La Réunion (Mundl et al. 2017; Rizo et al. 2019; Mundl-Petermeier et al. 2019; Jackson et al. 2020; Peters et al. 2021).

Despite the valuable information conveyed by isotopic systems, geochemical data provide indirect constraints on mantle heterogeneity, because they rely on the interpretation of time-integrated data collected at specific locations at the Earth's surface. A consistent interpretation of geochemical data requires accounting for the evolving convective dynamics of the mantle. This is because convection generates deformations that progressively homogenize heterogeneities, via the repeated action of stretching and folding combined with diffusion at smaller length scales (Ottino 1989; Olson et al. 1984; Gurnis 1986; Hoffman & McKenzie 1985). In this framework, the isotopic signatures from long-lived radioactive parent elements have raised the question of the ability of heterogeneities to survive in a solid, yet vigorously convecting mantle. This has been the subject of a number of studies since the 1980s (e.g. Tackley 2015, and references therein).
However, the isotopic anomalies recently measured for the short-lived radiogenic systems pose a new, fundamental question of the ability of ancient (Hadean) heterogeneities to survive through most of the Earth's history. The physical conditions that prevailed during the early Earth were considerably different from that of the solid-state convecting mantle. Indeed, this time period is characterized by extreme events that inevitably resulted in extensive melting episodes. Impacts of planetesimals and planetary embryos generated enormous amounts of heat during the end of the accretion phase. In particular, the giant Moon-forming impact likely melted the entire Earth and formed a global silicate magma ocean (Stevenson 1987; Solomatov & Stevenson 1993; Abe 1997; Canup 2008; Rubie et al. 2015). Smaller and repeated impacts may also have caused multiple events of shallow magma oceans (Elkins-Tanton 2012). These substantial sources of heat were supplemented by intense early radiogenic heating and the gravitational energy released by core formation, implying that the Earth most likely experienced at least one global magma ocean stage (Senshu et al. 2002; Monteux et al. 2009; Samuel et al. 2010).

The presence of a magma ocean leads to radically different convective dynamics compared to that of the solid-state mantle. The key quantity responsible for these differences is the silicate dynamic viscosity, $\eta$. At solid-state, experiments (Karato & Wu 1993) and postglacial rebound studies (Mitrovica & Forte 2004) infer $\eta \sim 10^{20}-10^{22}$ Pa s. While these values are relatively high, they nonetheless lead to Rayleigh numbers $Ra \propto 1/\eta$ (characterizing the convective vigour) ranging between $10^6$ and $10^8$. This enables the mantle to convect vigorously, but in the laminar (creep) regime where the Prandtl number $Pr \propto \eta$ (whose inverse expresses the importance of inertia) is larger than $10^{21}$.

For a fully molten silicate mantle, experiments (Rubie et al. 2003; Liebske et al. 2005) and ab initio calculations (Karki & Stixrud 2010) indicate that the viscosity drops by 20 orders of magnitude or more ($\eta \sim 10^{-3}-10^2$ Pa s) compared to that of the solid-state mantle. Such low viscosities lead to the following distinct dynamics: (i) The convective vigour in a magma ocean at $Ra \sim 10^{30}$ (Solomatov 2015) results in convective motions that are orders of magnitude faster than in the solid-state case. (ii) With $Pr \sim 10^{-1}-10^4$, about 20 orders of magnitude smaller than for the solid-state mantle, inertia can no longer be neglected. This leads to a different style of convective dynamics (Grossmann & Lohse 2000) where turbulence, otherwise absent in a solid-state mantle, now plays a dominant role. (iii) The low magma ocean viscosity implies that rotational effects may affect the dynamics, with Ekman number $Ek \propto \eta$ (expressing the relative importance of viscous forces to the Coriolis force) as low as $Ek \sim 10^{-13}-10^{-14}$ (Maas et al. 2021).

The extremely turbulent flow of the magma ocean transports heat very efficiently compared to the laminar flow and the slow convective cooling of the solid mantle. Consequently, the lifetime of a magma ocean is short. 1-D thermal evolution models have shown that a magma ocean could last from 1000 yr to a million years in the presence of a thick atmosphere (Solomatov 2000; Lebrun et al. 2013; Nikolaou et al. 2019; Salvador et al. 2017; Lichtenberg et al. 2021).

Previous works focusing on convective mixing processes in the solid mantle and in the magma ocean, have addressed different questions. The solid mantle stage is by far the most studied because it is more easily accessible to numerical and analogue modelling. Pioneering 2-D simulations (Hoffman & McKenzie 1985; Olson et al. 1984; Christensen 1989) showed that the size of the heterogeneities decreases exponentially with time, leading to efficient mixing, in agreement with conclusions in Gurnis (1986) and in Kellogg &

Turcotte (1987) that a few kilometres thick subducted crust gets reduced to a few centimetres by convective stirring. 3-D kinematic models, accounting for both poloidal and toroidal flow components, revealed that well-mixed domains, with exponential stretching, can coexist next to poorly mixed domains, where laminar stretching prevails (Ferrachat & Ricard 1998). However, there is a general agreement (van Keken et al. 2002; Coltice & Schmalzl 2006; Tackley 2015) that passive heterogeneities (i.e. heterogeneities that have similar properties than the surrounding fluid, hence not affecting the dynamics) have a short lifetime with respect to the age of the Earth. Studies that considered rheological heterogeneities found that more viscous mantle domains are less stretched and can be preserved for longer timescales (Manga 1996; Becker et al. 1999; Ballmer et al. 2017; Gülcher et al. 2020). Similarly, compositionally denser heterogeneities can remain unmixed for longer periods, as shown in Farnetani & Samuel (2003) who relied on Lyapunov exponents to measure mixing efficiency. The aforementioned studies only account for large-scale convective motions, nonetheless small-scale convection can significantly affect the mixing time by reducing it from billions of years to hundreds of millions of years (Samuel & King 2014). Other studies have investigated the influence of viscosity reduction induced by either temperature increase from the presence of insulating continental lids (Samuel et al. 2011) and hotter early evolutionary stages (Maurice et al. 2017), or associated with lowermost mantle phase changes (Samuel & Tosi 2012). They found that small effective viscosity strongly enhances stirring efficiency by increasing the magnitude of convective velocities and the time dependence of the flow. Earlier studies (e.g. Olson et al. 1984; Gurnis 1986) that underlined the positive correlation between the magnitude of convective velocities and stirring efficiency are in agreement with this result.

In the context of magma oceans, analogue experiments of impact-induced mixing (Deguen et al. 2011, 2014; Landeau et al. 2021), geodynamic considerations and numerical modelling (Rubie et al. 2003; Ricard et al. 2009; Samuel et al. 2010) focused on two main issues: The first one is to quantify the degree of chemical and isotopic equilibration between the metallic core of the impactor and the silicate magma ocean, which affects the partitioning of chemical species between the core and the mantle (see Carlson et al. 2014, and references therein). Deguen et al. (2014) showed that if the impactor size is small compared to the magma ocean depth, the equilibration can be fast, because turbulent mixing increases the metal-silicate interfacial area. Instead, for large impactor size, the flow evolves towards metallic diapirs that sink rapidly, with limited chemical interactions (Samuel et al. 2010; Dahl & Stevenson 2010; Wacheul et al. 2014). Recently, Landeau et al. (2021) derived scaling laws from impact-induced mixing experiments and found that the mass transfer between metal and silicates is larger (up to 20 times) than previous estimates. Although the above-mentioned works investigated mixing processes in terrestrial magma oceans, they did not account for natural convection or rotation. By accounting for both convective motions and rotational effects in spherical geometry Maas et al. (2021) revealed that the extent of metal-silicate equilibration of incoming impactors in a magma ocean also depends on the latitude of the impact. This could provide a mechanism to generate zonal heterogeneities, because at high latitudes the impact-delivered metal droplets have a limited dispersion in the magma ocean and their settling is fast, contrary to what is observed for low latitude impacts.

Similarly, other studies that considered convective dynamics in a magma ocean focused on the interplay between newly formed crystals and the convecting magma ocean. Numerical simulations in
Cartesian domains (Maas & Hansen 2015) or in spherical geometry (Maas & Hansen 2019) revealed that the depth and the location of crystal accumulation are affected by the rotational strength of the planet. These analogue and numerical experiments have generally showed that it is possible to create heterogeneities in a silicate magma ocean.

To the best of our knowledge, no studies have characterized the mixing of passive heterogeneities associated with the convective dynamics of a fully molten terrestrial magma ocean. Therefore, in this paper we conducted numerical experiments of convection in a fully molten medium in a spherical shell to characterize the stirring efficiency. We focus on the regime where the Coriolis force (rotation) is weak compared to buoyancy forcing (convection), which is plausible for a terrestrial magma ocean (Solomatov 2015; Hoink et al. 2006; Moeller & Hansen 2013). Uncertainties on the values of the physical quantities associated with magma oceans (e.g. viscosity, rotation period) imply that regimes where rotation dominates the dynamics may also be relevant to terrestrial magma oceans (e.g. Moeller & Hansen 2013; Maas & Hansen 2015, 2019) however, we will defer these cases to a future study. We derived scaling laws to predict the mixing times for passive heterogeneities in a convecting fluid with homogeneous physical properties (e.g. viscosity or density).

We first describe the numerical model used for liquid-state convection, our computational approach to measure the stirring efficiency, and the benchmark we performed to validate our method. Then, we discuss the structures of the different regimes of convection observed in our simulations. Focusing on the regime of convection dominated by the buoyancy force, we derive scaling laws for the characteristics of the flow and verify its relevance to a magma ocean's conditions. Following this, we define the processes which control stirring and develop a scaling law for the stirring efficiency. Finally, we apply our scaling laws to predict the mixing time for passive heterogeneities of various initial sizes in a vigorously convecting terrestrial magma ocean, and discuss the implications on the origin and the storage of early heterogeneities imprinted in present-day lavas.

## 2 NUMERICAL MODEL AND METHODS

We model finite Prandtl number convection under the Boussinesq approximation in a rotating spherical shell. The dimensionless inner radius $r_i = 1$ and outer radius $r_o = 2$, leading to an aspect ratio $r_i/r_o = 0.5$. Both boundaries are kept at a constant temperature: $T_o = 0$ at the top and $T_i = \Delta T$ at the bottom, where $\Delta T$ is the superadiabatic temperature difference that drives convective motions. We use free-slip boundary conditions that are appropriate for a fully molten mantle with no stagnant lid at the top, and with a liquid core at the bottom. To normalize the governing equations we use the length-scale $L = r_o - r_i$, which corresponds to the thickness of the convective shell, the temperature scale $\Delta T$ and the time scale $L^2/\nu$, where the kinematic viscosity $\nu$ is the ratio of dynamic viscosity over density ($\nu = \eta/\rho$).

The dimensionless equations that describe this system are:

The conservation of mass:

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

the conservation of energy:

$$ \partial_t T + \mathbf{u} \cdot \nabla T = \frac{1}{Pr} \nabla^2 T, \qquad (2) $$

and the conservation of momentum:

$$ \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} - \nabla^2 \mathbf{u} + \frac{2}{Ek} \mathbf{e}_z \times \mathbf{u} + \nabla P = \frac{Ra}{Pr} \frac{\mathbf{r}}{r_o} T, \qquad (3) $$

where $\mathbf{u}$ is the velocity vector, $T$ is the temperature, $P$ is the dynamic pressure, $\mathbf{e}_z$ is a unit vector parallel to the direction of the rotation axis and $\mathbf{r}$ the unit vector parallel to gravity. From left to right, the terms of the momentum equation correspond to inertial effects (acceleration and advection), momentum diffusion, the Coriolis force, the driving of the flow via pressure gradients, and the buoyancy force. Even though the focus of this study is the regime where rotational effects are negligible, we purposely kept the Coriolis term in the equation to ensure a posteriori that we did not artificially bias the dynamics by not including the rotational effect in the mathematical formulation of the problem.

The dimensionless governing parameters that appear in the above equations are the Rayleigh number $Ra$, the Prandtl number $Pr$ and the Ekman number $Ek$:

$$ Ra = \frac{\alpha g \Delta T L^3}{\nu \kappa}, \qquad Pr = \frac{\nu}{\kappa}, \qquad Ek = \frac{\nu}{\Omega L^2}, \qquad (4) $$

where $\alpha$ is the thermal expansion coefficient, $g$ the gravity acceleration, $\kappa$ the thermal diffusivity, which are all assumed to be constant. Similarly, the viscosity $\eta$ is assumed to be constant. $Ra$ expresses the vigour of convection and corresponds to the ratio of the driving thermal buoyancy force over resisting effects (thermal and momentum diffusion). $Pr$ is the ratio of momentum diffusivity over heat diffusivity, and its inverse characterizes the importance of inertial effects. Finally, $Ek$ is the ratio of the viscous force over the Coriolis force, where $\Omega (\text{rad s}^{-1})$ is the rotational speed about $\mathbf{e}_z$. Hence, rotational effects increase with decreasing the Ekman number.

The spherical geometry we considered is motivated by several reasons: (1) global terrestrial magma oceans are confined to spherical shells, (2) this geometry is the most suitable to include the effect of the Coriolis force on the global scale, (3) at high Reynolds number values expected for a fully molten magma ocean turbulence is a fully 3-D phenomenon and (4) even in absence of rotation, Cartesian geometries with periodic side boundaries generate strong zonal flows, which prevent the development of convective motions (Wang 2020).

We used the code PARODY-JA (Dormy 1997; Aubert 2019; Schaeffer 2013) initially developed to model Boussinesq convection and magnetic induction at finite Prandtl number in a rotating spherical shell. The code relies on a poloidal-toroidal decomposition to ensure a solenoidal velocity field and was successfully benchmarked (Christensen et al. 2001). The governing equations are solved using a finite-difference scheme parallelized with the MPI library along the radial direction, and a spectral method parallelized using OpenMP directives for the angular coordinates. A standard grid resolution for our simulations would be $n_r = 350$ points in the radial direction and spherical harmonics of degree $l$ and order $m$ of 200. $n_r$ is chosen as to ensure that at least 10 points are present in each thermal boundary layer, which is obtained via refining the grid close to the boundaries of the domain. We made several modifications to adapt PARODY-JA for the modelling of a silicate magma ocean: we removed the magnetic induction terms and equations, and we set a constant gravity profile with radius $r$, which is more appropriate for the Earth's mantle (Dziewonski & Anderson 1981). Our main addition to the code was the implementation of a routine to measure the efficiency of convective stirring, as described below.
## 2.1 Computational approach to measure the stirring efficiency

We estimated the efficiency of convective stirring using the Finite-Time Lyapunov Exponents (FTLEs), which provide a measure of the strain rate along Lagrangian particle trajectories. FTLEs characterize the deformation speed of a small, initially spherical, fluid parcel (McKenzie 1979). This method has already been used to study mixing properties and to derive mixing times in the mantle (Ferrachat & Ricard 2001; Farnetani & Samuel 2003; Coltice & Schmalzl 2006; Tackley 2015) and in geophysical flows.

The FTLE computation follows the approach described in Farnetani & Samuel (2003). It considers the deformation of a fluid parcel of size $\xi$ around a particle advected by the flow. Between two instants $t_0$ and $t_1$ the parcel will travel from $\mathbf{x}_0$ to $\mathbf{x}_1$ and will deform and rotate, leading to a change in $\xi$ from $\xi_0$ to $\xi_1$ as follows:

$$ \xi_1 = M(\mathbf{x}_0, t_0; \mathbf{x}_1, t_1)\xi_0, \qquad (5) $$

where $M$ is a linear operator that contains the information about the deformation history in the vicinity of the Lagrangian particle. $M$ is obtained by integrating the following equation, with the velocity gradients around a particle, along its trajectory:

$$ \frac{dM}{dt} = D_x u(\mathbf{x}(t), t)M. \qquad (6) $$

We use $M$ to compute the right Cauchy-Green strain tensor $M^T M$, whose eigenvalues $\sigma_i$ (with $i = 1, 2, 3$) represent deformation along the principal perpendicular axes and allow us to compute the FTLEs using:

$$ \lambda_i = \frac{1}{2} \left( \frac{1}{t_1 - t_0} \right) \ln \sigma_i . \qquad (7) $$

In 3-D space there are three eigenvalues corresponding to three different FTLE values: $\lambda_1 > \lambda_2 > \lambda_3$. The incompressibility constraint imposes that their sum is equal to zero at all times. The maximum (positive) FTLE, $\lambda_1$, characterizes the stretching along one principal direction, while the minimum (and negative) FTLE, $\lambda_3$, represents the maximum shrinking of the fluid parcel.

Once the simulations have reached a statistically steady-state stage, the system has forgotten its initial condition, and we randomly seed Lagrangian tracers in the spherical shell domain. We use passive tracers to track heterogeneities that have the same physical properties (e.g. density, viscosity) as the surrounding fluid. These passive tracers are therefore advected by the flow without affecting it. We then compute the FTLEs along the Lagrangian trajectories described by each tracer. This approach results in large sets of coupled ordinary differential equations (ODE), which are numerically integrated using a total variation diminishing (TVD) Runge-Kutta scheme (Harten 1984). The number of tracers is of the same order as the number of nodes used for the spatial resolution of the governing flow equations. For a standard case we typically use $n_r = 350$ points in the radial direction, spherical harmonics of degree $l$ and order $m$ of 200, leading to $n_r \times l \times m = 14 \times 10^6$ tracers. We tested the influence of the number of tracers on our results by considering up to eight times more tracers (e.g. eight tracers per grid node instead of one) and found that the difference in the final values of the Lyapunov exponents is very small (<0.1 per cent), thereby confirming the robustness of our results using one tracer per node.

## 2.2 Numerical benchmark of the FTLE computation

We benchmarked our implementation of the FTLE computation against an analytical solution of the Navier-Stokes equations, the ABC flow (Arnold 1965) in a Cartesian periodic domain $(x, y, z)$:

$$ \begin{aligned} u_x &= A \sin z + C \cos y, & u_y &= B \sin x + A \cos z, \\ u_z &= C \sin y + B \cos x. & & \end{aligned} \qquad (8) $$

Following Babaee et al. (2017) we set the value of the parameters to $A = \sqrt{3}$, $B = \sqrt{2}$ and $C = 1$. This generates chaotic Lagrangian trajectories (Dombre et al. 1986; Haller 2001). We compute the FTLE for one particle initially located at $x = 4$, $y = 0.6$ and $z = 0$ in the domain. Fig. 1(a) shows the trajectory of a single tracer within a periodic square domain of length $2\pi$ travelling along a streamline in the chaotic flow. We compute the FTLEs both analytically and using our routine with a second-order TVD Runge-Kutta integration scheme. Fig. 1(b) shows time evolution of the positive, intermediate, and negative FTLEs obtained using numerical and analytical solutions. Both solutions give similar results in good agreement with Babaee et al. (2017). We also verified that the sum of the three exponents remains close to zero (as expected for a solenoidal flow). This is shown in Fig. 1(c) where the sum of the FTLE is smaller than $7 \times 10^{-7}$ at all times. These results validate our implementation of the FTLE approach. In the following, we refer to $\lambda_1$ as $\lambda$ since we exclusively focus on the value of the maximum exponent to characterize the stirring efficiency.

The computation of the Lyapunov exponents requires the integration of two systems of coupled ODEs associated with each Lagrangian tracer (Samuel 2018). One corresponds to the advection of the tracer, and the other to the computation of the right Cauchy-Green deformation tensor (eq. 6).

Since we use passive Lagrangian tracers, the FTLE computation is not subject to a Courant-Friedrichs-Lewy stability criterion along each spatial direction: $u \Delta t < \Delta h$, with $u$ the corresponding velocity component, $\Delta t$ the time step and $\Delta h$ the grid spacing, unlike the time integration performed in PARODY-JA. Therefore, we can decouple the two systems and use a larger integration time step for the FTLEs. However, increasing the time step associated with the computation of the FTLEs and the advection of the tracers, $\Delta t_{FTLE}$, will degrade the precision of the numerical solution. We therefore performed multiple tests to determine the optimal value for $\Delta t_{FTLE}$ yielding both a good precision and a reduced computational time. Fig. 2(a) shows the total computational time associated with PARODY alone (circles) or with both PARODY and the FTLE computation (squares) as a function of $\Delta t_{FTLE}/\Delta t_{PARODY}$, where $\Delta t_{PARODY}$ is the time step used for solving the equations of motion. Note that the total computational time for PARODY alone depends on $\Delta t_{FTLE}/\Delta t_{PARODY}$ because the time associated with the output of the velocity field increases with smaller values of $\Delta t_{FTLE}/\Delta t_{PARODY}$ (i.e. more frequent output of the velocity fields). Using a second-order Runge-Kutta scheme with $\Delta t_{FTLE} = 50\Delta t_{PARODY}$ results in a negligible extra-cost associated with the FTLEs for a typical problem size composed of 350 radial points and spherical harmonics of degree and order $l = m = 200$.

We define $\lambda_{mean}$, the arithmetic spatial mean of the FTLEs averaged over 50 time steps. Fig. 2(b) shows the error on $\lambda_{mean}$ for first-, second- and third-order Runge-Kutta TVD integration schemes as a function of $\Delta t_{FTLE}/\Delta t_{PARODY}$. The error is computed with respect to the solution obtained using a third-order Runge-Kutta and $\Delta t_{FTLE} = \Delta t_{PARODY}$. Using a second-order Runge-Kutta method with $\Delta t_{FTLE} = 50\Delta t_{PARODY}$ results in an error of 1.9 per cent for the computation of $\lambda_{mean}$. Therefore, this larger time step for the FTLE computation
<figure><img src="image_3.png" /><figcaption>**Figure 1.** (a) Trajectory of a tracer in an ABC flow where $A = \sqrt{3}$, $B = \sqrt{2}$ and $C = 1$ inside a periodic square domain of length $2\pi$. The colour of the trajectory represents time, blue and red crosses indicate the initial and final position, respectively. (b) Analytical (lines) and numerical (squares) computation of the maximum $\lambda_1$ (blue), intermediate $\lambda_2$ (green) and minimum $\lambda_3$ (red) Lyapunov exponents for a tracer in the same ABC flow. The numerical solution is computed using a second-order TVD Runge-Kutta integration scheme. (c) Time evolution of the sum of the FTLEs.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 2.** (a) Total computation time as a function of $\Delta t_{FTLE}/\Delta t_{PARODY}$, the ratio of the time step associated with the FTLEs over the time step associated with the computation of the temperature and velocity fields. The green dots correspond to the computation of the flow alone (PARODY-JA), the red squares represent the computation time for PARODY-JA together with the FTLE computation. (b) Error associated with the estimate of the arithmetic average finite-time Lyapunov exponent for different Runge-Kutta (RK) TVD schemes as a function of $\Delta t_{FTLE}/\Delta t_{PARODY}$. The value used for the FTLE is the arithmetic mean of the FTLE over 50 time steps after the computation converged asymptotically. The error is based on simulations using a 3rd order Runge-Kutta and $\Delta t_{FTLE} = \Delta t_{PARODY}$.</figcaption></figure>

effectively removes the extra-cost at the expense of only a small error on $\lambda_{mean}$.

## 3 RESULTS

We performed simulations with $Ra = [10^3 - 10^9]$, $Ek = [10^{-5} - 10^{-1}]$ and a constant Prandtl number $Pr = 1$. By exploring a broad range of Rayleigh and Ekman numbers we identified different flow regimes, which are dominated either by the thermal buoyancy force or by the Coriolis force. Our focus is to establish a proper understanding of stirring processes in a non-linear (turbulent) regime of convection dominated by the buoyancy force with inertial effects. This regime is relevant to a terrestrial magma ocean and is the closest to the classical Rayleigh-Bénard convection.

### 3.1 Non-linear regimes of rotating convection

The dynamic regimes are determined by the balance between the buoyancy and the Coriolis forces. To describe the competition between the two forces, we rely on the convective Rossby number:

$$Ro_c = \frac{Ra^{1/2} Ek}{Pr^{1/2}}, \qquad (9)$$

defined in Gastine et al. (2016) as a proxy for the ratio of the buoyancy force to the Coriolis force. Low values of $Ro_c$ correspond to regimes dominated by the Coriolis force and high values to regimes dominated by the buoyancy force. For example, for $Ra = 10^8$ and within the range $Ek = [2.5 \times 10^{-5} - 10^{-2}]$, we identified flows with different structures and characterized the driving force of the flow using $Ro_c$. Fig. 3 shows the temperature field for three different flow regimes with varying rotational effects. In Fig. 3(a) we observe turbulent convective structures rising from the base of the domain towards the surface of the shell, with no visible impact of rotation. The flow appears to be dominated by the buoyancy force, with $Ro_c \gg 1$, and is compatible with the regime of non-rotating convection where rotational effects are weak. At lower Ekman number (Fig. 3b) and $Ro_c = 2.5$, both convective and rotational effects have an impact on the flow structure. In this case the convective upwellings start to be deviated by the zonal rotational flow around the equatorial axis, and columns parallel to the equatorial axis develop. This corresponds to a transitional regime, where both the rotational zonal flow and convective structures coexist. In this regime small changes in the balance between the buoyancy and the Coriolis force can lead to a flow dominated by either force. At even lower Ekman number (Fig. 3c), the Coriolis force dominates and rotational
<figure><img src="image_5.png" /><figcaption>**Figure 3.** Dimensionless temperature field represented in a quarter of the spherical shell at fixed Ra = 10⁸ and Pr = 1. (a) At Ek = 10⁻² the regime is dominated by the buoyancy force. (b) At Ek = 2.5 × 10⁻⁴ the regime is transitional. (c) At Ek = 2.5 × 10⁻⁵ the regime is rotation-dominated. For each case we also provide the convective Rossby number Ro<sub>c</sub>. The difference observed in flow regime from (a) to (c) is an effect of lowering the supercriticality R = Ra Ek¹/³ (Gastine et al. 2016) from R ~ 2 × 10⁵ to R ~ 73.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 4.** Dimensionless radial profiles of the horizontally averaged temperature for numerical experiments at Ra = 10⁸ and for various Ekman numbers.</figcaption></figure>

effects essentially control the flow. This corresponds to a rotation-dominated regime where a zonal flow strongly damps the convective structures.

Fig. 4 shows the radial temperature profile for the cases displayed in Fig. 3. For Ek = 10⁻², the profile is purely convective, with constant temperature in the bulk and well-developed thermal boundary layers, as expected when the flow is dominated by the buoyancy force. The relatively cold average bulk temperature observed occurs because of the spherical curvature and the resulting associated asymmetry at the inner and outer boundary layers due to the aspect ratio. This results in a larger temperature drop at the base of the domain (Gastine et al. 2015). For Ek = 10⁻⁴, the impact of rotational effects creates a temperature gradient in the bulk, but a non-negligible part of the temperature gradient still occurs in the thermal boundary layers, indicating the coexistence of rotational and convective effects of comparable magnitude (transitional regime). For Ek = 2.5 × 10⁻⁵, the temperature profile is essentially diffusive, such that most of the temperature variation occurs across the entire shell. This results from the dominant Coriolis force that induces a strong zonal flow.

Consequently, by using both Ro<sub>c</sub> (Fig. 3) and the temperature profiles (Fig. 4) we can discriminate among the different regimes of non-linear rotating convection.

## 3.2 Buoyancy-dominated regime

To better describe the convective dynamics of the buoyancy-dominated regime, we derived scaling laws for the Nusselt number and for the convective Reynolds number. The Nusselt number, that is the ratio of convective heat flux to the conductive heat flux taken at the outer radius, is defined as:

$$Nu = \frac{r_o}{r_i} \frac{Q L}{\rho C_p \kappa \Delta T}, \quad (10)$$

where Q is the convective heat flux per unit surface at the outer boundary, ρ the fluid density and C<sub>p</sub> the specific heat at constant pressure. The convective Reynolds number characterizes the root mean squared (RMS) velocity:

$$Re_c = \sqrt{2 \sum_{l=1}^{l_{max}} \sum_{m=1}^{l} e_l^m}, \quad (11)$$

where e is the dimensionless kinetic energy density (i.e. the kinetic energy per unit volume) and the overbar indicates arithmetic averaging over the entire spherical shell domain.

In this regime, the flow is mainly controlled by the intensity of convection and is similar to the classical Rayleigh-Bénard convection without rotation. Therefore, at constant Pr, the quantities associated with the flow scale only with the Rayleigh number. Fig. 5(a) shows the Nusselt number obtained for cases with different Ra values and for Ek = [10⁻³ – 10⁻¹] in the buoyancy-dominated regime, leading to the following scaling:

$$Nu = (0.081 \pm 0.006) Ra^{0.336 \pm 0.004}. \quad (12)$$

Gastine et al. (2016) inferred an asymptotic scaling Nu<sup>NR</sup> ∝ Ra<sup>1/3</sup>, with Nu<sup>NR</sup> the Nusselt number in the non-rotating regime of non-linear convection in a spherical shell. Analogue experiments by Cheng et al. (2015) and numerical simulations by
<figure><img src="image_7.png" /><figcaption>**Figure 5.** (a) Nusselt number as a function of the Rayleigh number for buoyancy-dominated cases. The black line $Nu = 0.081 Ra^{0.336}$ is the fit to the numerical simulations. (b) Convective Reynolds number as a function of the Rayleigh number for the cases displayed in (a) along with the corresponding fit $Re_c = 0.17 Ra^{0.497}$.</figcaption></figure>

Iyer et al. (2020) found, respectively, $Nu = 0.075 Ra^{0.322}$ and $Nu = 0.0525 \pm 0.006 Ra^{0.331 \pm 0.002}$. The Nusselt number scaling we obtained is therefore consistent with theoretical, numerical and analogue experimental results. The differences are essentially in the pre-factors and likely result from differences in boundary conditions, geometry and gravity profile. Following a similar approach for the convective Reynolds number, we compared our scaling :

$$Re_c = (0.17 \pm 0.03) Ra^{0.497 \pm 0.009}, \qquad (13)$$

(see Fig. 5b) with the asymptotic scaling $Re_c^{NR} \sim Ra^{1/2}$ by Gastine et al. (2016), and with $Re = (0.1555 \pm 0.006) Ra^{0.458 \pm 0.006}$ by Iyer et al. (2020). The good agreement between our results and previous studies indicates that in this regime, where convective motions account for approximately 90 per cent of the kinetic energy, the convective Reynolds number behaves in the same way as the classical Reynolds number. In summary, the quantities characterizing the flow in our numerical experiments behave in the same way as experimental, numerical and theoretical results from the literature for non-rotating convection.

## 3.3 Boundary of the buoyancy-dominated regime of convection

The scaling laws we derived above for $Nu$ and $Re_c$ are only valid within the buoyancy-dominated regime and cease to work as rotational effects significantly affect the flow. Simulations where the output parameters (i.e. $Nu$, $Re_c$) start deviating from the scaling laws displayed in Figs 5(a) and (b) mark the boundary of the buoyancy-dominated regime. Therefore, we considered that the buoyancy-dominated regime of convection corresponds to the region in the ($Ra$, $Nu$) space where the Nusselt number lies within $\pm 10$ per cent of the prediction obtained from eq. (12). Fig. 6 shows the Nusselt number for various cases associated with different values of $Ra$ and $Ek$. The grey area for simulations with $Nu < 1.4$ is the weakly non-linear regime (i.e. the onset of convection) and is irrelevant to the study of the magma ocean. The buoyancy-dominated regime as defined above is displayed in red, and simulations that lie within this region are considered for the characterization of the stirring efficiency. For cases where rotational effects become significant compared to convective motions, the Nusselt number deviates from the prediction of eq. (12). Hence, the lower the $Ek$, the higher the $Nu$ (or $Ra$) necessary to reach the buoyancy-dominated regime. For multiple ($Ra$, $Ek$) values, this corresponds to the criterion $Ra \ge 10 Ek^{-2}$, defined using the limit of the buoyancy-dominated regime described above and in Fig. 6.

<figure><img src="image_8.png" /><figcaption>**Figure 6.** Nusselt number as a function of the Rayleigh number for $Ek = [10^{-2} - 10^{-4}]$. The grey region represents the weakly non-linear regime (i.e. the onset of convection). The red region is the buoyancy-dominated regime, defined as the domain where the value of $Nu$ lies within $\pm 10$ per cent of the prediction from eq. (12) (as displayed by the black dotted line).</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 7.** Regime diagram for our numerical experiments in the ($Ra$, $Ek$) space with $Ro_c$ indicated via the colour scale. Within the non-linear regime of convection we defined a buoyancy-dominated regime (red area), a transitional regime and a rotation-dominated regime (white area). The weakly non-linear regime is indicated by a grey area, and cases that belong to this regime are not shown as they are irrelevant to the study of terrestrial a magma ocean.</figcaption></figure>

This criterion can be explained using the convective Rossby number, $Ro_c$, which is a proxy for the relative importance of the buoyancy force compared to the Coriolis force (Gastine et al. 2016). Using eq. (9), the convective Rossby number for $Pr = 1$ is $Ro_c = Ra^{1/2} Ek$, or $Ra = Ro_c^2 Ek^{-2}$. Hence, our criterion $Ra \ge 10 Ek^{-2}$ corresponds to $Ro_c \ge \sqrt{10} \sim 3.16$. $Ro_c \sim O(\infty)$ is characteristic of flows where the influence of both forces is comparable. Therefore, below this criterion the force balance shifts towards regimes where the dynamics is no longer dominated by convective motions and rotational effects become relevant.

Fig. 7 shows a regime diagram of our numerical experiments in the ($Ra$, $Ek$) space with $Ro_c$ indicated via the colour scale. The buoyancy-dominated regime is depicted in red and bounded by $Ra > 10 Ek^{-2}$ ($Ro_c \gtrsim 3$), and the weakly non-linear regime in grey is bounded by $Nu > 1.4$. In between both previous regimes, the white region delimits where non-linear convective dynamics are impacted by rotational effects.
<figure><img src="image_10.png" /><figcaption>Figure 8. Average positive finite-time Lyapunov exponent $\lambda_{\text{mean}}$ as a function of the Rayleigh number in the buoyancy-dominated regime. The black line corresponds to the fit $\lambda_{\text{mean}} = 0.0084 Ra^{0.734}$.</figcaption></figure>

## 3.4 Stirring efficiency in the buoyancy-dominated regime of convection

To characterize the stirring efficiency, we derived a scaling law for the FTLE distribution. For each case pertaining to the buoyancy-dominated regime, we found that the FTLE distribution is Gaussian and spatially homogeneous. Therefore, the average FTLE value for the entire domain, $\lambda_{\text{mean}}$, is a relevant proxy to characterize the stirring efficiency. Additionally, we found that the associated standard deviation of the FTLEs Gaussian distribution normalized by $\lambda_{\text{mean}}$ is roughly constant with a value of ~0.1. This allows one to predict the entire distribution of the FTLEs, rather than just its average value. We averaged $\lambda_{\text{mean}}$ over 50 time steps after the computation of the FTLEs has converged to an asymptotic distribution. This time averaging is performed to account for temporal fluctuations of the FTLEs. Fig. 8 displays $\lambda_{\text{mean}}$ as a function of Ra in the buoyancy-dominated regime. The observed dependence of $\lambda_{\text{mean}}$ with Ra can be fitted by the following expression:

$$ \lambda_{\text{mean}} = (0.0084 \pm 0.0002) Ra^{0.734 \pm 0.001}. \quad (14) $$

For the solid mantle, previous studies established that the Lagrangian strain rate is proportional to the Reynolds number or the RMS velocity (Coltice & Schmalzl 2006; Samuel et al. 2011). However, in our case, the Lagrangian strain rate scales as ~$Ra^{0.73}$ while the Reynolds number scales as ~$Ra^{0.5}$. A major difference between the solid mantle case and our simulations is the Prandtl number, considered infinite for the first, while $Pr = 1$ for the latter. Therefore, inertial effects could impact convective stirring inside the domain and lead to the observed differences.

To explain eq. (14) we computed the Eulerian strain rate in the model domain and analysed its spatial distribution. In practice, this corresponds to the computation of the FTLEs for one time step, without taking into account the entire trajectory of fluid parcels. By doing so we identified domains where the deformation is maximum, as shown in Fig. 9(a), an equatorial cut for a numerical experiment with $Ra = 10^8$, $Ek = 10^{-2}$. Regions where the Eulerian strain rate is particularly high are clearly visible (see also black ellipse) and they correspond to high convective heat transport (Fig. 9b). Moreover, maxima of Eulerian strain rate are found within a thin boundary layer at the base, and at the top, of the model domain, as shown in Fig. 10(a). To explain the occurrence and the location of these maxima, we focused on the horizontal velocity gradients and located the maximum horizontal velocity gradients at the beginning of a thin layer at the boundary in Fig. 10(b). We refer to these layers as ‘velocity gradient layers’. This velocity gradient (and associated FTLE increase) is a consequence of upwelling and

downwelling convective structures. To satisfy the conservation of mass that requires a velocity divergence equal to zero, the radial velocities associated with upwelling or downwelling currents are deviated laterally as they approach the boundaries of the domain, thereby creating pure shear. This is visible in Fig. 10(c) with roughly constant radial and horizontal velocity in the bulk, but closer to the boundaries the radial velocity decreases rapidly until it reaches 0, whereas the horizontal velocity increases simultaneously. Regions characterized by pure shear, also called hyperbolic regions (Gurnis 1986; Ottino 1989; Haller 2001), are known to yield the most efficient stirring efficiency because pure shear leads to an exponential rate of deformation. The presence of hyperbolic regions embedded in the velocity gradient layers explains the sharp increase in stirring efficiency in these areas, which control the global (averaged) efficiency depicted in Fig. 8. In our numerical experiments, this velocity gradient is reinforced by the free-slip boundary conditions. No-slip boundary conditions would create drag resistance at the boundaries, slowing the horizontal velocities and possibly mitigating the velocity gradient responsible for the intense mixing we observed in these regions.

With this knowledge, we can attempt to replicate our scaling law for $\lambda_{\text{mean}}$ (eq. 14) using dimensional analysis. Since the FTLEs are a measure of the Lagrangian strain rate their dimension is $t^{-1}$. We therefore make the reasonable assumption that the strain rate scales as the ratio of a convective velocity over a convective length scale. For the velocity scale, we use the convective Reynolds number, which is akin to the dimensionless convective RMS velocity and scales with ~$Ra^{1/2}$ (Fig. 5b and eq. 13). We use the dimensionless velocity gradient layer thickness $\delta^L$ or $\delta^U$ for the length scale, with the superscripts L and U indicating the lower and upper boundary, respectively. As discussed above, this choice for the length scale is motivated by the observation of regions with maximum stirring efficiency located in the velocity gradient layers (Fig. 10a). This leads to:

$$ \dot{\lambda} \sim \frac{Re_c}{\delta}. \quad (15) $$

To measure the thickness of the velocity gradient layers, we use the radial profile of the laterally-averaged horizontal velocity, limiting ourselves to cases in the buoyancy-dominated regime (e.g. cases 1–2, 4–13, 16–19, 24–28 in Table A1). The velocity gradient layer thicknesses, determined using the depth where the horizontal velocity gradient is maximal, scale as $\delta^U = (3.33 \pm 0.75)Ra^{-0.301 \pm 0.0181}$ and $\delta^L = (2.93 \pm 0.44)Ra^{-0.297 \pm 0.020}$, as seen in Fig. 11(a), where the uncertainties indicated correspond to one standard deviation. Using eq. (15) we obtain the two following scalings:

$$ \lambda \sim \frac{Re_c}{\delta^U} = 0.054 Ra^{0.798} \qquad \lambda \sim \frac{Re_c}{\delta^L} = 0.059 Ra^{0.794}. \quad (16) $$

Considering the uncertainty determined for both $\delta$ and $Re_c$, we obtain the 95 per cent confidence intervals (corresponding to two standard deviation uncertainty):

$$ \begin{aligned} \lambda &\sim \frac{Re_c}{\delta^U} = (0.054 \pm 0.031) Ra^{0.798 \pm 0.10} \\ \lambda &\sim \frac{Re_c}{\delta^L} = (0.059 \pm 0.027) Ra^{0.794 \pm 0.11}. \end{aligned} \quad (17) $$

Fig. 11(b) compares the above theoretical expressions, with the scaling derived from our numerical experiments (eq. 14, dashed black line). The dependency of the FTLEs on the Rayleigh number ($\lambda_{\text{mean}} \sim Ra^{0.73}$) lies within the confidence intervals for both velocity
<figure><img src="image_11.png" /><figcaption>**Figure 9.** (a) Equatorial cut of the Eulerian strain rate for Ra = $10^8$, Ek = $10^{-2}$ and Pr = 1. The black ellipse highlights a zone of high stirring efficiency. (b) Equatorial cut of the temperature field for the same numerical experiment. The black ellipse highlights a relevant convective structure.</figcaption></figure>

<figure><img src="image_12.png" /><figcaption>**Figure 10.** Results of a numerical experiment at Ra = $10^8$ and Ek = $10^{-2}$. (a) Dimensionless radial profile of the horizontally averaged Eulerian strain rate. (b) Dimensionless radial profile of the horizontally averaged horizontal RMS velocity gradient for the same numerical experiment. The thicknesses of the velocity gradient layers are indicated by the red lines. (c) Dimensionless radial profile of the horizontally averaged horizontal and radial RMS velocities.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Velocity gradient layer thicknesses as a function of the Rayleigh number in the buoyancy-dominated regime. The solid lines are fits $\delta^L = 2.93Ra^{-0.297}$ and $\delta^U = 3.33Ra^{-0.301}$. (b) Ratio of the convective Reynolds number over the upper and lower boundary layer thickness. The green and red areas represent the 95 per cent confidence intervals for both ratios. The black dashed line is the scaling for the average FTLE derived from the numerical experiments.</figcaption></figure>

gradient layers (eq. 17). Therefore, the scaling law derived from the numerical experiments is consistent with stirring being controlled by convection, with higher FTLE values located in the vicinity of active convective structures and maximum stirring efficiency

located in the velocity gradient layers. The above analysis provides an explanation for the differences between the scaling laws obtained here for the FTLEs and those previously derived for solid-state convection.
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Regime diagram in the (*Ra*, *Ek*) space. The black rectangle indicates the range explored in our numerical experiments which can be extrapolated to terrestrial magma ocean's conditions, that is the purple box (corresponding to *Ek* = 10⁻¹⁶ – 10⁻¹⁰ and *Ra* = 10²⁶ – 10³²).</figcaption></figure>

## 3.5 Application to a terrestrial magma ocean

To reach the conditions of a terrestrial magma ocean we need to extrapolate our results over ~15 orders of magnitude for the Rayleigh number and ~10 orders of magnitude for the Ekman number. Fig. 12 shows the regime diagram for *Ra* up to 10³⁴ and *Ek* as low as 10⁻¹⁸ and *Pr* = 1. The most plausible values of *Ra* and *Ek* for a magma ocean are indicated by a purple region, which includes *Ra* ~ 10²⁸ (Solomatov 2015) and *Ek* ~ 10⁻¹³ (Maas et al. 2021). The purple area clearly shows that the weakly non-linear regime is not relevant, whereas the buoyancy-dominated regime can correspond to (*Ra*, *Ek*) pairs appropriate for a terrestrial magma ocean.

Similar methods of extrapolation have been used to study planetary geodynamo in magneto-hydrodynamics numerical experiments (Christensen & Aubert 2006). Since we considered the buoyancy-dominated regime of non-linear convection that is relevant for a terrestrial magma ocean, we can extrapolate the scaling laws derived for the FTLEs. This enables us to deduce a mixing time for passive heterogeneities in a terrestrial magma ocean.

Vigorously convective flows with turbulent structures generate time-dependent and chaotic stirring, which is associated with regions where the deformation is dominated by pure shear (Khakhar et al. 1984; Haller 2001). Under these conditions, a passive, and initially spherical heterogeneity of size $d_0$ will stretch and shrink at an exponential rate (Gurnis 1986; Kellogg & Turcotte 1987). Considering the direction of maximum shrinking, the time evolution of the size of a passive heterogeneity, *d*, writes:

$$d(t) = d_0 e^{-\lambda t}. \qquad (18)$$

When the size of the shrinking heterogeneity becomes comparable to, or smaller than, the chemical diffusion length $\sqrt{\kappa_c t}$, where $\kappa_c$ is the chemical diffusivity, chemical diffusion homogenizes the heterogeneity with its surroundings, leading to complete mixing. Following Gurnis (1986) we combine eq. (18) with the diffusion length, and obtain a criterion for a heterogeneity to be mixed by diffusion:

$$\sqrt{\kappa_c t} \ge d_0 e^{-\lambda t}, \qquad (19)$$

which can be used to derive the mixing time $\tau$, corresponding to the moment where the minimum size of shrinking heterogeneity equals the diffusion length: $\sqrt{\kappa_c \tau} = d_0 e^{-\lambda \tau}$.

Initially spherical heterogeneities of 10 or 1000 km, evolving in our numerical experiment at *Ra* = 10⁸, *Ek* = 10⁻² and *Pr* = 1, would

<figure><img src="image_15.png" /><figcaption>**Figure 13.** Time evolution of the minimum length of heterogeneities in a terrestrial magma ocean. The green area represents the length below which a heterogeneity is completely homogenized by diffusion processes. The yellow areas represent the range of thicknesses for a passive heterogeneity, using a 95 per cent confidence interval corresponding to two standard deviations for the FTLE distribution.</figcaption></figure>

be mixed in ~890 Ma and ~2.9 Ga, respectively. This is considerably larger than the life time of a magma ocean that can range from hundreds of years to a few tens of million years in the presence of a thick atmosphere (Abe 1993; Lebrun et al. 2013; Solomatov 2015). However, in a thick terrestrial magma ocean the Rayleigh number is approximately 20 orders of magnitude higher than that of the solid-state mantle. This would result in a stirring efficiency ~15 orders of magnitude larger (eq. 14), thereby reducing the mixing time of heterogeneities by approximately the same amount.

In Fig. 13, we consider a magma ocean depth of 2000 km, *Ra* = 10²⁸, *Ek* = 10⁻¹¹, *Pr* = 1 and $\kappa_c$=10⁻¹⁰ m² s⁻¹ (Rubie et al. 2003) and show the evolution of the minimum size distribution for two heterogeneities of initial length 10 and 1000 km. As mentioned above the knowledge of $\lambda$ (from eq. 18) and the fact that the standard deviation for the FTLE distribution corresponds to ~0.1 $\lambda$ allows one to reconstruct the entire Gaussian distribution of the FTLEs, which is then used to predict the evolution of the minimum size distribution of passive mantle heterogeneities.

The smaller heterogeneity is mixed within only 6 × 10⁻⁶ to 10⁻⁵ yr, or approximately 4 min. Similarly, the larger heterogeneity is mixed in ~5 min. Consequently, passive heterogeneities of any size are rapidly homogenized by convective motions and cannot survive in a convecting terrestrial magma ocean in this regime. Our results provide a first approximation for the time required to homogenize passive heterogeneities in a vigorously convecting terrestrial magma ocean where rotation does not dominate the dynamics. This has strong implications for the interpretation of geochemical data, since we predict mixing times that are ~7–8 orders of magnitudes smaller than the life time of a terrestrial magma ocean in the fully liquid stage (i.e. ~10–1000 yr; Solomatov 2000; Lebrun et al. 2013).

# 4 DISCUSSION

Modelling explicitly the extremely high Rayleigh number [O (10²⁸)] and low Ekman number [O (10⁻¹⁵)] of a planetary magma ocean is currently out of reach. However, fluid dynamics simulations can attain the flow regime relevant to a magma ocean even though the pair of (*Ra*, *Ek*) do not have the aforementioned values. In particular,
once the system is in the buoyancy-dominated regime, the depen-
dence of the results on *Ek* vanishes such that this parameter can be
discarded, and our study has characterized the main dependence of
mixing properties on *Ra*.

The mixing time derived from our simulations implies passive
heterogeneities had to form late enough in the crystallization se-
quence of the magma ocean to ensure that they only underwent
solid-state convection. Moreover, no further magma ocean stage
had to occur after the formation of the heterogeneities. If the en-
ergy from a subsequent impact is large enough to remelt entirely
the mantle, heterogeneities would evolve through a purely liquid
convection phase, causing their rapid mixing.

In the light of Hf–W systematics, our findings imply that the last
major magma ocean stage occurred while $^{182}$Hf was extant. The
possibility that heterogeneous domains could form during crystal-
lization of the magma ocean is supported by molecular dynam-
ics simulations (Deng & Stixrude 2021) showing that Hf is more
compatible than W. There is still debate on how proceeded the
solidification of the magma ocean, namely ‘bottom up’ (Andrault
et al. 2011) or from mid-depth (Stixrude et al. 2009; Boukaré et al.
2015), however the presence of a basal magma ocean (Labrosse et al.
2007) and its crystallization would lead to heterogeneous tungsten
anomalies (Deng & Stixrude 2021). Alternatively, core formation
has been invoked to explain Hf–W observations. Hf is a lithophile
element, while W is moderately siderophile, therefore the metal-
lic core (Hf-poor, W-rich) is a reservoir with a strongly negative
$\mu^{182}$W (Touboul et al. 2012). A key issue is then how to transfer
a ‘core signature’ to the lowermost mantle. High pressure experi-
ments (Yoshino et al. 2020) suggest that it is possible to transfer W
from the core to the mantle because grain boundary diffusion of W
is a fast and strongly temperature-dependent process. This supports
the idea that a fraction of a lowermost mantle might be thermally and
chemically equilibrated with the core. Following Mundl-Petermeier
et al. (2020), the existence of a core–mantle equilibrated reservoir
could explain the $\mu^{182}$W anomalies found in modern OIB, since an
entrainment of only 0.3 per cent would be sufficient to create the
observed W anomalies in hotspot lavas.

Our models did not account for the presence of active hetero-
geneities (i.e. heterogeneities associated with different material
properties than the surrounding fluid) that could affect the flow
dynamics. Indeed, viscosity or density differences associated with
heterogeneities may alter the flow regime, resulting in different
mixing times. However, we expect that these effects will not appre-
ciably change our conclusions, because our mixing time estimates
are on the order of a few minutes, which is negligible with respect
to the lifetime of a fully liquid magma ocean (i.e. at least 1000 yr).
Therefore, active heterogeneities would have to affect the stirring ef-
ficiency by many orders of magnitude to invalidate our conclusions.
Another limitation is the fixed Prandtl number of 1 used in our study.
This value of *Pr* is plausible for terrestrial magma oceans and it sim-
plifies the comparison with other 3-D spherical convection studies
(e.g. Gastine *et al.* 2016). The large plausible range of magma ocean
viscosities, can result in either lower and larger *Pr* values. However,
this is out of the scope of this current work, and characterizing the
effect of the Prandtl number on our results will be deferred to a
forthcoming study. Lastly, the estimated mixing times are associ-
ated with a buoyancy-dominated regime of non-linear convection,
which represents only a part of the regimes relevant to a terrestrial
magma ocean, as seen in Fig. 12. Magma ocean's conditions could
also correspond to regimes where the Coriolis force strikes a bal-
ance with the buoyancy force or dominates over it, as depicted in
Figs 3(b) and (c), respectively. The effect of the Coriolis force is

to change the structure of the flow, because the rotating circulation
deviates the convective structures and creates convective columns
(Wang *et al.* 2021), possibly affecting the mixing efficiency. This
will be the focus of a future study.

5 CONCLUSIONS

We conducted numerical simulations, in a rotating spherical shell with free-slip boundary conditions, to characterize convective mixing for turbulent convection, and applied it to a terrestrial magma ocean. Our numerical experiments covered Rayleigh numbers from $10^3$ to $10^9$ and Ekman numbers from $10^{-5}$ to $10^{-1}$ with a fixed Prandtl number of 1, focusing on the buoyancy-dominated regime of a fully molten magma ocean with important inertial effects.

We observed that this regime is similar to classical Rayleigh-Bénard convection. The characteristics of the flow follow the same scaling as the one derived from previous theoretical, experimental and numerical studies, that is $Nu \sim Ra^{1/3}$ and $Re_c \sim Ra^{1/2}$. The regime is bounded by $Ra \ge 10 Ek^{-2}$ (Fig. 7) and beyond this limit, the influence of rotational effects becomes important.

For a buoyancy-dominated magma ocean, stirring is controlled
by the intensity of convection. We observed that stirring efficiency
is contained within velocity gradient layers located close to the inner
and outer boundaries, where the transition from radial to horizontal
velocities creates hyperbolic regions and large velocity gradients
that promote deformation. This effect is enhanced by the free-slip
boundary conditions we considered.

Away from the domain boundaries, the stirring efficiency is spa-
tially homogeneous and follows a Gaussian distribution. The latter
is well described by a power law $\lambda_{\text{mean}} = 0.008Ra^{0.734}$ with a constant
standard deviation of $\sim 0.1 \lambda_{\text{mean}}$ observed for all our experiments.
This relationship between the convective vigour and the stirring
efficiency is due to the fact that the average Lagrangian strain rate
$\lambda_{\text{mean}}$ scales as the ratio of convective velocities to the thickness of
the velocity gradient layers, where the deformation is maximum.

Our results indicate that for a liquid terrestrial magma ocean in
the buoyancy-dominated regime, stirring is homogeneous and ex-
tremely efficient. In this context, passive heterogeneities of 1000 km,
a length comparable to the thickness of a global magma ocean, are
mixed within ~5 min, which is considerably shorter than the life-
time of a terrestrial magma ocean in the fully molten stage (i.e. at
most ~1000 yr). Hence, our model favours scenarios where geo-
chemical heterogeneities do not survive global liquid magma ocean
stages, but may instead result from late accretion processes or from
core-mantle interactions.

ACKNOWLEDGMENTS

We thank N. Tosi and an anonymous reviewer for their constructive
comments, and G. Choblet for his editorial handling. We acknowl-
edge the funding the European Research Council (ERC SHRED,
Grant Agreement no. 833632: Survival of Hadean REmnants in a
Dynamic mantle) awarded to Catherine Chauvel. Numerical com-
putations were partly performed on the S-CAPAD/DANTE plat-
form, IPGP, France and the GENCI HPC resources of IDRIS under
allocations A0110412958 and A0130412958.

DATA AVAILABILITY

Data and code will be made available upon reasonable request to
the corresponding author.
REFERENCES

Abe, Y., 1993. Thermal evolution and chemical differentiation of the ter-
restrial magma ocean, in *Evolution of the Earth and Planets*, pp. 41–
54, eds Takahashi, E., Jeanloz, R. & Rubie, D., American Geophysical
Union.

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

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.

Arnold, V.I., 1965. Sur la topologie des ecoulements stationnaires des fluides parfaits, *C. R. Acad. Sci. Paris*, **261**, 17–20.

Aubert, J., 2019. Approaching Earth's core conditions in high-resolution geodynamo simulations, *Geophys. J. Int.*, **219**(Suppl 1), S137-S151.

Babaee, H., Farazmand, M., Haller, G. & Sapsis, T.P., 2017. Reduced-order description of transient instabilities and computation of finite-time Lyapunov exponents, *Chaos*, **27**(6), doi:10.1063/1.4984627.

Ballmer, M.D., Houser, C., Hernlund, J.W., Wentzcovitch, R.M. & Hirose, K., 2017. Persistence of strong silica-enriched domains in the Earth's lower mantle, *Nat. Geosci.*, **10**, 236–241.

Becker, T.W., Kellogg, J.B. & O'Connell, R.J., 1999. Thermal constraints on the survival of primitive blobs in the lower mantle, *Earth planet. Sci. Lett.*, **171**, 351–365.

Boukaré, 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.

Canup, R.M., 2008. Accretion of the Earth, *Phil. Trans. R. Soc. A*, **336**, 4061-4075.

Carlson, R.W. et al., 2014. How did Early earth become our modern world?,
*Annu. Rev. Earth planet. Sci.*, **42**, 151–178.

Cheng, J.S., Stellmach, S., Ribeiro, A., Grannan, A., King, E.M. & Aurnou, J.M., 2015. Laboratory-numerical models of rapidly rotating convection in planetary cores, *Geophys. J. Int.*, **201**(1), 1–17.

Christensen, U.R., 1989. Mixing by time-dependent convection, *Earth planet. Sci. Lett.*, **95**, 382–394.

Christensen, U.R. & Aubert, J., 2006. Scaling properties of convection-
driven dynamos in rotating spherical shells and application to planetary
magnetic fields, *Geophys. J. Int.*, **166**(1), 97–114.

Christensen, U.R. et al., 2001. A numerical dynamo benchmark, *Phys. Earth planet. Inter.*, **128**(1-4), 25-34.

Coltice, N. & Schmalzl, J., 2006. Mixing times in the mantle of the early Earth derived from 2-D and 3-D numerical simulations of convection, *Geophys. Res. Lett.*, **33**(23), doi:10.1029/2006GL027707.

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

Deguen, R., Olson, P. & Cardin, P., 2011. Experiments on turbulent metal-
silicate mixing in a magma ocean, *Earth planet. Sci. Lett.*, **310**(3), 303–
313.

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

Deng, J. & Stixrude, L., 2021. Deep fractionation of Hf in a solidifying magma ocean and its implications for tungsten isotopic heterogeneities in the mantle, *Earth planet. Sci. Lett.*, **562**, doi:10.1016/j.epsl.2021.116873.

Dombre, T., Frisch, U., Greene, J.M., Hénon, M., Mehr, A. & Soward, A.M., 1986. Chaotic streamlines in the ABC flows, *J. Fluid Mech.*, **167**, 353-391.

Dormy, E., 1997. *Modelisation numerique de la dynamo terrestre*, PhD thesis, Institut de Physique du Globe de Paris, doi:1997GLOB0003.

Dziewonski, A.M. & Anderson, D.L., 1981. Preliminary reference Earth model, *Phys. Earth planet. Inter.*, **25**(4), 297-356.

Elkins-Tanton, L.T., 2012. Magma oceans in the inner solar system, *Annu. Rev. Earth planet. Sci.*, **40**(1), 113–139.

Farnetani, C.G. & Samuel, H., 2003. Lagrangian structures and stirring in the Earth's mantle, *Earth planet. Sci. Lett.*, **206**(3), 335-348.

Ferrachat, S. & Ricard, Y., 1998. Regular vs. chaotic mantle mixing, *Earth planet. Sci. Lett.*, **155**(1), 75-86.

Ferrachat, S. & Ricard, Y., 2001. Mixing properties in the Earth's mantle: effects of the viscosity stratification and of oceanic crust segregation, *Geochem. Geophys. Geosyst.*, **2**, doi:10.1029/2000GC000092.

Gastine, T., Wicht, J. & Aurnou, J.M., 2015. Turbulent Rayleigh-Bénard convection in spherical shells, *J. Fluid Mech.*, **778**, 721-764.

Gastine, T., Wicht, J. & Aubert, J., 2016. Scaling regimes in spherical shell rotating convection, *J. Fluid Mech.*, **808**, 690-732.

Grossmann, S. & Lohse, D., 2000. Scaling in thermal convection: a unifying theory, *J. Fluid Mech.*, **407**, 27-56.

Gülcher, A.J.P., Gebhardt, D.J., Ballmer, M.D. & Tackley, P.J., 2020. Variable dynamic styles of primordial heterogeneity preservation in the Earth's lower mantle, *Earth planet. Sci. Lett.*, **536**, doi:10.1016/j.epsl.2020.116160.

Gurnis, M., 1986. Stirring and mixing in the mantle by plate-scale flow: large persistent blobs and long tendrils coexist, *Geophys. Res. Lett.*, **13**, 1474-1477.

Haller, G., 2001. Distinguished material surfaces and coherent structures in three-dimensional fluid flows, *Physica D*, **149**(4), 248-277.

Harten, A., 1984. On a class of high resolution total-variation-stable finite-
difference schemes, *SIAM J. Numer. Anal.*, **21**(1), 1–23.

Hoffman, N.R.A. & McKenzie, D.P., 1985. The destruction of geochemical heterogeneities by differential fluid motions during mantle convection, *82*, 163–206.

Hofmann, A.W. & White, W.M., 1982. Mantle plumes from ancient oceanic-
crust, *Earth planet. Sci. Lett.*, **57**, 421–436.

Hoïnk, T., Schmalzl, J. & Hansen, U., 2006. Dynamics of metal-silicate separation in a terrestrial magma ocean, *Geochem. Geophys. Geosyst.*, **7**, doi:10.1029/2006GC001268.

Horan, M., Carlson, R.W., Walker, R.J., Jackson, M., Garcon, M. & Norman, M., 2018. Tracking Hadean processes in modern basalts with 142-Neodymium, *Earth planet. Sci. Lett.*, **484**, 184-191.

Iyer, K.P., Scheel, J.D., Schumacher, J. & Sreenivasan, K.R., 2020. Classical 1/3 scaling of convection holds up to Ra = 10¹⁵, *Proc. Natl. Acad. Sci. U.S.A.* **117**(14), 7594-7598.

Jackson, M.G. et al., 2020. Ancient helium and tungsten isotopic signatures preserved in mantle domains least modified by crustal recycling, *Proc. Natl. Acad. Sci. U.S.A.*, **117**(49), 30 993–31 001.

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

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

Kellogg, L.H. & Turcotte, D.L., 1987. Homogenization of the mantle by convective mixing and diffusion, *Earth planet. Sci. Lett.*, **81**(4), 371-378.

Khakhar, D., Chella, R. & Ottino, J., 1984. *Stretching, Chaotic Motion, and Breakup of Elongated Droplets in Time Dependent Flows*, Mexico Univ. Nat. Autonom.

Labrosse, S., Hernlund, J.W. & Coltice, N., 2007. A crystallizing dense magma ocean at the base of the Earth's mantle, *Nature*, **450**(7171), 866-869.

Landeau, M., Deguen, R., Phillips, D., Neufeld, J.A., Lherm, V. & Dalziel, S.B., 2021. Metal-silicate mixing by large Earth-forming impacts, *Earth planet. Sci. Lett.*, **564**, doi:10.1016/j.epsl.2021.116888.

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.*, **118**(6), 1155-1176.

Lichtenberg, T., Bower, D.J., Hammond, M., Boukrouche, R., Sanan, P., Tsai, S. & Pierrehumbert, R.T., 2021. Vertically resolved magma ocean-protoatmosphere evolution: H₂, H₂O, CO₂, CH₄, CO, O₂, and N₂ as primary absorbers, *J. geophys. Res.*, **126**(2), 1-25.

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

Maas, C. & Hansen, U., 2015. Effects of Earth's rotation on the early differ-
entiation of a terrestrial magma ocean, *J. geophys. Res.*, **120**, 7508–7525.
Maas, C. & Hansen, U., 2019. Dynamics of a terrestrial magma ocean under planetary rotation: a study in spherical geometry, *Earth planet. Sci. Lett.*, **513**, 81–94.

Maas, C., Manske, L., Wünnemann, K. & Hansen, U., 2021. On the fate of impact-delivered metal in a terrestrial magma ocean, *Earth planet. Sci. Lett.*, **554**, doi:10.1016/j.epsl.2020.116680.

Manga, M., 1996. Mixing of heterogeneities in the mantle: effect of viscosity differences, *Geophys. Res. Lett.* **23**(4), 403–406.

Maurice, M., Tosi, N., Samuel, H., Plesa, A.-C., Hüttig, C. & Breuer, D., 2017. Onset of solid-state mantle convection and mixing during magma ocean solidification, *J. geophys. Res.*, **122**(3), 577–598.

McKenzie, D., 1979. Finite deformation during fluid flow, *Geophys. J. R. astr. Soc.*, **58**(3), 689–715.

Mitrovica, J. & Forte, A., 2004. A new inference of mantle viscosity based upon joint inversion of convection and glacial isostatic adjustment data, *Earth planet. Sci. Lett.*, **225**(1–2), 177–189.

Moeller, A. & Hansen, U., 2013. Influence of rotation on the metal rain in a Hadean magma ocean, *Geochem., Geophys., Geosyst.*, **14**(4), 1226–1244.

Monteux, J., Ricard, Y., Coltice, N., Dubuffet, F. & Ulvrovà, M., 2009. A model of metal-silicate separation on growing planets, **287**, 353–362.

Mundl, A., Touboul, M., Jackson, M.G., Day, J.M.D., Kurz, M.D., Lekic, V., Helz, R.T. & Walker, R.J., 2017. Tungsten-182 heterogeneity in modern ocean island basalts, *Science*, **356**(6333), 66–69.

Mundl-Petermeier, A., Walker, R.J., Jackson, M.G., Blichert-Toft, J., Kurz, M.D. & Halldórrsson, S.A., 2019. Temporal evolution of primordial tungsten-182 and 3He/4He signatures in the Iceland mantle plume, *Chem. Geol.*, **525**, 245–259.

Mundl-Petermeier, A., Walker, R.J., Fischer, R.A., Lekic, V., Jackson, M.G. & Kurz, M.D., 2020. Anomalous ¹²W in high 3He/4He ocean island basalts: fingerprints of Earth's core?, *Geochim. Cosmochim. Acta*, **271**, 194–211.

Nikolaou, A., Katyal, N., Tosi, N., Godolt, M., Grenfell, J.L. & Rauer, H., 2019. What factors affect the duration and outgassing of the terrestrial magma ocean?, *Astrophys. J.*, **875**(1), doi:10.3847/1538-4357/ab08ed.

Olson, P., Yuen, D.A. & Balsiger, D., 1984. Mixing of passive heterogeneities by mantle convection, *J. geophys. Res.*, **89**(B1), 425–436.

Ottino, J.M., 1989. *The Kinematics of Mixing: Stretching, Chaos and Transport*, Cambridge Univ. Press.

Peters, B.J., Mundl-Petermeier, A., Carlson, R.W., Walker, R.J. & Day, J. M.D., 2021. Combined lithophile-siderophile isotopic constraints on Hadean processes preserved in ocean island basalt sources, *Geochem. Geophys. Geosyst.*, **22**(3), doi:10.1029/2020GC009479.

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.

Rizo, H. et al., 2019. ¹²W evidence for core-mantle interaction in the source of mantle plumes, *Geochem. Perspect. Lett.*, **11**, 6–11.

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

Rubie, D.C., Nimmo, F. & Melosh, H.J., 2015. Formation or Earth's core, in *Treatise on Geophysics*, 2nd edn, Vol. 9, pp. 43–74, eds Schubert, G. & Stevenson, D.J., Elsevier.

Salvador, A., Massol, H., Davaille, A., Marcq, E., Sarda, P. & Chassefière, E., 2017. The relative influence of H₂O and CO₂ on the primitive surface conditions and evolution of rocky planets, *J. geophys. Res.*, **122**(7), 1458–1486.

Samuel, H., 2018. A deformable particle-in-cell method for advective transport in geodynamic modelling, *Geophys. J. Int.*, **214**, 1744–1773.

Samuel, H. & King, S.D., 2014. Mixing at mid-ocean ridges controlled by small-scale convection and plate motion, *Nat. Geosci.*, **7**, 603–606.

Samuel, H. & Tosi, N., 2012. The influence of post-perovskite strength on the Earth's mantle thermal and chemical evolution, *Earth planet. Sci. Lett.*, **323-324**, 50–59.

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.

Samuel, H., Aleksandrov, V. & Deo, B., 2011. The effect of contaminants on mantle convective stirring, *Geophys. Res. Lett.*, **38**(4), doi:10.1029/2010GL046056.

Schaeffer, N., 2013. Efficient spherical harmonic transforms aimed at pseudospectral numerical simulations, *Geochem. Geophys. Geosyst.*, **14**(3), 751–758.

Senshu, H., Kuramoto, K. & Matsui, T., 2002. Thermal evolution of a growing Mars, *J. geophys. Res.*, **107**, doi:10.1029/2001JE001819.

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

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

Solomatov, V.S. & Stevenson, D.J., 1993. Nonfractional crystallization of a terrestrial magma ocean, *J. geophys. Res.*, **98**(E3), 5391–5406.

Stevenson, D.J., 1987. Origin of the Moon - the collision hypothesis, *Annu. Rev. Earth planet. Sci.*, **15**, 271–315.

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.

Tackley, P.J., 2015. Mantle geochemical geodynamics, in *Treatise on Geophysics*, 2nd edn, pp. 521–585, ed. Schubert, G., Elsevier.

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

van Keken, P.E., Hauri, E.H. & Ballentine, C.J., 2002. Mantle mixing: the generation, preservation, and destruction of chemical heterogeneity, *Annu. Rev. Earth planet. Sci.*, **30**, 493–525.

Wachele, J.-B., Le Bars, M., Monteux, J. & Aurnou, J.M., 2014. Laboratory experiments on the breakup of liquid metal diapirs, **403**, 236–245.

Wang, Q., Chong, K.L., Stevens, R.J.A.M., Verzicco, R. & Lohse, D., 2020. From zonal flow to convection rolls in Rayleigh-BéNard convection with free-slip plates, *J. Fluid Mech.*, **905**, doi:10.1017/jfm.2020.793.

Wang, G., Santelli, L., Lohse, D., Verzicco, R. & Stevens, R.J.A.M., 2021. Diffusion-free scaling in rotating spherical Rayleigh-BéNard convection, *Geophys. Res. Lett.*, **48**(20), e2021GL095017.

Yoshino, T., Makino, Y., Suzuki, T. & Hirata, T., 2020. Grain boundary diffusion of W in lower mantle phase with implications for isotopic heterogeneity in oceanic island basalts by core-mantle interactions, *Earth planet. Sci. Lett.*, **530**, doi:10.1016/j.epsl.2019.115887.
APPENDIX: SIMULATIONS

<table><caption>**Table A1.** Quantities of interest for the numerical experiments conducted in this study. Measurements of the boundary layer thickness are limited to cases [1–2, 9–13, 17–19, 24–28]. Other cases were either not relevant (i.e. they did not pertain to the buoyancy-dominated regime) or lacked sufficient points at the same Ra values to allow for robust scalings. The computation of $\lambda_{mean}$ was not performed for several cases in the buoyantly dominated regime ([10, 14–15, 21–23, 32–33]) because of the associated prohibitive computational time. Additionally, $\lambda_{mean}$ was not computed for most cases outside of the buoyancy-dominated regime [29–31, 39–50].</caption><thead><tr><th>Case number</th><th>Ra</th><th>Ek</th><th>Nu</th><th>Re<sub>c</sub></th><th>δ<sup>L</sup></th><th>δ<sup>U</sup></th><th>V<sub>rms</sub></th><th>λ<sub>mean</sub></th><th>nr × l</th></tr></thead><tbody><tr><td>1</td><td>10<sup>7</sup></td><td>10<sup>-1</sup></td><td>17.91</td><td>587.95</td><td>0.0293</td><td>0.0256</td><td>623.1</td><td>1195.5</td><td>300 × 88</td></tr><tr><td>2</td><td>10<sup>8</sup></td><td>10<sup>-1</sup></td><td>40.03</td><td>1793.8</td><td>0.0123</td><td>0.0117</td><td>1874.5</td><td>6293.0</td><td>320 × 133</td></tr><tr><td>3</td><td>5 × 10<sup>4</sup></td><td>10<sup>-2</sup></td><td>2.79</td><td>41.16</td><td>—</td><td>—</td><td>41.44</td><td>24.3</td><td>200 × 66</td></tr><tr><td>4</td><td>10<sup>5</sup></td><td>10<sup>-2</sup></td><td>3.84</td><td>61.09</td><td>0.1073</td><td>0.1031</td><td>60.54</td><td>42.2</td><td>200 × 66</td></tr><tr><td>5</td><td>1.8 × 10<sup>5</sup></td><td>10<sup>-2</sup></td><td>4.94</td><td>85.67</td><td>0.0842</td><td>0.0832</td><td>82.76</td><td>66.1</td><td>200 × 66</td></tr><tr><td>6</td><td>2.3 × 10<sup>5</sup></td><td>10<sup>-2</sup></td><td>5.44</td><td>95.09</td><td>0.0839</td><td>0.0819</td><td>95.38</td><td>88.9</td><td>200 × 66</td></tr><tr><td>7</td><td>2.8 × 10<sup>5</sup></td><td>10<sup>-2</sup></td><td>5.80</td><td>105.36</td><td>0.0712</td><td>0.0674</td><td>103.15</td><td>88.4</td><td>200 × 66</td></tr><tr><td>8</td><td>4 × 10<sup>5</sup></td><td>10<sup>-2</sup></td><td>6.48</td><td>125.73</td><td>0.0691</td><td>0.0637</td><td>122.18</td><td>112.1</td><td>200 × 66</td></tr><tr><td>9</td><td>10<sup>6</sup></td><td>10<sup>-2</sup></td><td>8.64</td><td>201.10</td><td>0.0508</td><td>0.0498</td><td>202.03</td><td>207.4</td><td>200 × 88</td></tr><tr><td>10</td><td>5 × 10<sup>6</sup></td><td>10<sup>-2</sup></td><td>12.71</td><td>436.53</td><td>0.0346</td><td>0.0315</td><td>441.30</td><td>—</td><td>220 × 88</td></tr><tr><td>11</td><td>10<sup>7</sup></td><td>10<sup>-2</sup></td><td>17.68</td><td>606.63</td><td>0.0250</td><td>0.0237</td><td>613.17</td><td>1161.5</td><td>220 × 88</td></tr><tr><td>12</td><td>5 × 10<sup>7</sup></td><td>10<sup>-2</sup></td><td>30.11</td><td>1285.38</td><td>0.0182</td><td>0.0142</td><td>1267.04</td><td>—</td><td>300 × 133</td></tr><tr><td>13</td><td>10<sup>8</sup></td><td>10<sup>-2</sup></td><td>40.00</td><td>1788.85</td><td>0.0119</td><td>0.0116</td><td>1794.10</td><td>6541.3</td><td>350 × 133</td></tr><tr><td>14</td><td>5 × 10<sup>4</sup></td><td>5 × 10<sup>-3</sup></td><td>1.7</td><td>22.67</td><td>—</td><td>—</td><td>22.41</td><td>—</td><td>200 × 66</td></tr><tr><td>15</td><td>10<sup>5</sup></td><td>5 × 10<sup>-3</sup></td><td>2.8</td><td>47.54</td><td>—</td><td>—</td><td>48.09</td><td>—</td><td>200 × 66</td></tr><tr><td>16</td><td>10<sup>6</sup></td><td>5 × 10<sup>-3</sup></td><td>8.3</td><td>195.45</td><td>0.0471</td><td>0.0510</td><td>193.22</td><td>193.0</td><td>200 × 88</td></tr><tr><td>17</td><td>10<sup>7</sup></td><td>5 × 10<sup>-3</sup></td><td>17.30</td><td>591.61</td><td>0.0295</td><td>0.0241</td><td>601.23</td><td>1124.8</td><td>220 × 100</td></tr><tr><td>18</td><td>10<sup>8</sup></td><td>5 × 10<sup>-3</sup></td><td>40.35</td><td>1822.09</td><td>0.0118</td><td>0.0109</td><td>1835.67</td><td>6314.8</td><td>350 × 133</td></tr><tr><td>19</td><td>10<sup>8</sup></td><td>2.5 × 10<sup>-3</sup></td><td>39.80</td><td>1754.99</td><td>0.0131</td><td>0.0098</td><td>1741.04</td><td>6307.1</td><td>350 × 133</td></tr><tr><td>20</td><td>10<sup>6</sup></td><td>10<sup>-3</sup></td><td>4.02</td><td>145.60</td><td>—</td><td>—</td><td>145.60</td><td>125.8</td><td>220 × 88</td></tr><tr><td>21</td><td>2.5 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>—</td><td>101.33</td><td>—</td><td>—</td><td>104.28</td><td>—</td><td>220 × 88</td></tr><tr><td>22</td><td>5 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>11.00</td><td>422.60</td><td>—</td><td>—</td><td>421.75</td><td>—</td><td>220 × 88</td></tr><tr><td>23</td><td>7.5 × 10<sup>6</sup></td><td>10<sup>-3</sup></td><td>—</td><td>498.09</td><td>—</td><td>—</td><td>499.78</td><td>—</td><td>220 × 88</td></tr><tr><td>24</td><td>10<sup>7</sup></td><td>10<sup>-3</sup></td><td>15.70</td><td>600.04</td><td>0.0231</td><td>0.0260</td><td>600.19</td><td>1096.8</td><td>220 × 100</td></tr><tr><td>25</td><td>2.5 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>20.70</td><td>796.24</td><td>0.0192</td><td>0.0161</td><td>796.24</td><td>2203.1</td><td>300 × 133</td></tr><tr><td>26</td><td>5 × 10<sup>7</sup></td><td>10<sup>-3</sup></td><td>29.03</td><td>1170.57</td><td>0.0151</td><td>0.0165</td><td>1170.57</td><td>3545.7</td><td>300 × 133</td></tr><tr><td>27</td><td>10<sup>8</sup></td><td>10<sup>-3</sup></td><td>39.15</td><td>1666.13</td><td>0.0121</td><td>0.0129</td><td>1,666.13</td><td>6234.4</td><td>350 × 133</td></tr><tr><td>28</td><td>5 × 10<sup>8</sup></td><td>10<sup>-3</sup></td><td>72.50</td><td>3800.00</td><td>0.0075</td><td>0.0070</td><td>3800.00</td><td>20422.9</td><td>350 × 133</td></tr><tr><td>29</td><td>10<sup>6</sup></td><td>7.5 × 10<sup>-4</sup></td><td>2.90</td><td>129.54</td><td>—</td><td>—</td><td>128.19</td><td>—</td><td>220 × 88</td></tr><tr><td>30</td><td>5 × 10<sup>6</sup></td><td>7.5 × 10<sup>-4</sup></td><td>9.10</td><td>529.15</td><td>—</td><td>—</td><td>531.09</td><td>—</td><td>220 × 100</td></tr><tr><td>31</td><td>10<sup>7</sup></td><td>7.5 × 10<sup>-4</sup></td><td>13.60</td><td>744.31</td><td>—</td><td>—</td><td>728.92</td><td>—</td><td>300 × 133</td></tr><tr><td>32</td><td>5 × 10<sup>7</sup></td><td>7.5 × 10<sup>-4</sup></td><td>26.00</td><td>1140.18</td><td>—</td><td>—</td><td>1137.9</td><td>—</td><td>350 × 133</td></tr><tr><td>33</td><td>10<sup>8</sup></td><td>5 × 10<sup>-4</sup></td><td>36.41</td><td>1618.64</td><td>—</td><td>—</td><td>1621.01</td><td>—</td><td>350 × 133</td></tr><tr><td>34</td><td>10<sup>8</sup></td><td>2.5 × 10<sup>-4</sup></td><td>32.70</td><td>2353.72</td><td>—</td><td>—</td><td>2382.65</td><td>5521.6</td><td>350 × 133</td></tr><tr><td>35</td><td>10<sup>7</sup></td><td>10<sup>-4</sup></td><td>2.33</td><td>330.05</td><td>—</td><td>—</td><td>333.42</td><td>188.1</td><td>350 × 133</td></tr><tr><td>36</td><td>5 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>11.45</td><td>1281.03</td><td>—</td><td>—</td><td>1254.09</td><td>1622.3</td><td>350 × 133</td></tr><tr><td>37</td><td>7.5 × 10<sup>7</sup></td><td>10<sup>-4</sup></td><td>15.70</td><td>1570.99</td><td>—</td><td>—</td><td>1601.54</td><td>2477.0</td><td>350 × 133</td></tr><tr><td>38</td><td>10<sup>8</sup></td><td>10<sup>-4</sup></td><td>20.20</td><td>1775.39</td><td>—</td><td>—</td><td>1775.97</td><td>3350.6</td><td>400 × 150</td></tr><tr><td>39</td><td>2.5 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>39.70</td><td>3233.57</td><td>—</td><td>—</td><td>3282.39</td><td>—</td><td>420 × 150</td></tr><tr><td>40</td><td>5 × 10<sup>8</sup></td><td>10<sup>-4</sup></td><td>51.00</td><td>4483.30</td><td>—</td><td>—</td><td>4498.24</td><td>—</td><td>420 × 180</td></tr><tr><td>41</td><td>10<sup>9</sup></td><td>10<sup>-4</sup></td><td>77.50</td><td>5422.18</td><td>—</td><td>—</td><td>5484.67</td><td>—</td><td>450 × 180</td></tr><tr><td>42</td><td>10<sup>8</sup></td><td>5 × 10<sup>-5</sup></td><td>10.44</td><td>1788.85</td><td>—</td><td>—</td><td>1718.44</td><td>—</td><td>400 × 150</td></tr><tr><td>43</td><td>2.5 × 10<sup>8</sup></td><td>5 × 10<sup>-5</sup></td><td>22.97</td><td>3046.31</td><td>—</td><td>—</td><td>3029.92</td><td>—</td><td>420 × 150</td></tr><tr><td>44</td><td>5 × 10<sup>8</sup></td><td>5 × 10<sup>-5</sup></td><td>38.18</td><td>3780.74</td><td>—</td><td>—</td><td>3804.07</td><td>—</td><td>420 × 180</td></tr><tr><td>45</td><td>10<sup>8</sup></td><td>2.5 × 10<sup>-5</sup></td><td>3.80</td><td>1513.94</td><td>—</td><td>—</td><td>1518.53</td><td>—</td><td>400 × 150</td></tr><tr><td>46</td><td>5 × 10<sup>8</sup></td><td>2.5 × 10<sup>-5</sup></td><td>20.58</td><td>4898.98</td><td>—</td><td>—</td><td>4836.79</td><td>—</td><td>450 × 150</td></tr><tr><td>47</td><td>10<sup>9</sup></td><td>2.5 × 10<sup>-5</sup></td><td>36.67</td><td>6693.28</td><td>—</td><td>—</td><td>6664.02</td><td>—</td><td>480 × 200</td></tr><tr><td>48</td><td>2.5 × 10<sup>8</sup></td><td>10<sup>-5</sup></td><td>3.03</td><td>2144.76</td><td>—</td><td>—</td><td>2103.11</td><td>—</td><td>450 × 250</td></tr><tr><td>49</td><td>5 × 10<sup>8</sup></td><td>10<sup>-5</sup></td><td>6.33</td><td>3898.72</td><td>—</td><td>—</td><td>3894.77</td><td>—</td><td>450 × 250</td></tr><tr><td>50</td><td>10<sup>9</sup></td><td>10<sup>-5</sup></td><td>14.45</td><td>6403.12</td><td>—</td><td>—</td><td>6410.42</td><td>—</td><td>450 × 300</td></tr></tbody></table>