Chemical Convection and Stratification in the Earth's Outer Core

1. Introduction

Convection in the Earth's outer core is presently driven by buoyancy sources of both thermal and compositional origin. The thermal source results from the combination of the heat extracted by the mantle, the latent heat released at the bottom of the fluid core by the crystallization of the inner core and possibly the decay of radioactive elements. The chemical source is the consequence of the crystallization of the inner core which releases light elements into the fluid core (Braginsky, 1963; Braginsky and Roberts, 1995) and may also originate from chemical interactions with the mantle through magnesium exsolution (Badro et al., 2016; O'Rourke and Stevenson, 2016) or crystallization of SiO₂ (Hirose et al., 2017). A major difference between these buoyancy sources is tied to their molecular diffusivities which differ by several orders of magnitude, the chemical diffusivity being much weaker than its thermal counterpart. Different values have been proposed for the thermal and chemical molecular diffusivities, respectively κ^T and κ^ξ, and for the kinematic viscosity ν of liquid iron in core conditions, with typical ranges: κ^T ~ 5 × 10⁻⁶ − 10⁻⁵ m²s⁻¹ (Poirier, 1988; de Wijs et al., 1998; Roberts and King, 2013), κ^ξ ~ 5 × 10⁻⁹ − 2 × 10⁻⁸ m²s⁻¹ (Dobson, 2002; Pozzo et al., 2013; Ichikawa and Tsuchiya, 2015; Posner et al., 2017), and ν ~ 5 × 10⁻⁷ − 5 × 10⁻⁶ m²s⁻¹ (de Wijs et al., 1998; Vočadlo et al., 2000; Perrillat et al., 2010). These estimates lead to the following ranges for the Prandtl, Schmidt, and Lewis numbers (respectively, Pr = ν/κ^T, Sc = ν/κ^ξ, and Le = κ^T/κ^ξ = Sc/Pr) which are formed from the ratios of these diffusivities:

$\begin{array}{ll}Pr\simeq 0.05-1;Sc\simeq 25-1,000;Le\simeq 250-5,000,& (1)\end{array}$

the Schmidt number Sc being the chemical equivalent to the Prandtl number. Convection in the Earth's outer core is therefore a combination of thermal convection, characterized by a Prandtl number close to 1 or smaller, and chemical convection associated with a much larger Schmidt number.

Classically, the diffusivity difference is ignored in numerical simulations. The thermal and compositional fields are combined into one single variable, the so-called “codensity” (Braginsky and Roberts, 1995; Lister and Buffett, 1995), under the assumption that the mixing action of the unresolved sub-grid-scale turbulence would lead to comparable effective diffusivities for temperature and composition,

$\begin{array}{ll}{\kappa }^T={\kappa }^{\xi }={\kappa }^{\text{turb}},& (2)\end{array}$

where κ^turb is a homogeneous and isotropic “turbulent” or “eddy” diffusivity. Since turbulent mixing would affect all quantities in a similar way, the Prandtl, Schmidt, and Lewis numbers should be of order one:

$\begin{array}{ll}Pr=Sc=Le=1.& (3)\end{array}$

This approximation is convenient to implement in numerical codes but relies on debatable approximations. First, the turbulent regime of the core is not known, though under the effect of rotation it is probably anisotropic (Braginsky and Roberts, 1995; Matsushima, 2004) as suggested by the low value of the Rossby number Ro = U/DΩ ~ 10⁻⁶, where U, D, and Ω are, respectively the typical velocity of the fluid, the thickness of the liquid core and the Earth's rotation rate. Rotating MHD turbulence is extremely complex and some of the results of the classic theory of turbulence (established for non-rotating and non-magnetic cases) may not apply in this context (Nataf and Gagniere, 2008; Nataf and Schaeffer, 2015). Turbulence can also be studied through different prisms: while some models favor a spectral description, other authors (Loper and Moffatt, 1993; Moffatt and Loper, 1994; Loper, 2007) consider buoyant plumes or blobs as the elementary structures at the origin of turbulent motions in the core. It is unlikely that the “turbulent” diffusivity parameterization accounts for all these complexities, even though it could provide a reasonable first-order description in some specific turbulent regimes. Second, the codensity approximation is particularly problematic where the net contribution of both fields produces a stable density gradient. In such context, turbulence is probably partially or totally suppressed and any notion of “turbulent” diffusivity becomes highly questionable (Braginsky and Roberts, 1995). The existence of stably stratified layers has been evoked for several planetary cores (Christensen, 2006, 2015; Manglik et al., 2010) including that of the Earth (Gomi et al., 2013). These arguments provide motivation for treating the temperature and composition fields separately without invoking any uncertain parametrization of the turbulence. Furthermore, although certainly still a matter of debate (Konôpková et al., 2016), recent high estimates of the thermal conductivity of the core (de Koker et al., 2012; Pozzo et al., 2012; Gomi et al., 2013; Ohta et al., 2016) suggest that, at present, a large fraction of the heat leaving the Earth's core is simply conducted down the adiabat and is not available for driving thermal convection. Compositional convection may thus play a crucial role in powering Earth's dynamo (Lister and Buffett, 1995; Labrosse, 2015) so that, to first-order, the dynamics of the fluid core can be modeled by considering only compositional convection with a large Schmidt number. Very different properties can be expected for the dynamics of rapidly-rotating convection at high Prandtl/Schmidt numbers compared to more moderate values (Busse, 2002; Simitev and Busse, 2003, 2005). However, the large Prandtl/Schmidt number case has not yet been systematically explored in numerical simulations.

In addition to the high value of the Schmidt number, another fundamental specificity of compositional convection arises from the complexity of the chemical processes occurring at the “mushy layer” that may exist at the interface between the fluid outer core and the crystallizing inner core. Experiments conducted on model systems involving the freezing of ammonium chloride in solution (Chen and Chen, 1991; Huguet et al., 2016) showed that the gravitational instability manifests itself through eruptions of buoyant fluid in the form of chemical plumes emerging through localized “chimneys” that are spontaneously created in the mushy layer. Fluid from the bulk core percolates through the mush and becomes gradually lighter as it converges horizontally into the chimneys. Contrary to thermal plumes, which are more rapidly diffused, chemical plumes emanating from the mushy layer can be expected to keep a coherent shape for a longer time due to their lower molecular diffusivity. This led several authors to speculate about the potential role played by such chemical plumes on the turbulence, convective dynamics, and dynamo action (Gubbins et al., 1982; Moffatt, 1989; Loper and Moffatt, 1993; Bergman and Fearn, 1994; Braginsky, 1994; Moffatt and Loper, 1994; Loper, 2007). Moffatt and Loper (1994) conjectured that light fluid erupts from the mushy layer rather in the form of chemical “blobs” which would keep a coherent shape throughout their path across the core controlled by a magnetostrophic force balance, a prediction challenged by St Pierre (1996) who argued, using numerical simulations, that such blobs would rapidly break-up into plate-like structures elongated in the direction of the rotation axis. The existence of chemical blobs was then reported a few years later in laboratory experiments as the result of instabilities causing the chemical plumes to bend horizontally and break-up (Claßen et al., 1999). However, the authors did not observe the blobs elongation anticipated by St Pierre (1996), possibly because these experiments are good analogs of only the polar regions on Earth. The predictions for the existence and dynamics of such chemical structures have not yet been verified in global-scale numerical simulations.

A second interesting phenomenon envisioned by Moffatt and Loper (1994) is that the blobs may not completely mix with the bulk of the fluid core but accumulate below the core mantle boundary, thereby gradually forming a chemically stratified layer at the top of the core via a mechanism analogous to the “filling-box” models (Baines and Turner, 1969). This prediction is particularly noteworthy since the existence of a stably stratified layer in the outermost region of the core has been suggested by numerous seismic observations for several decades (Lay and Young, 1990; Tanaka, 2007; Helffrich and Kaneshima, 2010; Kaneshima, 2018). Buffett (2014) also invokes the presence of a stratified layer below the CMB to explain the 60-year period in the geomagnetic secular variation as well as fluctuations in the dipole field. A stably stratified layer at the top of the core is expected to have important implications for the geodynamo, the dynamics, and evolution of the core (Christensen and Wicht, 2008; Nakagawa, 2011; Helffrich and Kaneshima, 2013; Labrosse, 2015; Olson et al., 2017). Though several mechanisms have been proposed to explain the formation of such a layer, its origin remains unknown. The layer could be chemically stratified and produced by the accumulation of “chemical blobs” as proposed by Moffatt and Loper (1994), chemical exchanges with the mantle (Buffett and Seagle, 2010) or result from the incomplete mixing between the primitive core of the Earth and that of a giant impactor (Landeau et al., 2016). Alternatively, the layer could be thermally stratified if the CMB heat flow is sub-isentropic (Gubbins et al., 1982; Labrosse et al., 1997; Lister and Buffett, 1998; Gomi et al., 2013; Labrosse, 2015; Olson et al., 2017, 2018; Christensen, 2018). So far, no numerical simulation has been able to self-consistently produce a stable chemically stratified layer at the top of the core that would emerge from the “chemical blobs” scenario advocated by Moffatt and Loper (1994). The reason may be partly inherent to the codensity framework in which the large majority of the simulations have been conducted up to the present day.

Here, we propose to explore pure chemical convection at high Schmidt number by means of numerical simulations and test whether the accumulation of chemically buoyant parcels of fluid can lead to the formation of a stably stratified layer at the top of the core. Focusing on compositional effects only allows us to gain insight into the specific phenomena involved, independently of an additional forcing of thermal buoyancy. However, though moderately large Prandtl numbers can be reached with the classical codes (for instance, Simitev and Busse, 2005 use up to Pr = 20), the study of larger Prandtl or Schmidt numbers requires prohibitive resolutions. Here we use a recently developed version of the PARODY code including a particle-in-cell method (Bouffard et al., 2017) to conduct a first exploratory study of rapidly-rotating pure chemical convection in the end-member case of infinite Schmidt number. Note that we consider Sc = +∞ by setting κ^ξ = 0 and that our approach thus differs from that of Zhang and Busse (1990), Zhang (1991, 1992) who achieved infinite Prandtl number by using an infinite viscosity (ν = +∞). As will be shown, the phenomena involved in rotating pure compositional convection are rich and numerous. In this exploratory study, we aim at giving a first descriptive and qualitative insight into this particular physics. We also assess whether the stably stratified layer at the top of the core can be formed by the accumulation of chemically buoyant parcels of fluids emanating from the inner core boundary.

This paper is organized as follows. In section 2, the mathematical model and the numerical method are described. The results of a series of numerical simulations of pure compositional convection are then presented. The chemical structures and their dynamics are first studied in section 3. In section 4 we show that, when no light elements are allowed to cross the core-mantle boundary, a chemically stratified layer systematically forms at the top of the spherical shell in all simulations. The potential implications for the Earth's outer core are discussed in section 5. The final section 6 gives a general conclusion.

