# Geochemistry, Geophysics, Geosystems®

## RESEARCH ARTICLE

10.1029/2024GC011891

**Key Points:**

* Rotation can structure a magma ocean in different domains with distinct stirring efficiency and limited mass exchange

* Stirring efficiency in each domain leads to homogenization of heterogeneities in timescales short compared to solidification onset

* The limited mass exchange between the domains could allow to preserve large-scale heterogeneous reservoirs over long timescales

**Correspondence to:**

H. Samuel,
samuel@ipgp.fr

**Citation:**

Thomas, B., Samuel, H., Farnetani, C. G., Aubert, J., & Chauvel, C. (2024). The influence of rotation on the preservation of heterogeneities in magma oceans.

*Geochemistry, Geophysics, Geosystems*, 25, e2024GC011891. https://doi.org/10.1029/2024GC011891

Received 19 SEP 2024

Accepted 7 OCT 2024

**Author Contributions:**

**Conceptualization:** B. Thomas,
H. Samuel, C. G. Farnetani, J. Aubert,
C. Chauvel

**Formal analysis:** B. Thomas, H. Samuel,
C. G. Farnetani, J. Aubert

**Funding acquisition:** C. Chauvel

**Investigation:** B. Thomas

**Methodology:** B. Thomas, H. Samuel,
C. G. Farnetani, J. Aubert

**Project administration:** C. Chauvel

**Resources:** H. Samuel, C. G. Farnetani

**Software:** B. Thomas, H. Samuel,
J. Aubert

**Supervision:** H. Samuel, C. G. Farnetani,
C. Chauvel

**Validation:** B. Thomas, H. Samuel,
C. G. Farnetani, J. Aubert

© 2024 The Author(s). Geochemistry, Geophysics, Geosystems published by Wiley Periodicals LLC on behalf of American Geophysical Union. This is an open access article under the terms of the Creative Commons Attribution-NonCommercial-NoDerivs License, which permits use and distribution in any medium, provided the original work is properly cited, the use is non-commercial and no modifications or adaptations are made.

## The Influence of Rotation on the Preservation of Heterogeneities in Magma Oceans

B. Thomas¹, H. Samuel¹, C. G. Farnetani¹, J. Aubert¹, and C. Chauvel¹

¹Université Paris Cité, Institut de Physique du Globe de Paris, CNRS, Paris, France

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

**Abstract** Understanding the composition of lavas erupted at the surface of the Earth is key to reconstruct the long-term history of our planet. Recent geochemical analyses of ocean island basalt samples indicate the preservation of ancient mantle heterogeneities dating from the earliest stages of Earth's evolution (Péron & Moreira, 2018, https://doi.org/10.7185/geochemlet.1833), when a global magma ocean was present. Such observations contrast with fluid dynamics studies which demonstrated that in a magma ocean the convective motions, primarily driven by buoyancy, are extremely vigorous (Gastine et al., 2016, https://doi.org/10.1017/jfm.2016.659) and are therefore expected to mix heterogeneities within just a few minutes (Thomas et al., 2023, https://doi.org/10.1093/gji/ggad452). To elucidate this paradox we explored the effects of the Earth's rapid rotation on the stirring efficiency of a magma ocean, by performing state-of-the-art fluid dynamics simulations of low-viscosity, turbulent convective dynamics in a spherical shell. We found that rotational effects drastically affect the convective structure and the associated stirring efficiency. Rotation leads to the emergence of three domains with limited mass exchanges, and distinct stirring and cooling efficiencies. Still, efficient convective stirring within each region likely results in homogenization within each domain on timescales that are short compared with the solidification timescales of a magma ocean. However, the lack of mass exchange between these regions could lead to three or four large-scale domains with internally homogeneous, but distinct compositions. The existence of these separate regions in a terrestrial magma ocean suggests a new mechanism to preserve distinct geochemical signatures dating from the earliest stages of Earth's evolution.

**Plain Language Summary** Geochemical heterogeneities from short-lived radionuclide parent-daughter systems date back to the very beginning of Earth's history, because the parent element became extinct a few tens of millions years after the solar system formation. Yet, geochemical heterogeneities from short-lived radioactive elements are observed in some present-day ocean island basalts. These geochemical observations prompted us to explore the mixing efficiency of the young Earth, when its silicate mantle was at a liquid state. We conducted numerical simulations in three-dimensional spherical geometry to study turbulent convection and stirring efficiency of magma ocean. For the flow regime with important rotational effects our results show that the magma ocean is separated in different domains (i.e., a polar, a columnar and an outer domain). Within each domain the convective stirring is very efficient, such that the mixing time of heterogeneities is short compared with the solidification timescales. Interestingly, we find that the mass exchange between domains is limited, leading to the possible preservation of large-scale (domain size) heterogeneous reservoirs. This mechanism could explain how geochemical heterogeneities from the early Earth might be observed in modern ocean island basalts.

## 1. Introduction

With a long-lasting plate tectonics history and its ongoing volcanic activity, the Earth is a dynamic planet. Erupted volcanic products convey key information on its dynamic interior, opening a rich time-integrated window into its 4.5 billion years long evolution. Indeed, the composition of lava samples collected at the surface of the Earth results from processes that shaped the planet throughout its history. In particular, the geochemical diversity of lava samples suggests that the Earth's mantle is heterogeneous at various scales, and the identification of distinct isotopic signatures in the geochemical record indicates that the observed mantle heterogeneity results from distinct source types of different ages, such as episodically or continuously recycled continental and oceanic crustal material, or ancient primordial matter (Ballentine et al., 2010; Hofmann, 1997; Weis et al., 2023; Zindler & Hart, 1986). Thanks to high-precision measurements of isotopic ratios from short-lived radionuclide parent-daughter pair system of Iodine (I)/Xenon (Xe), where $^{129}$I decays into $^{129}$Xe, and $^{130}$Xe is not radiogenic, low $^{129}$Xe/$^{130}$Xe ratios measured in OIB samples compared to Mid-Ocean Ridge Basalts indicate a source with an
**Visualization:** B. Thomas

**Writing – original draft:** B. Thomas, H. Samuel, C. G. Farnetani, J. Aubert, C. Chauvel

**Writing – review & editing:** B. Thomas, H. Samuel, C. G. Farnetani, J. Aubert, C. Chauvel

initially low I/Xe (Parai & Mukhopadhyay, 2021; Péron & Moreira, 2018). The half-life of ¹²⁹I (16 Myr) is so short that this parent element became extinct soon after the formation of the solar system, indicating that the low ¹²⁹Xe/¹³⁰Xe signatures found in OIBs are ancient and were created during an early magma ocean stage. Various interpretations have been suggested to explain the presence of the low ¹²⁹Xe/¹³⁰Xe values in plume sources. These include the delivery of exogenic material, the inefficient degassing of parts of the mantle, or a relatively dry mantle that would not mix efficiently (Caracausi et al., 2016; Pahlevan & Stevenson, 2007; Parai, 2022; Yuan et al., 2023). The subject is a matter of debate, and its timing relative to the Moon formation event remains unclear.

A silicate magma ocean is the likely result of the extreme events experienced by the young Earth, which caused extensive heating: the gravitational energy released by core formation (Monteux et al., 2009; Samuel et al., 2010), the intense early radiogenic heating (Urey, 1955), and the considerable amounts of energy generated by impacts of planetesimals and planetary embryos (Abe, 1997; Canup, 2008; V. S. Solomatov & Stevenson, 1993; Stevenson, 1987), including the giant Moon-forming impact that melted the entire Earth, resulting in one or several episodes of global silicate magma oceans (Rubie et al., 2015).

Contrary to the signatures of other short-lived isotopic systems that is, Samarium - Neodymium (¹⁴⁶Sm-¹⁴²Nd) that can be caused by early melting and solidification processes (Bouvier & Boyet, 2016; Boyet & Carlson, 2005; Caro et al., 2005), or Hafnium - Tungsten (¹⁸²Hf-¹⁸²W) which can result from chemical diffusion of W from the core into the mantle (Mundl-Petermeier et al., 2020; Rizo et al., 2019), leading to ambiguous and debated interpretations on their origins, the detection of ancient fingerprints in the I/Xe system suggests unequivocally the preservation of ancient mantle isotopic heterogeneities (Moreira et al., 2007; Parai & Mukhopadhyay, 2021; Parai et al., 2012; Tucker et al., 2012; Weis et al., 2023) dating from a time period concomitant with the presence of a terrestrial magma ocean (Abe, 1993, 1997; Lebrun et al., 2013; Nikolaou et al., 2019; Salvador et al., 2017). Indeed, regassing processes would not be sufficient to account for the observed ancient ¹²⁹Xe/¹³⁰Xe signatures (Mukhopadhyay, 2012; Parai & Mukhopadhyay, 2021; Parai et al., 2012; Weis et al., 2023). Addressing the key question of the origin of observed ancient isotopic fingerprints requires understanding the conditions under which heterogeneities present in a magma ocean can survive or be homogenized by vigorous convective motions.

The survival of ancient heterogeneities in a magma ocean is directly related to the ability of turbulent convective motions to deform and to progressively reduce the size of compositional variations that may eventually be homogenized by chemical diffusion on small length scales (Gurnis, 1986; Kellogg & Turcotte, 1987). The convective dynamics of a magma ocean strongly depends on the dynamic viscosity, η, of molten silicates, which is by far the most variable and least constrained physical quantity that affects magma ocean's dynamics. Importantly, magma ocean viscosity (η ~ 10⁻³–10² Pa s) is orders of magnitude lower than the viscosity of solid-state mantle, being similar to seawater viscosity, as shown by laboratory experiments (Liebske et al., 2005) and ab initio calculations (Karki & Stixrude, 2010). Such a low-viscosity magma ocean leads to an extremely vigorous and turbulent convection (Grossmann & Lohse, 2000), characterized by a considerably high Rayleigh number, Ra ∝ 1/η ~ 10³⁰ (V. Solomatov, 2015). The liquid viscosity also leads to a low Prandtl number (which expresses the ratio of momentum to thermal diffusion) Pr ∝ η, indicating that inertial effects are not negligible, contrary to solid-state convection.

Both high Ra and low Pr favor vigorous motion associated with large magnitudes of convective velocities (Gastine et al., 2016; Grossmann & Lohse, 2000), leading to a common assumption that convective mixing in silicate magma oceans is efficient and homogenization of heterogeneities is extremely rapid (Pahlevan & Stevenson, 2007). This assumption was tested in Thomas et al. (2023), which investigated the efficiency of stirring in magma oceans dominated by buoyancy forces, and showed that the strong deformation induced by the vigorous convective motion results in the homogenization of passive heterogeneities in just a few tens of hours. However, these findings contrast with the measurements of ancient isotopic signatures in the I/Xe system that indicate the preservation of heterogeneities in an early magma ocean.

Yet, while the low viscosity promotes vigorous convecting motion and efficient stirring, it also exposes the magma ocean to the influence of the planet's rotation rate Ω, by enhancing the Coriolis force. Rotational effects can be expressed by the Ekman number, Ek ∝ η/Ω. After the Moon-forming impact, the proto-Earth could have had a rotation period Ω of only 2–5 hr (Agnor et al., 1999; Ćuk & Stewart, 2012), corresponding to Ekman number
values as low as Ek ~ 10⁻¹⁵ for a magma ocean (Maas et al., 2021). This indicates the possibility of rotation exerting a significant role on magma ocean convective dynamics.

