Original Research ARTICLE
Chemical Convection and Stratification in the Earth's Outer Core
- 1Department Planets and Comets, Max Planck Institute for Solar System Research, Göttingen, Germany
- 2Laboratoire de Planétologie et Géodynamique, Université de Nantes, Nantes, France
- 3Laboratoire de Géologie de Lyon, Université de Lyon, École Normale Supérieure de Lyon, Université Lyon-1, CNRS, UMR 5276, Lyon, France
Convection in the Earth's outer core is driven by buoyancy sources of both thermal and compositional origin. The thermal and compositional molecular diffusivities differ by several orders of magnitude, which can affect the dynamics in various ways. So far, the large majority of numerical simulations have been performed within the codensity framework that consists in combining temperature and composition, assuming artificially enhanced diffusivities for both variables. In this study, we use a particle-in-cell method implemented in a 3D dynamo code to conduct a first qualitative exploration of pure compositional convection in a rotating spherical shell. We focus on the end-member case of infinite Schmidt number by totally neglecting the compositional diffusivity. We show that compositional convection has a very rich physics that deserves several more focused and quantitative studies. We also report, for the first time in numerical simulations, the self-consistent formation of a chemically stratified layer at the top of the shell caused by the accumulation of chemical plumes and blobs emitted at the bottom boundary. When applied to likely numbers for the Earth's core, some (possibly simplistic) physical considerations suggest that a stratified layer formed in such a scenario would be probably weakly stratified and may be compatible with magnetic observations.
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 SiO2 (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−6 − 10−5 m2s−1 (Poirier, 1988; de Wijs et al., 1998; Roberts and King, 2013), κξ ~ 5 × 10−9 − 2 × 10−8 m2s−1 (Dobson, 2002; Pozzo et al., 2013; Ichikawa and Tsuchiya, 2015; Posner et al., 2017), and ν ~ 5 × 10−7 − 5 × 10−6 m2s−1 (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:
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,
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:
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−6, 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,
where ξ0 and ρ0 are the reference compositional mass fraction and density, respectively, and β denotes the coefficient of chemical expansivity.
In addition to the mass conservation,
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,
where is the fluid velocity, t the time, Ω the Earth's rotation rate around an axis oriented by , Π the dynamic pressure and go the gravitational acceleration at the top of the shell at radius r = ro. The Lorentz force is neglected. We adopt the viscous time scale D2/ν in which the thickness of the shell D = ro − ri is taken as the typical length scale, ri being the radius of the inner boundary. Hence, using (ρ0νΩ) as a scaling for pressure, one can write the dimensionless momentum equation,
in which E is the Ekman number, defined by
and Raξ a “modified” compositional Rayleigh (or Roberts-Rayleigh) number,
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 based on the chemical diffusivity,
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
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
In the end-member case Sc = +∞ (κξ = 0) considered here, Equation (12) becomes a first-order hyperbolic equation:
2.2. Boundary Conditions
We adopt non-penetration and no-slip conditions at both boundaries which implies
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 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 via a source term distributed within a thin layer of thickness h above the ICB, so that the integrated added mass fraction matches the ICB chemical flux:
The choice of an adequate function for 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:
The constants in the exponential function are empirically set to σ = h/5 and h = 0.4E1/3. These values will be justified later in section 3.1. The scaling factor σ0 is computed to satisfy Equation (15). Mathematically, this formulation is equivalent to solving the transport equation
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 , where 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:
Using the corresponding compositional buoyancy flow injected at the inner boundary,
one can derive the following scale for composition:
The expression of the modified compositional Rayleigh number then becomes
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 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π/mc-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.
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. (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 ϕ0 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−5 − 10−3. 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,
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. 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. (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. 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−3, Ro = 1.9 × 10−3 and Ro = 4.5 × 10−4. 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. 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. (A) Equatorial section of the compositional field in case 10 at time t = 0.066 and Ro = 4.4 × 10−4. (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 , where H is the height of water tank and B0 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−4 and Ro = 1.9 × 10−3, 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−5 and Ro = 4.4 × 10−4 (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 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. (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 rb above which the mean Brunt-Väisälä frequency N is positive:
where denotes the average density at radius r . For pure compositional convection, Equation (23) is equivalent to
in which is the mean composition at radius r. The bottom of the stratified layer rb, defined by Equation (24), is indicated by dashed lines on Figures 11A,B, which show radial profiles of the mean composition and (N/Ω)2 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 rb which require time averaging. In addition, the definition of rb 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. (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/Ω)2 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−3 and only 20% for case 6 with E = 3 × 10−4 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−5 s−1) 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−3 s−1) (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. 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
Assuming the scaling law (22) can be used as a good proxy for the size of chemical plumes in core conditions, we can estimate
Taking E = 10−15 and Re = 108 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
Taking Ro = U/DΩ ~ 10−6 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 δ,
where τvisc denotes the viscous time D2/ν. For example, E = 10−4 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−4 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
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 , 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
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−11 m s−1 (Labrosse, 2015), Ω ~ 7 × 10−5 s−1 and g ~ 3 m s−2 leads to
The parameter f was found to be close to 10−3 − 10−2 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−4 m s−1 below the CMB (Finlay and Amit, 2011). Estimations from numerical simulations predict a slightly higher root-mean-square (RMS) velocity of order 10−3 m s−1 (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−2 m s−1 for the chemical blobs gives
Figure 12 shows the mean Brunt-Väisälä frequency N obtained when spreading linearly a density anomaly Δρ ~ 10−7ρ over a stratified layer of thickness H,
For instance, g ~ 10 m s−2 and H = 100 km would give N ~ 3 × 10−6 s−1, 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−5 s−1 compatible with magnetic observations would require either a thinner stratified layer (H ~ 10 km) or a higher density anomaly such that Δρ/ρ ~ 10−6. This density ratio ratio would then give f ~ 10−5 and w ~ 0.1 m s−1. 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−3 s−1) called for by some seismic studies can be achieved only by invoking unrealistically high blobs velocities (w ≳ 102 m s−1). 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−5 s−1) and, very likely, seismically invisible, we cannot exclude stronger stratifications depending on the precise processes occurring in the ICB region.
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−7.
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.
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.
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.
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.
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
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
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
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
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
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
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
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
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
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
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
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
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
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.
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
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.
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
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
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
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
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
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
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
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
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).
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
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
Keywords: compositional convection, stratification, core dynamics, particle-in-cell, infinite Schmidt number
Citation: Bouffard M, Choblet G, Labrosse S and Wicht J (2019) Chemical Convection and Stratification in the Earth's Outer Core. Front. Earth Sci. 7:99. doi: 10.3389/feart.2019.00099
Received: 15 July 2018; Accepted: 23 April 2019;
Published: 10 May 2019.
Edited by:Peter E. Driscoll, Carnegie Institution for Science (CIS), United States
Reviewed by:Binod Sreenivasan, Indian Institute of Science (IISc), India
Kenneth Philip Kodama, Lehigh University, United States
Copyright © 2019 Bouffard, Choblet, Labrosse and Wicht. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Mathieu Bouffard, firstname.lastname@example.org