2. Mathematical and Numerical Modeling

2.1. Principal Equations

The fluid considered is a mixture of a heavy component, representing iron and nickel, and a light constituent, whose compositional mass fraction is denoted by ξ. As thermal effects are neglected here, the density ρ can be modeled as a function of the compositional mass fraction only,

$\begin{array}{ll}\rho ={\rho }_0\left[1-\beta (\xi -{\xi }_0)\right],& (4)\end{array}$

where ξ₀ and ρ₀ are the reference compositional mass fraction and density, respectively, and β denotes the coefficient of chemical expansivity.

In addition to the mass conservation,

$\begin{array}{ll}\overrightarrow{\nabla }·\overrightarrow{\text{u}}=0,& (5)\end{array}$

we solve for the momentum conservation, expressed in the form of the well-known Navier-Stokes equation, for a Newtonian rotating fluid in the Boussinesq approximation,

$\begin{array}{ll}\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}·\overline{\overline{\nabla }}\overrightarrow{u}+2\text{Ω}({\overrightarrow{e}}_{\text{z}}\times \overrightarrow{u})=-\frac{1}{{\rho }_0}\overrightarrow{\nabla }\Pi +\nu {\overrightarrow{\nabla }}^2\overrightarrow{u}+\beta \xi g_o\frac{\overrightarrow{r}}{r_o},& (6)\end{array}$

where $\overrightarrow{u}$ is the fluid velocity, t the time, Ω the Earth's rotation rate around an axis oriented by ${\overrightarrow{e}}_z$, Π the dynamic pressure and g_o the gravitational acceleration at the top of the shell at radius r = r_o. The Lorentz force is neglected. We adopt the viscous time scale D²/ν in which the thickness of the shell D = r_o − r_i is taken as the typical length scale, r_i being the radius of the inner boundary. Hence, using (ρ₀νΩ) as a scaling for pressure, one can write the dimensionless momentum equation,

$\begin{array}{ll}E\left(\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}·\overline{\overline{\nabla }}\overrightarrow{u}-{\overrightarrow{\nabla }}^2\overrightarrow{u}\right)+2{\overrightarrow{e}}_{\text{z}}\times \overrightarrow{u}+\overrightarrow{\nabla }\text{Π}=R{a}_{\xi }\xi \frac{\overrightarrow{r}}{r_o},& (7)\end{array}$

in which E is the Ekman number, defined by

$\begin{array}{ll}E=\frac{\nu }{\Omega D^2},& (8)\end{array}$

and Ra_ξ a “modified” compositional Rayleigh (or Roberts-Rayleigh) number,

$\begin{array}{ll}R{a}_{\xi }=\frac{g_0\beta C^{*}D}{\nu \Omega },& (9)\end{array}$

where C*/D is a measure for the global compositional variation driving convection. C* depends on the prescribed boundary conditions which are discussed in section 2.2. Note that, in the infinite Prandtl number case, the classical Rayleigh number constructed using the time scale $D^2/{\kappa }_{\xi }$ based on the chemical diffusivity,

$\begin{array}{ll}R{a}_{\xi }^{{\kappa }^{\xi }}=\frac{g_0\beta C^{*}D}{{\kappa }^{\xi }\Omega },& (10)\end{array}$

cannot be used as it is infinite, which means that even a very small value of Ra_ξ would lead to buoyancy-driven fluid motion and the notion of super-criticality is no longer relevant in this case.

In the Earth's fluid core, chemicals can be transported by advection and diffusion. The corresponding equation yields

$\begin{array}{ll}\frac{\partial \xi }{\partial t}+\overrightarrow{\text{u}}·\overrightarrow{\nabla }\xi ={\kappa }^{\xi }{\overrightarrow{\nabla }}^2\xi +{\sigma }^{\xi },& (11)\end{array}$

in which σ^ξ is a numerical sink/source term ensuring that the mean composition remains constant throughout the simulation (see section 2.2 for details). The equivalent dimensionless equation reads

$\begin{array}{ll}\frac{\partial \xi }{\partial t}+\overrightarrow{\text{u}}·\overrightarrow{\nabla }\xi =\frac{1}{Sc}{\overrightarrow{\nabla }}^2\xi +{\sigma }^{\xi }.& (12)\end{array}$

In the end-member case Sc = +∞ (κ^ξ = 0) considered here, Equation (12) becomes a first-order hyperbolic equation:

$\begin{array}{ll}\frac{\partial \xi }{\partial t}+\overrightarrow{\text{u}}·\overrightarrow{\nabla }\xi ={\sigma }^{\xi }.& (13)\end{array}$

2.2. Boundary Conditions

We adopt non-penetration and no-slip conditions at both boundaries which implies

$\begin{array}{ll}\overrightarrow{\text{u}}(r=r_o)=\overrightarrow{\text{u}}(r=r_i)=\overrightarrow{0}.& (14)\end{array}$

The first-order nature of Equation (13) combined with the condition (14) on the velocity implies that no chemical “signal” imposed at a boundary can propagate inside the shell, neither by advection nor diffusion. The value of the composition at a boundary has therefore no effect on the compositional field inside the shell and will in turn not be affected by the latter. Consequently, once the initial condition has been specified (for instance ξ = 0 everywhere), there is mathematically no need to impose any boundary condition for the compositional field. At the top boundary, we do not consider chemical interactions with the mantle so that the compositional flux equals zero. This physical condition is automatically satisfied by the mathematical formulation of the problem as explained above. At the bottom boundary, however, we wish to model a total flux of light elements $I_i^{\xi }$ resulting from the crystallization of the inner core. Since this cannot be accommodated in the framework of our mathematical formulation, we inject the mass flux $I_i^{\xi }$ via a source term ${\sigma }_i^{\xi }$ distributed within a thin layer of thickness h above the ICB, so that the integrated added mass fraction matches the ICB chemical flux:

$\begin{array}{ll}4\pi {\displaystyle {\int }_{r=r_i}^{r_i+h}{r}^2}{\sigma }_i^{\xi }(r)\text{d}r=I_i^{\xi }.& (15)\end{array}$

The choice of an adequate function for ${\sigma }_i^{\xi }(r)$ requires some care. Indeed, since the velocity vanishes in the Ekman layer to accommodate condition (14), any buoyant fluid injected at the bottom of this layer tends to “stick” to the boundary in the absence of chemical diffusion, hardly moves upward and does not participate in the convective dynamics. Adding the light elements slightly above in a region that extends beyond the Ekman layer (see Figure 1) allows some destabilization and prevents a detrimental accumulation of light elements in the layer over time. This can be performed smoothly using a Gaussian-like function:

$\begin{array}{ll}{\sigma }_i^{\xi }(r)=\frac{{\sigma }_0}{r^2}exp\left\{-{\left[r-\left(r_i+\frac{h}{2}\right)\right]}^2/{\sigma }^2\right\}.& (16)\end{array}$

The constants in the exponential function are empirically set to σ = h/5 and h = 0.4E^1/3. These values will be justified later in section 3.1. The scaling factor σ₀ is computed to satisfy Equation (15). Mathematically, this formulation is equivalent to solving the transport equation

$\begin{array}{ll}\frac{\partial \xi }{\partial t}+\overrightarrow{u}·\overrightarrow{\nabla }\xi ={\sigma }^{\xi }+{\sigma }_i^{\xi }H\left[(r_i+h)-r)\right]& (17)\end{array}$

and imposing no boundary conditions for the composition (at least, “mathematically” speaking). The injection of light elements occurring in the ICB region is accounted for by the local source term ${\sigma }_i^{\xi }{H}\left[(r_i+h)-r)\right]$, where ${H}$ is a Heaviside function, and is balanced by a uniform sink term σ^ξ in the entire shell so that the mean composition remains constant with time:

$\begin{array}{ll}{\displaystyle {\int }_{r_i}^{r_o}\left({\sigma }^{\xi }+{\sigma }_i^{\xi }H\left[(r_i+h)-r)\right]\right)}{r}^2\text{d}r=0.& (18)\end{array}$

Using the corresponding compositional buoyancy flow injected at the inner boundary,

$\begin{array}{ll}{f}_i^{\xi }=\beta g{I}_i^{\xi },& (19)\end{array}$

one can derive the following scale for composition:

$\begin{array}{ll}{C}^{*}=\frac{f_i^{\xi }}{4\pi D\beta {\rho }_{0}g\nu }.& (20)\end{array}$

The expression of the modified compositional Rayleigh number then becomes

$\begin{array}{ll}R{a}_{\xi }=\frac{f_i^{\xi }}{4\pi {\rho }_0{\nu }^2\Omega }.& (21)\end{array}$

FIGURE 1

Figure 1. Illustration of the shape of the source term ${\sigma }_i^{\xi }$ applied above the inner core boundary.

2.3. Numerical Method

Though moderate Schmidt numbers (< 50) may be numerically accessible with most classic dynamo codes, the infinite Schmidt number case requires a specific numerical method. Indeed, the time integration schemes classically used in the spherical dynamo codes (Crank-Nicolson for diffusion and Adams-Bashforth for the non-linear terms) are conditionally stable and produce spurious oscillations when the chemical diffusivity is set to zero. Other numerical methods such as TVD schemes have been specifically designed for hyperbolic equations but still introduce some amount of numerical diffusion and/or “clipping” effects (see Munz, 1988 for comparison of several TVD schemes). An interesting alternative is the “particle-in-cell” (PIC) method which has been extensively used for a wide range of physical problems (see Brackbill et al., 1988; Huang et al., 1992; Tackley and King, 2003; Tskhakaya et al., 2007; Sironi and Spitkovsky, 2009 for a non-exhaustive list of examples).

In this work, we use a recently developed version of the dynamo code PARODY-JA (Aubert et al., 2008) which includes a PIC method to treat the compositional field. Details on the numerical method can be found in Bouffard et al. (2017), together with the results of two benchmarks that validate the implementation of the PIC method. We performed a total of 11 numerical simulations summarized in Table 1 to conduct a first exploratory study of non-magnetic pure chemical convection in a rapidly-rotating spherical shell in the end-member case of infinite Schmidt number (κ^ξ = 0). For one simulation (case 9*), the homogeneous sink term σ^ξ in Equation (17) is replaced by a sink term near the top boundary equivalent to ${\sigma }_i^{\xi }$ that equilibrates the light elements injected at the bottom. As will be emphasized later, we limited ourselves to a small number of cases in this first preliminary study due to the high computational cost, particularly at low Ekman numbers. In order to decrease the computation time, we considered flows which are 2π/m_c-periodic along the longitude (see Table 1). In all cases, we ensured that the width of the computational domain remained large compared to the size of the individual chemical structures (see section 3), so that no physical property of chemical convection is lost by assuming these symmetries. The addition of two cases in full spherical shell (cases 1 and 2) for which no qualitative difference was observed with the other cases gives further support to this statement. Also, for reasons that will be explained in section 4.2, no simulation was run until a statistically stationary state was reached. As the mean values of diagnostic quantities (for instance, the kinetic energy) fluctuate throughout one simulation, we did not report output values such as Reynolds and Rossby numbers in Table 1. The Rossby numbers given in the figures throughout the manuscript correspond to instantaneous values associated with the snapshots displayed.