Previous studies of convection influenced by the Coriolis force (Aurnou et al., 2015; Kunnen, 2021; Plumley & Julien, 2019) have shown how convective motions are deviated into columns parallel to the rotational axis (Busse & Carrigan, 1976; Zhang, 1992). When the azimuthal and cylindrically radial components of the velocity field are correlated, these columnar structures create Reynolds stresses (i.e., stresses arising from the turbulent fluctuations of velocities) through inertia. These stresses generate large-scale zonal circulations which can dominate over natural convective motions (Christensen, 2002; Gastine et al., 2013). While the physical processes governing the transition from non-rotating to rotating convection have been described in Gastine et al. (2016) in spherical geometry, only a few studies investigated the effect of rotation on the compositional evolution of a magma ocean. Maas and Hansen (2015) studied the outcome of crystals in a rotating magma ocean and revealed that the location of crystal accumulation is influenced by the planet's rotational strength. Maas et al. (2021) studied the behavior of metal droplets resulting from impacts in a silicate magma ocean and showed that the extent of metal-silicate chemical equilibration depends on the latitude of the impact. These works underlined the potential influence of rotation on the fate of heterogeneities in a terrestrial magma ocean. To date, no systematic study has explicitly quantified the efficiency of convecting mixing in magma oceans strongly influenced by rotation.

## 2. Modeling

### 2.1. Fluid Dynamics Model

To model convection in a rotating terrestrial magma ocean in spherical geometry, we solve for Boussinesq convection in a spherical shell with inertial and rotational effects. To conform the model to the context of a silicate magma ocean, we use a constant gravitational acceleration throughout the shell. The inner and outer boundaries of the model domain are located at dimensionless radii $r_i = 1$ and $r_o = 2 r_i$, where free-slip conditions and a constant temperature difference $\Delta T$ are applied with $T_i = \Delta T$ and $T_o = 0$.

The characteristic scales used to normalize the governing equations are: $\Delta T$ for the temperature, $L = r_o - r_i$ (i.e., the thickness of the magma ocean shell) for the distances, $L^2 \rho / \eta$ for time where $\rho$ is the density. These result in the following set of dimensionless conservation equations for mass:

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

for energy:

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

and for momentum:

$$ \partial_t \mathbf{u} + \mathbf{u} \cdot \nabla \mathbf{u} - \nabla^2 \mathbf{u} + \frac{2}{\text{Ek}} \mathbf{e}_z \times \mathbf{u} + \nabla p = \frac{\text{Ra}}{\text{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 the unit vector parallel to the rotational axis, and $\mathbf{r}$ the unit vector parallel to the gravity vector.

Here we introduce the dimensionless numbers used in this set of equations: the Rayleigh number $Ra = \alpha g \Delta T L^3 \rho / (\eta \kappa_T)$ expresses the vigor of convection and corresponds to the ratio of the driving thermal buoyancy forces over the diffusive terms. The Prandtl number $Pr = \eta / (\rho \kappa_T)$ that expresses the ratio of momentum diffusivity to heat diffusivity (i.e., the ratio of viscous forces to inertia). And lastly, the Ekman number $Ek = \eta / (\rho \Omega L^2)$, whose inverse characterizes the relative importance of rotation over viscous effects. The physical parameters used here are listed in Table 1 with values characteristic of a magma ocean.
<table><caption>**Table 1**
Typical Values of Physical Parameters in a Terrestrial Magma Ocean</caption><thead><tr><th>Name</th><th>Symbol</th><th>Typical value</th></tr></thead><tbody><tr><td>Thermal expansion coefficient</td><td>α</td><td>3 10<sup>-5</sup> K<sup>-1</sup></td></tr><tr><td>Gravity</td><td>g</td><td>9.81 m s<sup>-2</sup></td></tr><tr><td>Temperature</td><td>T</td><td>4000 K</td></tr><tr><td>Magma ocean thickness</td><td>L</td><td>2.885 10<sup>6</sup> m</td></tr><tr><td>Dynamic viscosity</td><td>η</td><td>10<sup>-3</sup> – 10<sup>2</sup> Pa s</td></tr><tr><td>Density</td><td>ρ</td><td>4000 kg m<sup>-3</sup></td></tr><tr><td>Thermal diffusivity</td><td>κ<sub>T</sub></td><td>10<sup>-6</sup> m<sup>2</sup> s<sup>-1</sup></td></tr><tr><td>Angular velocity</td><td>Ω</td><td>7.3 10<sup>-5</sup> – 10<sup>-3</sup> rad s<sup>-1</sup></td></tr></tbody></table>

## 2.2. Stirring Model

To estimate the stirring efficiency we used Finite-Time Lyapunov Exponents, which provide a measure of the Lagrangian strain rate. Finite-time Lyapunov Exponents characterize the rate of deformation of fluid parcels along their trajectories (McKenzie, 1979). The approach to measure the Lagrangian strain rate is described in Thomas et al. (2023) and references therein. During the time period $t_1 - t_0$, a given fluid parcel (represented by passive Lagrangian tracers) is advected from $\mathbf{x}_0$ to $\mathbf{x}_1$ and its size changes form $\xi_0$ to $\xi_1$. This deformation is characterized by:

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

where $M$ is a linear operator containing the deformation history of the fluid parcel along the trajectory of the Lagrangian particle. $M$ is obtained by integrating:

$$ \frac{DM}{Dt} = D_x \mathbf{u}(\mathbf{x}(t), t)M, \quad (5) $$

using the velocity gradients determined (via tri-linear interpolations of grid values), along the particle trajectory. 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 principal orthogonal axes. The FTLEs are then obtained from $\sigma_i$ as follows:

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

In 3D space there are three eigenvalues corresponding to three different FTLEs: $\lambda_1 > \lambda_2 > \lambda_3$. The maximum (positive) FTLE, $\lambda_1$, characterizes the maximum stretching, this is the eigenvalue of interest for our problem and we refer below as $\lambda$.

To prevent bias from the initial condition of a numerical experiment, the computation of the Lagrangian strain rate described above is performed only once the flow has reached a statistically steady-state. We then randomly seed our domain with passive tracers (i.e., particles that do not affect the flow dynamics). The number of used tracers is of the same order of magnitude as the number of grid nodes used to discretize the governing conservation equations.

## 2.3. Chaotic Mixing Model

Mixing results from the combination of two processes (Ricard, 2015; Villermaux, 2019): (a) stirring induced by fluid motion that reduces the size of the heterogeneities through deformation, and (b) chemical diffusion that homogenizes the compositional heterogeneities. The latter process is only efficient on small scales and therefore comes into play once the heterogeneity size has been significantly reduced by stirring processes. We derive below a simple model that accounts for both aforementioned processes. The turbulent flow (and even steady laminar flow in 3D geometry (Dombre et al., 1986; Ferrachat & Ricard, 1998)) produces chaotic stirring via pure shear deformation, leading to the exponential stretching of fluid parcels with time (Gurnis, 1986; Kellogg & Turcotte, 1987). Consequently, the size of a heterogeneity $d(t)$ shrinks in one preferential direction (Gurnis, 1986; Samuel & King, 2014) as:

$$ d(t) = d_0 \exp(-\lambda t), \quad (7) $$

where $d_0$ is the initial size of the heterogeneity (assumed to be spherical) and $\lambda$ is the Lagrangian strain rate. Under the effect of pure shear, an heterogeneity is considered mixed once the effect from diffusion become relevant in front of the stretching induced by convective stirring. This is indicated by the Batchelor scale which characterize the length at which stretching from the stirring is perfectly compensated by the spreading from diffusion (Batchelor, 1959):
<table><caption>**Table 2**
Simulations in the Rotation-Dominated Regime Used for the Study of the Stirring Efficiency</caption><thead><tr><th>Experiment</th><th>Ra</th><th>Ek</th><th>Nu</th><th>Re<sub>c</sub></th><th>δ<sup>U</sup></th><th>δ<sup>L</sup></th><th>λ<sub>columnar</sub></th><th>λ<sub>outer</sub></th><th>λ<sub>polar</sub></th><th>nr × l</th></tr></thead><tbody><tr><td>1</td><td>5.00E+07</td><td>1.00E-04</td><td>11.45</td><td>4.86E+02</td><td>1.85E-02</td><td>1.55E-02</td><td>1414.00</td><td>1860.00</td><td>1906.00</td><td>400 × 200</td></tr><tr><td>2</td><td>7.50E+07</td><td>1.00E-04</td><td>15.70</td><td>6.69E+02</td><td>1.70E-02</td><td>1.44E-02</td><td>2115.00</td><td>2918.00</td><td>2759.00</td><td>400 × 200</td></tr><tr><td>3</td><td>1.00E+08</td><td>1.00E-04</td><td>20.20</td><td>8.38E+02</td><td>1.50E-02</td><td>1.43E-02</td><td>3129.00</td><td>4054.00</td><td>3929.00</td><td>400 × 200</td></tr><tr><td>4</td><td>1.00E+08</td><td>5.00E-05</td><td>10.44</td><td>5.10E+02</td><td>1.80E-02</td><td>1.10E-02</td><td>1633.00</td><td>2360.00</td><td>2359.00</td><td>400 × 200</td></tr><tr><td>5</td><td>2.50E+08</td><td>5.00E-05</td><td>22.97</td><td>1.20E+03</td><td>8.00E-03</td><td>1.10E-02</td><td>4443.00</td><td>6659.00</td><td>6808.00</td><td>400 × 200</td></tr><tr><td>6</td><td>2.50E+08</td><td>2.50E-05</td><td>11.21</td><td>5.48E+02</td><td>1.55E-02</td><td>7.60E-03</td><td>2730.00</td><td>3926.00</td><td>3757.00</td><td>400 × 200</td></tr><tr><td>7</td><td>5.00E+08</td><td>2.50E-05</td><td>20.58</td><td>1.30E+03</td><td>1.50E-02</td><td>7.60E-03</td><td>4101.00</td><td>9333.00</td><td>8598.00</td><td>420 × 300</td></tr><tr><td>8</td><td>1.00E+09</td><td>2.50E-05</td><td>36.67</td><td>2.28E+03</td><td>6.00E-03</td><td>8.00E-03</td><td>8009.00</td><td>17995.00</td><td>15600.00</td><td>420 × 300</td></tr><tr><td>9</td><td>5.00E+08</td><td>1.00E-05</td><td>6.33</td><td>5.66E+02</td><td>1.20E-02</td><td>5.00E-03</td><td>1662.00</td><td>3446.00</td><td>2861.00</td><td>420 × 300</td></tr><tr><td>10</td><td>7.50E+08</td><td>1.00E-05</td><td>10.50</td><td>7.75E+02</td><td>1.15E-02</td><td>5.00E-03</td><td>2500.00</td><td>6161.00</td><td>5720.00</td><td>420 × 300</td></tr><tr><td>11</td><td>1.00E+09</td><td>1.00E-05</td><td>14.45</td><td>1.00E+03</td><td>1.10E-02</td><td>4.50E-03</td><td>3128.00</td><td>8922.00</td><td>8493.00</td><td>420 × 300</td></tr><tr><td>12</td><td>7.50E+08</td><td>7.50E-06</td><td>6.70</td><td>6.32E+02</td><td>1.15E-02</td><td>3.50E-03</td><td>1741.00</td><td>4292.00</td><td>3998.00</td><td>420 × 300</td></tr><tr><td>13</td><td>5.00E+09</td><td>2.25E-06</td><td>9.88</td><td>1.18E+03</td><td>6.00E-03</td><td>1.80E-03</td><td>3149.00</td><td>13981.00</td><td>11400.00</td><td>460 × 340</td></tr></tbody></table>

Note. The first table lists the cases in the rotation-dominated regimes for which we conducted our study. The quantities of interest are listed: Ra, Ek, Nu, Re<sub>c</sub>, δ<sup>U</sup>, δ<sup>L</sup> and λ for all three regions and the resolution nr × l, the radial resolution nr and l controlling the truncation of the spherical harmonics.

$$l_B = \sqrt{\frac{\kappa_c}{\lambda}}, \qquad (8)$$

where $κ_c = 10^{-9} \text{ m}^2 \text{ s}^{-1}$ is the diffusion coefficient. Using the expression above for $d(t)$ confronted with $l_B$, we obtain the mixing time $τ$:

$$\tau = \frac{1}{2} \frac{\ln(\lambda d_0^2)}{\lambda_c}. \qquad (9)$$

As seen in the equation above, the dependence of the mixing time on the diffusivity $κ_c$ is weak (logarithmic). This indicates that the mixing behavior of distinct chemical species (e.g., noble gases like Xenon and trace elements like Uranium) even with significantly different values of chemical diffusivities will be essentially similar.

### 3. Characteristics of Rapidly Rotating Convective Flow

We conducted 27 numerical experiments (see Tables 2 and 3) where we explored the parameter range $Ek = [8.75 \times 10^{-7} - 10^{-4}]$ and $Ra = [10^7 - 6.25 \times 10^9]$ at constant $Pr = 1$. The value of $Pr$ was chosen because it corresponds to the case where momentum diffusion and thermal diffusion are comparable. This $Ra - Ek$ range enables us to reach regimes of convection, with the effects of rotation, which lie between the buoyancy-dominated regime and the onset of convection (i.e., weakly rotating convection), as defined in Thomas et al. (2023). The convective Rossby number $Ro_c = Ra^{1/2}Ek Pr^{-1/2}$ (Gilman, 1977), which characterizes the relative importance of the buoyancy forcing to the Coriolis force, was probably between ~0.03 and ~100 (Maas & Hansen, 2015; Moeller & Hansen, 2013). According to the boundary defined in Thomas et al. (2023) for the buoyancy-dominated regime, a convective Rossby number inferior to $O(10)$ would lead to a flow with significant rotational effects. In this study, the $Ro_c$ of our simulations ranges from 0.061 to 1.58 to focus on the regime where rotation is important. The other end-member case was studied in Thomas et al. (2023). Hence, in our numerical experiments the balance between rotation and convection is similar to what could be expected for a terrestrial magma ocean, and it satisfies the condition $Ro_c < 10$ of rotating convection.
<table><caption>**Table 3**
Simulations in the Transitional Regime or With Unstable Circulations Which We Did Not Use for the Study of the Stirring Efficiency</caption><thead><tr><th>Experiment</th><th>Ra</th><th>Ek</th><th>Nu</th><th>nr × l</th></tr></thead><tbody><tr><td>14</td><td>1.00E+07</td><td>1.00E-04</td><td>2.33</td><td>350 × 180</td></tr><tr><td>15</td><td>2.50E+08</td><td>1.00E-04</td><td>39.70</td><td>400 × 200</td></tr><tr><td>16</td><td>1.00E+08</td><td>2.50E-05</td><td>3.80</td><td>400 × 200</td></tr><tr><td>17</td><td>1.00E+08</td><td>1.00E-05</td><td>1.42</td><td>400 × 200</td></tr><tr><td>18</td><td>2.50E+08</td><td>1.00E-05</td><td>3.03</td><td>400 × 200</td></tr><tr><td>19</td><td>5.00E+08</td><td>7.50E-06</td><td>4.1</td><td>420 × 300</td></tr><tr><td>20</td><td>5.00E+08</td><td>5.00E-06</td><td>2.61</td><td>420 × 300</td></tr><tr><td>21</td><td>7.50E+08</td><td>5.00E-06</td><td>3.68</td><td>420 × 300</td></tr><tr><td>22</td><td>5.00E+08</td><td>3.75E-06</td><td>1.90</td><td>420 × 300</td></tr><tr><td>23</td><td>7.50E+08</td><td>3.75E-06</td><td>2.68</td><td>420 × 300</td></tr><tr><td>24</td><td>3.75E+09</td><td>1.00E-06</td><td>2.59</td><td>460 × 340</td></tr><tr><td>25</td><td>5.00E+09</td><td>1.00E-06</td><td>3.20</td><td>460 × 340</td></tr><tr><td>26</td><td>6.25E+09</td><td>1.00E-06</td><td>3.91</td><td>460 × 340</td></tr><tr><td>27</td><td>5.00E+09</td><td>8.75E-07</td><td>3.02</td><td>460 × 380</td></tr></tbody></table>

*Note.* The following table lists the simulations outside the rotation-dominated regime, either in the transitional regime or with unstable circulations. Only *Ra*, *Ek*, and *Nu* are specified since we did not focus on the other quantities.

### 3.1. Flow Structure Under the Effect of Rotation

To describe liquid convection influenced by rotation, we first focus on one simulation at $Ra = 2.5 \times 10^8$ and $Ek = 2.5 \times 10^{-5}$, thereafter called the reference case. The temperature field (Figure 1a) shows that the convective structures are deviated from the radial direction and form columns parallel to the rotational axis, as observed in previous studies (Busse & Carrigan, 1976; Christensen, 2002; Gastine et al., 2013; Wang et al., 2021; Zhang, 1992). Most importantly, the deviation induced by the rotational effects prevents hot material from reaching the outer edge of the shell around the equator. Consequently, in the equatorial region convective motions are inhibited by rotation (Gilman, 1977; Yadav et al., 2015). Such temperature field profoundly differs from the thermal structure observed in non-rotating convection (i.e., for higher values of *Ek*). This is because rotation organizes the flow in three main regions, as shown by the zonal velocity field, $V_{zonal}$, (i.e., rotating circulations) displayed in Figure 1b. The zonal velocities define three distinct large-scale circulations: (a) The innermost polar region, characterized by a prograde circulation with respect to the direction of rotation, located at the North and South poles on each side of the core; (b) a columnar region, characterized by a retrograde circulation, where the flow deviates the convective structures into columns; (c) the outermost-equatorial region, characterized by a prograde circulation. The boundaries between the distinct regions are defined by $V_{zonal} = 0$, which marks the transition between prograde/retrograde circulations. Such flow structure is also observed in other geophysical contexts (e.g., convection in liquid iron cores of terrestrial

planets) (Christensen, 2002; Zhang, 1992). In a flow dominated by rotation, the curvature at the boundary of the domain tilts the convective columns and leads to positive angular momentum, resulting in strong prograde circulations compensated by retrograde circulations. This leads to the structure observed in Figure 1b (Gastine et al., 2013). Similar structures, but with a different number of circulations, are observed in models of wind patterns for gas giants like Jupiter or Saturn (Giampieri & Dougherty, 2004; Heimpel & Aurnou, 2007; Lemasquerier et al., 2023), or with inverted zonal velocity directions in icy planets like Uranus (Fry et al., 2012) and Neptune (Martin et al., 2012), or even icy moons like Titan (Amit et al., 2020).

<figure><img src="image_2.png" /><figcaption>**Figure 1.** Results of the numerical experiments for the reference model with $Ra = 2.5 \times 10^8$ and $Ek = 2.5 \times 10^{-5}$.</figcaption></figure>

(a) Dimensionless temperature field. The convective heat transport is deviated from the radial direction and forms columns parallel to the rotational axis. This limits higher convective intensity at the outer edge of the equator far from the columns.

(b) Dimensionless zonal velocity field. The flow structure is organized in prograde (blue) and retrograde (orange) circulations. $V_{zonal} \sim 0$ (white) defines three regions, thereafter called polar, columnar and outer.
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Lateral snapshots of the zonal velocity fields for the case with inverted prograde/retrograde circulations. (a) Case with $Ra = 2.5 \ 10^8$, $Ek = 10^{-4}$. (b) The reference case at $Ra = 2.5 \ 10^8$, $Ek = 2.5 \ 10^{-5}$. (c) Case with no circulation at the poles at $Ra = 10^8$, $Ek = 2.5 \ 10^{-5}$. The values written inside refer to the maximum of the prograde (+) and retrograde (−) velocities in the snapshot.</figcaption></figure>

The observed flow structure is adequately described in terms of $r_h$, the non-dimensional horizontal radial distance from the rotational axis, with values ranging between $r_h = 0$ along the rotational axis, to $r_h = r_o = 2 r_i$ at the outermost edge of the spherical shell domain. For our reference case ($Ra = 2.5 \ 10^8$, $Ek = 2.5 \ 10^{-5}$), the boundaries between the three distinct regions inferred from $V_{zonal} = 0$ are located at $r_h \sim 0.7 r_i$ and at $r_h \sim 1.5 r_i$. For the Earth with a global magma ocean, the above dimensionless values for $r_h$ would correspond to 2,450 km and 5,250 km, respectively.

## 3.2. Limits of the Flow Regime

The flow regime observed for the reference case ($Ra = 2.5 \ 10^8$, $Ek = 2.5 \ 10^{-5}$, Figure 2b) is not the only distinct regime we observed in our numerical experiments. An inversion of the prograde and retrograde circulations is observed for only one case at the same $Ra = 2.5 \ 10^8$ but for a higher value of the Ekman number ($Ek = 10^{-4}$). Figure 2a shows a flow with an outermost, farthest from the rotation axis, retrograde circulation. This circulation is flanked by a large prograde circulation around the mid latitudes up to a small retrograde circulation at the poles (only visible at the South Pole for the snapshot displayed in Figure 2a). This configuration contrasts with the reference case displayed in Figure 2b. Note that the polar region is small and close to the rotation axis, thus its boundary no longer corresponds to the cylinder tangent to the core. Gastine et al. (2013) studied this transition of regime, which represents the change from a flow dominated by the rotation to a regime where the effect of the Coriolis force is of similar magnitude than the buoyancy forcing. This leads to a flow with a force balance closer to a quasi-geostrophic equilibrium, resulting in vigorous 3D convection. This 3D convection generates efficient stirring, which homogenizes the flow and its angular momentum, resulting in this new observed structure (Aurnou et al., 2007). We observe this transition by increasing $Ek$, at constant $Ra$, which corresponds to a reduction of the relative importance of the Coriolis force with respect to the buoyancy force. Since the relative importance of rotation and convection seems to control the regime transition (Aurnou et al., 2007; Gastine et al., 2013), we can use the convective Rossby number to characterize the boundaries of the regimes. The only simulation showing this regime transition has $Ro_c = 1.58$. Previous numerical studies, with similar conditions, observed the same transition at $Ro_c \ge \mathcal{O}(1-4)$ (Maas & Hansen, 2019; Yadav et al., 2015). In our case the transition occurs between $Ro_c = 1$ and $Ro_c = 1.58$. The reasonable assumption that the transition occurs at $Ro_c > \mathcal{O}(1-1.5)$, defines a new limit for the regime studied here, that we represent using the condition
<figure><img src="image_4.png" /><figcaption>**Figure 3.** Comparison between the reference case $Ra = 2.5 \ 10^8$, $Ek = 2.5 \ 10^{-5}$ (black line) and a case with the same $Ek$ but lower $Ra = 10^8$ (gray line). (a) The laterally averaged zonal velocity versus the distance to the rotation axis $r_h$. (b) The laterally averaged $Nu$ versus the distance to the rotation axis $r_h$. (c) $Nu$ for each region as a function of the supercriticality $R$. Below $R = 90$, the rotating circulation does not start in the polar region. This also implies that the polar region cannot be defined using the inversion of zonal circulation, and is instead approximated using the smallest polar region from our cases at $R > 90$. See text for further details.</figcaption></figure>

$Ro_c < 1.25$. This implies that the conditions for a rotation-dominated are now bounded within $0.03 < Ro_c < 1.25$. Therefore, when studying the stirring efficiency we will exclude the simulation with $Ro_c = 1.58$.

We observed another regime transition, this time caused by a reduced convective turbulence. Figure 2c shows the lateral cross section of numerical experiments at constant $Ek = 2.5 \ 10^{-5}$, but at $Ra = 10^8$. In this case, the prograde circulation in the polar region is no longer visible, indicating that a transition regime can also be induced by decreasing the influence of the buoyancy force relative to the Coriolis force. To better characterize this transition, we analyze the zonal velocity and the convective heat transfer as a function of $r_h$ (Figure 3).

Figure 3a shows the magnitude of zonal velocities for the reference case (black line) and for the case at the same $Ek$ but with a lower Rayleigh number $Ra = 10^8$ (gray line). At lower $Ra$, the prograde circulation completely stops in the polar region, resulting in the disappearance of the boundary between the polar and the columnar regions. Figure 3b shows the dimensionless convective heat transfer $Nu$, laterally averaged as a function of $r_h$. For the case at lower $Ra$ (gray line), the flow suffers a drastic decrease in convective vigor, both in the polar region and
<figure><img src="image_5.png" /><figcaption>**Figure 4.** *Ra* versus *Ek* regime diagram for the simulations conducted at *Pr* = 1. All the different regimes and boundaries are indicated, see text for further details. The (*Ra*,*Ek*) pairs in the white region correspond to the “rotation-dominated” regime. The black rectangle corresponds to our reference case displayed in Figure 3. Boundaries for the buoyancy-dominated regime and the weakly nonlinear regime were defined in Thomas et al. (2023).</figcaption></figure>

in the columnar region, compared to the reference case (black line). The presence of the circulations is due to the
coupling of convective velocities with the rotational effects (Christensen, 2002) and we can assume that the
transition results from the weak convective vigor that cannot initiate a stable circulation in the polar region. This
transition is not directly caused by the relative importance of the convective forcing to the Coriolis force, since the
transition is observed at distinct values of Ro_c and for different Ekman numbers. We rather conclude that the
regime transition is caused by a weakening of convection that prevents the start of the circulation in the polar
region.

To ensure that our simulations are in the regime where the circulation in the polar region is stable, we relied on the
supercriticality $R = Ra Ek^{4/3}$, which expresses the strength of the convective turbulence of a flow (Chan-
drasekhar, 1961). More specifically, reaching a turbulent regime of convection requires $R$ beyond a given
minimum threshold value (Plumley & Julien, 2019). We find that below $R = 90$ the convective turbulence is too
weak to create and sustain a circulation in the polar region. The zonal circulation and a stable convection at the
poles start only when $R \gtrsim 90$. Figure 3c shows *Nu* for the three regions as a function of *R*, for our complete set of
simulations.

Below the threshold $R = 90$ the values of $Nu$ are represented using the global average because the polar region is not defined. The values are small and increase following a similar trend to what is observed for the polar region at $R > 90$. Once $R > 90$, the polar circulation starts and the regime observed for the reference case is attained. Simulations at high supercriticality reach high $Nu$ values, characteristic of vigorous convection. Figure 3c provides us with another means to constrain the $Ra - Ek$ limit values appropriate to study steady rotating circulations: simulations with $R < 90$ should be discarded.

The complete regime diagram is shown in Figure 4, where the discarded simulations lie in the gray hatched area, labeled “Unstable circulations.” The constraint of steady rotating circulations reduces the number of simulations from 26 to 13, corresponding to those lying in the white domain in Figure 4, labeled “Strong rotation.” In the white domain, the (*Ra*, *Ek*) pairs ensure to be in the “rotation-dominated” regime we choose to study.
<figure><img src="image_6.png" /><figcaption>**Figure 5.** (a) Nusselt number versus colatitude at fixed Ekman number ($Ek = 2.5 \times 10^{-5}$) and variable Rayleigh number. (b) Nusselt number as a function of the supercriticality $R$. The error bars represent one standard deviation on the measurement of the time-averaged Nusselt number. The dashed line represents the scaling law $Nu = 0.12 R^{0.86}$. The cases from (a) are highlighted using their color and full symbols.</figcaption></figure>

In summary: to study the rotation-dominated regime of convection, we derived two main constraints. First, the constraint on stable polar circulation and turbulent convection requires $R > 90$. The regime below this boundary is irrelevant for a terrestrial magma ocean, because a fully molten magma ocean at $Ra = 10^{30}$ and $Ek = 10^{-16}$ would be associated with a supercriticality $R = 4.64 \times 10^8$ that is several orders of magnitude above 90. The second constraint requires that $Ro_c < 1.25$ to insure the Coriolis force dominates over the buoyancy forcing. The buoyancy-dominated regime has $Ro_c > 10$ (pink area in Figure 4), the rotation-dominated regime has $Ro_c < 1.25$ (white area in Figure 4), and for the intermediate $Ro_c$ (dark blue area) we found a flow regime with inverted circulations, which is out of the purpose of this study. Our numerical experiments with (Ra, Ek) pairs within the white area (Figure 4) enable us to study the end-member regime for a magma ocean with rotational effects. Characterizing the effect of rotation on the convective stirring efficiency is important because it allows comparing the mixing times of the rotation-dominated regime with those of the buoyancy-dominated regime.

# 4. Scaling Laws for Liquid Convection

## 4.1. Convective Heat Transfer

Here we focus on the behavior of the convective heat transfer in the rotation-dominated regime. Since the flow has lost its spherical symmetry, we have to consider how the Nusselt number varies as a function of the colatitude. Figure 5a shows $Nu$ versus the colatitude $\theta$ for three cases at fixed Ekman number $Ek = 2.5 \times 10^{-5}$ and variable Rayleigh number $Ra = 2.5 \times 10^8, 5 \times 10^8, 10^9$. The heat transfer is maximum at mid-colatitudes (polar and columnar regions) with two peaks at $\theta \sim 30^{\circ}$, and $\theta \sim 150^{\circ}$, while it is minimum in the equatorial region ($\theta = 90^{\circ}$). Heat transfer at mid-colatitudes increases strongly with increasing values of $Ra$, whereas at the equator, it increases only marginally. Previous studies of 3D spherical convection with rotation, such as Wang et al. (2021), also found that $Nu$ varies with latitude. However, free-slip and no-slip simulations display different characteristics, with a sharp increase of $Nu$ at the equator for no-slip boundary conditions (Yadav et al., 2015). For our free-slip simulations the small increase of $Nu$ at the equator is marginal compared to the increase observed at mid-colatitudes. In the rotation-dominated regime, the reduced convective heat flux at the equator could reduce the global heat transfer. This may increase the solidification time of the magma ocean compared to the buoyancy-dominated regime.

To assess the average cooling rate of a terrestrial magma ocean, we focus on the global Nusselt number averaged over the entire spherical shell surface. Figure 3c showed that in each region $Nu$ follows a comparable dependency with the supercriticality $R$. By performing a least-square fit of $Nu(R)$, as displayed in Figure 5b, we obtain the following scaling law:

$$Nu = 0.12 \pm 0.04 R^{0.86 \pm 0.06}, \qquad (10)$$
<figure><img src="image_7.png" /><figcaption>**Figure 6.** (a) Vertical cross-section with Lagrangian strain rate values calculated for the reference case ($Ra = 2.5 \ 10^8$ and $Ek = 2.5 \ 10^{-5}$). Black dashed lines indicate the boundary between the three regions. Lagrangian strain rate values are maximum in the columnar region, remain high in the polar region, but decrease in the outer region. (b) Histograms of Lagrangian strain rate distributions for the polar (purple), the columnar (red) and the outer (green) regions.</figcaption></figure>

showing that the intensity of the convective heat transfer is controlled by the supercriticality. It is important to note that this expression was obtained for a fixed Prandtl number, and therefore a dependency on Pr is expected for *Nu*. Constraining the dependence of *Nu* on *Pr* would require conducting simulations at various *Pr*, *Ra* and *Ek* numbers, which is out of the scope of this study. A similar result can be found in the literature, but with different values for the exponent of *R*. King et al. (2012) found a scaling law with a proportionality $Nu \propto R^3$ and Cheng et al. (2015) found $Nu \propto R^{3.6}$. These exponents are significantly higher than what we observe, probably because both studies were conducted in Cartesian box domains, which would represent only the region close to the rotation axis (similar to polar regions in spherical geometry). It is well known that such box geometry tends to increase the heat transfer compared to the spherical geometry (Cheng et al., 2015; King et al., 2012; Pueschel et al., 2013). The spherical configuration associated with the free-slip boundary conditions significantly reduces the convective vigor. Even compared to no-slip spherical convection, the values of *Ra* and *Pr* required to reach the same vigor of convection are considerably more extreme. This is the reason why we observe a trend for *Nu* distinct from the asymptotic scalings from the literature.

As shown in Figure 5a, the heat transfer is high in the polar and columnar regions, but low at the equator. This important difference could cause a reduction of the exponent in the scaling law $Nu(R)$ for the spherical geometry (Yadav et al., 2015), a scaling law only for the polar region yields $Nu \propto R$, and is closer to the scaling laws found in the literature. However, Gastine et al. (2016) found the scaling law $Nu \propto Ra^{3/2} Ek^2 = R^{3/2}$ in spherical geometry for no-slip boundary conditions, which favors efficient heat transfer (Julien et al., 2016). For our configuration, reaching a proportionality $Nu \propto R^{3/2}$ may require exploring lower values of the Ekman number, but we cannot provide a definitive conclusion at this stage.

## 4.2. Convective Stirring in Rotating Liquid Convection

Convective stirring is directly linked to the ability of the flow to produce pure shear deformation (Haller, 2001; Khakhar et al., 1984). We quantified the stirring efficiency by computing the Lagrangian strain rate (Khakhar et al., 1984) using Finite-Time Lyapunov exponents with the approach described in Farnetani and Samuel (2003) and in Thomas et al. (2023). The Lagrangian strain rate characterizes the rate of deformation of small, initially spherical, fluid parcels along Lagrangian particle trajectories (Farnetani & Samuel, 2003; McKenzie, 1979).
The measured Lagrangian strain rates for the reference case (Figure 6a) are high in the polar region, indicating efficient stirring, and reach a maximum in the columnar region, where the high deformations are induced by the zonal flow and by the turbulent motions occurring in the convective columns. With $r_h$ increasing toward the outer region, a sharp drop in strain rate is observed. This is due to the fact that in the outer region, turbulent motions are weaker and the convective motions are strongly deviated by rotation, resulting in a reduced stirring efficiency, because extreme rotation-induced motion promotes simple shear flow, which is associated with considerably weaker stirring efficiency than pure shear deformation triggered by buoyancy forces.

The obtained Lagrangian strain rate distribution (Figure 6b) is Gaussian for the polar region, indicating a homogeneous stirring efficiency in this domain (Farnetani & Samuel, 2003; Samuel & King, 2014), while the outer and the columnar regions show non-symmetrical deviations toward higher and lower values, respectively. These non-Gaussian deviations result from two effects: (a) convective columns create turbulent motions, causing interactions between the columnar and the outer regions, and (b) the distinct strain rate values in these two regions differ significantly, implying that even limited inter-region mass exchanges generate noticeable deviations in the corresponding distributions. The polar region instead is less affected by these processes because the neighboring columnar region exhibits similar Lagrangian strain rates. Therefore, the non-Gaussian distributions result essentially from limited inter-region exchanges, rather than reflecting the intrinsic stirring efficiency in each region. Consequently, in the following we neglect deviations from Gaussian distributions for all regions. This reasonable assumption enables us to describe the entire strain rate distribution in each region using only the arithmetic average dimensionless Lagrangian strain rate, $\lambda$, and the corresponding standard deviation, $\sigma_{\lambda}$. The distribution in each region is characterized by a distinct $\sigma_{\lambda}$, with no clear dependency on Ra and Ek. The size of the regions is specific to each simulation and is defined according to $r_h$, where $V_{zonal} = 0$ on the profile of $V_{zonal}(r_h)$. We report values of $\sigma_{\lambda}/\lambda \approx 0.17, 0.23$, and $0.27$ for the polar, the columnar, and the outer regions, respectively.

The distinct distributions of $\lambda$ suggest that, in this regime of convection, the efficiency of stirring is local to each region. We first attempted two-parameter fits within each region to obtain

$$ \lambda_{\text{polar}} = (0.022 \pm 0.007) Ra^{1.06 \pm 0.044} Ek^{0.80 \pm 0.077}, \quad (11a) $$

$$ \lambda_{\text{columnar}} = (0.0082 \pm 0.0014) Ra^{1.14 \pm 0.050} Ek^{0.84 \pm 0.076}, \quad (11b) $$

$$ \lambda_{\text{outer}} = (3.21 \pm 0.33) Ra^{0.88 \pm 0.039} Ek^{0.99 \pm 0.06}. \quad (11c) $$

We do however anticipate that stirring should be due to the same physical mechanism in the three regions, as also suggested by the similarity between the exponents obtained in Equations 11a–11c. The similar values of exponents obtained for dependencies in Ra and Ek also suggest a common dependence on the product of these two parameters. The result of the corresponding fit is

$$ \lambda_{\text{polar}} = (0.58 \pm 0.24) (Ra Ek)^{1.01 \pm 0.12}, \quad (12a) $$

$$ \lambda_{\text{columnar}} = (0.40 \pm 0.14) (Ra Ek)^{1.06 \pm 0.1}, \quad (12b) $$

$$ \lambda_{\text{outer}} = (0.35 \pm 0.14) (Ra Ek)^{0.99 \pm 0.05}. \quad (12c) $$

This confirms similar scaling laws in the three regions as well as a simple linear dependence on Ra Ek. By forcing this dependence we finally obtain

$$ \lambda_{\text{polar}} = (0.64 \pm 0.09) Ra Ek, \quad (13a) $$

$$ \lambda_{\text{columnar}} = (0.72 \pm 0.12) Ra Ek, \quad (13b) $$

$$ \lambda_{\text{outer}} = (0.32 \pm 0.05) Ra Ek, \quad (13c) $$

which may be expressed as $\lambda \propto \eta/(\rho L^2) Ra Ek$ using a viscous timescale. A subsequent normalization by $\Omega$: $\lambda/\Omega \propto \alpha g \Delta T L/(\Omega L)^2$ Pr further highlights the relative importance of the buoyancy forcing to the Coriolis
<figure><img src="image_8.png" /><figcaption>**Figure 7.** (a) Average Lagrangian strain rate $\lambda$ in each region as a function of the product Ra Ek. The scaling law for each region is represented by the dashed lines. The error bars represent one standard deviation uncertainty on the data. (b) Convective Reynolds number as a function of the product Ra Ek.</figcaption></figure>

force in this regime through the term $\alpha g \Delta T / (\Omega L)^2$. Note that the residues of the different fits for $\lambda$ are comparable: 11% for the three parameter, 17% for the two parameters, and 18% for the single parameter fit. Therefore, in the following we resort to the fit with a single prefactor parameter, which is more suitable for extrapolations over multiple orders of magnitude from our experimental range to magma ocean conditions.

Here differences in stirring efficiency between the three regions reduce to the scaling prefactors, with efficiency in the polar/columnar region being 1.5–4 times higher than in the outer region. The scaling laws from this set of equations are displayed in Figure 7a, along with the data points.

Another underlying implication of Equation 13 may be understood by using in dimensional form (i.e., $s^{-1}$), which leads to

$$ \lambda \propto \frac{\alpha g \Delta T}{\Omega L} \left( \frac{\eta}{\rho \kappa_T} \right). \qquad (14) $$

We note that the dependence in $Pr = \eta/(\rho\kappa_T)$ within this equation is not testable with our data set containing only cases at $Pr = 1$, such that normalizing time using $L^2/\kappa_T$, instead of $L^2\rho/\eta$ as fundamental time scale, leads to a scaling law where this dependence in $Pr$ also disappears. This suggests a secondary or vanishing control of diffusive effects on the mixing efficiency.

Following Thomas et al. (2023) we introduce the convective (or non-zonal) Reynolds number:

$$ Re_c = \frac{\rho U_{\text{conv}} L}{\eta} = \sqrt{2 \sum_{l=1}^{l_{\text{max}}} \sum_{m=1}^{l} e_l^m}. \qquad (15) $$

Here $e_l^m$ is the dimensionless kinetic energy density (i.e., the kinetic energy per unit volume) at spherical harmonics degree $l$ and order $m$, while the overbar indicates temporal averaging over the entire spherical shell domain. Starting from $m = 1$, the contribution from the axis-symmetric flows has been excluded. This provides a dimensionless measure of the Root Mean Squared (RMS) non-zonal convective velocity $U_{\text{conv}}$. Figure 7b also supports a dependency of $Re_c$ with the product Ra Ek:

$$ Re_c \approx 0.1 Ra Ek, \qquad (16) $$

with a possible $Pr$ dependency of $Re_c$ that we cannot investigate here. This result obtained with models where $Pr = 1$ is compatible with that obtained in Gastine et al. (2016), $Re_c \propto (Ra Ek/Pr)^{1.25}$ to $(Ra Ek/Pr)^1$. Assuming $Re_c \propto RaEk/Pr$ these results again suggest a simple dimensional law of the form
<figure><img src="image_9.png" /><figcaption>**Figure 8.** Regime diagram as a function of Ra, Ek for extreme values characteristic of a magma ocean. The point cloud represents magma ocean conditions. The color bar shows the convective Rossby number. Computations are conducted for the rotation-dominated regime and for the buoyancy-dominated regime and exclude the transitional regime, which we did not explore in this study.</figcaption></figure>

$$U_{\text{conv}} \propto \frac{\alpha g \Delta T}{\Omega}, \qquad (17)$$

and that, similarly to Thomas et al. (2023), the dimensional stirring efficiency is simply proportional to the shear of the convective flow:

$$\lambda \propto \frac{U_{\text{conv}}}{L}. \qquad (18)$$

We however note that the shear length scale in Thomas et al. (2023) was that of a velocity gradient layer near the boundaries, while in this case the relevant length scale appears to be commensurate with the size of the system. In a turbulent flow strongly constrained by rotation, we indeed expect an inverse cascade of energy leading to the formation of structures at the Rhines scale $\ell \sim L\sqrt{Re_cEk} \propto L Ro_c$, (e.g., Gastine et al., 2016) which is typically large in our data set. Because the Rhines scale is approached at levels of hydrodynamic turbulence significantly higher than those present in our models, our data set is unable to discriminate between this scale and the system size in Equation 18, but it remains possible that stirring is indeed controlled by shear at this turbulent macroscale. If this is the case, for parameter values typical of a rotation-dominated magma ocean (see Figure 8 below), the mixing times predicted with the system scale (as done in the next section) may be further decreased by a factor up to 100 to account for this.

## 5. Mixing Time of Heterogeneities in a Rotating Terrestrial Magma Ocean

To apply our results, we compute the values of Ra and Ek for physical parameters characteristic of a terrestrial magma ocean. The other governing parameter at play is the Prandtl number, Pr. In the context of terrestrial magma oceans all physical quantities but the viscosity appearing in Pr are well constrained. It is therefore reasonable to consider Pr primarily fixes the viscosity. Since our study was conducted at Pr = 1 we cannot freely vary the viscosity, which must fulfill the condition Pr = 1, leading to $\eta/\rho = \kappa_T$. However, the angular velocity can vary from $\Omega = 7.3 \ 10^{-5} \text{ rad s}^{-1}$, characteristic of modern Earth, to $\Omega = 10^{-3} \text{ rad s}^{-1}$, for a faster rotating Earth following giant impacts (Agnor et al., 1999; Ćuk & Stewart, 2012). Moreover, the thickness of the magma ocean on Earth can range from a fully molten mantle to a mantle molten only up to half its depth (Rubie et al., 2015). Finally, we choose to vary the temperature difference between the surface temperature and the potential temperature of the magma ocean from 100 to 2000 K. Temperature differences of 1 K would still lead to Ra values characteristic of a terrestrial magma ocean (Lebrun et al., 2013). The other relevant physical parameters are listed in Table 1.

We compute Ra, Ek, and $Ro_c$ using the parameters range described above, and we obtain the regime diagram shown in Figure 8. The point cloud is obtained by sampling randomly in the physical parameters range. We can see that at Pr = 1, the magma ocean has a high probability to be in the rotation-dominated regime. We speculate that higher values of Pr would tend to favor the buoyancy-dominated regime. Figure 8 also shows that a terrestrial magma ocean could be in the “transitional regime,” where the Coriolis force and the buoyancy force are close to equilibrium. Since our study does not characterize this regime we do not consider the simulations at $Ro_c > 1.25$. Therefore, the mixing times are computed only for simulations with $0.01 < Ro_c < 1.25$ and $10 < Ro_c < 100$.

Using the scaling laws for $\lambda$ (Equation 13) for the rotation-dominated regime, and the scaling law of $\lambda$ for the buoyancy-dominated regime from Thomas et al. (2023), we can compute the evolution of the size of an heterogeneity $d(t)$ for these two “end-member” flow regimes. Convective stirring reduces $d(t)$ exponentially (Hoffman & McKenzie, 1985; Khakhar et al., 1984) until it reaches the Batchelor scale (Batchelor, 1959). The mixing time of heterogeneities, $\tau$, is obtained by solving Equation 9 and is plotted versus $Ro_c$ in Figure 9. We considered a distribution of stirring efficiency sampled randomly within two standard deviations $\sigma_\lambda$ of the distribution of $\lambda$ for each case. In this way, in addition to randomly sampling the values of the physical quantities involved in the governing parameters, we also account for the Gaussian distribution of the Lyapunov exponents
<figure><img src="image_10.png" /><figcaption>**Figure 9.** Mixing time of heterogeneities as a function of the convective Rossby number for simulations at a terrestrial magma ocean's conditions. The mixing for the three regions in the rotation-dominated regime is represented as one single point cloud due to the similarity in mixing time. The mixing time in rotation-dominated convection is approximately six orders of magnitude longer than in the buoyancy-dominated regime but remains small (a few days), compared to Lebrun et al. (2013) estimates of solidification timescales for a magma ocean (yellow area).</figcaption></figure>

observed and quantified in our numerical experiments. For the rotation-dominated regime, we have different stirring efficiencies in each region (Figure 6). However, using our scaling laws in Equation 13, the discrepancies is only of approximately a factor 2, which is negligible compared to the difference in timescales shown in Figure 9. Hence, the mixing time in the rotation-dominated regime for the three regions is represented as only one point cloud.

We obtain mixing times ranging from ~50 days to ~1 hr in the rotation-dominated-regime. This is very short compared to timescales for the solidification of the magma ocean. The solidification time estimated in Lebrun et al. (2013) is shown in Figure 9 as the yellow region. It ranges from 10 years (e.g., Lebrun et al., 2013), corresponding to the formation of the first crystals (i.e., determined as the time at which the adiabat crosses the liquidus curve), to 1 Myr for the complete solidification of the magma ocean. There is a general consensus that a magma ocean started solidifying within a few tens of years and finished solidifying as late as a few tens of million years (Abe, 1993; Elkins-Tanton, 2008; Lebrun et al., 2013; Nikolaou et al., 2019). Magma ocean solidification is not expected to start at *t* = 0 because it is likely that the thermal state of the magma ocean was initially hotter than the liquidus of silicates (Lebrun et al., 2013; Nikolaou et al., 2019), in particular after a giant impact event (Rubie et al., 2015). However, these solidification timescales were derived from scaling laws on *Nu* for the buoyancy-dominated regime and so far, no study has quantified the solidification time of a magma ocean in the rotation-dominated regime. In addition, previous studies

assumed instantaneous exsolution of super saturated volatiles. This process increases the greenhouse effect, therefore it reduces the cooling and delays the onset of solidification in terrestrial magma oceans (Salvador & Samuel, 2023). On the other hand, in Section 4.1 we hinted at the possibility that the convective heat transfer could be reduced for a rotation-dominated magma ocean. Consequently, the magma solidification times inferred from coupled magma ocean-atmosphere studies remains associated with large uncertainties. Our mixing time estimates represent a lower bound because they do not explicitly account for the dispersion of the heterogeneities within each domain. While a good correspondence is found between mixing times inferred from the Lagrangian strain rate and the dispersion time of heterogeneities at infinite Prandtl number (Ferrachat & Ricard, 2001; Samuel et al., 2011), the dispersion time may be longer in finite Prandtl number flows, and may be governed by the overturn time. However, in a vigorously convecting magma ocean that is characterized by small overturn times these dispersion times would remain small compared to the life time of the magma ocean in the fully liquid state. Figure 9 also shows that in the buoyancy-dominated regime the mixing time can be orders of magnitude smaller than the mixing time in rotation dominated flow. As mentioned above, such small mixing times may overestimate stirring efficiency and may instead be bounded by the overturn time, which would still result in a mixing time considerably smaller than the lifetime of a magma ocean. Hence, for conditions of a fully molten terrestrial magma ocean, even if the rotation-dominated regime leads to longer mixing times than for the buoyancy-dominated regime, the convective mixing seems to homogenize passive heterogeneities in each of the different regions before the onset of solidification.

## 6. Mass Exchange Between Regions

Even though convective mixing tends to erase heterogeneities in the rotation-dominated regime, the flow structure may play an important role in preserving heterogeneities. The separation of the flow into three distinct regions (or four domains if the polar regions on the North and South hemispheres are disconnected on each side of the core) might be of interest for the survival of heterogeneities, because of the limited mass exchanges across each region. The distribution of heterogeneities in a rotating magma ocean has been studied in Maas et al. (2021), focusing on the settling of particles following impacts. The authors found differences in the settling time of heterogeneities according to the latitude of the impact, moreover they showed that particles tend to remain in the same region for extended periods of time. In our case, we do not account for active heterogeneities, therefore we
<figure><img src="image_11.png" /><figcaption>**Figure 10.** Percentage of tracers that escaped their starting region versus the number of transits, for the outer region (green), the columnar region (red) and the polar region (blue) for a case with $Ra = 2.5 \times 10^8$, $Ek = 2.5 \times 10^{-5}$, and $Ro_c = 0.39$.</figcaption></figure>

do not consider settling, but we can measure the mass exchanges induced by the flow between the different regions.

Using passive Lagrangian tracers, we measured the fraction of fluid parcels that leaves/enters each of the three identified regions. To account for small time perturbations of the boundary's position between the regions we added an error of $\pm 0.05 r_i$ on each side of the boundaries. For the reference case ($Ra = 2.5 \times 10^8$, $Ek = 2.5 \times 10^{-5}$, $Ro_c = 0.39$) we initially seeded $10^6$ tracers in each region. Figure 10 shows the percentage of tracers escaping each region over a time period corresponding to ~20 transits. A transit timescale refers to the characteristic travel time for fluid parcels within a region, and is computed as the ratio of the thickness of the spherical shell to the global RMS average of the convective velocities. If we extrapolate to magma ocean conditions $Ra \sim 10^{30}$ and $Ek \sim 10^{-15}$ using Equation 16, 20 transits are reached for a timescale of approximately 100 years, which may be comparable or larger than the lifetime of a magma ocean in the fully liquid stage (Lebrun et al., 2013), despite large uncertainties, as mentioned above.

The velocity which contributes to the escaping of tracers from one region to another is the convective velocities. However, part of the convective velocities are parallel to the rotation axis, transporting heat in the convective columns, and do not participate in the spreading of tracers to the different regions. Hence, we might overestimate the velocities contributing to the transit in this timescale. The percentage of escaping tracers reaches a constant value after ~15 transit time for the outer and polar regions, and after ~13 transit time for the columnar region.

We note that the evolution of the inter-region mass exchanges displays an asymptotic behavior, which is also observed in tests we conducted for other simulations, indicating that the obtained results would remain valid for longer periods of time. Overall, the estimates of inter-region mass exchange suggest that an important fraction of tracers would remain in the starting region. Consequently, heterogeneities may remain in the same region for long periods of time during the magma ocean stage. A stricter criteria for the dispersion of tracers would be to look at the exchanges between the North and South poles due to the presence of the core and the size of the polar regions. Indeed, the exchanges of heterogeneities between the North and South polar regions are expected to be even smaller than those displayed in Figure 10, mainly because the core here acts as an impermeable boundary between these antipodal regions.

A terrestrial magma ocean would be at similar $Ro_c$ as our simulations (Maas & Hansen, 2015, 2019; Moeller & Hansen, 2013) and we can assume that the asymptotic behavior would also be observed for the transfer of heterogeneities in a terrestrial magma ocean. A reduced mass exchange between the distinct regions could lead to the preservation of compositional differences among the regions, even if the convective stirring efficiency is high (Figure 9). For example, if during the accretion period a region was subject to several impacts with distinct compositions, then this region may develop and preserve a distinct geochemical signature with respect to the other regions. In other words, the limited mass exchange between the regions of a magma ocean would favor the preservation of distinct geochemical fingerprints within each regions. Such a global magma ocean would then be composed of large scale regions that are well mixed, which would however retain distinct geochemical signatures due to the reduced inter-region mixing. Upon cooling of the magma ocean the solidification would freeze the distinct signatures, and part of these could be preserved during the subsequent less vigorous convective evolution at solid-state stage.

These results have first-order implications for the interpretation of geochemical observations that indicate the preservation of ancient heterogeneities in Earth's present-day mantle. The limited inter-region mass exchanges implies that passive heterogeneities could survive a liquid magma ocean in a regime of rotation-influenced convection, at odds with previous studies (Pahlevan & Stevenson, 2007; Thomas et al., 2023) where the influence of rotation was not considered or had a minor effect on the dynamics. Moreover, our findings might solve the paradox of the survival of low $^{129}$Xe/$^{130}$Xe isotopic anomalies dating from the Hadean era observed in present-day OIBs (Moreira et al., 2007; Mundl et al., 2017; Parai & Mukhopadhyay, 2021; Rizo et al., 2019). Since the $^{129}$Xe/$^{130}$Xe anomalies cannot be created by regassing processes (Mukhopadhyay, 2012; Parai & Mukhopadhyay, 2021; Parai
et al., 2012; Weis et al., 2023) they indicate the preservation of mantle heterogeneities formed during the magma ocean period.

Interestingly, as the magma ocean solidifies, its viscosity increases resulting in changes in Ra, Ek and Pr numbers. This implies that the dynamical regime of the solidifying magma ocean will progressively change. Eventually the effects of the rotation would vanish, possibly leading to conditions corresponding to the transitional or even buoyancy-dominated regime. Hence, there could be scenarios of solidification where the magma ocean goes from rotation-dominated to buoyancy-dominated as the fraction of crystals and the viscosity increase, leading to an even more complex stirring history of the magma ocean. In particular, the flow regime of a magma ocean would be affected by a progressive viscosity increase, before reaching the transition to solid-state convection. This would require to study three regimes (i.e., buoyancy-dominated, transitional and rotation-dominated) at variable Pr, fully coupled with the crystallization processes. Accounting for such coupling remains challenging due to the coexistence of various times scales involved in these different regimes. The case of a crystallizing magma ocean that is not explicitly considered in our experiments would also require investigating the role of rigid boundary conditions at its base rather than free slip conditions in future studies.

## 7. Conclusions

We investigated the dynamical behavior of a magma ocean subject to important rotational effects using computational fluid dynamics experiments in a spherical shell geometry. We first characterized the structure of the flow, then constrained the conditions leading to the rotation-dominated regime. We then focused on the dependency between the heat transfer ($Nu$) and the supercriticality $R$, the convective turbulence of the flow. We found that $Nu$ is best described by a power law dependence, with $Nu \sim R^{0.8}$. This dependence of $Nu$ on the governing parameters is lower than what is observed in previous studies on rotating convection (Cheng et al., 2015; Gastine et al., 2016; King et al., 2012).

Rotational effects split the flow into three (or four) distinct regions: one (or two) polar region(s) centered along the rotation axis, next to a columnar region located further away from the rotation axis, and a last outermost region. Each region displays distinct stirring properties, notably, the stirring efficiency in the outer region is approximately half the stirring efficiency of the polar and columnar regions. For all three regions, the stirring efficiency depends on the product Ra Ek. The convective velocities also display a similar dependence ($Re_c \sim Ra Ek$) in agreement with Gastine et al. (2016). Even though rotation inhibits convection in the rotation-dominated regime, the stirring efficiency remains controlled by the intensity of the buoyancy-driven convective motions. This observation may be due to the chaotic nature of convective motions, as opposed to the organized structure of rotating circulations.

Convective stirring at magma ocean conditions leads to mixing times ranging from approximately 5 hr to 50 days for passive heterogeneities in the regime we considered with Pr = 1 and $0.01 < Ro_c < 1.25$. This is orders of magnitude longer than what is observed in the buoyancy-dominated regime (Figure 9), but remains short compared to the solidification timescales of a magma ocean.

Despite its efficient convective stirring, the structure of the rotation-dominated regime creates a new way to preserve heterogeneities. The mass exchange between the regions is limited and displays an asymptotic behavior in time. This flow structure creates three or four impermeable large-scale domains with internally homogeneous, yet possibly distinct compositions. This constitutes a novel mechanism to possibly preserve heterogeneous compositions in a vigorously convecting magma ocean.

In this study, and in a previous companion paper (Thomas et al., 2023), we characterized stirring efficiency for the two end-members regimes of a fully molten magma ocean, in the rotation-dominated and in the buoyancy-dominated regimes. However, a systematic study of the transitional regime is necessary to obtain a complete view of the mixing history in the magma ocean. This transitional regime is observed between the buoyancy-dominated and rotation-dominated regimes (Gastine et al., 2016). While we could not study this regime with our set of simulations, it may be relevant to a terrestrial magma ocean conditions (Figure 8) (Maas & Hansen, 2015, 2019; Moeller & Hansen, 2013). Even though it displays a roughly similar structure to the rotation-dominated regime, apart from the inverted circulations (Figure 2a), the balance between the buoyancy forces and the Coriolis force hints toward stronger convective motions, and could lead to distinct dynamics, stirring processes, and mass exchange between the regions that should be investigated in future studies.
Acknowledgments

This study benefited from discussions with Renaud Deguen and Thomas Gastine. We thank Yanick Ricard and an anonymous reviewer for their thoughtful comments that improved the manuscript, and the Editor Boris Kaus. We acknowledge 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 computations were performed on the S-CAPAD/DANTE platform, IPGP, France, and using the GENCI HPC resources of IDRIS under allocations A0110412958 and A0130412958.

Data Availability Statement

The data and codes used to produce the figures in this study are available at Zenodo via https://doi.org/10.5281/zenodo.12624417 (Thomas et al., 2024).

References

Abe, Y. (1993). Thermal evolution and chemical differentiation of the terrestrial magma ocean. In E. Takahashi, R. Jeanloz, & D. Rubie (Eds.), *Evolution of the Earth and planets* (pp. 41–54). American Geophysical Union AGU. https://doi.org/10.1016/B978-0-444-53802-4.00155-X

Abe, Y. (1997). Thermal and chemical evolution of the terrestrial magma ocean. *Physics of the Earth and Planetary Interiors, 100*(1), 27–39. https://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**(1), 219–237. https://doi.org/10.1006/icar.1999.6201

Amit, H., Choblet, G., Tobie, G., Terra-Nova, F., Čadek, O., & Bouffard, M. (2020). Cooling patterns in rotating thin spherical shells— Application to Titan’s subsurface ocean. *Icarus*, **338**, 113509. https://doi.org/10.1016/j.icarus.2019.113509

Aurnou, J., Calkins, M., Cheng, J., Julien, K., King, E., Nieves, D., et al. (2015). Rotating convective turbulence in Earth and planetary cores. *Physics of the Earth and Planetary Interiors*, **246**, 52–71. https://doi.org/10.1016/j.pepi.2015.07.001

Aurnou, J., Heimpel, M., & Wicht, J. (2007). The effects of vigorous mixing in a convective model of zonal flow on the ice giants. *Icarus*, **190**(1), 110–126. https://doi.org/10.1016/j.icarus.2007.02.024

Ballentine, C. J., Keken, P. E. V., Porcelli, D., & Hauri, E. H. (2010). Numerical models, geochemistry and the zero-paradox noble-gas mantle. *Philosophical Transactions of the Royal Society of London. Series A*, **360**(1800), 2611–2631.

Batchelor, G. K. (1959). Small-scale variation of convected quantities like temperature in turbulent fluid Part 1. General discussion and the case of small conductivity. *Journal of Fluid Mechanics*, 5(1), 113–133. https://doi.org/10.1017/S002211205900009X

Bouvier, A., & Boyet, M. (2016). Primitive solar system materials and Earth share a common initial $^{142}$Nd abundance. *Nature*, **537**(7620), 399–402. https://doi.org/10.1038/nature19351

Boyet, M., & Carlson, R. W. (2005). $^{142}$Nd evidence for early global differentiation of the silicate Earth. *Science*, **309**, 576–581. https://doi.org/10.1126/science.1113634

Busse, F. H., & Carrigan, C. R. (1976). Laboratory simulation of thermal convection in rotating planets and stars. *Science*, **191**(4222), 81–83. https://doi.org/10.1126/science.191.4222.81

Canup, R. M. (2008). Accretion of the Earth. *Philosophical Transactions of the Royal Society A*, **336**(1883), 4061–4075. https://doi.org/10.1098/rsta.2008.0101

Caracausi, A., Avice, G., Burnard, P. G., Füri, E., & Marty, B. (2016). Chondritic xenon in the Earth’s mantle. *Nature*, **533**(7601), 82–85. https://doi.org/10.1038/nature17434

Caro, G., Bourdon, B., Wood, B. J., & Corgne, A. (2005). Trace-element fractionation in Hadean mantle generated by melt segregation from a magma ocean. *Nature*, **436**(7048), 246–249. https://doi.org/10.1038/nature03827

Chandrasekhar, S. (1961). *Hydrodynamic and hydrodynamic stability*. Oxford University Press.

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. *Geophysical Journal International*, **201**(1), 1–17. https://doi.org/10.1093/gji/ggu480

Christensen, U. R. (2002). Zonal flow driven by strongly supercritical convection in rotating spherical shells. *Journal of Fluid Mechanics*, **470**, 115–133. https://doi.org/10.1017/S0022112002002008

Čuk, M., & Stewart, S. T. (2012). Making the Moon from a fast-spinning Earth: A giant impact followed by resonant despinning. *Science*, **338**(6110), 1047–1052. https://doi.org/10.1126/science.1225542

Dombre, T., Frisch, U., Greene, J. M., Hénon, M., Mehr, A., & Soward, A. M. (1986). Chaotic streamlines in the ABC flows. *Journal of Fluid Mechanics*, **167**(-1), 353–391. https://doi.org/10.1017/S0022112086002859

Elkins-Tanton, L. T. (2008). Linked magma ocean solidification and atmospheric growth for Earth and Mars. *Earth and Planetary Science Letters*, **271**(1–4), 181–191. https://doi.org/10.1016/j.epsl.2008.03.062

Farnetani, C. G., & Samuel, H. (2003). Lagrangian structures and stirring in the Earth’s mantle. *Earth and Planetary Science Letters*, **206**(3), 335–348. https://doi.org/10.1016/S0012-821X(02)01085-3

Ferrachat, S., & Ricard, Y. (1998). Regular vs. chaotic mantle mixing. *Earth and Planetary Science Letters*, **155**(1), 75–86. https://doi.org/10.1016/S0012-821X(97)00200-8

Ferrachat, S., & Ricard, Y. (2001). Mixing properties in the Earth’s mantle: Effects of the viscosity stratification and of oceanic crust segregation. *Geochemistry, Geophysics, Geosystems*, **2**(4), 1013. https://doi.org/10.1029/2000GC000092

Fry, P. M., Sromovsky, L. A., de Pater, I., Hammel, H. B., & Rages, K. A. (2012). Detection and tracking of subtle cloud features on Uranus. *The Astronomical Journal*, **143**(6), 150. https://doi.org/10.1088/0004-6256/143/6/150

Gastine, T., Wicht, J., & Aubert, J. (2016). Scaling regimes in spherical shell rotating convection. *Journal of Fluid Mechanics*, **808**, 690–732. https://doi.org/10.1017/jfm.2016.659

Gastine, T., Wicht, J., & Aurnou, J. M. (2013). Zonal flow regimes in rotating anelastic spherical shells: An application to giant planets. *Icarus*, **225**(1), 156–172. https://doi.org/10.1016/j.icarus.2013.02.031

Giampieri, G., & Dougherty, M. K. (2004). Rotation rate of Saturn’s interior from magnetic field observations. *Geophysical Research Letters*, **31**(16), L16701. https://doi.org/10.1029/2004GL020194

Gilman, P. A. (1977). Nonlinear dynamics of Boussinesq convection in a deep rotating spherical shell. *Geophysical & Astrophysical Fluid Dynamics*, **8**(1), 93–135. https://doi.org/10.1080/03091927708240373

Grossmann, S., & Lohse, D. (2000). Scaling in thermal convection: A unifying theory. *Journal of Fluid Mechanics*, **407**, 27–56. https://doi.org/10.1017/S0022112099007545

Gurnis, M. (1986). Stirring and mixing in the mantle by plate-scale flow: Large persistent blobs and long tendrils coexist. *Geophysical Research Letters*, **13**, 1474–1477. https://doi.org/10.1029/GL013i013p01474

Haller, G. (2001). Distinguished material surfaces and coherent structures in three-dimensional fluid flows. *Physica D: Nonlinear Phenomena*, **149**(4), 248–277. https://doi.org/10.1016/S0167-2789(00)00199-8

Heimpel, M., & Aurnou, J. (2007). Turbulent convection in rapidly rotating spherical shells: A model for equatorial and high latitude jets on Jupiter and Saturn. *Icarus*, **187**(2), 540–557. https://doi.org/10.1016/j.icarus.2006.10.023
Hoffman, N. R. A., & McKenzie, D. P. (1985). The destruction of geochemical heterogeneities by differential fluid motions during mantle convection. *Geophysical Journal International*, 82(2), 163–206. https://doi.org/10.1111/j.1365-246x.1985.tb05134.x

Hofmann, A. W. (1997). Mantle geochemistry: The message from oceanic volcanism. *Nature*, 385(6613), 219–229. https://doi.org/10.1038/385219a0

Julien, K., Aurnou, J. M., Calkins, M. A., Knobloch, E., Marti, P., Stellmach, S., & Vasil, G. M. (2016). A nonlinear model for rotationally constrained convection with Ekman pumping. *Journal of Fluid Mechanics*, 798, 50–87. https://doi.org/10.1017/jfm.2016.225

Karki, B., & Stixrude, L. (2010). Viscosity of MgSiO₃ liquid at Earth's mantle conditions: Implications for an early magma ocean. *Science*, 328(5979), 740–742. https://doi.org/10.1126/science.1188327

Kellogg, L. H., & Turcotte, D. L. (1987). Homogenization of the mantle by convective mixing and diffusion. *Earth and Planetary Science Letters*, 81(4), 371–378. https://doi.org/10.1016/0012-821X(87)90124-5

Khakhar, D., Chella, R., & Ottino, J. (1984). Stretching, chaotic motion, and breakup of elongated droplets in time dependent flows. In *Unknown host publication title*. National Autonomous University of Mexico.

King, E. M., Stellmach, S., & Aurnou, J. M. (2012). Heat transfer by rapidly rotating Rayleigh-Bénard convection. *Journal of Fluid Mechanics*, 691, 568–582. https://doi.org/10.1017/jfm.2011.493

Kunnen, R. P. J. (2021). The geostrophic regime of rapidly rotating turbulent convection. *Journal of Turbulence*, 22(4–5), 267–296. https://doi.org/10.1080/14685248.2021.1876877

Lebrun, T., Massol, H., Chassefière, E., Davaille, A., Marcq, E., Sarda, P., et al. (2013). Thermal evolution of an early magma ocean in interaction with the atmosphere. *Journal of Geophysical Research: Planets*, 118(6), 1155–1176. https://doi.org/10.1002/jgre.20068

Lemasquerier, D., Favier, B., & Le Bars, M. (2023). Zonal jets experiments in the gas giants? Zonostrophic regime. *Icarus*, 390, 115292. https://doi.org/10.1016/j.icarus.2022.115292

Liebske, C., Schmickler, B., Terasaki, H., Poe, B., Suzuki, A., Funakoshi, K., et al. (2005). Viscosity of peridotite liquid up to 13 GPa: Implications for magma ocean viscosities. *Earth and Planetary Science Letters*, 240(3–4), 589–604. https://doi.org/10.1016/j.epsl.2005.10.004

Maas, C., & Hansen, U. (2015). Effects of Earth's rotation on the early differentiation of a terrestrial magma ocean. *Journal of Geophysical Research: Solid Earth*, 120(11), 7508–7525. https://doi.org/10.1002/2015jb012053

Maas, C., & Hansen, U. (2019). Dynamics of a terrestrial magma ocean under planetary rotation: A study in spherical geometry. *Earth and Planetary Science Letters*, 513, 81–94. https://doi.org/10.1016/j.epsl.2019.02.016

Maas, C., Manske, L., Wünnemann, K., & Hansen, U. (2021). On the fate of impact-delivered metal in a terrestrial magma ocean. *Earth and Planetary Science Letters*, 554, 116680. https://doi.org/10.1016/j.epsl.2020.116680

Martin, S. C., de Pater, I., & Marcus, P. (2012). Neptune's zonal winds from near-IR Keck adaptive optics imaging in August 2001. *Astrophysics and Space Science*, 337(1), 65–78. https://doi.org/10.1007/s10509-011-0847-y

McKenzie, D. (1979). Finite deformation during fluid flow. *Geophysical Journal of the Royal Astronomical Society*, 58(3), 689–715. https://doi.org/10.1111/j.1365-246X.1979.tb04803.x

Moeller, A., & Hansen, U. (2013). Influence of rotation on the metal rain in a Hadean magma ocean. *Geochemistry, Geophysics, Geosystems*, 14(4), 1226–1244. https://doi.org/10.1002/ggge.20087

Monteux, J., Ricard, Y., Coltice, N., Dubuffet, F., & Ulvrovà, M. (2009). A model of metal-silicate separation on growing planets. *Earth and Planetary Science Letters*, 287(3–4), 353–362. https://doi.org/10.1016/j.epsl.2009.08.020

Moreira, M., Kunz, J., & Allègre, C. (2007). Rare gas systematics in popping rock: Isotopic and elemental compositions in the upper mantle rare gas systematics in popping rock. *Science*, 279(5354), 1178–1181. https://doi.org/10.1126/science.279.5354.1178

Mukhopadhyay, S. (2012). Early differentiation and volatile accretion recorded in deep-mantle neon and xenon. *Nature*, 486(7401), 101–104. https://doi.org/10.1038/nature11141

Mundl, A., Touboul, M., Jackson, M. G., Day, J. M. D., Kurz, M. D., Lekic, V., et al. (2017). Tungsten-182 heterogeneity in modern ocean island basalts. *Science*, 356(6333), 66–69. https://doi.org/10.1126/science.aal4179

Mundl-Petermeier, A., Walker, R. J., Fischer, R. A., Lekic, V., Jackson, M. G., & Kurz, M. D. (2020). Anomalous ¹⁸²W in high ³He/⁴He ocean island basalts: Fingerprints of Earth's core? *Geochimica et Cosmochimica Acta*, 271, 194–211. https://doi.org/10.1016/j.gca.2019.12.020

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? *The Astrophysical Journal*, 875(1), 11. https://doi.org/10.3847/1538-4357/ab08ed

Pahlevan, K., & Stevenson, D. J. (2007). Equilibration in the aftermath of the lunar-forming giant impact. *Earth and Planetary Science Letters*, 262(3–4), 438–449. https://doi.org/10.1016/j.epsl.2007.07.055

Parai, R. (2022). A dry ancient plume mantle from noble gas isotopes. *Proceedings of the National Academy of Sciences of the United States of America*, 119(29), https://doi.org/10.1073/pnas.2201815119

Parai, R., & Mukhopadhyay, S. (2021). Heavy noble gas signatures of the North Atlantic Popping Rock 2iID43: Implications for mantle noble gas heterogeneity. *Geochimica et Cosmochimica Acta*, 294, 89–105. https://doi.org/10.1016/j.gca.2020.11.011

Parai, R., Mukhopadhyay, S., & Standish, J. J. (2012). Heterogeneous upper mantle Ne, Ar and Xe isotopic compositions and a possible Dupal noble gas signature recorded in basalts from the Southwest Indian Ridge. *Earth and Planetary Science Letters*, 359, 227–239. https://doi.org/10.1016/j.epsl.2012.10.017

Péron, S., & Moreira, M. (2018). Onset of volatile recycling into the mantle determined by xenon anomalies. *Geochemical Perspectives Letters*, 9, 21–25. https://doi.org/10.7185/geochemlet.1833

Plumley, M., & Julien, K. (2019). Scaling laws in Rayleigh-Bénard convection. *Earth and Space Science*, 6(9), 1580–1592. https://doi.org/10.1029/2019EA000583

Pueschel, M. J., Terry, P. W., Jenko, F., Hatch, D. R., Nevins, W. M., Görler, T., & Told, D. (2013). Extreme heat fluxes in gyrokinetic simulations: A new critical *β*. *Physical Review Letters*, 110(15), 155005. https://doi.org/10.1103/PhysRevLett.110.155005

Ricard, Y. (2015). Physics of mantle convection. *Treatise on geophysics*, 7, 23–71. https://doi.org/10.1016/B978-0-444-53802-4.00127-5

Rizo, H., Andrault, D., Bennett, N. R., Humayun, M., Brandon, A., Vlastélic, I., et al. (2019). 182W evidence for core-mantle interaction in the source of mantle plumes. *Geochemical Perspectives Letters*, 11, 6–11. https://doi.org/10.7185/geochemlet.1917

Rubie, D. C., Nimmo, F., & Melosh, H. J. (2015). Formation of Earth's core. In G. Schubert & D. J. Stevenson (Eds.), *Treatise on geophysics* (2nd ed., Vol. 9, pp. 43–74). 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. *Journal of Geophysical Research: Planets*, 122(7), 1458–1486. https://doi.org/10.1002/2017JE005286

Salvador, A., & Samuel, H. (2023). Convective outgassing efficiency in planetary magma oceans: Insights from computational fluid dynamics. *Icarus*, 390, 115265. https://doi.org/10.1016/j.icarus.2022.115265
Samuel, H., Aleksandrov, V., & Deo, B. (2011). The effect of continents on mantle convective stirring. *Geophysical Research Letters*, 38(4), L04307. https://doi.org/10.1029/2010GL046056

Samuel, H., & King, S. D. (2014). Mixing at mid-ocean ridges controlled by small-scale convection and plate motion. *Nature Geoscience*, 7(8), 603–606. https://doi.org/10.1038/ngeo2208

Samuel, H., Tackley, P. J., & Evonuk, M. (2010). Heat partitioning in terrestrial planets during core formation by negative diapirism. *Earth and Planetary Science Letters*, 290(1–2), 13–19. https://doi.org/10.1016/j.epsl.2009.11.050

Solomatov, V. (2015). Magma oceans and primordial mantle differentiation. In G. Schubert (Ed.), *Treatise on geophysics* (2nd ed., pp. 81–104). Elsevier. https://doi.org/10.1016/B978-0-444-53802-4-00155-X

Solomatov, V. S., & Stevenson, D. J. (1993). Nonfractional crystallization of a terrestrial magma ocean. *Journal of Geophysical Research*, 98(E3), 5391–5406. https://doi.org/10.1029/92JE02579

Stevenson, D. J. (1987). Origin of the Moon - The collision hypothesis. *Annual Review of Earth and Planetary Sciences*, 15(1), 271–315. https://doi.org/10.1146/annurev.ea.15.050187.001415

Thomas, B., Samuel, H., Farnetani, C. G., Aubert, J., & Chauvel, C. (2023). Mixing time of heterogeneities in a buoyancy-dominated magma ocean. *Geophysical Journal International*, 236(2), 764–777. https://doi.org/10.1093/gji/ggad452

Thomas, B., Samuel, H., Farnetani, C. G., Aubert, J., & Chauvel, C. (2024). The influence of rotation on the preservation of heterogeneities in magma oceans [Dataset]. *Zenodo*. https://doi.org/10.5281/zenodo.12624417

Tucker, J. M., Mukhopadhyay, S., & Schilling, J.-G. (2012). The heavy noble gas composition of the depleted MORB mantle (DMM) and its implications for the preservation of heterogeneities in the mantle. *Earth and Planetary Science Letters*, 355–356, 244–254. https://doi.org/10.1016/j.epsl.2012.08.025

Urey, H. C. (1955). The cosmic abundances of potassium, uranium and thorium and the heat balances of the Earth, the Moon, and Mars. *Proceedings of the National Academy of Sciences*, 41(3), 127–144. https://doi.org/10.1073/pnas.41.3.127

Villermaux, E. (2019). Mixing versus stirring. *Annual Review of Fluid Mechanics*, 51(1), 245–273. https://doi.org/10.1146/annurev-fluid-010518-040306

Wang, G., Santelli, L., Lohse, D., Verzicco, R., & Stevens, R. J. A. M. (2021). Diffusion-free scaling in rotating spherical Rayleigh–Bénard convection. *Geophysical Research Letters*, 48(20), e2021GL095017. https://doi.org/10.1029/2021GL095017

Weis, D., Harpp, K. S., Harrison, L. N., Boyet, M., Chauvel, C., Farnetani, C. G., et al. (2023). Earth’s mantle composition revealed by mantle plumes. *Nature Reviews Earth & Environment*, 4(9), 604–625. https://doi.org/10.1038/s43017-023-00467-0

Yadav, R., Gastine, T., Christensen, U., Duarte, L., & Reiners, A. (2015). Effect of shear and magnetic field on the heat-transfer efficiency of convection in rotating spherical shells. *Geophysical Journal International*, 204(2), 1120–1133. https://doi.org/10.1093/gji/ggv506

Yuan, Q., Li, M., Desch, S. J., Ko, B., Deng, H., Garnero, E., et al. (2023). Moon-forming impactor as a source of LLVPs in Earth’s deep mantle. *Nature*, 623(7985), 95–99. https://doi.org/10.1038/s41586-023-06589-1

Zhang, K. (1992). Spiralling columnar convection in rapidly rotating spherical fluid shells. *Journal of Fluid Mechanics*, 236, 535–556. https://doi.org/10.1017/S0022112092001526

Zindler, A., & Hart, S. (1986). Chemical geodynamics. *Annual Review of Earth and Planetary Sciences*, 14(1), 493–571. https://doi.org/10.1146/annurev.ea.14.050186.002425