TABLE 1

Table 1. Summary of the parameters used for the different simulations.

3. Properties and Dynamics of Chemical Plumes and Blobs

3.1. Destabilization of the Bottom Light Layer

The continuous injection of light elements creates a lighter layer immediately above the bottom boundary. This configuration is prone to the Rayleigh-Taylor instability which, after some time, causes the layer to destabilize, forming thin rising filamentary chemical plumes with a “sheet-like” or “curtain” shape elongated in the direction of the rotation axis owing to the Taylor-Proudman constraint (see Figure 2A). Such structures are qualitatively in agreement with previous laboratory experiments by Cardin and Olson (1992) and are also typically observed for thermal plumes in simulations of rotating convection in a spherical shell (see, for instance, Figure 2 from Gastine et al., 2016). However, there seems to be a transition from “sheet-like” structures to more isolated plumes when crossing the tangent cylinder (see, for example, Figures 2A,B).

FIGURE 2

Figure 2. (A) Azimuthal section of the compositional field corresponding to case 8 shortly after the destabilization of the initial bottom light layer. To better visualize the curtain shape, an azimuthal average was performed over a smaller interval Δϕ = 0.01 around the longitude ϕ₀ centered on a chemical “curtain.” (B) Isosurface corresponding to a constant composition for case 3. The vertical dotted line shows the position and orientation of the rotation axis.

When they rise, the first chemical plumes carry away light elements and make the bottom light layer more heterogeneous. As a consequence, the following generations of compositional “sheets” assume a more complex and discontinuous structure since light elements are not always locally available to supply the rising curtains (see Figure 2B). Plumes may also merge and their location evolves over time.

We measured the mean thickness δ of these sheet-like or filamentary rising plumes for five different values of the Ekman number in the range 10⁻⁵ − 10⁻³. We estimated δ for the first generation of plumes by plotting zonal profiles of the composition in the equatorial plane, at a short distance above the bottom light layer. For each “peak” corresponding to a plume in the compositional profile, we define the width of the plume as the distance over which the composition is positive. This is facilitated by the fact that, at the beginning of the simulation, the background composition is homogeneous and negative due to the volumetric sink term σ^ξ. The results are displayed on Figure 3 and support a scaling law with a 1/3 exponent,

$\begin{array}{ll}\delta ~0.2E^{0.328},& (22)\end{array}$

which is the same exponent as the onset of classical thermal convection in a rotating spherical shell (Jones et al., 2000). In order to rule out any extra dependence on the Rayleigh number, we varied this parameter over two orders of magnitude while keeping a fixed Ekman number. No measurable difference was found for the average thickness of the chemical plumes. The thickness of the injection layer was also varied. The size of the chemical plumes was found to be affected as soon as the thickness of the injection layer becomes smaller than the size of the chemical plumes called for by Equation (22). The constant h in Equation (16) was therefore adjusted to be slightly larger than the predicted size of the plumes, but kept as small as possible in order to model local processes occurring close to the bottom boundary.

FIGURE 3

Figure 3. Mean thickness of the chemical plumes measured in the equatorial plane for different values of the Ekman number.

The exact value of the pre-factor in the scaling law (22) may be slightly sensitive to the parameters σ and h characterizing the injection layer in Equation (16) and depends on the definition of δ but remains small anyway compared to the typical pre-factors for thermal convection. Consequently, the size of chemical plumes is always a very small fraction of the shell's thickness. This has an important numerical implication for the grid (and tracers) resolution which must be very high to correctly resolve chemical plumes. This requirement is visible on the chemical spectrum shown on Figure 4 left, which is remarkably flat before dropping abruptly. The small size of the chemical structures also produces smaller scales in the velocity field. For comparable Rayleigh numbers, the kinetic energy spectra contain significantly less energy at low degrees and drop much more slowly for compositional convection compared to thermal (codensity) convection (Figure 4B). However, in the simulations listed in Table 1, resolving the chemical structures sets a stronger constraint on the resolution than the kinetic energy spectra. The 1/3 exponent in the scaling law (22) implies that the resolution should be approximately doubled for each decade dropped for the Ekman number, which justifies the high spherical harmonic degrees employed in the simulations (see Table 1). At fixed Ekman number, the numerical demand of compositional convection is significantly higher than an equivalent codensity run. We therefore restricted ourselves to a small number of simulations in this exploratory study and considered longitudinal symmetries to decrease the computation time (see Table 1). Note that the fraction of the computation time spent in the PIC method is dominant but decreases when lower Ekman numbers are employed.

FIGURE 4

Figure 4. (A) m spectrum of the chemical field for case 3. (B) kinetic energy spectra corresponding to the benchmark proposed by Breuer et al. (2010) and extended to the infinite Lewis number case (see Bouffard et al., 2017). δ_ξ denotes the fraction of chemical forcing.

3.2. Dynamics of Chemical Plumes

3.2.1. Effect of the Taylor-Proudman Constraint on the Ascension of Plumes

When rotational effects remain moderate, chemical plumes rise almost radially (see Figure 5A). When the rotational constraint is increased and therefore the Rossby number Ro = U/DΩ (with U the root-mean-squared velocity) is lowered, the tangent cylinder (TC), an imaginary cylinder parallel to the rotation axis that touches the inner core at the equator, becomes less penetrable (Aurnou et al., 2003; Jones, 2015). The velocity field outside the TC is mostly invariant along the rotation axis by assuming a columnar form (see Figure 6). The light elements released at mid to high latitudes tend to stay inside the TC and rise more vertically along the z coordinate (see Figures 5A–C, 7B). The light elements released in the vicinity where the TC touches the inner boundary still have a chance to cross this imaginary boundary. As a result, the mass fraction of light elements is significantly higher inside the TC (which only represents 15% of the volume of the spherical shell) than outside. This is visible on Figures 5A–C and we checked that the mass fraction of light elements inside (resp. outside) the TC is systematically larger (smaller) than the mean composition in the entire shell.

FIGURE 5

Figure 5. Azimuthal average of the chemical field for cases 3 (t = 1; A), 6 (t = 1; B), and 8 (t = 0.4; C). Rossby numbers are, respectively: Ro = 9.4 × 10⁻³, Ro = 1.9 × 10⁻³ and Ro = 4.5 × 10⁻⁴. The color scales have been deliberately modified and saturated in order to better visualize the structure of the compositional field in the convecting region.

FIGURE 6

Figure 6. Azimuthal section of the velocity perpendicular to (Vs) and parallel to (Vz) the rotation axis, for case 8 at time t = 0.55.

FIGURE 7

Figure 7. (A) Equatorial section of the compositional field in case 10 at time t = 0.066 and Ro = 4.4 × 10⁻⁴. (B) Azimuthal average of the chemical field corresponding to the same snapshot. In both figures, color scales have been slightly saturated to improve the contrast of the compositional structures.

3.2.2. Secondary “Blobs” Instabilities

Due to the absence of chemical diffusion, rising chemical plumes can be affected by secondary instabilities, some of which are usually not observed for thermal plumes (Kerr and Mériaux, 2004; Kerr et al., 2008). These instabilities may have distinct origins in the polar and equatorial portions of the shell since rotation acts differently in these regions.

Several laboratory experiments have been dedicated to the study of the dynamics of buoyant plumes in a rotating environment. These were conducted in rotating cylinders in which the rotation axis is parallel to the gravity vector, a situation that is analogous to polar regions in the Earth's outer core. When discharging a single point turbulent plume at the top of a water-filled rotating tank, Frank et al. (2017) reported a systematic anticyclonic precession of the plume in the whole range spanned for the Rossby number $Ro={(B_0{\Omega }^{-3})}^{1/4}{H}^{-1}=0.02-1.3$, where H is the height of water tank and B₀ the source buoyancy flux. The vertical extent of the precession region and the time necessary for the development of precession were found to decrease with the Rossby number. Helical motions were also reported for turbulent thermal plumes in laboratory experiments of rotating convection modeling the tangent cylinder region (Aurnou et al., 2003). In another series of experiments, in which laminar chemical plumes are produced by the crystallization of a mush at the bottom of a rotating cylinder cooled from below, Claßen et al. (1999) observed three-dimensional helical undulations of the chemical plumes, both in the rotating and non-rotating cases. In this experiment, increasing rotation exacerbates these undulations and deflects the chemical plumes horizontally. Past a critical angle, the plumes become susceptible to a Rayleigh-Taylor instability which causes chemical “blobs” to detach from the initial plumes and rise vertically. This so-called “blob-instability” was found to occur closer to the mush when the Ekman number decreases. It is not clear whether the phenomena reported by Claßen et al. (1999), Aurnou et al. (2003), and Frank et al. (2017) for the dynamics of chemical plumes have a common physical origin or result from distinct instabilities.

We observed similar behaviors for chemical plumes in our simulations. On Figure 5B with E = 3 × 10⁻⁴ and Ro = 1.9 × 10⁻³, the pronounced plume at the north pole shows some helical undulation before a chemical blob detaches and continues to rise, a phenomenon similar to that reported by Claßen et al. (1999). This phenomenon is less easily visible for case 10 at a lower Ekman number E = 3 × 10⁻⁵ and Ro = 4.4 × 10⁻⁴ (see Figure 7). In this case, blobs still detach from the initial plumes but undulations and horizontal bending are often not as manifest. However, the plume denoted by an arrow on Figure 7 experiences a helical motion within a short distance from the plume's source which bears some resemblance with the precessions and undulations reported, respectively by Aurnou et al. (2003) and Frank et al. (2017) (for a comparison, see respectively Figures 2, 6 in these papers). The local Rossby number Ro_δ ~ U/δΩ computed at the scale of the plume δ based on Equation (22) is close to 0.07 in the case corresponding to Figure 7 and falls in the range explored by Frank et al. (2017), though the authors use a slightly different definition of the Rossby number.

We also observed other “blob-instabilities” similar to that reported by Claßen et al. (1999) in the equatorial plane, though the underlying physical mechanisms seem to contrast with the polar region. An illustration of a “blob-instability” in the equatorial plane is displayed on Figure 8 where the plume indicated by an arrow bends horizontally before a secondary “blob-like” plume detaches and continues to rise. The bending occurs closer to the boundary when the Ekman number is decreased. However, a fraction of the plumes bend in a direction that is not compatible with a deflection by the Coriolis force. Furthermore, it typically takes a few tens of rotation times (for instance about 20 in the case corresponding to Figure 9B) before the plumes bend and break-up, which suggests that the bending is not primarily caused by the Coriolis force. In addition, blobs may still detach from the plumes even in the absence of bending (see, for instance, Figure 9A). Different instabilities may therefore affect the dynamics of chemical plumes and should be clarified in the future. We did not observe any obvious and significant elongation of chemical blobs in the direction of the rotation axis as predicted by St Pierre (1996), but this may be because monitoring the shape of chemical blobs over time is quite arduous. However, in the long term, the fragmentation and stirring of the chemical plumes and blobs combined with the absence of chemical diffusion produces a granular “flaky” aspect in the compositional field (see, for example, Figure 9A), a feature that contrasts with the smoothness of the temperature (codensity) field in low Prandtl number simulations.

FIGURE 8

Figure 8. Equatorial snapshot of the compositional field for case 9*.

FIGURE 9

Figure 9. (A) Equatorial section of the compositional field at t = 6 for case 1. The color scale has been saturated so that the top stratified layer becomes more apparent. The white dashed lines delimit qualitatively three regions in the spherical shell: a bottom plumes region, an intermediate more or less well mixed blobs region and a stratified layer at the top. (B) Equatorial section of the radial velocity field (in color) for the same snapshot. Arrows denote local horizontal velocity vectors.

3.2.3. Jet-Instability and Mixing of the Chemical Plumes

Though plumes have a rather laminar behavior in our simulations (which operate at low Reynolds numbers), Figure 8 shows that some laminar chemical plumes evolve to a more “turbulent” state characterized by entrainment of ambient fluid that causes the gradual dilution of the plumes. This is also visible for a lower Ekman number on Figure 7A, which shows very thin chemical plumes transforming into large blobs that reach the top boundary. This transition to a more “turbulent” state is qualitatively similar to the so-called “jet” instability, in which a straight laminar flow destabilizes after some time and leads to a more complex jet flow (Michalke, 1984). The initial jet can fragment into smaller drops or break-up into blobs with a size comparable to the width of the initial jet. Eddies may also form at the interface and start entraining ambient fluid, causing the jet to widen. Figure 10 shows that there is a qualitative agreement between the dynamics of some plumes from our simulations and laboratory experiments of the jet instability. We observe that, at the beginning of the simulations, when the background composition is homogeneous, the first plumes quickly start entraining ambient fluid. In the later generations of plumes however, when the background composition is more heterogeneous, a significant fraction of plumes bend and break-up into blobs before they start entraining surrounding fluid. This suggests that the background composition plays a role in the dynamics of plumes, though its precise influence remains to be investigated.

FIGURE 10

Figure 10. (Left) Examples of chemical plumes for case 9* after the initial bottom light layer starts destabilizing at the beginning of the simulation. (Right) Examples of some regimes of jet instability in laboratory experiments. Taken from Tang et al. (2003).

The efficiency of plumes mixing caused by the entrainment of heavier fluid differs between the rotating and non-rotating cases. Many studies have pointed out that turbulent mixing and lateral entrainment of fluid are inhibited in rotating flows (Helfrich and Battisti, 1991; Fernando et al., 1998; Levy and Fernando, 2002; Goodman et al., 2004; Gastine et al., 2016). Mixing processes may also be different inside and outside the TC. In our simulations, lateral entrainment of ambient fluid seems partially inhibited in polar regions where the plumes and blobs are little diluted during their vertical ascension. In contrast, in the equatorial plane, plumes experience more significant stirring, though lateral entrainment may still be partially inhibited compared to a non-rotating case. The different mixing efficiency of chemical plumes between polar and equatorial regions is visible on Figures 5B,C. A more detailed and quantitative study involving more turbulent flows will be necessary to better understand the mixing dynamics of chemical plumes.

4. Formation of a Chemically Stratified Layer at the Top of the Shell

4.1. Accumulation of Chemical Plumes and Blobs at the Top of the Spherical Shell

In all our simulations, a fraction of the rising chemical plumes and blobs reach the top of the spherical shell where they accumulate to form a stably stratified layer. This layer is clearly visible on equatorial sections of the chemical field (Figure 9A), azimuthal projections of the ϕ-average compositional field (Figure 5), and on radial profiles of the mean composition (Figure 11A) which all show that the mass fraction of light elements increases radially in the top portion of the shell. Defining the extent of the stably stratified layer is a delicate task. One possibility is to define the bottom of the stratified region as the radius r_b above which the mean Brunt-Väisälä frequency N is positive:

$\begin{array}{ll}{N}^2=-\frac{g}{\rho }\frac{\partial \overline{\rho }}{\partial r}>0\text{for}r>r_b,& (23)\end{array}$

where $\overline{\rho }(r)$ denotes the average density at radius r . For pure compositional convection, Equation (23) is equivalent to

$\begin{array}{ll}\frac{\partial \overline{\xi }}{\partial r}>0\text{for}r>r_b,& (24)\end{array}$

in which $\overline{\xi }(r)$ is the mean composition at radius r. The bottom of the stratified layer r_b, defined by Equation (24), is indicated by dashed lines on Figures 11A,B, which show radial profiles of the mean composition and (N/Ω)² ratio, respectively. Similar shapes are observed for these profiles in the other simulations. Figure 11B illustrates that the Brunt-Väisälä frequency increases very rapidly in the top half portion of the shell where the N/Ω ratio reaches a value close to 0.7 in this case. In contrast, the bottom half of the layer is weakly stratified with N/Ω < 0.1. A first disadvantage of the definition (23) is that the weak slope of N in the lower part of the stable layer causes strong fluctuations in the value of r_b which require time averaging. In addition, the definition of r_b remains somehow arbitrary since there is dynamically no neat boundary delimiting the convective and stratified regions. Indeed, the low values of the N/Ω ratio across the layer do not suppress radial movements caused by penetrating plumes and blobs and allow for the penetration of columnar convection as expected from the linear theory (Takehiro and Lister, 2001) (see, for instance, Figure 6). Simpler definitions of the layer based, for example, on the intersection of tangents, are also arbitrary and difficult to implement in practice due to the complexity of the shape of the chemical profiles. In fact, these profiles are self-similar and do not seem to fit any classical mathematical function. Determining a relevant criterion for the definition of the top light layer will be the object of future investigations.

FIGURE 11

Figure 11. (A) Radial profile of the mean composition at different times (in viscous time units) for case 1. The dashed lines indicate the bottom of the layer as defined by Equation (24). (B) (N/Ω)² ratio at time t = 6 for the same case.

Figure 9 suggests that the spherical shell can be roughly divided into three regions as anticipated by Loper (2007): a bottom part dominated by plumes, a more or less “well-mixed” intermediate blobs region and a stably stratified layer at the top. The height of the plumes region is also imprecise and fluctuates spatially due to different plume dynamics inside and outside the TC. Moreover, as we will explain below, the relative extent of these regions evolves with time and depends on the Ekman number: the stratified layer thickens over time, but shrinks when the Ekman number is decreased.

4.2. Dynamics and Evolution of the Layer

4.2.1. Latitudinal Topography

Figures 5A–C show that the thickness of the stably stratified layer is lower at the equator and increases toward the poles, creating a latitudinal topography. We attribute this to a direct consequence of the tendency of high latitudes plumes to rise almost vertically along the rotation axis and remain inside the TC, which results in a differential growth rate of the layer between the low and high latitudes.

Furthermore, the mass fraction of light elements at the very top of the stratified layer also increases from the equator to the poles. We make the hypothesis that this is due to the different mixing efficiency of the plumes inside and outside the TC. In polar regions, plumes rise almost vertically and experience less mixing and dilution than in the equatorial plane. Consequently, they reach the stratified region with a higher mass anomaly than equatorial plumes which experience more fragmentation and mixing.

4.2.2. Time Evolution

After a first phase of rapid growth of the layer caused by the destabilization of the initial bottom light layer, the thickness of the layer [based on definition (24)] increases very slowly. Figure 11 shows that the stratified layer extends to about 35% of the spherical shell's thickness after 2 viscous times, whereas it takes 6 viscous times for the layer to fill 45% of the shell. Several tens of viscous times are in fact necessary to reach a statistically steady-state in the simulations. The resulting computational demand is enormous, especially for low Ekman numbers, and we ran none of our simulations until the statistically steady-state was reached in the present study. Even though running the simulations for very long times is probably not necessary (since the Earth is not in a statistically steady state), fitting the time evolution of the layer's thickness with a particular mathematical function may still require costly simulations. Comparison of the different runs suggests that the growth rate is dependent on the Ekman number and, to a lesser extent, the Rayleigh number. For instance, after one viscous time, the stratified layer extends to about 25% of the shell for case 1 with E = 10⁻³ and only 20% for case 6 with E = 3 × 10⁻⁴ and comparable Rayleigh numbers (see Figures 5A,B). Deriving scaling laws will require a proper definition of the layer and a more systematic (and cumbersome) exploration of the parameter space that we will perform in the future.

4.3. Analogy With the “Filling-Box” Models

This situation bears some interesting resemblance to the “filling-box” models studied experimentally by Baines and Turner (1969). In these experiments, a buoyant turbulent plume is injected at the bottom of a box of finite volume. As it moves upwards, the first plume entrains ambient fluid that is then mixed internally by the efficient turbulence of the plume. When it reaches the top boundary, the plume spreads horizontally and forms a homogenous light layer at the top of the box separated from the rest of the fluid by a density jump. When the second rising plume enters this light region, it entrains lighter fluid and therefore reaches the top of the box even lighter than the first plume. It then spreads horizontally, pushing downwards the first light layer. The process repeats, gradually forming a stably stratified layer at the top of the box that slowly grows until it fills the entire box (see Figure 1 of the paper by Worster and Huppert, 1983 for a schematic illustration).

In the idealized “filling-box” model described above, parametrization of the turbulent plume allows to derive analytically the time evolution of the density profile and the position of the first front (Baines and Turner, 1969). However, a series of differences prevents a direct extrapolation of these results to our simulations of chemical convection. In addition to the spherical geometry and a higher number of plumes, two major specificities of our simulations are the presence of rotation and a different plumes dynamics (undulations, fragmentation into blobs, jet instability, inhibition of lateral mixing…). This makes the classic parametrization of “turbulent plumes” dubious in this context. The filling-box models also assume the total absence of turbulence and mixing in the stratified region which may not be the case for our simulations in which no clear separation in the dynamics can be drawn between the bottom of the layer and the convecting region.

The compositional profiles shown on Figure 11 also contrast with the density profiles in the “filling-box” models notably by the absence of a density jump at the bottom, and by the profile of the Brunt-Väisälä frequency which increases with the radial position whereas it decreases with height in the filling-box models. This confirms that different mechanisms occur in our simulations and that the adaptation of the original filling-box model to the case of rotating chemical convection in a spherical shell requires extra investigations.

5. Implications for the Earth's Core

5.1. Limitations of the Present Study

Determining the implication of the results shown above for the Earth is a delicate task that is impeded by several factors. First, assessing the “plausibility” of a given scenario for the formation of a stably stratified layer below the CMB is complicated by the apparent inconsistency between studies based on magnetic and seismic observations (see Table 2 for a summary). Buffett et al. (2016) show that a rather thin (< 140 km) and weakly stratified (N ≲ 5.4 − 6.1 × 10⁻⁵ s⁻¹) layer can explain both the 60-year period in the geomagnetic secular variations together with fluctuations in the dipole field and, possibly, length of the day. Dynamo simulations suggest that the morphology of the magnetic field becomes too little Earth-like when the stratified layer is thicker than 400 km (Olson et al., 2017; Christensen, 2018). Analysis of the appearance of inverse flux patches in the southern hemisphere suggests an upper bound of only 100 km (Gubbins, 2007). Another study of geomagnetic secular variations remains unconclusive (Lesur et al., 2015). On the side of seismology, some studies based on the analysis of SmKS waves traveling in the outermost core are in favor of a thick layer (up to 450 km) with a strong stratification (N ~ 10⁻³ s⁻¹) (Helffrich and Kaneshima, 2010; Kaneshima and Matsuzawa, 2015; Kaneshima, 2018). However, two seismic studies found no evidence for a stratification (Alexandrakis and Eaton, 2010; Irving et al., 2018), though this may simply indicate that the layer is too weakly stratified to leave a seismic signature. Second, the limited set of simulations presented in this study does not allow for the derivation of scaling laws, which prevents any direct extrapolation to the Earth's core. Third, several physical ingredients are not considered in the present study such as the magnetic field, thermal convection, and turbulence.

TABLE 2

Table 2. List of the main studies that proposed some constraints on the stratified layer at the top of the Earth's core based on either magnetic or seismic observations.

The presence of a magnetic field can influence the dynamics of buoyant parcels of fluid in a non-trivial way (Moffatt and Loper, 1994; Sreenivasan and Narasimhan, 2017) and should be incorporated in future studies. One may also argue that adding the destabilizing effect of thermal convection (omitted in this study) could destroy the chemically stratified layer. Although some preliminary simulations of thermochemical convection (not shown in this paper) suggest this is not the case—at least when both fields contribute equally to the buoyancy flux—the dynamics of thermochemical convection remains to be explored. Note that thermal effects may not necessarily be destabilizing. The recently revised value of the thermal conductivity (de Koker et al., 2012; Pozzo et al., 2012; Gomi et al., 2013; Ohta et al., 2016) seems to imply that the top of the core is thermally stratified (Labrosse, 2015). The presence of a thermally stratified layer at the top of the core may ease the formation of chemical stratification as well. The present study also ignores any possible double-diffusive phenomena that may affect the dynamics of a stably stratified layer in various ways.

On the other hand, the method of light elements injection in our simulations may not adequately reflect the complex chemical processes occurring in the inner core boundary region. Indeed, injecting light elements via a homogeneous source term above the bottom boundary tends to generate plumes with a sheet-like structure instead of the cylindrical shape that would be expected if plumes erupt from isolated chimneys. A step forward toward a more realistic modeling of the mushy layer could consist of injecting light elements within randomly distributed “patches” mimicking local chimneys, although this brings extra complexities since the number of mush chimneys in the ICB conditions is not constrained. Furthermore, the scaling law (22) may not be relevant to predict the size of chemical plumes in core conditions if these plumes erupt from localized chimneys formed in the mushy layer. The initial width of the plumes would in that case likely be controlled by the dynamics of the mushy layer rather than by the Rayleigh-Taylor destabilization of a homogeneous light layer. Unfortunately, the size of the chimneys formed in the mush is presently unknown in core conditions. An alternative scenario is the presence of a slurry layer above the inner core boundary (Loper and Roberts, 1981; Gubbins et al., 2008; Wong et al., 2018) which may not be adequately modeled by a simple compositional volumetric source term. Also, note that the term “compositional convection” referred to in this study could be subject to a semantic debate: Rayleigh-Taylor instability and Rayleigh-Bénard convection are both buoyancy-driven instabilities but the former is for two layers of different densities while the latter seems to imply a more continuous density variation associated with diffusion. It is not clear whether the nature of the first instability affects what we observe in the shell away from the boundary.

Another crucial question is whether chemical plumes would experience a complete turbulent mixing in core conditions. Laboratory experiments by Jellinek et al. (1999), in which a buoyant fluid is released at the bottom of a non-rotating box, show that plumes experience efficient turbulent mixing as soon as their associated plumes Reynolds number exceeds 100, which can prevent the formation of a stratified layer at the top of the tank. In our simulations, the local Reynolds number at the scale δ of the chemical plumes yields

$\begin{array}{ll}R{e}_{\delta }~U\delta /\nu .& (25)\end{array}$

Assuming the scaling law (22) can be used as a good proxy for the size of chemical plumes in core conditions, we can estimate

$\begin{array}{ll}R{e}_{\delta }~0.2Re{E}^{1/3}.& (26)\end{array}$

Taking E = 10⁻¹⁵ and Re = 10⁸ for the core, one obtains Re_δ ~ 200, which suggests that plumes would be totally mixed during their ascension in the core, a phenomenon that cannot be captured in our simulations that operate at low global-scale Reynolds numbers (< 100) and therefore, even lower plumes' Reynolds numbers. However, as rotation is known to partially inhibit lateral entrainment of fluid and subsequent mixing (Helfrich and Battisti, 1991; Fernando et al., 1998; Levy and Fernando, 2002; Goodman et al., 2004; Gastine et al., 2016), it is not clear whether chemical plumes would experience significant mixing and dilution in core conditions. The local Rossby number based on the plumes length scale called for by Equation (22) writes

$\begin{array}{ll}R{o}_{\delta }~\frac{U}{\delta \Omega }~5Ro{E}^{-1/3}.& (27)\end{array}$

Taking Ro = U/DΩ ~ 10⁻⁶ in core conditions, we get Ro_δ ~ 0.5 which suggests that rotational effects may still be influencial for chemical plumes. Future investigations will be necessary to quantify the mixing of chemical plumes and the role played by rotation in more turbulent flow regimes.

One may also question the relevance of the infinite Schmidt number approximation used in the simulations. The dynamics of a rising chemical plume is susceptible to be affected by chemical diffusion if the latter acts at a length scale comparable with the size of the plume (or larger) during its typical rising time. Using the scaling law (22) for the average width of the plumes, one can estimate the typical time τ necessary for chemical diffusion to operate at the scale of a chemical plume or blob of size δ,

$\begin{array}{ll}\tau ~\frac{{\delta }^2}{{\kappa }^{\xi }}~0.04Sc{E}^{2/3}\times {\tau }_{visc},& (28)\end{array}$

where τ_visc denotes the viscous time D²/ν. For example, E = 10⁻⁴ and Sc ~ 100 give τ ~ 0.01τ_visc. This time must be compared with the typical time it takes for a chemical plume or blob to reach the top boundary, which likely depends on the convective regime. For instance, in case 8 with E = 10⁻⁴ and Sc = +∞, chemical plumes typically reach the top stratified region within 0.1τ_visc. For a finite value of Sc, the typical diffusion length ℓ during this time period would be

$\begin{array}{ll}\ell ~D\sqrt{0.1/Sc}~D/30~3\delta .& (29)\end{array}$

For intermediate Schmidt numbers, chemical diffusion would therefore not be negligible for the dynamics of chemical plumes. A higher value Sc = 1, 000, corresponding to the upper part of the plausible range expected for the Schmidt number (cf. Equation 1), would give $\ell ~D\sqrt{0.1/Sc}~D/100~\delta$, a more moderate but still non negligible value. Chemical diffusion could also be a significant contributor to the long-term dynamics of the chemically stratified layer, since the latter builds over long time scales in the simulations. Diffusion would probably “smooth” the granular “flaky” aspect of the compositional field (visible for instance on Figure 9A) over time and cause the stratified layer to grow downwards, like in the models of Buffett and Seagle (2010) and Gubbins and Davies (2013). In the case of turbulent flows, chemical diffusion may also act at small scales and allow compositional mixing. Only stirring occurs when chemical diffusion is neglected. The infinite Schmidt number approximation may thus be debatable for the modeling of chemical convection, but this end-member case is certainly a better approximation than the codensity approach in the sense that it can capture physical aspects that are specific to the weakly-diffusive compositional field.

5.2. Physical Constraints on the Properties of a Chemically Stratified Layer Formed by Accumulation of Chemical Blobs

Despite the numerous limitations of the present study listed above that prevent any direct extrapolation of our results to the Earth's core, some physical considerations can be invoked to constrain the properties of a chemically stratified layer built by the accumulation of light chemical blobs emitted at the ICB. In such a scenario, the total density anomaly Δρ_L across the stably stratified layer is bounded by that Δρ of the chemical blobs erupting from the mushy layer (Δρ_L ≤ Δρ). The density anomaly Δρ of such chemical blobs is not known and may have varied throughout the history of the inner core. Here we try to obtain an order of magnitude for this density anomaly using very simple hypotheses and calculations. Following Moffatt (1989) (but see Moffatt and Loper, 1994 for a more elaborate study), we assume that the dynamics of erupting chemical blobs is governed by a force balance between the Coriolis and buoyancy forces. Using the approximate relation ρ_s ~ 1.05ρ, where ρ_s is the density of the inner core and ρ that of the bulk fluid core, one can write

$\begin{array}{ll}{f}^{1/2}\frac{\Delta \rho }{\rho }~{\left(\frac{\Omega Ṙ}{20g}\right)}^{1/2};f^{1/2}w~\left(\frac{gṘ}{20\Omega }\right),& (30)\end{array}$

where f denotes the surface fraction of the mushy layer occupied by chimneys, w the upward velocity of the chemical blobs and Ṙ the growh rate of the inner core. Taking Ṙ ~ 10⁻¹¹ m s⁻¹ (Labrosse, 2015), Ω ~ 7 × 10⁻⁵ s⁻¹ and g ~ 3 m s⁻² leads to

$\begin{array}{ll}{f}^{1/2}\frac{\Delta \rho }{\rho }~3\times 10^{-9};f^{1/2}w~2\times 10^{-4}{\text{m s}}^{-1}.& (31)\end{array}$

The parameter f was found to be close to 10⁻³ − 10⁻² in some laboratory experiments (Copley et al., 1970), but is unknown for the mushy layer in core conditions. The range of acceptable values is however constrained by the typical velocity expected for the fluid core. The latter can be inferred from the inversion of magnetic observations and was estimated to be of order 4 × 10⁻⁴ m s⁻¹ below the CMB (Finlay and Amit, 2011). Estimations from numerical simulations predict a slightly higher root-mean-square (RMS) velocity of order 10⁻³ m s⁻¹ (Christensen and Aubert, 2006). Since velocities have a statistical distribution, the velocities of erupting blobs may be significantly larger than the RMS velocity in the convecting region. Assuming a typical velocity 10⁻² m s⁻¹ for the chemical blobs gives

$\begin{array}{ll}f~4\times 10^{-4};\frac{\Delta \rho }{\rho }~1.5\times 10^{-7}.& (32)\end{array}$

Figure 12 shows the mean Brunt-Väisälä frequency N obtained when spreading linearly a density anomaly Δρ ~ 10⁻⁷ρ over a stratified layer of thickness H,

$\begin{array}{ll}N~\sqrt{\frac{g}{H}\frac{\Delta \rho }{\rho }}.& (33)\end{array}$

For instance, g ~ 10 m s⁻² and H = 100 km would give N ~ 3 × 10⁻⁶ s⁻¹, which is just slightly smaller than the estimation required by Buffett (2014) and Buffett et al. (2016) to explain some properties of the geomagnetic field, and more than two orders of magnitude smaller than the value suggested by Helffrich and Kaneshima (2010) based on seismic observations. Obtaining a mean buoyancy frequency of order N ~ 10⁻⁵ s⁻¹ compatible with magnetic observations would require either a thinner stratified layer (H ~ 10 km) or a higher density anomaly such that Δρ/ρ ~ 10⁻⁶. This density ratio ratio would then give f ~ 10⁻⁵ and w ~ 0.1 m s⁻¹. It should be noted that the mean Brunt–Väisälä frequency of the layer estimated using Equation (33) constitutes an upper bound corresponding to the ideal case in which there is no dilution of the chemical blobs along their path to the top of the core (Δρ_L ~ Δρ). Provided the force balance proposed by Moffatt (1989) is reasonable, it is therefore clear that matching the strong stratification (N ≲ 10⁻³ s⁻¹) called for by some seismic studies can be achieved only by invoking unrealistically high blobs velocities (w ≳ 10² m s⁻¹). One may however criticize the simplicity of this assumption by arguing that the properties of the blobs erupting from the mush chimneys would rather be controlled by the dynamics of the mushy layer. Though it seems more likely that a stratified layer at the top of the core formed by the accumulation of chemical blobs emitted at the ICB would be weakly stratified (most likely N < 10⁻⁵ s⁻¹) and, very likely, seismically invisible, we cannot exclude stronger stratifications depending on the precise processes occurring in the ICB region.

FIGURE 12

Figure 12. Mean Brunt-Väisälä frequency (orange) and N/Ω ratio (blue) obtained when spreading linearly a density difference Δρ across a stratified layer of given thickness (x-axis), for Δρ/ρ = 10⁻⁷.

Should the accumulation of chemical blobs alone be unable to explain the presence of a thick stably stratified layer at the top of the core, it might still play an important role. Gubbins and Davies (2013) studied the possible formation of a chemically stratified layer by barodiffusion. Though their model produces a layer that may be compatible with observations, the authors mention that a crucial prerequisite in their scenario is the presence of an initial “embryo” of layer that is not mixed by convection. The latter could well be produced by the accumulation of incompletely mixed chemical blobs emitted at the inner core boundary. Though the absence of scaling laws prevents us from making a prediction on the thickness of such a layer, the stratified layer rapidly extends well beyond the Ekman layer in all our simulations. We therefore suggest that the accumulation of incompletely mixed chemical plumes may form a stratified layer of sufficient thickness to allow a thicker layer to grow according to the scenario of Gubbins and Davies (2013).

6. Discussion and Conclusion

Using a new version of the PARODY code including a particle-in-cell method, we ran a first series of exploratory simulations of pure chemical convection at infinite Schmidt number in a rotating spherical shell. In these simulations, light elements are injected within a thin layer above the bottom boundary. This preliminary work remains descriptive and qualitative, but reveals that chemical convection in the Earth's outer core is a rich and vast topic that warrants several more focused and quantitative studies. We showed that, principally as a result of their lower diffusivity, chemical plumes have a complex dynamics that may be influenced by several distinct instabilities in the different regions of the shell. We observed secondary instabilities such as undulations, plumes precession, fragmentation into chemical “blobs” that comply qualitatively well with some previous predictions and laboratory experiments. Plumes may also experience jet-instability and start widening and entraining ambient fluid. A particularly interesting result is the systematic formation of a chemically stratified layer at the top of the shell in all simulations. The layer has a latitudinal topography and slowly grows owing to the accumulation of incompletely mixed chemical plumes and blobs emitted at the bottom boundary.

The existence of a stably stratified layer at the top of the core has been proposed for several decades based on seismic and magnetic observations (Lay and Young, 1990; Tanaka, 2007; Helffrich and Kaneshima, 2010; Buffett, 2014; Kaneshima, 2018) as well as core evolution models (Labrosse, 2015). The debate concerning the existence and properties of this layer (thermal and/or chemical stratification, thickness, and strength of stratification) remains open and our exploratory study does not bring definitive answers to any of these points. However, to our knowledge, this is the first time that the self-consistent formation of a stratified layer is observed in numerical simulations. The layer forms through mechanisms that bear some similarity with the scenario envisioned by Moffatt and Loper (1994) and the “filling-box” models (Baines and Turner, 1969).

As our simulations operate with control parameters that are very remote from core conditions and miss some important ingredients like thermal convection, magnetic field and turbulence, it is not clear whether the results shown in this study would hold when approaching more realistic conditions. Future studies should aim at incorporating these essential ingredients and systematically exploring the parameter space in order to understand more quantitatively the physics of compositional convection and make predictions for the Earth. Our preliminary study suggests that a chemically stratified layer formed by the accumulation of chemically buoyant parcels of fluid released at the bottom of the core would probably be weakly stratified and thus, likely seismically invisible. However, this relies on a force balance assumption that may be irrelevant. Constraints on the stratification that can be produced by such a scenario may rather come from the precise dynamics of the mushy layer which is presently poorly understood and deserves future investigations.

Author Contributions

MB ran the simulations, analyzed the data, and wrote the manuscript. GC, SL, and JW significantly contributed to the physical reflection, the interpretation of the results, and to the elaboration of the paper.

Funding

This study is partly funded by the French Agence Nationale de la Recherche (grant 446 number ANR-15-CE31-0018-01, MaCoMaOc). The simulations were performed on the ADA machine of the IDRIS center (Orsay, Paris), with grant number i2016047436 and on the Hydra cluster of the Max Planck Society (Garching, Germany).

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank David Gubbins, Binod Sreenivasan, and Kenneth P. Kodama for their valuable comments which were very helpful to improve the manuscript. We also thank Philippe Cardin and Hagay Amit for fruitful discussions on this topic.

References

Alexandrakis, C., and Eaton, D. W. (2010). Precise seismic-wave velocity atop Earth's core: No evidence for outer-core stratification. Phys. Earth Planet. Inter. 180, 59–65. doi: 10.1016/j.pepi.2010.02.011

CrossRef Full Text | Google Scholar

Aubert, J., Aurnou, J., and Wicht, J. (2008). The magnetic structure of convection-driven numerical dynamos. Geophys. J. Int. 172, 945–956. doi: 10.1111/j.1365-246X.2007.03693.x

CrossRef Full Text | Google Scholar

Aurnou, J., Andreadis, S., Zhu, L., and Olson, P. (2003). Experiments on convection in Earth's core tangent cylinder. Earth Planet. Sci. Lett. 212, 119–134. doi: 10.1016/S0012-821X(03)00237-1

CrossRef Full Text | Google Scholar

Badro, J., Siebert, J., and Nimmo, F. (2016). An early geodynamo driven by exsolution of mantle components from Earth's core. Nature 536:326. doi: 10.1038/nature18594

PubMed Abstract | CrossRef Full Text | Google Scholar

Baines, W., and Turner, J. (1969). Turbulent buoyant convection from a source in a confined region. J. Fluid Mech. 37, 51–80. doi: 10.1017/S0022112069000413

CrossRef Full Text | Google Scholar

Bergman, M. I., and Fearn, D. R. (1994). Chimneys on the Earth's inner-outer core boundary? Geophys. Res. Lett. 21, 477–480. doi: 10.1029/94GL00028

CrossRef Full Text | Google Scholar

Bouffard, M., Labrosse, S., Choblet, G., Fournier, A., Aubert, J., and Tackley, P. J. (2017). A particle-in-cell method for studying double-diffusive convection in the liquid layers of planetary interiors. J. Comput. Phys. 346, 552–571. doi: 10.1016/j.jcp.2017.06.028

CrossRef Full Text | Google Scholar

Brackbill, J. U., Kothe, D. B., and Ruppel, H. M. (1988). FLIP: a low-dissipation, particle-in-cell method for fluid flow. Comput. Phys. Commun. 48, 25–38. doi: 10.1016/0010-4655(88)90020-3

CrossRef Full Text | Google Scholar

Braginsky, S. (1963). Structure of the F layer and reasons for convection in the Earth's core. Soviet Phys. Dokl. 149, 8–10.

Google Scholar

Braginsky, S. (1994). MAC-oscillations of the hidden ocean of the core. J. Geomagn. Geoelectr. 45, 1517–1538. doi: 10.5636/jgg.45.1517

CrossRef Full Text | Google Scholar

Braginsky, S. I., and Roberts, P. H. (1995). Equations governing convection in Earth's core and the geodynamo. Geophys. Astrophys. Fluid Dyn. 79, 1–97. doi: 10.1080/03091929508228992

CrossRef Full Text | Google Scholar

Breuer, M., Manglik, A., Wicht, J., Trümper, T., Harder, H., and Hansen, U. (2010). Thermochemically driven convection in a rotating spherical shell. Geophys. J. Int. 183, 150–162. doi: 10.1111/j.1365-246X.2010.04722.x

CrossRef Full Text | Google Scholar

Buffett, B. (2014). Geomagnetic fluctuations reveal stable stratification at the top of the Earth's core. Nature 507:484. doi: 10.1038/nature13122

PubMed Abstract | CrossRef Full Text | Google Scholar

Buffett, B., Knezek, N., and Holme, R. (2016). Evidence for MAC waves at the top of Earth's core and implications for variations in length of day. Geophys. J. Int. 204, 1789–1800. doi: 10.1093/gji/ggv552

CrossRef Full Text | Google Scholar

Buffett, B. A., and Seagle, C. T. (2010). Stratification of the top of the core due to chemical interactions with the mantle. J. Geophys. Res. 115:B04407. doi: 10.1029/2009JB006751

CrossRef Full Text | Google Scholar

Busse, F. (2002). Is low Rayleigh number convection possible in the Earth's core? Geophys. Res. Lett. 29, 9–1. doi: 10.1029/2001GL014597

CrossRef Full Text | Google Scholar

Cardin, P., and Olson, P. (1992). An experimental approach to thermochemical convection in the Earth's core. Geophys. Res. Lett. 19, 1995–1998. doi: 10.1029/92GL01883

CrossRef Full Text | Google Scholar

Chen, C., and Chen, F. (1991). Experimental study of directional solidification of aqueous ammonium chloride solution. J. Fluid Mech. 227, 567–586. doi: 10.1017/S0022112091000253

CrossRef Full Text | Google Scholar

Christensen, U. R. (2006). A deep dynamo generating Mercury's magnetic field. Nature 444, 1056–1058. doi: 10.1038/nature05342

PubMed Abstract | CrossRef Full Text | Google Scholar

Christensen, U. R. (2015). Iron snow dynamo models for Ganymede. Icarus 247, 248–259. doi: 10.1016/j.icarus.2014.10.024259

CrossRef Full Text | Google Scholar

Christensen, U. R. (2018). Geodynamo models with a stable layer and heterogeneous heat flow at the top of the core. Geophys. J. Int. 215, 1338–1351. doi: 10.1093/gji/ggy352

CrossRef Full Text | Google Scholar

Christensen, U. R., and Aubert, J. (2006). Scaling properties of convection-driven dynamos in rotating spherical shells and application to planetary magnetic fields. Geophys. J. Int. 166, 97–114. doi: 10.1111/j.1365-246X.2006.03009.x

CrossRef Full Text | Google Scholar

Christensen, U. R., and Wicht, J. (2008). Models of magnetic field generation in partly stable planetary cores: applications to Mercury and Saturn. Icarus 196, 16–34. doi: 10.1016/j.icarus.2008.02.013

CrossRef Full Text | Google Scholar

Claßen, S., Heimpel, M., and Christensen, U. (1999). Blob instability in rotating compositional convection. Geophys. Res. Lett. 26, 135–138. doi: 10.1029/1998GL900227

CrossRef Full Text | Google Scholar

Copley, S., Giamei, A. F., Johnson, S., and Hornbecker, M. (1970). The origin of freckles in unidirectionally solidified castings. Metall. Mater. Trans. B 1, 3455–3455. doi: 10.1007/BF02643435

CrossRef Full Text | Google Scholar

de Koker, N., Steinle-Neumann, G., and Vlček, V. (2012). Electrical resistivity and thermal conductivity of liquid Fe alloys at high P and T, and heat flux in Earth's core. Proc. Nat. Acad. Sci. U.S.A. 109, 4070–4073. doi: 10.1073/pnas.1111841109

PubMed Abstract | CrossRef Full Text | Google Scholar

de Wijs, G. A., Kresse, G., Vočadlo, L., Dobson, D., Alfe, D., and Gillan, M. J., et al. (1998). The viscosity of liquid iron at the physical conditions of the Earth's core. Nature 392, 805–807. doi: 10.1038/33905

CrossRef Full Text | Google Scholar

Dobson, D. P. (2002). Self-diffusion in liquid Fe at high pressure. Phys. Earth Planet. Inter. 130, 271–284. doi: 10.1016/S0031-9201(02)00011-0

CrossRef Full Text | Google Scholar

Fernando, H., Chen, R., and Ayotte, B. (1998). Development of a point plume in the presence of background rotation. Phys. Fluids 10, 2369–2383. doi: 10.1063/1.869754

CrossRef Full Text | Google Scholar

Finlay, C. C., and Amit, H. (2011). On flow magnitude and field-flow alignment at Earth's core surface. Geophys. J. Int. 186, 175–192. doi: 10.1111/j.1365-246X.2011.05032.x

CrossRef Full Text | Google Scholar

Frank, D., Landel, J., Dalziel, S., and Linden, P. (2017). Anticyclonic precession of a plume in a rotating environment. Geophys. Res. Lett. 44, 9400–9407. doi: 10.1002/2017GL074191

CrossRef Full Text | Google Scholar

Gastine, T., Wicht, J., and Aubert, J. (2016). Scaling regimes in spherical shell rotating convection. J. Fluid Mech. 808, 690–732. doi: 10.1017/jfm.2016.659

CrossRef Full Text | Google Scholar

Gomi, H., Ohta, K., Hirose, K., Labrosse, S., Caracas, R., and Verstraete, M., et al. (2013). The high conductivity of iron and thermal evolution of the Earth's core. Phys. Earth Planet. Inter. 224, 88–103. doi: 10.1016/j.pepi.2013.07.010

CrossRef Full Text | Google Scholar

Goodman, J. C., Collins, G. C., Marshall, J., and Pierrehumbert, R. T. (2004). Hydrothermal plume dynamics on Europa: implications for chaos formation. J. Geophys. Res. 109:E03008. doi: 10.1029/2003JE002073

CrossRef Full Text | Google Scholar

Gubbins, D. (2007). Geomagnetic constraints on stratification at the top of Earth's core. Earth Planets Space 59, 661–664. doi: 10.1186/BF03352728

CrossRef Full Text | Google Scholar

Gubbins, D., and Davies, C. (2013). The stratified layer at the core–mantle boundary caused by barodiffusion of oxygen, sulphur and silicon. Phys. Earth Planet. Inter. 215, 21–28. doi: 10.1016/j.pepi.2012.11.001

CrossRef Full Text | Google Scholar

Gubbins, D., Masters, G., and Nimmo, F. (2008). A thermochemical boundary layer at the base of Earth's outer core and independent estimate of core heat flux. Geophys. J. Int. 174, 1007–1018. doi: 10.1111/j.1365-246X.2008.03879.x

CrossRef Full Text | Google Scholar

Gubbins, D., Thomson, C., and Whaler, K. (1982). Stable regions in the Earth's liquid core. Geophys. J. Int. 68, 241–251. doi: 10.1111/j.1365-246X.1982.tb06972.x

CrossRef Full Text | Google Scholar

Helffrich, G., and Kaneshima, S. (2010). Outer core compositional stratification from observed core wave speed profiles. Nature 468, 807–810. doi: 10.1038/nature09636

PubMed Abstract | CrossRef Full Text | Google Scholar

Helffrich, G., and Kaneshima, S. (2013). Causes and consequences of outer core stratification. Phys. Earth Planet. Inter. 223, 2–7. doi: 10.1016/j.pepi.2013.07.005

CrossRef Full Text | Google Scholar

Helfrich, K. R., and Battisti, T. M. (1991). Experiments on baroclinic vortex shedding from hydrothermal plumes. J. Geophys. Res. 96, 12511–12518. doi: 10.1029/90JC02643

CrossRef Full Text | Google Scholar

Hirose, K., Morard, G., Sinmyo, R., Umemoto, K., Hernlund, J., and Helffrich, G., et al. (2017). Crystallization of silicon dioxide and compositional evolution of the Earth's core. Nature 543, 99–102. doi: 10.1038/nature21367

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, K., Zhang, R., van Genuchten, M. T., Ewing, R., Brebbia, C., and Gray, W., et al. (1992). A simple particle tracking technique for solving the convection-dispersion equation. Comput. Methods Water Resour. 1, 87–96.

Google Scholar

Huguet, L., Alboussière, T., Bergman, M., Deguen, R., Labrosse, S., and Lesoeur, G. (2016). Structure of a mushy layer under hypergravity with implications for Earth's inner core. Geophys. J. Int. 204, 1729–1755. doi: 10.1093/gji/ggv554

CrossRef Full Text | Google Scholar

Ichikawa, H., and Tsuchiya, T. (2015). Atomic transport property of Fe–O liquid alloys in the Earth's outer core P, T condition. Phys. Earth Planet. Inter. 247, 27–35. doi: 10.1016/j.pepi.2015.03.006

CrossRef Full Text | Google Scholar

Irving, J. C., Cottaar, S., and Lekić, V. (2018). Seismically determined elastic parameters for Earth's outer core. Sci. Adv. 4:eaar2538. doi: 10.1126/sciadv.aar2538

PubMed Abstract | CrossRef Full Text | Google Scholar

Jellinek, A. M., Kerr, R. C., and Griffiths, R. W. (1999). Mixing and compositional stratification produced by natural convection: 1. experiments and their application to Earth's core and mantle. J. Geophys. Res. 104, 7183–7201.

Google Scholar

Jones, C. (2015). “Thermal and compositional convection in the outer core,” in Treatise on Geophysics: Core Dynamics, 2nd Edn, Vol. 8, eds P. Olson and G. Schubert (Amsterdam: Elsevier), 131–185.

Google Scholar

Jones, C. A., Soward, A. M., and Mussa, A. I. (2000). The onset of thermal convection in a rapidly rotating sphere. J. Fluid Mech. 405, 157–179. doi: 10.1017/S0022112099007235

CrossRef Full Text | Google Scholar

Kaneshima, S. (2018). Array analyses of SmKS waves and the stratification of Earth's outermost core. Phys. Earth Planet. Inter. 276, 234–246. doi: 10.1016/j.pepi.2017.03.006

CrossRef Full Text | Google Scholar

Kaneshima, S., and Helffrich, G. (2013). Vp structure of the outermost core derived from analysing large-scale array data of SmKS waves. Geophys. J. Int. 193, 1537–1555. doi: 10.1093/gji/ggt042

CrossRef Full Text | Google Scholar

Kaneshima, S., and Matsuzawa, T. (2015). Stratification of Earth's outermost core inferred from SmKS array data. Prog. Earth Planet. Sc. 2:15. doi: 10.1186/s40645-015-0046-5

CrossRef Full Text | Google Scholar

Kerr, R. C., and Mériaux, C. (2004). Structure and dynamics of sheared mantle plumes. Geochem. Geophys. 5:Q12009. doi: 10.1029/2004GC000749

CrossRef Full Text | Google Scholar

Kerr, R. C., Mériaux, C., and Lister, J. R. (2008). Effect of thermal diffusion on the stability of strongly tilted mantle plume tails. J. Geophys. Res. 113:B09401. doi: 10.1029/2007JB005510

CrossRef Full Text | Google Scholar

Konôpková, Z., McWilliams, R. S., Gómez-Pérez, N., and Goncharov, A. F. (2016). Direct measurement of thermal conductivity in solid iron at planetary core conditions. Nature 534:99. doi: 10.1038/nature18009

PubMed Abstract | CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Labrosse, S., Poirier, J.-P., and Mouël, J.-L. L. (1997). On cooling of the Earth's core. Phys. Earth Planet. Inter. 99, 1–17. doi: 10.1016/S0031-9201(96)03207-4

CrossRef Full Text | Google Scholar

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

CrossRef Full Text | Google Scholar

Lay, T., and Young, C. J. (1990). The stably-stratified outermost core revisited. Geophys. Res. Lett. 17, 2001–2004. doi: 10.1029/GL017i011p02001

CrossRef Full Text | Google Scholar

Lesur, V., Whaler, K., and Wardinski, I. (2015). Are geomagnetic data consistent with stably stratified flow at the core–mantle boundary? Geophys. J. Int. 201, 929–946. doi: 10.1093/gji/ggv031

CrossRef Full Text | Google Scholar

Levy, M., and Fernando, H. (2002). Turbulent thermal convection in a rotating stratified fluid. J. Fluid Mech. 467, 19–40. doi: 10.1017/S0022112002001350

CrossRef Full Text | Google Scholar

Lister, J. R., and Buffett, B. A. (1995). The strength and efficiency of thermal and compositional convection in the geodynamo. Phys. Earth Planet. Inter. 91, 17–30. doi: 10.1016/0031-9201(95)03042-U

CrossRef Full Text | Google Scholar

Lister, J. R., and Buffett, B. A. (1998). Stratification of the outer core at the core-mantle boundary. Phys. Earth Planet. Inter. 105, 5–19. doi: 10.1016/S0031-9201(97)00082-4

CrossRef Full Text | Google Scholar

Loper, D. (2007). Turbulence and small-scale dynamics in the core. Treat. Geophys. 8, 187–206. doi: 10.1016/B978-044452748-6/00131-0

CrossRef Full Text | Google Scholar

Loper, D. E., and Moffatt, H. K. (1993). Small-scale hydromagnetic flow in the Earth's core: rise of a vertical buoyant plume. Geophys. Astrophys. Fluid Dyn. 68, 177–202. doi: 10.1080/03091929308203567

CrossRef Full Text | Google Scholar

Loper, D. E., and Roberts, P. H. (1981). A study of conditions at the inner core boundary of the Earth. Phys. Earth Planet. Inter. 24, 302–307. doi: 10.1016/0031-9201(81)90117-5

CrossRef Full Text | Google Scholar

Manglik, A., Wicht, J., and Christensen, U. R. (2010). A dynamo model with double diffusive convection for Mercury's core. Earth Planet. Sci. Lett. 289, 619–628. doi: 10.1016/j.epsl.2009.12.007

CrossRef Full Text | Google Scholar

Matsushima, M. (2004). Scale similarity of mhd turbulence in the earth's core. Earth Planets Space 56, 599–605. doi: 10.1186/BF03352520

CrossRef Full Text | Google Scholar

Michalke, A. (1984). Survey on jet instability theory. Prog. Aerosp. Sci. 21, 159–199. doi: 10.1016/0376-0421(84)90005-8

CrossRef Full Text | Google Scholar

Moffatt, H. (1989). “Liquid metal MHD and the geodynamo,” Proceedings of IUTAM Symposium on Liquid Metal Magnetohydrodynamics (Riga), 403–412.

Google Scholar

Moffatt, H., and Loper, D. (1994). The magnetostrophic rise of a buoyant parcel in the Earth's core. Geophys. J. Int. 117, 394–402. doi: 10.1111/j.1365-246X.1994.tb03939.x

CrossRef Full Text | Google Scholar

Munz, C.-D. (1988). On the numerical dissipation of high resolution schemes for hyperbolic conservation laws. J. Comput. Phys. 77, 18–39. doi: 10.1016/0021-9991(88)90154-4

CrossRef Full Text | Google Scholar

Nakagawa, T. (2011). Effect of a stably stratified layer near the outer boundary in numerical simulations of a magnetohydrodynamic dynamo in a rotating spherical shell and its implications for Earth's core. Phys. Earth Planet. Inter. 187, 342–352. doi: 10.1016/j.pepi.2011.06.001

CrossRef Full Text | Google Scholar

Nataf, H.-C., and Gagniere, N. (2008). On the peculiar nature of turbulence in planetary dynamos. C R Phys. 9, 702–710. doi: 10.1016/j.crhy.2008.07.009

CrossRef Full Text | Google Scholar

Nataf, H.-C., and Schaeffer, N. (2015). “Turbulence in the core,” in Treatise on Geophysics,16 2nd Edn. (Amsterdam: Elsevier), 161–181.

Google Scholar

Ohta, K., Kuwayama, Y., Hirose, K., Shimizu, K., and Ohishi, Y. (2016). Experimental determination of the electrical resistivity of iron at Earth's core conditions. Nature 534, 95–98. doi: 10.1038/nature17957

PubMed Abstract | CrossRef Full Text | Google Scholar

Olson, P., Landeau, M., and Reynolds, E. (2017). Dynamo tests for stratification below the core-mantle boundary. Phys. Earth Planet. Inter. 271, 1–18. doi: 10.1016/j.pepi.2017.07.003

CrossRef Full Text | Google Scholar

Olson, P., Landeau, M., and Reynolds, E. (2018). Outer core stratification from the high latitude structure of the geomagnetic field. Front. Earth Sci. 6:140. doi: 10.3389/feart.2018.00140

CrossRef Full Text | Google Scholar

O'Rourke, J. G., and Stevenson, D. J. (2016). Powering Earth's dynamo with magnesium precipitation from the core. Nature 529, 387–389. doi: 10.1038/nature16495

PubMed Abstract | CrossRef Full Text | Google Scholar

Perrillat, J.-P., Mezouar, M., Garbarino, G., and Bauchau, S. (2010). In situ viscometry of high-pressure melts in the Paris–Edinburgh cell: application to liquid FeS. High Press. Res. 30, 415–423. doi: 10.1080/08957959.2010.494844

CrossRef Full Text | Google Scholar

Poirier, J. (1988). Transport properties of liquid metals and viscosity of the Earth's core. Geophys. J. Int. 92, 99–105. doi: 10.1111/j.1365-246X.1988.tb01124.x

CrossRef Full Text | Google Scholar

Posner, E. S., Rubie, D. C., Frost, D. J., and Steinle-Neumann, G. (2017). Experimental determination of oxygen diffusion in liquid iron at high pressure. Earth Planet. Sci. Lett. 464, 116–123. doi: 10.1016/j.epsl.2017.02.020

CrossRef Full Text | Google Scholar

Pozzo, M., Davies, C., Gubbins, D., and Alfe, D. (2012). Thermal and electrical conductivity of iron at Earth's core conditions. Nature 485, 355–358. doi: 10.1038/nature11031

PubMed Abstract | CrossRef Full Text | Google Scholar

Pozzo, M., Davies, C., Gubbins, D., and Alfè, D. (2013). Transport properties for liquid silicon-oxygen-iron mixtures at Earth's core conditions. Phys. Rev. B 87:014110. doi: 10.1103/PhysRevB.87.014110

CrossRef Full Text | Google Scholar

Roberts, P. H., and King, E. M. (2013). On the genesis of the Earth's magnetism. Rep. Prog. Phys. 76:096801. doi: 10.1088/0034-4885/76/9/096801

PubMed Abstract | CrossRef Full Text | Google Scholar

Simitev, R., and Busse, F. (2003). Patterns of convection in rotating spherical shells. New J. Phys. 5:97. doi: 10.1088/1367-2630/5/1/397

CrossRef Full Text | Google Scholar

Simitev, R., and Busse, F. (2005). Prandtl-number dependence of convection-driven dynamos in rotating spherical fluid shells. J. Fluid Mech. 532, 365–388. doi: 10.1017/S0022112005004398

CrossRef Full Text | Google Scholar

Sironi, L., and Spitkovsky, A. (2009). Synthetic spectra from particle-in-cell simulations of relativistic collisionless shocks. Astrophys. J. Lett. 707:L92. doi: 10.1088/0004-637X/707/1/L92

CrossRef Full Text | Google Scholar

Sreenivasan, B., and Narasimhan, G. (2017). Damping of magnetohydrodynamic waves in a rotating fluid. J. Fluid Mech. 828, 867–905. doi: 10.1017/jfm.2017.550

CrossRef Full Text | Google Scholar

St Pierre, M. G. (1996). On the local nature of turbulence in Earth's outer core. Geophys. Astrophys. Fluid Dyn. 83, 293–306. doi: 10.1080/03091929608208969

CrossRef Full Text | Google Scholar

Tackley, P. J., and King, S. D. (2003). Testing the tracer ratio method for modeling active compositional fields in mantle convection simulations. Geochem. Geophys. 4. doi: 10.1029/2001GC000214

CrossRef Full Text | Google Scholar

Takehiro, S.-I., and Lister, J. R. (2001). Penetration of columnar convection into an outer stably stratified layer in rapidly rotating spherical fluid shells. Earth Planet. Sci. Lett. 187, 357–366. doi: 10.1016/S0012-821X(01)00283-7

CrossRef Full Text | Google Scholar

Tanaka, S. (2007). Possibility of a low P-wave velocity layer in the outermost core from global SmKS waveforms. Earth Planet. Sci. Lett. 259, 486–499. doi: 10.1016/j.epsl.2007.05.007

CrossRef Full Text | Google Scholar

Tang, L.Masutani, S. M., et al. (2003). “Laminar to turbulent flow liquid-liquid jet instability and breakup,” in The Thirteenth International Offshore and Polar Engineering Conference (Honolulu, HI).

Google Scholar

Tskhakaya, D., Matyash, K., Schneider, R., and Taccogna, F. (2007). The Particle-In-Cell Method. Contrib. Plasm. Phys. 47, 563–594. doi: 10.1002/ctpp.200710072

CrossRef Full Text | Google Scholar

Vočadlo, L., Alfè, D., Price, G. D., and Gillan, M. J. (2000). First principles calculations on the diffusivity and viscosity of liquid Fe–S at experimentally accessible conditions. Phys. Earth Planet. Inter. 120, 145–152. doi: 10.1016/S0031-9201(00)00151-5

CrossRef Full Text | Google Scholar

Wong, J., Davies, C. J., and Jones, C. A. (2018). A Boussinesq slurry model of the F–layer at the base of Earth's outer core. Geophys. J. Int. 214, 2236–2249. doi: 10.1093/gji/ggy245

CrossRef Full Text | Google Scholar

Worster, M. G., and Huppert, H. E. (1983). Time-dependent density profiles in a filling box. J. Fluid Mech. 132, 457–466. doi: 10.1017/S002211208300172X

CrossRef Full Text | Google Scholar

Zhang, K.-K. (1991). Convection in a rapidly rotating spherical shell at infinite Prandtl number: steadily drifting rolls. Phys. Earth Planet. Inter. 68, 156–169. doi: 10.1016/0031-9201(91)90015-A

CrossRef Full Text | Google Scholar

Zhang, K.-K. (1992). Convection in a rapidly rotating spherical shell at infinite Prandtl number: transition to vacillating flows. Phys. Earth Planet. Inter. 72, 236–248. doi: 10.1016/0031-9201(92)90204-9

CrossRef Full Text | Google Scholar

Zhang, K.-K., and Busse, F. (1990). Generation of magnetic fields by convection in a rotating spherical fluid shell of infinite Prandtl number. Phys. Earth Planet. Inter. 59, 208–222. doi: 10.1016/0031-9201(90)90226-N

CrossRef Full Text | Google Scholar