Rotating double-diffusive convection in stably stratified planetary cores

TL;DR

This paper investigates rotating double-diffusive convection in planetary cores, revealing large-scale instabilities and flow regimes influenced by stratification and rotation, with implications for Earth's early core dynamics.

Contribution

It provides a comprehensive analysis of RDDC in spheres, highlighting the inviscid nature of instabilities, new scaling laws, and the impact of stratification on planetary core convection.

Findings

Large-scale convective motions occur at lower Rayleigh numbers in stably stratified spheres.

The critical Rayleigh number scales as Ra ~ Ek^{-1} in stably stratified regimes.

Double diffusion significantly reduces the critical Rayleigh number for Earth's core conditions.

Abstract

In planetary fluid cores, the density depends on temperature and chemical composition, which diffuse at very different rates. This leads to various instabilities, bearing the name of double-diffusive convection. We investigate rotating double-diffusive convection (RDDC) in fluid spheres. We use the Boussinesq approximation with homogeneous internal thermal and compositional source terms. We focus on the finger regime, in which the thermal gradient is stabilising whereas the compositional one is destabilising. First, we perform a global linear stability analysis in spheres. The critical Rayleigh numbers drastically drop for stably stratified fluids, yielding large-scale convective motions where local analyses predict stability. We evidence the inviscid nature of this large-scale double-diffusive instability, enabling the determination of the marginal stability curve at realistic planetary regimes. In particular, we show that in stably stratified spheres, the Rayleigh numbers $Ra$ at the onset evolve like $Ra \sim Ek^{-1}$, where $Ek$ is the Ekman number. This differs from rotating convection in unstably stratified spheres, for which $Ra \sim Ek^{-4/3}$. The domain of existence of inviscid convection thus increases as $Ek^{-1/3}$. Second, we perform nonlinear simulations. We find a transition between two regimes of RDDC, controlled by the strength of the stratification. Furthermore, far from the RDDC onset, we find a dominating equatorially anti-symmetric, large-scale zonal flow slightly above the associated linear onset. Unexpectedly, a purely linear mechanism can explain this phenomenon, even far from the instability onset, yielding a symmetry breaking of the nonlinear flow at saturation. For even stronger stable stratification, the flow becomes mainly equatorially-symmetric and intense zonal jets develop. Finally, we apply our results to the early Earth core. Double diffusion can reduce the critical Rayleigh number by four decades for realistic core conditions. We suggest that the early Earth core was prone to turbulent RDDC, with large-scale zonal flows.

Tables

Table 1. Table 1 : Dimensionless numbers characterising diffusive effects and typical values in the Earth’s liquid core (Braginsky & Roberts, 1995 ; Labrosse, 2015 ) and stably stratified stellar envelopes (Garaud et al., 2015 ) . Kinematic viscosity ν 𝜈 \nu , thermal diffusivity κ T subscript 𝜅 𝑇 \kappa_{T} , compositional diffusivity κ C subscript 𝜅 𝐶 \kappa_{C} , planetary angular velocity Ω s subscript Ω 𝑠 \Omega_{s} and radius R 𝑅 R .

Symbol	Name	Definition	Earth (current)	Stars
$L$	Lewis	${\kappa }_T/{\kappa }_C$	${10}^4$	${10}^3-{10}^7$
$Pr$	Prandtl	$\nu /{\kappa }_T$	$0.01-0.1$	${10}^{-6}$
$Sc$	Schmidt	$\nu /{\kappa }_C$	${10}^2-{10}^3$	${10}^{-3}-{10}^1$
$Ek$	Ekman	$\nu /\left({\Omega }_s{R}^2\right)$	${10}^{-15}$	${10}^{-18}$

Table 2. Table 2 : Critical azimuthal wave number m c superscript 𝑚 𝑐 m^{c} , compositional Rayleigh number R ​ a C c 𝑅 superscript subscript 𝑎 𝐶 𝑐 Ra_{C}^{c} and angular frequency ω c superscript 𝜔 𝑐 \omega^{c} at the linear onset ( σ = 0 𝜎 0 \sigma=0 ) of compositional overturning convection (i.e. for R ​ a T = 0 𝑅 subscript 𝑎 𝑇 0 Ra_{T}=0 ). Computations with SINGE at S ​ c = 3 𝑆 𝑐 3 Sc=3 .

$Ek$	$m^c$	$R{a}_C^c$	$\omega$
${10}^{-4}$	10	$1.59\times {10}^6$	$-7.74\times {10}^1$
${10}^{-5}$	20	$3.40\times {10}^7$	$-4.63\times {10}^2$
${10}^{-6}$	55	$7.60\times {10}^8$	$-2.32\times {10}^3$

Table 3. Table 3 : Parameters for the onset of convection at core conditions ( E ​ k = 10 − 15 𝐸 𝑘 superscript 10 15 Ek=10^{-15} ): critical thermal R ​ a T c 𝑅 superscript subscript 𝑎 𝑇 𝑐 Ra_{T}^{c} or compositional R ​ a C c 𝑅 superscript subscript 𝑎 𝐶 𝑐 Ra_{C}^{c} Rayleigh numbers, critical wave number m c superscript 𝑚 𝑐 m^{c} computed from table 2 of Jones et al. ( 2000 ) . Note that our dimensionless numbers E ​ k 𝐸 𝑘 Ek and R ​ a 𝑅 𝑎 Ra differ from theirs.

	Thermal		Compositional
$Pr$	$0.01$	$0.1$
$Sc$			${10}^2$
$R{a}_T^c$	$6.00\times {10}^{19}$	$1.42\times {10}^{20}$	$0$
$R{a}_C^c$	$0$	$0$	$1.08\times {10}^{21}$
$m^c$	$10755$	$22301$	$67570$

Table 4. Table 4 : Critical wave number m c superscript 𝑚 𝑐 m^{c} , thermal Rayleigh number R ​ a T c 𝑅 superscript subscript 𝑎 𝑇 𝑐 Ra_{T}^{c} and angular frequency ω c superscript 𝜔 𝑐 \omega^{c} at the marginal onset ( σ = 0 𝜎 0 \sigma=0 ) of thermal overturning convection (i.e. for R ​ a C = 0 𝑅 subscript 𝑎 𝐶 0 Ra_{C}=0 ). Computations at S ​ c = 3 𝑆 𝑐 3 Sc=3 and P ​ r = 0.3 𝑃 𝑟 0.3 Pr=0.3 . The first row is obtained at E ​ k = 10 − 4 𝐸 𝑘 superscript 10 4 Ek=10^{-4} and the second one at E ​ k = 10 − 5 𝐸 𝑘 superscript 10 5 Ek=10^{-5} .

$Ek$	$m^c$	$R{a}_T^c$	$\omega$
${10}^{-4}$	12	$9.86\times {10}^5$	$-5.48\times {10}^2$
${10}^{-5}$	40	$3.60\times {10}^7$	$-2.08\times {10}^3$

Equations

$${\rho }^{\ast }/{\rho }_m=1-{\alpha }_T(T^{\ast }-T_m)-{\alpha }_C(C^{\ast }-C_m)$$

$${\kappa }_T{\nabla }^2{T}_0^{\ast }=-{\mathcal{Q}}_T,{\kappa }_C{\nabla }^2{C}_0^{\ast }=-{\mathcal{Q}}_C,$$

$$T_0^{\ast }(r)$$

$$C_0^{\ast }(r)$$

$$T_0(r)=\frac{1-r^2}{Pr},C_0(r)=\frac{1-r^2}{Sc},$$

$$Pr=\frac{\nu }{{\kappa }_T},Sc=\frac{\nu }{{\kappa }_C},L=\frac{Sc}{Pr},$$

$$\frac{\partial \mathbf{u}}{\partial t}+(\mathbf{u}\cdot \nabla )\mathbf{u}$$

$$+\left(R{a}_T\Theta +R{a}_C\xi \right)r{1}_r,$$

$$\frac{\partial \Theta }{\partial t}+(\mathbf{u}\cdot \nabla )\Theta$$

$$\frac{\partial \xi }{\partial t}+(\mathbf{u}\cdot \nabla )\xi$$

$$\nabla \cdot \mathbf{u}$$

$$Ek=\frac{\nu }{{\Omega }_s{R}^2},$$

$$R{a}_T=\frac{{\alpha }_T{g}_0{\mathcal{Q}}_T{R}^6}{6\nu {{\kappa }_T}^2},R{a}_C=\frac{{\alpha }_C{g}_0{\mathcal{Q}}_C{R}^6}{6\nu {{\kappa }_C}^2}$$

$$\mathbf{u}\cdot 1_r=0,\mathbf{u}\times 1_r=0\text{at}r=1.$$

$$\frac{\partial \Theta }{\partial r}=\frac{\partial \xi }{\partial r}=0\text{at}r=1.$$

$$N^2=\mathbf{g}\cdot \nabla \left(\frac{\rho }{{\rho }_m}\right).$$

$$N_{0,T}^2=-{\alpha }_T\mathbf{g}\cdot \nabla T_0^{\ast }\text{and}{N}_{0,C}^2=-{\alpha }_C\mathbf{g}\cdot \nabla C_0^{\ast }$$

$$\frac{N_0^2}{{\Omega }_s^2}(r)=-2r^{2}E{k}^2\left(\frac{R{a}_T}{Pr}+\frac{R{a}_C}{Sc}\right)=-2r^{2}R{o}_c^2,$$

$$E_{\{u,t,c\}}=\frac{1}{2}\int \left\{|\mathbf{u}{|}^2,{\Theta }^2,{\xi }^2\right\}\mathrm{d}V.$$

$$N_0^2=N_{0,T}^2+N_{0,C}^2>0,\text{i.e.}R{a}_{T}L<-R{a}_C.$$

$$N_{0,T}^2>0,\text{i.e.}R{a}_T>0.$$

$$R_0=\frac{{\alpha }_T}{{\alpha }_C}\cdot \frac{|\nabla T_0^{\ast }|}{|\nabla C_0^{\ast }|}\sim \left|\frac{R{a}_T}{R{a}_C}\right|L,$$

$$1\le R_0<L,\text{i.e.}|R{a}_T|\le R{a}_C.$$

$$1\le R_0^{-1}\le \frac{Pr+1}{Pr+1/L}.$$

$$\left[\mathbf{u},\Theta ,\xi \right]\propto \exp \left(\mathrm{i}\mathbf{k}\cdot {\mathbf{r}}_0+\lambda t\right),$$

$$\left[\mathbf{u},\Theta ,\xi \right](\mathbf{r},t)=\left[\hat{\mathbf{u}},\hat{\Theta },\hat{\xi }\right](r,\theta )\exp \left[\mathrm{i}m\varphi +\lambda t\right],$$

$$\mathcal{A}\mathbf{X}=\lambda \mathcal{B}\mathbf{X},$$

$$\lambda =\lambda (m,s,R{a}_T,R{a}_C,Pr,Sc,Ek).$$

$$R{a}_C=-R{a}_T\text{for}R{a}_C,|R{a}_T|\to \infty .$$

$$l_{\perp }$$

$${\tilde{R}}_0=\frac{R_0-1}{L-1}.$$

$$Re=\sqrt{\frac{2}{V}{E}_u},Ro=ReEk,$$

Full text

1 Introduction

1.1 Geophysical context

Thermo-compositional convection stirs motions in the Earth’s core (Jones, 2015), that sustain large-scale magnetic fields via dynamo action. The thermal part is generated by the super-adiabatic thermal gradient. It mainly comes from the secular cooling of the core, driven by the heat extracted at the core-mantle boundary (CMB). Additionally, because of this cooling, latent heat is released by the crystallisation of the inner core (Verhoogen, 1961). Radioactive heat sources can also participate, although their contribution is debated (e.g. Hirao et al., 2006; Bouhifd et al., 2007; Chidester et al., 2017). The compositional part is sustained by the ejection of light elements into the fluid core, mainly due to the solidification of the inner core (e.g. Fearn & Loper, 1981). Currently, compositional buoyancy is expected to dominate over thermal buoyancy (Braginsky & Roberts, 1995; Lister & Buffett, 1995; Buffett et al., 1996). Few models have considered individual contributions of thermal and compositional buoyancies for the present dynamics of the core, by using experiments (Cardin & Olson, 1992), asymptotic models (Busse, 2002; Simitev, 2011) or numerical simulations (e.g. Glatzmaier & Roberts, 1996; Kutzner & Christensen, 2000; Hori et al., 2012; Bouffard, 2017).

The crystallisation of the inner core is a rather recent geophysical feature, initiated 1 Ga or 2 Ga ago (Labrosse, 2015). However, the geodynamo is active since at least 3.45 Ga (Usui et al., 2009; Tarduno et al., 2010), despite the absence of the main buoyancy source (crystallization of the inner core). Moreover, driving the early geodynamo by thermal buoyancy alone requires large secular cooling rates (Gubbins et al., 2003). Such fast cooling rates are problematic for most thermal histories (Labrosse, 2015), although allowed by the large remaining uncertainties (e.g. Williams, 2018). Prior the inner core crystallization, a large fraction of the core is expected to present a sub-adiabatic temperature (Nimmo, 2015; Labrosse, 2015), inhibiting (thermal) convective motions. Therefore, determining the origin of the fluid motions sustaining the early geodynamo is elusive.

It has been suggested that light elements, dissolved during the core formation (e.g. Badro et al., 2015), may have been exsolved due to the secular cooling (Buffett et al., 2000). The exsolution of buoyant magnesium oxide would provide compositional buoyancy, notably prior to the nucleation of the inner core (O’Rourke & Stevenson, 2016; Badro et al., 2016). This mechanism has been criticised, e.g. because the magnesium solubility in the core depends not only on the temperature but also strongly on the oxygen content (Du et al., 2017). Moreover, this scenario requires a core formation at extremely high temperature to incorporate a sufficient amount of magnesium. Instead, Hirose et al. (2017) advocated for top-down crystallisation of silicon oxides, incorporated in the core via the metal-segregation processes in a deep magma ocean at moderate temperatures. These non-standard mechanisms put forward the possibility to drive the early geodynamo by double-diffusive convection.

1.2 Double-diffusive convection

Double-diffusive convection (DDC) refers to various buoyancy-driven instabilities, generated by two different components of buoyancy. For planetary cores, we refer to thermal and chemical buoyancies. The two sources diffuse at different rates, with the thermal (fast) diffusivity $\kappa_{T}$ and the chemical (slow) one $\kappa_{C}$. Their ratio defines the dimensionless Lewis number $L=\kappa_{T}/\kappa_{C}$, which is expected to be at least $10^{3}$ (Braginsky & Roberts, 1995) in planetary cores (see table 1).

DDC takes different forms, depending on the value of $L$ and on the sign of the mean gradients of each individual component of the density. Classical convection occurs when both thermal and compositional gradients are destabilising. Then, we distinguish (i) the finger regime (Stern, 1960), when the chemical gradient is unstable and the thermal one stable, and (ii) the semi-convection quadrant (Spiegel, 1969) with a stabilising compositional gradient and a destabilising thermal one. Recently, double-diffusive effects have been evidenced even with slightly stabilising thermal gradients, leading to finger convection for unstable stratification (e.g. Kellner & Tilgner, 2014).

DDC has been mainly studied for oceanographic purposes (e.g. Schmitt, 1994; Radko, 2013). Applications has become also apparent in astrophysics (e.g. Garaud, 2018) or mantle physics (Hansen & Yuen, 1988, 1989, 1990). Rotational effects have been largely neglected in these works. Only a few studies investigated rotating double-diffusive convection (RDDC), usually by considering rotational effects in local Cartesian models. Under this assumption, rotation has essentially a stabilising effect (Acheson, 1980; Pearlstein, 1981; Moll et al., 2017; Sengupta & Garaud, 2018). Yet, the relevance of these local models remains elusive for rapidly rotating planetary cores. Indeed, a subtle interplay between the rapid rotation and the bounded spherical geometry is expected for RDDC. Notably, Busse (2002) predicted asymptotically the existence of double-diffusive convection at low Rayleigh numbers in rapidly rotating fluids cores, by extending his reduced annulus model (Busse, 1970). Simitev (2011) did confirm these predictions numerically in the annulus geometry. Finally, only few studies tackled RDDC in spherical geometries with both unstable buoyancies (Glatzmaier & Roberts, 1996; Breuer et al., 2010; Trümper et al., 2012; Takahashi, 2014), and even fewer with antagonist gradients (Manglik et al., 2010; Net et al., 2012).

1.3 Computational methods

Simulations of RDDC in spherical geometry are computationally challenging. A major difficulty is to use small enough values of $\kappa_{C}$ for fixed values of $\kappa_{T}$, to probe the regime $L\gg 1$. This means that the spatial resolution must be adequate, for simulating both the fine-scale compositional structures and the thermal ones. In addition, planetary cores are generally rapidly rotating, as measured by the dimensionless Ekman number $Ek\ll 1$ (table 1). Thus, RDDC must be investigated in the regime $Ek\ll 1$ simultaneously with $L\gg 1$. Eulerian numerical methods cannot presently encompass this broad range of length (and time) scales properly. Hence, computations are always performed for dimensionless parameters orders of magnitude away from core values.

To circumvent these issues, a ”particle-in-cell” (PIC) method has been developed (Bouffard, 2017; Bouffard et al., 2019). It models the compositional field in the limit $L\gg 1$ as a collection of advected particles, while keeping an Eulerian description for velocity and temperature fields. While PIC methods excel in the diffusionless limit $\kappa_{C}=0$, they suffer from several drawbacks at finite values of $L$. For instance, Bouffard et al. (2017) showed that the PIC approach currently does not compare well with proposed benchmarks of RDDC in spherical geometry (Breuer et al., 2010), obtained at finite values of $L$. Finally, even if mixing Eulerian and PIC methods may be desirable for initial value problems, this approach prevents from efficiently finding the instability onset. In contrast, the determination of the onset with Eulerian methods reduces to eigenvalue problems, which can be solved efficiently (e.g. for convection Net et al., 2012; Kaplan et al., 2017).

1.4 Outline

In this study, we aim at investigating numerically RDDC in spherical bodies. We are motivated by explaining the origin of the early geodynamo and by the potential importance of the double-diffusive effects highlighted by Busse (2002) and Simitev (2011). We will focus on rotating full spheres, without inner cores. Beyond the geophysical motivation, a full sphere geometry is the simpler configuration to illustrate the intricate influence of rotation and global geometry on RDDC. Moreover, we will employ the classical Eulerian description, for which efficient codes are available.

The paper is organised as follows. The formulation of the problem is described in §2, together with our numerical method of choice. In §3, we draw physical insights from existing local stability analyses. Then, we conduct a global stability analysis in spheres in §4, and we compare it with the asymptotic theory of RDDC in cylindrical geometry of Busse (2002) In §5, we perform nonlinear simulations to study the rotating finger convection (i.e. for a destabilising compositional gradient and a stabilising thermal one). In §6, we predict the onset of RDDC for core conditions and discuss the geophysical implications. Finally, we end the paper in §7 with a conclusion and outline several perspectives for geophysical and astrophysical bodies.

2 Description of the problem

2.1 Dimensional background state

We model RDDC in planetary cores by studying thermal and compositional Boussinesq convection in a rotating sphere. We consider a full sphere of radius $R$, filled with an homogeneous incompressible Newtonian fluid of density $\rho$, molecular kinematic viscosity $\nu$, thermal diffusivity $\kappa_{T}$ and compositional diffusivity $\kappa_{C}\ll\kappa_{T}$. The fluid is co-rotating with the sphere at the angular velocity $\boldsymbol{\Omega}=\Omega_{s}\boldsymbol{1}_{z}$ in the inertial frame. The fluid is also stratified in density under the (dimensional) imposed gravitational field $\boldsymbol{g}=-g_{0}r\,\boldsymbol{1}_{r}$, where $g_{0}R$ is the dimensional value of the gravity field at the outer spherical boundary $r=R$ and $\boldsymbol{1}_{r}$ is the unit radial vector in spherical coordinates $(r,\theta,\phi)$.

Within the Boussinesq approximation (Spiegel & Veronis, 1960), variations of the density $\rho^{*}$ due to the (dimensional) temperature $T^{*}$ and concentration of light elements $C^{*}$ are only taken in the buoyancy force. We use the following linear equation of state

[TABLE]

by assuming $|\rho^{*}-\rho_{m}|/\rho_{m}\ll 1$, where $(T_{m},C_{m},\rho_{m})$ are the mean reference values at $r=R$ and $(\alpha_{T},\alpha_{C})$ are the thermal and compositional expansion coefficients. In equation of state (1), $T^{*}$ is actually the departure from the adiabatic reference temperature profile. Similarly, $C^{*}$ is the departure from the compositional reference barodiffusive profile (Davies & Gubbins, 2011), which is rather small compared to the adiabatic density profile (Gubbins et al., 1979, 2004).

We work in the co-rotating reference frame. We study slight departures from a motionless, hydrostatic background state for the temperature $T_{0}^{*}$ and composition $C_{0}^{*}$. The latter profiles are governed by the dimensional temperature and composition equations in the Boussinesq approximation

[TABLE]

with $\mathcal{Q}_{T}$ and $\mathcal{Q}_{C}$ the thermal and compositional source (or sink) terms.

Thermo-compositional convection is sustained by the thermal and compositional gradients ($\nabla T_{0}^{*},\nabla C_{0}^{*}$). They can be maintained by (i) non-zero internal sources/sinks $(\mathcal{Q}_{T},\mathcal{Q}_{C})$, (ii) thermal or compositional fields externally imposed at the boundary or (iii) flux conditions. In the Earth’s core, the thermal gradient is mainly imposed by heat extracted at core-mantle boundary (CMB), yielding flux conditions. The compositional gradient is presently mainly driven by the crystallisation of the solid inner core (e.g. Loper & Roberts, 1981), while, in the early Earth, it may have been driven by the precipitation of light elements at the top of the core (O’Rourke & Stevenson, 2016; Badro et al., 2016). Hence, flux-type conditions are more relevant for compositional effects. Actually, the proper boundary condition ties the heat flux and the compositional flux to the local core dynamics (Braginsky & Roberts, 1995). This intricate condition has only been implemented in the anelastic simulations of Glatzmaier & Roberts (1996), who also treated separately thermal and chemical buoyancies. Yet, they assumed identical turbulent diffusivities, which discards double-diffusive effects.

However, the choice of the boundary conditions is less crucial for the dynamics in the full sphere geometry (investigated here) than in spherical shells (Kutzner & Christensen, 2000; Hori et al., 2012). To ensure stationary solutions, we assume that thermal and compositional background profiles are sustained by spatially homogeneous sources $(\mathcal{Q}_{T},\mathcal{Q}_{C})$. Hence, the dimensional solutions of equations (2) are

[TABLE]

Without loss of generality, we set $(T_{m},C_{m})=(0,0)$, because they do not play any dynamical role (only the gradients do have a role).

2.2 Dimensionless governing equations

For numerical convenience, we work with dimensionless quantities. We use the length scale $R$, the viscous time scale ${R^{2}}/{\nu}$, the temperature scale $(\nu\mathcal{Q}_{T}R^{2})/(6{\kappa_{T}}^{2})$ and the composition scale $(\nu\mathcal{Q}_{C}R^{2})/(6{\kappa_{C}}^{2})$. Note that temperature and composition scales can be either positive or negative, depending on the signs of $(\mathcal{Q}_{T},\mathcal{Q}_{C})$. In the following, we write the dimensionless velocity, temperature and composition without asterisk to differentiate them from their dimensional counterparts. In dimensionless form, dimensional background state (3) yields

[TABLE]

with

[TABLE]

the Prandtl number, the Schmidt number and the Lewis number.

In the co-rotating frame, we assume centrifugal effects to be small compared to the self-gravity of the fluid sphere $\boldsymbol{g}$. This condition is typically met in planetary cores, such that we neglect the centrifugal buoyancy in the Boussinesq equations (Lopez et al., 2013). We denote $\boldsymbol{u}$ the dimensionless velocity field, $\Theta$ the dimensionless temperature and $\xi$ the dimensionless concentration departing from motionless background state (2). The governing dimensionless equations are

[TABLE]

with $\boldsymbol{u}$ the dimensionless velocity field, $p$ the dimensionless reduced pressure (including the centrifugal force). In equations (6), we have introduced the Ekman number

[TABLE]

the thermal and compositional Rayleigh numbers

[TABLE]

which can be either positive or negative, depending on the signs of $(\mathcal{Q}_{T},\mathcal{Q}_{C})$, Typical values of numbers $[Pr,Sc,L,Ek]$ are given in table 1.

Equations (6) are supplemented by boundary conditions (BC). At the outer spherical boundary modeling the CMB, the velocity field satisfies the non-penetration and no-slip boundary conditions in the co-rotating frame, i.e.

[TABLE]

For the thermal and compositional perturbations $(\Theta,\xi)$, we impose zero radial fluxes

[TABLE]

The above boundary conditions (10) are relevant for a planetary core, in which the CMB controls heat and compositional fluxes through the profiles $T_{0}(r)$ and $C_{0}(r)$. The fixed temperature flux models the heat flux exctracted by the mantle, while the fixed compositional flux models exsolution of light elements by the mantle (e.g. O’Rourke & Stevenson, 2016).

2.3 Brunt-Väisälä frequency

To compare heat and composition gradients, we introduce the total dimensional Brunt-Väisälä frequency $N$. The latter is defined in the Boussinesq approximation by (e.g. Bullen, 1975)

[TABLE]

The fluid is stably stratified in density when $N^{2}>0$, neutral when $N^{2}=0$ and unstably stratified when $N^{2}<0$. The total dimensional Brunt-Väisälä frequency characterising the background state, denoted $N_{0}$ in the following, is such that $N^{2}_{0}=N^{2}_{0,T}+N^{2}_{0,C}$ where

[TABLE]

are the thermal and compositional contributions. Solutions (3) show that positive values of $(\mathcal{Q}_{T},\mathcal{Q}_{C})$ (respectively negative) give destabilising (respectively stabilising) thermal and compositional gradients.

To compare the rotational effects with the stratification, a relevant quantity is the square of the Brunt-Väisälä frequency normalised by the fluid angular velocity $\Omega_{s}$. In dimensionless variables, it reads for the background state

[TABLE]

where $Ro_{c}=Ek\,(Ra_{T}/Pr+Ra_{C}/Sc)^{1/2}$ is the double-diffusive convective Rossby number. Formula (13) is illustrated in figure 1. Because of the quadratic radial dependence in the background state (3), the background Brunt-Väisälä frequency is linear in $r$ in our model. In pure thermal convection ($Ra_{C}=0,Ra_{T}>0$), $Ro_{c}$ is often employed as a proxy of the ratio between buoyancy and Coriolis forces (Gastine et al., 2016). In the strongly stratified regime, characterised by $Ro_{c}\gg 1$, the scaling properties become reminiscent to non-rotating convection, whereas turbulent rotating convection is expected when $Ro_{c}\ll 1$. Hence, we can expect a similar distinction between a strongly stratified regime of RDDC, when $|{N_{0}^{2}}/{\Omega_{s}^{2}}|\gg 1$ (i.e. $|Ro_{c}|\gg 1$), and a weakly stratified regime when $|{N_{0}^{2}}/{\Omega_{s}^{2}}|\ll 1$ (i.e. $|Ro_{c}|\ll 1$).

2.4 Numerics in spheres

We will employ the classical Eulerian description in spherical geometry to solve equations (6). So far, most Eulerian simulations of RDDC have neglected the distinction between thermal and compositional buoyancies. This lead to the co-density approach, first proposed by Lister & Buffett (1995) and Braginsky & Roberts (1995), in which the two components have the same diffusivities $\kappa_{T}=\kappa_{C}$. This assumption is widely used (e.g. Schaeffer et al., 2017) and is mostly motivated by simplicity and numerical convenience, reducing by one both the number of parameters and equations. The proposed justification is that these molecular diffusivities should be replaced by a turbulent one, accounting for the mixing by unresolved small-scale eddies. However, this assumption is highly questionable and only possibly valid for highly turbulent flows, as found for overturning convection (Nataf & Schaeffer, 2015). Additionally, it filters out double-diffusive effects.

Only few Eulerian codes have treated separately the two buoyant components in spherical geometry, by using pseudo-spectral methods (Glatzmaier & Roberts, 1996; Manglik et al., 2010; Net et al., 2012; Takahashi, 2014) or finite volumes (Breuer et al., 2010). Here, we use the linear code SINGE (https://bitbucket.org/vidalje/singe) and the nonlinear code XSHELLS (https://nschaeff.bitbucket.io/xshells/), which are both open-source codes. We have implemented in both codes the composition equation (6c) to account for double-diffusive effects. The SINGE code has been used for linear computations of waves (Vidal & Schaeffer, 2015) and convection onsets (Gastine et al., 2016; Kaplan et al., 2017) in spherical geometry. On the other hand, XSHELLS can simulate turbulent flows in several contexts (Schaeffer et al., 2017; Kaplan et al., 2017, 2018), scaling on thousands of cores, by using a domain decomposition in the radial direction (MPI and OpenMP standards). XSHELLS solves the dynamical equations with a second order time-stepping scheme and treats the diffusive terms implicitly, while the nonlinear terms are handled explicitly.

Both codes use a pseudo-spectral method, by describing the velocity field $\boldsymbol{u}$ with poloidal and toroidal scalars (e.g. Backus, 1986). Then, poloidal and toroidal scalars are expanded onto spherical harmonics $Y_{l}^{m}(\theta,\phi)$ of degree $l$ and azimuthal wave number $m$, truncated at $(l_{\max},m_{\max})$ in the simulations. Similarly, temperature $\Theta$ and composition $\xi$ are also expanded onto spherical harmonics. The two codes use second order finite differences in radius with $N_{r}$ points and spherical harmonic expansions provided by the fast SHTns library (Schaeffer, 2013). At the origin ($r=0$), geometric conditions are applied: scalar fields ($\Theta$, $\xi$) can have a non-zero value at the origin. Since it must be independent of $\theta$ and $\phi$, only shperical harmonic $l=0$ is allowed. Similary for vector fields that have a non-zero vector at the origin, only $l=1$ is allowed. These all translate into appropriate boundary conditions, that distinguish $l=0$, $l=1$ and $l>1$ and which are used in both codes. For nonlinear simulations, numerical instabilites can arise because of the clustering of points near the origin. In the XSHELLS code, these instabilities are suppressed by truncating the spherical harmonic degree at $l_{tr}(r)=1+\sqrt{r/r_{s}}\,l_{max}$ with $r_{s}=0.5$. The XSHELLS code passes benchmarks designed to highlight issues arising at the origin (Marti et al., 2014).

The typical spatial resolution at $Ek=10^{-5}$ is $N_{r}=192$, $l_{\max}=120$, $m_{\max}=110$. For the most demanding nonlinear simulations (at large $Ra_{C},Ra_{T}$), the numerical resolution is $N_{r}=384$, $l_{\max}=320$ and $m_{\max}=300$. For such simulations, we show in figure 2 typical instantaneous spectra of the volume average of kinetic, thermal and compositional energies defined by

[TABLE]

Spectra are numerically well converged. We have also integrated the dynamics over several viscous time units (to skip any possible transient) to ensure reliable numerical results.

3 Insights from local stability analyses

Composition and heat do not play a symmetrical role when $L\neq 1$. Several canonical situations occur and various local stability criteria have been devised for non-rotating fluids (e.g. Garaud, 2018). Although the spherical geometry is natural for planetary cores, ruling out boundary effects yields physical insights for the stability. We briefly apply them for the background state (3).

Pioneering stability criteria have been inferred for non-rotating, diffusionless stellar interiors. Ledoux (1947) obtained the stability criterion (in dimensional and dimensionless forms)

[TABLE]

Note that in the absence of compositional effects, Ledoux criteria (15) reduces to the Schwarzschild criterion (Schwarzschild & Härm, 1959)

[TABLE]

When the background state is both Schwarzschild ($N_{0,T}^{2}<0$) and Ledoux unstable ($N_{0}^{2}<0$), the fluid is prone to overturning convection driven by thermal and compositional buoyancies.

However, Ledoux and Schwarzschild criteria (15)-(16) are not sufficient when heat (rapid diffuser) and composition (slow diffuser) have opposite destabilising/stabilising effects. Actually, the stability of the system depends on the density ratio $R_{0}$ (Stern, 1960), given by

[TABLE]

in which the last estimate holds for our background state (3) at the outer boundary. When the fluid is Ledoux unstable, i.e. $R_{0}\leq 1$, the system is usually prone to overturning convection, but also sometimes to finger convection (Schmitt, 2011). When the fluid is stable according to Ledoux criterion (15) the situation depends on the values of $(Pr,Sc)$. On the one hand, the situation $N^{2}_{0,T}>0$ (i.e. $Ra_{T}<0$) and $N^{2}_{0,C}<0$ (i.e. $Ra_{C}>0$) refers to the finger regime. In addition to overturning convection for $R_{0}\leq 1$, the finger configuration is prone to double-diffusive instabilities when (Baines & Gill, 1969)

[TABLE]

In that case, several finger DDC patterns can develop. On the other hand, the situation $N^{2}_{T,0}<0$ and $N^{2}_{C,0}>0$ refers to the semi-convection regime (Spiegel, 1969). The fluid is prone to double-diffusive instabilities when (e.g. Radko, 2013)

[TABLE]

Based on typical values of dimensionless Lewis and Prandtl numbers given in table 1, we can expect many celestial fluid bodies to be unstable against double-diffusive convection according to criteria (17)-(19).

The aforementioned local criteria do not account for rotational effects. Because the background state (3) is spatially varying, we cannot directly use plane-wave perturbations usually employed in local analyses (e.g. Cébron et al., 2013). However, the spatial extent of a local model is much smaller than the size of the global domain. Hence, we can linearise (3) in first approximation around a given position $\boldsymbol{r}_{0}$ to use WKB-type perturbations of the form

[TABLE]

with $\boldsymbol{k}$ the local wave vector and $\lambda=\sigma+\mathrm{i}\omega$ the eigenvalue, where $\sigma\geq 0$ is the growth rate (or damping rate if $\sigma\leq 0$) and $\omega$ is the angular frequency. Note that perturbations (20) differ from WKB-type perturbations considered by Yano (1992) and Jones et al. (2000) for thermal convection. Indeed, the latter perturbations are exponentially decaying in the cylindrical direction around a given cylindrical radius (to fulfill the boundary conditions). After some algebra, this yields a polynomial equation for the eigenvalue $\lambda$, similar to the one obtained by Sengupta & Garaud (2018), valid at the local position $\boldsymbol{r}_{0}$ of colatitude angle $\theta$. Note that Braginsky (2006) obtained a similar polynomial but considered a truncated version of the Coriolis force. As first obtained by Sengupta & Garaud (2018), the local analysis shows that the aforementioned non-rotating criteria are asymptotically valid for weakly rotating RDDC. Moreover, it indubitably shows that fastest-growing unstable waves for local rotating finger convection are largely unaffected by rotation. The unstable waves span the height of the local domain, with typical wave numbers $\boldsymbol{k}\boldsymbol{\cdot}\boldsymbol{g}=0$, called elevator waves (e.g. Sengupta & Garaud, 2018). All other waves are merely stabilised by rotation. Moreover, the range of density ratios $R_{0}$ for which RDDC takes place is unchanged compared to non-rotating DDC, given by (18) in the finger regime.

However, this local behaviour may be misleading. Indeed, it is known that WKB-type local solutions do not necessarily provide approximations to the complete three-dimensional global solutions. For instance, they can severely differ for thermal convection in the limit $Ek\ll 1$ (Busse, 1970; Soward, 1977; Yano, 1992; Jones et al., 2000). Therefore, the local analysis, predicting unavoidably elevator modes as the fastest-growing modes, is likely inaccurate to describe the onset of RDDC for rapidly rotating cores, and we turn to a global stability analysis.

4 Global stability analysis

4.1 Generalised eigenvalue problem

In this section, we perform a global linear stability analysis of background state (3). To do so, we discard the nonlinear terms ($\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\boldsymbol{u},\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\Theta,\boldsymbol{u}\boldsymbol{\cdot}\boldsymbol{\nabla}\xi)$ in equations (6). The symmetries of the background state and the linearised equations leads to uncoupled families of modes. The axisymmetry implies all azimuthal wave numbers $m$ are uncoupled and can be considered separately. Similarly, the reflexion symmetry about the equatorial plane implies the same for symmetric ($s=1$) and anti-symmetric ($s=-1$) modes with respect to that plane. Thus, for a given $m$ and symmetry $s$, we expand the linear perturbations in spherical coordinates $(r,\theta,\phi)$ as

[TABLE]

where $\lambda=\sigma+\mathrm{i}\omega$ is the complex eigenvalue with the growth rate $\Re_{e}(\lambda)=\sigma$ and the angular frequency $\Im_{m}(\lambda)=\omega$. Substituting expansions (21) into equations (6) yields the generalised eigenvalue problem (in symbolic form)

[TABLE]

with $\boldsymbol{X}=[\widehat{\boldsymbol{u}},\widehat{\Theta},\widehat{\xi}]$ the state vector and $(\mathcal{A},\mathcal{B})$ two linear operators, associated with the left and right hand sides of equations (6) and taking into account boundary conditions (9)-(10). Problem (22) is a boundary value problem, giving the dispersion relation for the complex eigenvalue

[TABLE]

From relation (23), the linear onset of convection is defined by the marginal state $\sigma=0$, realized for a set of Rayleigh numbers $(Ra_{T},Ra_{C})$ for given values of $(m,s,Ek,Pr,Sc)$.

We use the SINGE code (Vidal & Schaeffer, 2015) to solve the generalised eigenvalue problem (22), by using an efficient sparse eigenvalue solver provided by the SLEPC library (Hernández et al., 2005). At the parameters of our study, we found that the onset of RDDC is systematically governed by equatorially symmetric ($s=1$) perturbations (i.e. they have a lower onset than anti-symmetric perturbations). This is similar to purely thermal convection in spheres (e.g. Busse, 1970; Jones et al., 2000) and RDDC with an inner-core at similar parameters (Net et al., 2012). Nevertheless, the antisymmetric modes may still play a role (Landeau & Aubert, 2011; Net et al., 2012), see here §5.

We survey dispersion relation (23) by fixing all parameters except one of the Rayleigh numbers $Ra_{X}$ (where $X$ can be $T$ or $C$), that we vary until the growth rate $\sigma=0$ is bracketed within a small tolerance. This is done automatically by the SINGE code using an optimization procedure based on Brent’s method. Having computed a collection of Rayleigh numbers $Ra_{X}$ at the onset for various azimuthal wave numbers $m$, we can usually define the critical number $Ra_{X}^{c}$ obtained for the critical wave number $m^{c}$, yielding the minimum Rayleigh number over all computed azimuthal numbers.

4.2 Marginal stability

4.2.1 Convection for unstably stratified fluids

We set $Sc=3$ and $Pr=0.3$, giving a Lewis number $L=10$, and report in table 2 the critical parameters at the onset of pure compositional (overturning) convection for $Ra_{T}=0$. As already noticed for pure thermal convection (Zhang, 1992; Jones et al., 2000), we report only a broad agreement between our global numerical results and local predictions at the onset (e.g. Busse, 1970). The critical Rayleigh numbers $Ra_{C}^{c}$ are typically under-estimated by a factor two in the local theories (compared to the numerics), whereas the critical wave number $m^{c}$ and the angular frequency $\omega$ are over-estimated (not shown).

Then, we investigate the stability in the presence of an additional stabilising thermal background, which we refer to as the finger regime ($Ra_{C}\geq 0,Ra_{T}\leq 0$). For many fixed $Ra_{T}<0$, we determine the critical value of the compositional Rayleigh number $Ra_{C}^{c}$, reported in figure 3 for three Ekman numbers $Ek=\{10^{-6},10^{-5},10^{-4}\}$. When $|Ra_{T}|\ll|Ra_{C}|$, the preferred modes of convection are almost that of a pure compositional convection, with an onset almost unchanged. Indeed, double-diffusive effects become significant only when $|Ra_{T}|\sim|Ra_{C}|$. This behaviour has also been observed in thick shells (Net et al., 2012).

4.2.2 Inviscid convection for stably stratified fluids

When $Ra_{T}<0$ is further decreased, double-diffusive effects start playing a significant role when the fluid is stably stratified in density. For some values of $Ra_{T}$, there are now three values of $Ra_{C}$ that give $\sigma=0$, and $\sigma$ does not evolve monotonically with $Ra_{C}$. The marginal stability curve $\sigma=0$ takes schematically the form of a tongue in the $Ra_{c}-Ra_{T}$ diagram (figure 3), stretching towards lower $Ra_{C}$ within the stably stratified domain. Within this tongue, convection occurs at $Ra_{C}$ and $m$ much lower than for $Ra_{T}=0$ (typically $m\leq 10$), down to $m=1$ near the edges. This effect gets more important as $Ek$ is reduced, as observed in figure 3. At $Ek=10^{-6}$, $Ra_{C}$ in the tongue can go down to 10 times lower than the minimum $Ra_{C}^{c}$ of pure chemical convection. Furthermore, the smaller the $m$, the lower $Ra_{C}$ is at the onset. Hence, the critical wave number $m^{c}$ severely drops, e.g. from $m^{c}=20$ to $m^{c}=1$ at $Ek=10^{-5}$. This contradicts local stability analyses (e.g. Sengupta & Garaud, 2018), which do not capture this puzzling double-diffusive behaviour. Indeed, the existence of the double-diffusive tongue is due to the interplay between global rotation and the bounded geometry, as outlined by Busse (2002). However, note that the onset of modes with large azimuthal wave numbers $m$ is almost unaffected by these effects, in agreement with the asymptotic limit of short-wavelength perturbations.

When $Ra_{T}$ is still further reduced, the critical $Ra_{C}^{c}$ increases again for all wave numbers. Ultimately, the stability curves for all azimuthal numbers $m$ collapse onto the asymptotic regime of non-rotating finger convection (18), i.e.

[TABLE]

However, we show in appendix A that limit (18) is not always valid in the sphere, depending on the thermal and compositional boundary conditions.

Because the edge of the tongue consists of a large-scale $m=1$ mode, we can expect being able to compute the onset with SINGE at the parameters of the Earth’s core. We remark that the tongue is stunningly invariant when plotted using inviscid dimensionless numbers, as shown in figure 4. We have also checked that $Pr$ and $Sc$ play only a role through the Lewis number $L$. These two observations prove the inviscid nature of the low Rayleigh number double diffusive convection. To our knowledge, this behaviour has not been noticed by previous authors, although it can be inferred from the theory of Busse (2002), see appendix C.

Furthermore, the tongue only weakly depends on the Lewis number when $L\gg 1$. Hence, the black curve displayed in figure 4, computed at $Ek=10^{-15}$, fully characterises the convection onset within a stably stratified sphere, for any Ekman number $Ek\leq 10^{-4}$. In particular, the lowest value of $Ra_{C}$ in this regime is given by $Ra_{C}\simeq 52\,Ek^{-1}$ for $Ra_{T}\simeq-26\,Ek^{-1}$. Because the viscous convection onsets at $Ra_{C}\sim Ek^{-4/3}$ (Busse, 1970; Jones et al., 2000), the Ekman number controls the transition between inviscid low-Rayleigh number convection and the standard viscous convection. Thus, the domain of existence of inviscid convection increases as $Ek^{-1/3}$.

This behaviour supports the possibility of convection in planetary cores at low Rayleigh numbers (compared to the ones for pure compositional rotating convection). However, unlike the suggestion of Busse (2002) who mistakenly considered the non-rotating limit, the unstable Rayleigh numbers are not reduced down to non-rotating values, but rather to $Ra_{C}\simeq 52\,Ek^{-1}$. Note that the correct behaviour $Ra_{C}\sim Ek^{-1}$ is actually present in the annulus model (see appendix C).

We also remark that these effects subsist with other boundary conditions, but the shape of the unstable tongue varies as shown in appendix A. Interestingly, the asymptotic limit from local theory $Ra_{T}=-Ra_{C}$ is not always relevant (as pointed out above). Finally, note that for the semi-convection quadrant ($Ra_{C}\leq 0,Ra_{T}\geq 0$ – reported in appendix B), we find a similar behaviour with almost no effect of small stabilising compositional gradients. However, for stably stratified fluids, the marginal curves $\sigma=0$ are significantly different, and should be studied in future work.

4.3 Eigenmodes at the onset

The rapid rotation does provide constraints on the velocity structure, not taken into account in local (unbounded) analyses. For instance in convective rotating spheres with the no-slip condition, flows approximately obey the Taylor-Proudman theorem (Greenspan, 1968). This constraint yields quasi-geostrophic (QG) columnar motions, almost invariant along the rotation axis $\boldsymbol{1}_{z}$, as recovered numerically by SINGE (Kaplan et al., 2017). Then, we show in figure 5 and 6 the spatial pattern of several eigenmodes at the onset of finger convection. They are representative of our linear numerical results, and do not depend much on the viscosity.

The eigenmode at the onset of almost pure compositional convection is shown in figure 5a. The flow is in the form of spiraling columnar rolls (Zhang, 1992), extending spirally from near latitude $60^{\circ}$ to the equatorial region. For this mode, the composition and temperature perturbations are phase-shifted by about $90^{\circ}$. Spiraling modes appear to be the preferred modes of convection for the moderate value $Sc=3$. However, in the limit $Sc\gg 1$, spiraling is expected to be small (Zhang, 1992; Guervilly, 2010). In figure 5b, we show a typical low-frequency mode ($m=60$) computed at $Ra_{C}=10^{9}=-Ra_{T}$. For this mode, the composition and temperature perturbations are indistinguishable. In that case, the critical Rayleigh number for all the modes are close, such that several modes are likely to be triggered in a slightly supercritical state.

Then, we show in figure 6 the $m=1$ mode at the onset within the double-diffusive tongue of figure 7. At the tip of this tongue ($Ek\,Ra_{C}\simeq 52$ in figure 6a) the mode is quite simple and spans the whole sphere and is almost stationary. Remarkably, the composition and temperature perturbations are phase-shifted by about $45^{\circ}$. The flow exhibits features reminiscent of quasi-geostrophy (columns aligned along the rotation axis). For stronger forcing ($Ek\,Ra_{C}\gtrsim 1000$ in figures 6b,c), the mode increases in complexity, with several zeros in the direction parallel to the rotation axis. There, it is no longer columnar and could not be captured by the quasi-geostrophic approach (Busse, 2002; Simitev, 2011).

5 Nonlinear simulations of RDDC

5.1 Nonlinear onset

As illustrated in figure 7a, the linear global analysis predicts the existence of alternating stabilising and destabilising double-diffusive effects when increasing $Ra_{C}$ for a fixed $Ra_{T}$ at the upper edge of the inviscid tongue. We compare in figure 7b computations performed with SINGE and XSHELLS at $Ek=10^{-5}$, along the profile $Ra_{T}=-1.7\times 10^{6}$ shown in figure 7a. The growth rate computed with XSHELLS (during the exponential growth) is in perfect agreement with the eigenvalue computations.

Then, we aim at determining if this effect survives against finite-amplitude perturbations in nonlinear simulations. To do this, we have run the simulations sequentially for increasing value of $Ra_{C}$, and using the output of the previous simulation as initial state. Starting from a linearly stable background state, increasing $Ra_{C}$ first destabilises the system, leading to RDDC within the unstable tongue. Then, further increasing $Ra_{C}$ from a previous nonlinear state (at smaller $Ra_{C}$) abruptly inhibits the previously established RDDC when $Ra_{C}$ gets out of the tongue. This is counter-intuitive as restabilisation occurs even though the compositional profile has a priori a stronger destabilising gradient. Finally, overturning convection sets up again in the system for larger values of $Ra_{C}>3.4\times 10^{7}$. Similarly, we also find that the double-diffusive tongue subsists nonlinearly by varying $Ra_{T}$ at a fixed $Ra_{C}$ (not shown).

Thus, we have shown that this double-diffusive tongue is a linear mechanism, that persists against nonlinear perturbations of finite amplitude. We have found no evidence from the numerics that RDDC may onsets through a subcritical bifurcation, as recently obtained in pure thermal convection at much lower $Ek$ and $Pr$ (Kaplan et al., 2017).

5.2 Double-diffusive structures

In the following, we have conducted nonlinear simulations for stably stratified fluids along the profile $Ra_{T}=-Ra_{C}/3$ shown in figure 7 as a diagonal line. Along this profile, the density ratio (17) is kept constant $R_{0}=L/3$ but the background Brunt-Väisälä frequency increases according to formula (13). Note that we have also performed non-linear simulations in the semi-convection quadrant, as briefly discussed in appendix B.

Within the double-diffusive tongue, for typical compositional Rayleigh numbers $Ra_{C}\leq 10^{8}$, the nonlinear solutions are reminiscent of the eigenmodes at the linear onset (not shown). However, for higher Rayleigh number ($Ra_{C}\geq 10^{8}$), many high-wavenumber modes are unstable (see figure 7), leading to extremely thin convection fingers, elongated in the direction of the rotation axis due to the rapid rotation (figure 8).

In non-rotating systems, finger DDC leads to spatial scales intrinsically governed by the fast (thermal) diffusion and viscosity (e.g. Radko, 2013). Recently, Bouffard (2017) proposed another empiric scaling law in the presence of rotation. These two scaling laws predict the typical length of density structures in the equatorial plane $l_{\perp}$. They read respectively in the non-rotating and rotating regimes (with our variables)

[TABLE]

in which the rightmost forms involving $Ra_{C}$ are only valid for profiles characterised by $Ra_{T}\propto Ra_{C}$. Note that scalings (25a)-(25b) are expected for large enough values of the Lewis number. In addition, the typical horizontal size of the fingers is reasonably well approximated by prediction (25a) in the non-rotating case, even for moderate values of $L$. Indeed, relation (25a) holds for local computations at $L=3$ (see figure 7a of Traxler et al., 2011).

We assess their relevance for RDDC against 3D simulations performed at the finite value of $L=10$ in figure 9. We have determined the approximate number of fingers in the equatorial plane to estimate $l_{\perp}$. We observe two regimes, with a transition between $Ra_{C}\simeq 5\times 10^{8}$ and $Ra_{C}\simeq 1.5\times 10^{9}$. Our measurements do not seem to be in obvious agreement with the previous scaling laws, but the decrease of $l_{\perp}$ with increasing $Ra_{C}$ slows down at the transition, as predicted. The transition occurs for the Brunt-Väisälä frequency $N_{0}/\Omega_{s}\simeq 0.5$, and will be seen in several other diagnostics in the following (see below). We did not test the dependence of $l_{\perp}$ with the Ekman number $E$, which is predicted by eq. 25b. This would require to reduce the Ekman number, and run several high-resolution simulations at the edge of what is feasible.

For the simulations in the strongly stratified regime at $Ra_{C}\geq 10^{9}$, we may look for density staircases (Stern & Turner, 1969). The latter are made of stacks of well-mixed convective layers, separated by stably stratified shells for the total density profile (e.g. Stellmach et al., 2011). However, we have not found any evidence of density staircases in our simulations. In the non-rotating regime, Brown et al. (2013) found that local simulations performed at low values of the reduced density ratio $\widetilde{R}_{0}\ll 0.01$ exhibit properties consistent with density layering, with

[TABLE]

The finger regime is mapped into $0\leq\widetilde{R}_{0}\leq 1$. We have $\widetilde{R}_{0}\sim 0.2$ in our simulations, such that the absence of staircases is expected even for non-rotating fluids. Thus, performing more turbulent simulations, at lower values of $\widetilde{R}_{0}\ll 1$, appears necessary to investigate the interplay between rotational effects and density staircases.

5.3 Turbulence and transport

We now focus on specific features of finger convection in the turbulent regime. To quantify the nonlinear outcome, we compute in figure 10 the root mean square (rms) Reynolds and Rossby numbers

[TABLE]

with $V=4\pi/3$ the dimensionless spherical volume and $E_{u}$ the kinetic energy defined by formula (14). We have used the time average of $E_{u}$ in the saturated regime to determine the rms velocity. We have also separated $Re$ and $Ro$ based on total and non-zonal poloidal energies, to illustrate several regimes of finger convection.

First, when $N_{0}/\Omega_{s}\lesssim 0.5$, the Reynolds numbers based on total and radial velocities both exhibit the same scaling $Re\propto Ra_{C}^{0.93}$. However, when $N_{0}/\Omega_{s}\gtrsim 0.5$, another regime appears. Although $Re$ based on the total velocity is still nearly proportional to $Ra_{C}$, the scaling of $Re$ based on the poloidal energy is suddenly altered for $Ra_{C}\geq 10^{9}$, yielding $Re\propto Ra_{C}^{0.24}$. Hence, for strong stratification, radial (poloidal) motions are inhibited, while toroidal ones are not. This behaviour is consistent with scaling arguments and simulations of sustained stratified turbulence (Billant & Chomaz, 2001; Brethouwer et al., 2007). Indeed, a transition is expected between two turbulent regimes, characterised by strong and weak radial (here poloidal) motions. Such a dichotomy has been also evidenced in pioneering global simulations of tidally driven stratified flows (Vidal et al., 2018).

We can compare our results with the unbounded RDDC recently studied by Sengupta & Garaud (2018). They find the local Reynolds number $Re_{l}$, based on the convective velocity (analog to our non-zonal poloidal energy), to scale as

[TABLE]

In our case, this formula gives a constant value of $\sim 1$, in apparent contrast with the evolution of $Re$ shown in figure 10a. However, $Re$ is based on the (global) radius of our sphere, which is not relevant to estimate a local Reynolds number. Using rather the local length $l_{\perp}$, we estimate $Re_{l}=l_{\perp}\,Re$ in our simulations. As shown in figure 10b, we then obtain a constant $Re_{l}\simeq 1$ in agreement with formula (28). In both regimes $N_{0}/\Omega_{s}\leq 0.5$ and $N_{0}/\Omega_{s}\geq 0.5$, we thus recover the behaviour observed in unbounded RDDC. This can be understood with the following physical argument. Finger convection works because during the motion of a fluid particle, the temperature can be exchanged with its surroundings. Hence, the thermal diffusion time-scale $l_{\perp}^{2}/\kappa_{T}$ must not be smaller than the advection time-scale $l_{\perp}/u$. This leads to the condition $u\,l_{\perp}/\kappa_{T}\lesssim 1$, that is the Péclet number is of order one. This also translates into $Re_{l}\lesssim 1/Pr$, which is consistent with our findings (fig. 10b, $Re_{l}\simeq 1$ independent of $Ra_{C}$). Note that, because we have set $Pr=0.3$, the two predictions cannot be distinguished.

We now turn to the efficiency of convective transport of temperature and composition, which are quantified by the Nusselt $Nu_{T}$ and Sherwood $Sh$ numbers respectively. Their value is 1 for pure diffusion, and increase with increasing convection strength. In a convective sphere with internal sources and fixed flux at the outer boundary, they are given by

[TABLE]

with $[T_{0},C_{0}](r)$ the dimensionless background profiles (4) and the rms values of temperature and compositional perturbations $[\Theta_{rms},\xi_{rms}](r)$, defined from thermal and compositional energies at the radius $r$. In figure 11, we observe that $Nu_{T}$ is only weakly affected by varying $Ra_{C}$, yielding $Nu_{T}-1\leq 10^{-2}$. This is in agreement with local models of non-rotating finger convection. Indeed, Brown et al. (2013) showed that $Nu_{T}$ is always low and drops to 1 as $L$ (or $R_{0}$) is increased. This shows that the significant thermal diffusion necessary for finger patterns to develop always dominates the heat transport.

The compositional Nusselt number exhibits more significant variations. When $Ra_{C}$ increases in the regime $N_{0}/\Omega_{s}\lesssim 0.5$ defined above, $Sh$ increases up to $Sh\sim 2$. Thus, the turbulent compositional flux is enhanced, for a fixed $R_{0}$ and an increasing strength of the background stratification along the profile. Then, in the second regime ($N_{0}/\Omega_{s}\gtrsim 0.5$), increasing further $Ra_{C}$ does not yield significant changes in $Sh$.

Note that the scaling of the Nusselt and Sherwood numbers in figure 11, are in agreement with the laws of non-rotating finger convection (e.g. Garaud, 2018). Indeed, they predict constant Nusselt and Sherwood numbers for constant buoyancy ratio $R_{0}$ and Prandtl number $Pr$.

By contrast, the regime $N_{0}/\Omega_{s}\lesssim 0.5$ is more puzzling. One one hand, no clear scaling was observed for $Sh$, $Nu_{T}$ or $l_{\perp}$, but on the other hand the $Re\sim Ra_{C}$ scaling found in this regime has also been put forward in rotating thermal convection (Guervilly et al., 2019), also obtained at low $Pr$ and low $Nu$.

Figure 12 shows the evolution of the convective power for the same set of simulations, as a function of $Ra_{C}=-3Ra_{T}$. The convective power – the work done by buoyancy forces – is the sum of the thermal buoyancy power $P_{T}=\alpha_{T}(T_{0}^{*}+\Theta)\boldsymbol{g}\boldsymbol{\cdot}\boldsymbol{u}$ and the solutal buoyancy power $P_{C}=\alpha_{C}(C_{0}^{*}+\xi)\boldsymbol{g}\boldsymbol{\cdot}\boldsymbol{u}$. Three regimes may be distinguished here. Within the anomalous inviscid tongue ($Ra_{C}\lesssim 4\times 10^{-7}$, see Fig. 7a), the compositional buoyancy almost balances the thermal buoyancy. Only a small amount of net power drives convection. For larger $Ra_{C}$, the solutal buoyancy overcomes more easily the thermal stabilizing gradient, leading to more efficient convection. This coincides with the identification of small-scales fingers (see Fig. 9). Interestingly, for the largest forcings, we find that the net convective power $P$ evolves like $P\simeq 0.013Ra_{C}$. A scaling of $P$ proportional to $Ra_{C}$ is also reported in standard convection and convective dynamos (e.g. Christensen & Aubert, 2006), but with a proportionality constant close to one.

5.4 Zonal flows

For strong forcings, two different symmetries of zonal flow emerge which are discussed below.

5.4.1 Equatorial anti-symmetry for moderate forcing

Within the regime $N_{0}/\Omega_{s}\lesssim 0.5$, an equatorially anti-symmetric, differential rotation emerges from saturated quasi-geostrophic motions. A typical temporal evolution is summarised by figure 13a. Initially, the saturated nonlinear flow is dominated by quasi-geostrophic vortices (see figure 13b). They are associated with density fingers (thicker than the ones illustrated in figure 8, obtained at larger $Ra_{C}$). Columnar motions are predominant as long as $t\leq 2$ in the simulation (figure 13a). Then, an anti-symmetric flow grows on few viscous time units (figure 13a). Meanwhile, the energy of the equatorially symmetric flow is significantly reduced, such that the total energy of the fluid remains roughly constant. This anti-symmetric flow is mainly toroidal, consisting of a strong differential rotation: the flow is prograde in the Northern hemisphere and retrograde in the the Southern hemisphere (see figure 13c), and is associated with a segregation of both compositional and temperature anomalies in one hemisphere (not shown). This is the first report of such flows in finger convection. Actually, as shown in figure 14, they appear just above the linear onset for the equatorially anti-symmetric and axisymmetric (EAA) mode. By contrast, Landeau & Aubert (2011) found the appearance of the EAA mode in purely thermal convection much further above its onset. In our case, the EAA flow appearing in the nonlinear regime is clearly linked to the crossing of the associated linear threshold. It is quite unexpected that, far from the instability onset, a purely linear mechanism can explain the symmetry breaking of a nonlinear flow at saturation. This highlights the potential importance of linear modes even far from the global stability threshold in this systems. Furthermore, it emphasizes the importance of long simulations spanning several diffusion times.

5.4.2 Equatorial symmetry for strong forcing

Figure 14 shows that the EAA is overtaken by equatorially symmetric zonal flows for $N_{0}/\Omega_{s}\gtrsim 0.5$. This contrasts with Landeau & Aubert (2011), where the EAA mode increasingly dominates with forcing. Being mainly toroidal, the zonal flow does not affect the Reynolds number based on poloidal non-zonal energy shown in figure 10. The more $Ra_{C}$ increases, the larger the amplitude of this zonal flow, which quickly dominates all other components, as shown in figure 14. A typical kinetic energy spectrum is shown as a function of $m$ in figure 2b. The zonal $m=0$ component has an amplitude up to several orders of magnitude larger than the non-zonal components.

This zonal flow has a strong radial dependence, as illustrated in figure 15a. In the bulk (here at $r=0.5$), the zonal flow is prograde. However, it naturally exhibits multiple alternating prograde and retrograde jets at the outer spherical boundary, with rich dynamics (figure 15b).

These zonal flows may be seen as the manifestation in spheres of large-scale vortices (LSV) found in (unbounded) local simulations. LSV are conspicuous in local simulations of rotating finger convection in the polar regions (Sengupta & Garaud, 2018). They also appear in rotating semi-convection (Moll et al., 2017) and in rotating pure-thermal convection (e.g. Guervilly et al., 2014; Guervilly & Hughes, 2017). Julien et al. (2018) argued that the formation of zonal flows and jets is a robust feature resulting from an inverse energy cascade, provided that the flow is strongly anisotropic. In our simulations it is the zonal flows that allow rapid velocities ($Re\sim 1000$) to be reached by convection in stably stratified fluids.

6 Towards planetary core conditions

6.1 Linear onset in the early Earth

We have highlighted, using linear and nonlinear simulations, that rotation has surprising effects for rotating finger convection. Now, a complete quantitative picture of the onset of rotating finger convection is emerging for typical core conditions ($Ek=10^{-15}$ and $L=10^{4}$). Indeed, we remind the reader that double-diffusive effect are negligible except for $|Ra_{T}|\sim|Ra_{C}|$. Hence, we gather in table 3 the parameters of pure thermal or compositional convection, as predicted by the global theory of Jones et al. (2000). Moreover, thanks to the inviscid nature of the instability in the stably stratified region ($N_{0}^{2}>0$), we have already determined the onset of finger convection at Earth’s core conditions (see figure 4).

The scenario is illustrated in figure 16. For core conditions, we clearly observe the possibility of convection at reduced Rayleigh number, immensely facilitated by the stably stratified thermal profile. First, the wave number at the onset is strongly reduced within the tongue, yielding typical values $m\leq 10$. The growth rate increases with $Ra_{C}$ from a few to several thousands per viscous time-scale within the tongue.

Then, as expected, the critical Rayleigh number at the onset of pure compositional convection $Ra_{C}^{c}\sim 10^{21}$ (table 3) is orders of magnitude larger than the critical value at the upper edge of the double-diffusive tongue, i.e. $Ra_{C}^{c}\simeq 10^{17}$. The composition Rayleigh number $Ra_{C}$ is thus reduced by four decades for the early Earth by adding a stabilising temperature gradient.

6.2 Speculative estimates for the early Earth

To investigate the relevance of RDDC in the early Earth, we can use orders of magnitude arguments. They are presently highly speculative, due e.g. to the large modeling uncertainties. They will be certainly revised by future additional constraints, provided by mineral physics and thermal models of the Earth.

A typical estimate of the compositional Rayleigh number is

[TABLE]

where $\Delta\rho_{C}$ is the typical density yielding the compositional buoyancy (due to light elements) and $\rho_{m}$ the typical density of the core. Following Jones (2015), typical values are $R=3500$ km for the radius of the core, $Rg_{0}=10\,\text{m.s}^{-2}$ for the gravitational acceleration, $\nu=10^{-6}$ m2.s*-1*for the (molecular) kinematic viscosity and $Sc=10^{2}$ for the Schmidt number. The amount of light elements attributable to compositional sources is highly debated. The compositional gradient was likely destabilising in the early Earth, due to the exsolution or crystallisation of light elements in the core. The equilibration at high temperatures in the aftermath of giant impacts would be responsible for a small amount of magnesium to partition into the core, yielding the exsolution of light magnesium oxides in the core (Badro et al., 2016; O’Rourke & Stevenson, 2016). This mechanism is energetically efficient, since precipitating a layer of magnesium-bearing material with a typical thickness of 10 km above the CMB would be equivalent to crystallising the entire inner core (O’Rourke & Stevenson, 2016). Instead of invoking such singular events, Hirose et al. (2017) advocated the crystallization of silicon dioxide. Nonetheless, in the two scenarios, roughly the same mass of light elements is precipitated/crystallised. Thus, typical (speculative) upper bounds are $0.2$ wt % of precipitated magnesium-bearing minerals (see figure 2 of O’Rourke & Stevenson, 2016) or $0.4$ wt % of crystallised silicon dioxide (Hirose et al., 2017). Based on these two scenarios, we may consider the typical value $\Delta\rho_{C}/\rho_{m}=10^{-3}$ as an upper bound. Then, formula (30) yields the estimate $Ra_{C}/Sc\leq 10^{29}$ in the early Earth. This upper estimate is much larger than the critical values required at the onset (figure 16), typically $Ra_{C}^{c}\sim 10^{17}-10^{18}$ in the inviscid tongue and $Ra_{C}^{c}\sim 10^{21}$ for compositional convection (even without stabilising thermal effects). This suggest that the early Earth did undergo highly supercritical RDDC, either for unstably or stably stratified fluids.

The properties of convection would be certainly different in the two regimes. Hence, to argue in favour of one regime, we have to estimate the square of the total background Brunt-Väisälä frequency $N_{0}^{2}$. On the one hand, the compositional part is

[TABLE]

This gives the speculative estimate $-10^{-1}\leq{N_{0,C}^{2}}/{\Omega_{s}^{2}}\leq 0$. On the other hand, the presence of a thick, thermally stratified layer seems probable prior to the formation of the inner core if there was a sub-adiabatic heat flux at the top of the core $Q_{\text{cmb}}$. To our knowledge, there is no reliable agreement between thermal models of the Earth (e.g. Labrosse et al., 1997; Nimmo, 2015; Nakagawa, 2018). Therefore, we estimate a speculative upper bound from the difference between the total heat flux $Q_{\text{cmb}}$ and the adiabatic flux $Q_{a}$ at the CMB

[TABLE]

with $\mathcal{S}$ the surface of the outer core, $\alpha_{T}\sim 10^{-5}$ K*-1* the thermal expansion coefficient (Nimmo, 2015) and $k$ the thermal conductivity. The latter quantity is badly constrained (Williams, 2018), so does the thermal history of the Earth. Possible values are $40-160$ W.m*-1*.K*-1*. A broad range of values appears possible for $Q_{a}-Q_{\text{cmb}}$. Upper-bound estimates are presently a few TW, yielding the (highly) speculative bounds for a thermal stratification $0\leq N_{0,T}^{2}/\Omega_{s}^{2}\ll 1-10$. The upper bound values are very close to the plausible geophysical estimates of the thermal Brunt-Väisälä frequency at the top of the Earth in the present time (e.g. Labrosse et al., 1997; Buffett, 2014; Helffrich & Kaneshima, 2013). Consequently, $N_{0,T}^{2}+N_{0,C}^{2}$ may have been either positive or negative.

To sum up, the Early Earth may have been prone to either overturning convection (for unstably stratified fluids) or finger RDDC (for stratified fluids).

7 Conclusion

7.1 Summary

We have revisited rotating double-diffusive convection (RDDC) in planetary cores, by considering flows driven by buoyancy forces of thermal and compositional origins. We have studied RDDC with a Boussinesq model in a full sphere, with internal source and sink terms. We have separated thermal and compositional effects, to go beyond the codensity approach (Braginsky & Roberts, 1995; Lister & Buffett, 1995) commonly used in planetary simulations. We have mainly focused on the finger regime ($Ra_{C}\geq 0,Ra_{T}<0$), by considering stabilising thermal effects and destabilising compositional effects.

First, we have performed the linear stability analysis of background diffusive state (3) in the finger quadrant, by using a global (spherical) method. A global picture is now emerging. A quantitative proxy of the strength of rotational and stratified effects is the absolute value of the square of the dimensionless background Brunt-Väisälä frequency, i.e. the ratio $|N_{0}^{2}/\Omega_{s}^{2}|$. Overturning convection occurs for unstably stratified fluids ($N_{0}^{2}/\Omega_{s}^{2}\leq 0$). When overturning convection is controlled by rotational effects ($-1\leq N_{0}^{2}/\Omega_{s}^{2}\leq 0$), the onset is largely unaffected by double-diffusive effects when $|Ra_{T}|\ll Ra_{C}$ in the finger regime. Then, the linear spherical analysis recovers asymptotically the onset of non-rotating DDC in strongly stratified regime ($N_{0}^{2}/\Omega_{s}^{2}\gg 1$). On the other hand, it strongly differs in the other regime $N_{0}^{2}/\Omega_{s}^{2}\ll 1$ with local analyses. Indeed, local analyses predict that rotation has a simple stabilising effect, merely increasing the critical Rayleigh numbers at the onset. However, rotational effects are more subtle in the presence of double diffusion. Indeed, the global analysis shows that the linear onset of RDDC can occur for lower Rayleigh numbers for stably stratified fluids than for unstably stratified fluids. This phenomenon, first outlined by Busse (2002), is intrinsically due to rotational effects in the bounded spherical geometry. Therefore, they are filtered out by local models. The associated flows at the linear onset do not always take the form of quasi-geostrophic motions (aligned with the rotation axis), unlike in standard rotating convection (e.g. Zhang et al., 2007; Kaplan et al., 2017). In addition, for a specific combination of boundary conditions (namely fixed temperature and imposed composition flux), rotating double-diffusive convection surprisingly occurs for density ratios $R_{0}>L$, which is beyond the limit of non-rotating double-diffusive convection. In the finger regime, double-diffusive effects become preponderant only for stably stratified fluids ($N_{0}^{2}/\Omega_{s}^{2}\geq 0$). On the contrary, as discussed in appendix B, double-diffusive effects start playing a role even for unstably stratified fluids ($N_{0}^{2}/\Omega_{s}^{2}\leq 0$) in the semi-convection quadrant ($Ra_{C}\leq 0,Ra_{T}>0$).

Second, we have conducted high-resolution, nonlinear simulations for rotating stratified fluids ($N_{0}^{2}/\Omega_{s}^{2}\geq 0$) in the finger regime. Several nonlinear features have been obtained. Outside the DD tongue for large enough $Ra_{C}$, the flow structures (fingers) strongly differ from the linearly unstable tongue modes at the upper edge of the DD tongue. Moreover, we have identified a sharp transition outside the tongue in the rapidly rotating finger regime. This transition empirically occurs at $N_{0}/\Omega_{s}\simeq 0.5$ in the simulations, for the fixed value $Ek=10^{-5}$. In the first regime, the nonlinear flows exhibit equatorially anti-symmetric, large-scale zonal flows, which appears when the associated linear onset is crossed. In the second regime, strong equatorially symmetric zonal flows are sustained. The latter flows are reminiscent of the large-scale vortices found in local models of finger convection (e.g. Sengupta & Garaud, 2018). The turbulent properties, e.g. the output Reynolds or Nusselt numbers, are also significantly different in the two regimes. Notably, we have found scalings for the second regime that appear in broad agreement with the scalings proposed for local DDC.

Finally, we have succeeded in predicting the onset of RDDC numerically at core conditions, after noticing the inviscid nature of finger convection in the weakly stratified regime. We have shown that the combination of rotation and double-diffusive effects is strongly destabilising in the inviscid tongue for stably stratified fluids. The critical Rayleigh number is reduced by four decades for realistic core conditions. Then, we have crudely estimated the thermal and compositional stratification in the Early Earth. We support that it may have undergone highly turbulent RDDC, either in the overturning compositional convection (unstably stratified) or in the finger regime associated with strong zonal flows.

7.2 Perspectives

7.2.1 Discussion and improvements

A considerable amount of work remains to be done, e.g. to expand the surveyed parameter space and to refine the model. Further simulations are required to understand the nonlinear saturation of finger convection (figure 10), e.g. by varying $L$, $Ra_{T}$ and $Ra_{C}$. On the one hand, we have found that local scalings of non-rotating finger convection (Garaud, 2018) may qualitatively hold in the second rotating regime. Nonetheless, a more exhaustive numerical survey of the parameter space is required to assess their quantitative validity. Moreover, it remains an open question whether regimes of rotating thermal convection (Gastine et al., 2016) apply for RDDC, both for destabilising and stabilising density profiles. Therefore, this calls for assessing and possibly improving the scaling laws describing rotating convection in the presence of significant double-diffusive effects.

For numerical reasons, we have considered moderate values for the Lewis $L=10$ and Ekman $Ek=10^{-5}$ numbers in the nonlinear simulations. The value of $L$ is about two orders of magnitude smaller than the expected values in planetary cores. Larger values of $L$ may facilitate the generation double-diffusive structures. In particular, we have not found any density staircases (Stern & Turner, 1969), resulting from secondary instabilities. Several theories have been proposed in the non-rotating case (Stern & Turner, 1969; Radko, 2013). For the moderate values of $Pr$ characterising planetary cores, their generation may rely on the mixing by nonlinear internal waves (Garaud et al., 2015). Yet, these mechanisms remain to be confirmed in the presence of rapid rotation. Their existence may strongly affect the turbulent regime. Indeed, it has been shown that density staircases can increase the turbulent heat and compositional fluxes by several orders of magnitude (e.g. in oceanography Schmitt et al., 2005). Thus, the conditions of existence for density staircases in rotating finger convection remain unanswered and studying them deserves future work.

We have outlined that we cannot rule out RDDC in the Early Earth. Now, investigating the dynamo capability is necessary to assess the validity of the proposed mechanisms for the origin of the Early geodynamo (Badro et al., 2016; O’Rourke & Stevenson, 2016; O’Rourke et al., 2017; Hirose et al., 2017). The dynamo capability of rotating finger convection remains an open question. Typically, dynamo action requires $Rm>100$, where $Rm=RePm$ is the magnetic Reynolds number with $Pm=\nu/\eta$ the magnetic Prandtl number ($Pm\ll 1$ for cores) and $\eta$ the magnetic diffusivity. With this first study, we cannot establish scaling laws that would allow us to infer $Rm$ at core conditions. However, Fig. 10 shows that $Re$ can be large, possibly allowing large $Rm$ too. For large $Ra_{C}$, the flow organizes itself into strong large-scale zonal shears and weak small-scale fingers. Even though the radial velocity of the small-scale finger is small, the large-scale zonal shear is large. This situation could in principle sustain an $\alpha\omega$ dynamo, in which the large-scale shear is responsible for a so-called $\omega$-effect while the small-scale convection produces an $\alpha$-effect (e.g. Roberts, 1972). We have also checked that our flow displays a significant amount of helicity, an ingredient thought to be important to obtain an important $\alpha$-effect. From a numerical point of view, we reach $Re\sim 10^{3}$ in our simulations. In an $\alpha\omega$ dynamo context, the relevant magnetic Reynolds number would be the geometric mean $Rm^{*}$ of the $Rm$ based on the large-scale zonal flow and the $Rm$ based on the small-scale one (e.g. Roberts, 1972). According to Fig. 10a, this leads to $Rm^{*}\sim 200Pm$, potentially allowing dynamos for $Pm\gtrsim 1$.

Beyond the question of the dynamo capability, we can wonder about the strength of the generated magnetic field. In the case of simple convective dynamos (ie without double-diffusive effects), the field strength scales as $P^{1/3}$, where $P$ is the convective power (e.g. Christensen & Aubert, 2006). Despite the small values of $Nu$ and $Sh$ (see Fig. 11), we find significant buoyancy power in our simulations, scaling like $P\simeq 0.013Ra_{C}$ (see Fig. 12). This scaling is similar to the one found in standard convective dynamos (see Christensen & Aubert, 2006), differing only by the constant factor which is about 100 times smaller here. Assuming this scaling holds, we can expect strong magnetic fields to be generated, provided that $Ra_{C}$ is large enough. Nevertheless, the saturation of a dynamo driven by double-diffusive convection may behave differently. In addition, the large-scale zonal flows we have found in these simulations, which may persist for core conditions, are known to be important for the dynamo process in stratified interiors (e.g. Spruit, 2002) whereas it does not change much the radial transport (and thus the Nusselt and Sherwood numbers). Indeed, such zonal flows can sustain various hydrodynamic and magnetic instabilities (e.g Knobloch, 1982; Jouve et al., 2015). Hence, dynamo onset, field strength at saturation, and extrapolation to core conditions all require a future study of dynamo driven by double-diffusive convection in the turbulent rotating regime.

Recently, Guervilly & Cardin (2016) and Kaplan et al. (2017) found that the smooth (linear) onset of rapidly rotating thermal convection is replaced by (nonlinear) hysteresis cycles and subcritical behaviours, at small enough Ekman numbers. These effects may survive with double-diffusive effects in the overturning regime. Finger convection may also occur though a subcritical bifurcation when $L\gg 1$, as proposed for non-rotating stratified fluids in planar models (Veronis, 1965; Proctor, 1981). This mathematical observation has not been confirmed yet numerically. Notably, we have not found evidence supporting this behaviour in the numerics. However, these nonlinear effects may only appear for $L$ larger than in our simulations. Therefore, studying finite-amplitude perturbations appears of special interest to investigate the transition towards turbulence in RDDC when $L\gg 1$.

Finally, we have neglected so far several double-diffusive effects occurring in a binary mixture. More relevant compositional boundary conditions may be implemented, e.g. the intricate boundary condition proposed by Braginsky & Roberts (1995); Glatzmaier & Roberts (1996). Investigating additional binary effects in the thermal and heat fluxes is also worthy of interest (still in the Boussinesq approximation). They are only responsible for second order effects at the linear onset (e.g. Hort et al., 1992; Net et al., 2012), when a background state state is imposed. However, they may play a dynamical role in nonlinear simulations. For instance, barodiffusion is the tendency of light material to migrate down the pressure gradient. Barodiffusion sustains the accumulation of light elements at the top of the core (Gubbins & Davies, 2013), to naturally increase the Brunt-Väisälä frequency. Handling barodiffusion is not demanding numerically, e.g. in shells by considering a system forced by the boundaries (i.e. no background state) but with an additional mass sink (e.g. Davies & Gubbins, 2011; Bouffard, 2017). These effects should be considered for consistent future nonlinear simulations.

7.2.2 Towards planetary applications and beyond

Beyond the origin of the early geodynamo, the (possible) outermost stable stratification in the Earth’s core is another long standing geophysical issue (e.g. Loper & Roberts, 1981; Braginsky, 1993; Lister & Buffett, 1998). The existence of such a layer has been outlined by seismological (Helffrich & Kaneshima, 2010, 2013; Irving et al., 2018), geodetic (Buffett & Seagle, 2010) and geomagnetic (Gubbins, 2007; Buffett, 2014) data. The density stratification may have a thermal and/or compositional origin (e.g. Buffett & Seagle, 2010; Davies et al., 2018; Nakagawa, 2018; Bouffard et al., 2019). Indeed, the thermal conductivity has been revised upward by ab-inito calculations (Pozzo et al., 2012; de Koker et al., 2012; Pozzo et al., 2013) and experiments (Gomi et al., 2013; Ohta et al., 2016; Konôpková et al., 2016). This may favour an outer sub-adiabatic thermal stratification, but large thermodynamical uncertainties remain (Williams, 2018). Moreover, Mound et al. (2019) pointed out that this outermost stratification may be regional (rather than global), being generated by the lateral variations in heat flux at the core–mantle boundary . Stratification may be also sustained by the accumulation of light elements (e.g. Loper & Roberts, 1981). This stratified layer may affect the geodynamo (e.g. Olson et al., 2017; Christensen, 2018), e.g. by filtering small-scale internal convective motions (Vidal & Schaeffer, 2015) or trapping waves (Knezek & Buffett, 2018). However, this hypothetical layer may be prone to either rotating finger convection or semi-convection (Braginsky, 2006), making the internal core dynamics more complex. In particular, intense zonal flows could develop, as we have found in this work. Partially stratified core layers may also exist in other planets, e.g. Mercury (Manglik et al., 2010; Takahashi et al., 2019) or Venus (Jacobson et al., 2017). Therefore, it is of special interest to determine whether thermally and/or compositionally stably stratified layers can survive dynamically against RDDC.

In addition, double-diffusive effects are also relevant for giant planets (Stevenson, 1982), such as Saturn (Stevenson & Salpeter, 1977; Leconte & Chabrier, 2013) and Jupiter (Moll et al., 2017). Stellar interiors may also undergo DDC (Garaud, 2018), e.g. low-mass hosting exoplanets (Vauclair, 2004) or massive stars (e.g. Merryfield, 1995; Woosley et al., 2002). Even though they were largely neglected, rotational effects may be significant in these objects, e.g. for the giant planets of our Solar system which are rapidly rotating (9.9 hr for Jupiter and 10.7 hr for Saturn) or for some radiative stars (e.g. Jouve et al., 2015).

The validity of the Boussinesq model for compressible interiors should be assessed. The scalings for the typical length scale of density structures, applied to planetary Earth-like parameters, yield (Bouffard, 2017) $l_{\perp}\sim 20$ cm for rapid rotations cm and $l_{\perp}=40$ cm in the non-rotating case. Spiegel & Veronis (1960) showed that the Boussinesq approximation is relevant for dynamical scales smaller than the pressure scale height, typically one-tenth of the radius of stars. Therefore, the compressible dynamics may be surprisingly well described by using the Boussinesq approximation, as advocated in the non-rotating regime (Radko, 2016). A comparison between Boussinesq and anelastic models of RDDC (e.g. Glatzmaier & Roberts, 1996) is certainly worthy of interest for astrophysical objects.

In addition, gaseous planets would require to consider stress-free conditions for the flow. Our results show that, in the limit $Ek\ll 1$, stress-free conditions do not affect the onset of inviscid RDDC, which remains symmetric with respect to the equatorial plane. However, these bodies are characterised by much smaller values of $Pr\ll 1$ (compared to planetary cores). In this regime, flows at the onset can be equatorially anti-symmetric torsional modes. They sometimes appear as the preferred unstable modes of (pure) thermal convection in spheres in the limit $Pr\ll 1$ (e.g. at $Pr/Ek=10$), but only for stress-free conditions (Sánchez et al., 2016; Zhang et al., 2017) as commonly used for giant planets and stars. Moreover, polar anti-symmetric modes have also been found at the onset when $Pr\ll 1$, for (pure) thermal convection in thick (Garcia et al., 2008) and thin (Garcia et al., 2018) spherical shells. The nonlinear regime in the low-$Pr$ regime is expected to differ from the high-$Pr$ regime (e.g. in the non-rotating regime Garaud, 2018). Therefore, studying RDDC in the low-$Pr$ regime with stress-free conditions may lead to different double-diffusive effects than those previously obtained in shells (e.g. Net et al., 2012).

Finally, we remark that the large-scale inviscid mode in the stably stratified regime is always $m=1$, with a net flow at the center within the equatorial plane. Such a mode could constrain the translation direction of a freshly-nucleated inner core to be perpendicular to the rotation axis, in agreement with seismological observation of the hemispherical dichotomy of the inner core (see e.g. Deguen, 2012).

acknowledgments

This project was funded by ANR-14-CE33-0012 (MagLune). JV acknowledges the support of STFC Grant ST/R00059X/1. We thank the geodynamo team (Univ. Grenoble Alpes), R. Deguen and T. Gastine for fruitful discussions and comments. Computations were performed on the Froggy platform of CIMENT (https://ciment.ujf-grenoble.fr), supported by the Rhône-Alpes region (CPER0713 CIRA), OSUG2020 LabEx (ANR10 LABX56) and EquipMeso (ANR10 EQPX-29-01). ISTerre is part of Labex OSUG2020 (ANR10 LABX56). All figures were produced using matplotlib (http://matplotlib.org/) or paraview (http://www.paraview.org/).

Appendix A Other boundary conditions at the linear onset

We investigate the effects of different mechanical, thermal and compositional boundary conditions (BC) on RDDC in spheres. We substitute no-slip conditions (9) by stress-free conditions for the velocity field

[TABLE]

with $\boldsymbol{\mathcal{E}}=\left[\boldsymbol{\nabla}\boldsymbol{u}+(\boldsymbol{\nabla}\boldsymbol{u})^{\top}\right]/2$ the strain rate tensor (incompressible Newtonian fluid). Instead of fixed flux conditions (10), we consider fixed temperature or composition at the boundary

[TABLE]

Numerical results, computed with SINGE, have been performed for $m=1$ and $m=6$ at $Ek=10^{-4}$ and $m=1$ at $Ek=10^{-11}$. Given that the results lead to the same conclusions, we only show the results for $m=1$ and $Ek=10^{-11}$ in figure 17. Within the stable double-diffusive tongue given by the Ledoux criterion (15), the linear onset is independent of the mechanical conditions. For the low Ekman number considered here, using stress-free (33) or no-slip condition (9) leads to the same marginal stability curve (not shown). However, changing the boundary condition on the temperature or composition field has important effects on the shape of the marginal stability curve, but the latter still remains independent of viscosity. Surprisingly, with a fixed temperature and imposed buoyancy flux, the double-diffusive convection extends to $Ra_{T}<-Ra_{C}$, which corresponds to density ratios $R_{0}>L$. This linear instability, located beyond the expected range of finger convection, has been confirmed by time-stepping nonlinear simulations with XSHELLS (at $Ek\,Ra_{T}=-10^{3}$, $Ek\,Ra_{C}=6\times 10^{2}$ and $Ek=10^{-5}$).

Appendix B Semi-convection

The onset of RDDC in the semi-convection quadrant ($Ra_{T}>0,Ra_{C}<0$) is represented in figure 18 the linear computations at the onset computed with SINGE, for two values of $Ek$. The critical parameters at the onset of pure thermal convection are given in table 4, for completeness with table 2 for pure compositional convection. The onset of convection is largely insensitive to double-diffusive effects as long as $|Ra_{C}|\ll Ra_{T}$. This refers to the overturning regime of thermal convection. For higher $|Ra_{C}|$, double-diffusive effects start to be important when $|Ra_{C}|\geq Ra_{T}$. As in the finger regime, the marginal stability curve $\sigma$ takes the form of a tongue in the $Ra_{c}-Ra_{T}$ diagram (figure 3). However, double-diffusive effects become significant even for unstably stratified fluids ($N_{0}^{2}<0$), as opposed to the finger quadrant in which only stably stratified fluids ($N_{0}^{2}\geq 0$) are strongly affected. Within this tongue, modes with small azimuthal wave number are triggered at the onset, which also occurs for smaller thermal Rayleigh number than in the overturning regime. In the limit $|Ra_{C}|\to\infty$, RDDC reaches asymptotically the non-rotating regime predicted by formula (19). Then, we show in 19b an illustrative nonlinear simulation of semi-convection at $Ek=10^{-5}$ and $Ra_{T}=10^{8}$. Density structures exhibit larger spatial scales than the ones obtained in simulations within the finger regime (for similar absolute values of the Rayleigh numbers).

Appendix C Revisiting the annulus geometry

C.1 Mathematical formulation

We revisit the model of RDDC in a cylindrical annulus. A few misprints are present in Busse (2002), which also used other dimensionless variables. Furthermore, Busse (2002) made wrong assumptions when drawing his conclusions, mistakenly considering the non-rotating limit. Before taking the annulus model further, we clearly explain the theory, going through the derivation of the equations in our formalism.

For the sake of tractable analytical developments, Busse (1970) pointed out that a simplified model of QG convection in spheres should consider a thin cylindrical annulus, with sloping top and bottom boundaries. Using this asymptotic model, he investigated the onset of thermal convection with $Ra_{C}=0$ (Busse, 1986), and extended it to RDDC (Busse, 2002). This model considers a thin-gap geometry centered on the QG columns at the onset. Moreover, this asymptotic theory can embrace core conditions in the limit $Ek\ll 1$ and $L\gg 1$. The annulus geometry is illustrated in figure 20. We consider the cylindrical annulus region, located at the cylindrical radius $s_{0}$ in a full sphere rotating at the angular velocity $\Omega_{s}$. We use the small-gap approximation, by assuming $d/s_{0}\ll 1$. Thus, the effects of the spherical curvature can be neglected and we use the Cartesian coordinate system of unit vectors $(\boldsymbol{1}_{x},\boldsymbol{1}_{y},\boldsymbol{1}_{z})$ centered at $s_{0}$. The annular channel is bounded at top and bottom by rigid conical caps with the angle of inclination $\chi$. We denote $h$ the half-height of the cylindrical annulus (with respect to the equatorial plane). In the background state, the fluid is stratified in temperature and composition under the inward gravity field $\boldsymbol{g}=-g_{0}s_{0}\boldsymbol{1}_{x}$, which is constant at the scale of the annulus. The inner wall (respectively the outer one) is kept at the constant temperature $T_{0,1}^{*}$ and composition $C_{0,1}^{*}$ (respectively $T_{0,2}^{*}$ and $C_{0,2}^{*}$).

We choose the gap $d$ as length scale, $d^{2}/\nu$ as time scale, $\Delta T^{*}\,Pr=(T_{0,1}^{*}-T_{0,2}^{*})\,Pr$ as thermal scale and $\Delta C^{*}\,Sc=(C_{0,1}^{*}-C_{0,2}^{*})\,Sc$ as compositional scale. These thermal and compositional scales are the local analogues of the global scales chosen in the main text. Dimensionless variables are denoted in the following without asterisk. We assume that the slope $\chi$ of the upper and lower caps shown in figure 20 is small ($\chi\ll 1$), such that the local conductive background state is close to the one in the annulus of uniform depth (e.g. Busse, 1970). Hence, the dimensionless background state is

[TABLE]

Then, the local form of equations (6) for the dimensionless perturbations $(\boldsymbol{u},\Theta,\xi)$ takes the form

[TABLE]

We have introduced in equations (36) the local Ekman and Rayleigh numbers

[TABLE]

Note that Rayleigh numbers (37) are the local versions of the spherical Rayleigh numbers (8) introduced in the main text.

We seek velocity solutions of equations (36) with small variations along the rotation axis $\boldsymbol{1}_{z}$. Hence, the velocity takes the form of QG flows

[TABLE]

with $u_{z}$ the small vertical velocity (at the order $\chi$) and $\Psi$ the velocity stream function in the equatorial plane $(z=0)$. The linear onset given by equations (36) can be solved by considering stress-free, iso-thermal and iso-compositional boundaries

[TABLE]

We also assume that the upper and lower conical caps at $z=\pm h/(2d)$ are rigid, with fixed vertical thermal and compositional fluxes (Busse, 1986). This yields

[TABLE]

Other conditions are irrelevant in the analysis. In particular, the neglected viscous boundary layer vanishes in the limit $Ek\ll 1$ in the annular geometry (Hunter, 1967). This is in agreement with the observation that the viscous boundary condition is of second order importance for the onset of convection in spheres (Zhang & Jones, 1993; Jones et al., 2000), at least for not too small values of $Pr$ at fixed $Ek$ (Zhang et al., 2017). Although these boundary conditions are not physically realistic (Braginsky & Roberts, 1995), they do not hinder from investigating the leading order double-diffusive effects. Then, following Busse (2002), we take the $z$-component of the curl of momentum equation (36a) and average it over $z$ (from bottom to upper caps). This yields at first order in $\chi$ (see Busse, 1986, for the derivation)

[TABLE]

with the two-dimensional horizontal Laplacian $\Delta_{\perp}=\partial^{2}/\partial x^{2}+\partial^{2}/\partial y^{2}$ and the parameter

[TABLE]

In the rapidly rotating limit $\widetilde{Ek}\ll 1$, $\beta$ is a leading order parameter containing the effects of the boundary curvature (the so-called $\beta$-effect).

We assume periodicity in the $\boldsymbol{1}_{x}$ direction to satisfy boundary conditions (39), yielding the form of the solutions (e.g. Busse, 1986)

[TABLE]

where $\left[\widehat{\Psi},\widehat{\Theta},\widehat{\xi}\right]$ are complex-valued amplitudes, $\lambda$ is the complex eigenvalue with $\Re_{e}(\lambda)=\sigma$ the growth rate, $m$ is the azimuthal wave number and $n$ is the degree of spatial complexity along the horizontal direction. We substitute solutions (41) into equations (36b)-(36c) and (41). They can be recast into a single equation for $\widehat{\Psi}$. This equation can be recast into the original form introduced by Busse (2002), i.e.

[TABLE]

with $a^{2}=m^{2}+n^{2}\pi^{2}$ and by introducing the thermal and compositional Rayleigh numbers in the Busse’s notation

[TABLE]

Because all azimuthal wave numbers $m$ are separated, the marginal stability curve $\sigma=0$ is obtained by minimising the critical Rayleigh number over all values of $m$. In the following, we will survey the properties of RDDC in the annulus geometry by varying Busse’s parameters $(R_{T},R_{C})$.

In our notations, the growth rate $\sigma=\Re_{e}(\lambda)$ of RDDC is predicted by the following polynomial equation (Busse, 2002)

[TABLE]

with $a^{2}=m^{2}+\pi^{2}$, $[Ra_{C},Ra_{T}]$ the Rayleigh numbers (8) and $\beta$ a geometrical parameter in the annulus geometry. When double-diffusive effects are negligible, the onsets of pure thermal or compositional rotating convection are naturally recovered, given by the critical values

[TABLE]

with the function

[TABLE]

C.2 New asymptotic predictions

In the limit $Ek\ll 1$, we have obtained an analytical expression for the double-diffusive onset from formula (46). This contradicts the prediction Busse (2002), made by mistakenly considering the non-rotating limit. Within the double-diffusive tongue, the onset is given by (see details in the supplementary material)

[TABLE]

with negative (respectively positive) values of $Ra_{T}^{c}$ in the finger (respectively semi-convection) regime. Predictions (49) agree very well with our numerical simulations in the sphere (see figure 4). In particular, we recover that the results do not depend on $Sc$ and $Pr$, but only on $L$ and $Ek$. Moreover, for each $m$, the minimum $|Ra_{C}^{c}|$ is located along the line

[TABLE]

and is given by

[TABLE]

This corrected expression of the reduced onset agrees with our numerical results in the sphere (figure 4). Note that we also recover that, near this point, the onset is independent on $L$, and thus only depends on $Ek$.

C.2.1 Matching the annulus to the sphere

Simitev (2011) showed numerically that $n=1$ is always the most unstable radial wave number in the annulus geometry. So, we have fixed $n=1$ in the following, as originally considered by Busse (2002). Then, parameters (37)-(42) are local parameters. Moreover, the latter parameter $\beta$ is constant in the thin-gap approximation. However, the spherical curvature, here measured by $\chi$, is spatially varying in the sphere. For a matching to the sphere, these local parameters should be adjusted at the location of the QG structure at the onset, as schematically illustrated in figure 20. Indeed, $\beta$ strongly depends on the critical cylindrical radius $s_{0}$ at which columnar QG motions first appear, which is known to vary in spheres (Jones et al., 2000). Similarly, $(\widetilde{Ra}_{T},\widetilde{Ra}_{C})$ depend not only on the global Rayleigh numbers $(Ra_{T},Ra_{C})$ introduced in the main text (8), but also on the local position $s_{0}$.

Therefore, $(\widetilde{Ra}_{T},\widetilde{Ra}_{C},\beta)$ are free parameters in the model. To heuristically link the local and global parameters, we introduce one adjustable parameters $\Gamma$ such that

[TABLE]

$\Gamma$ should depend on the dimensionless parameters at the onset, i.e. $\Gamma=\Gamma(Ek,Pr,Sc,Ra_{T},Ra_{C})$. Thus, this parameter is not a priori uniquely determined.

C.2.2 Benchmark with SINGE

We now compare the prediction of the previous model with the actual data given by SINGE. To do so, we have adjusted $\Gamma$ such that the marginal stability curve $\sigma=0$, predicted by (44), coincides with the critical Rayleigh numbers at the onset of pure compositional convection ($Ra_{T}=0$) as computed by SINGE. We show in figure 21 the superposition of the marginal stability curve $\sigma=0$ determined by SINGE and the stability map predicted by equation (44) in the finger quadrant.

Several points are worthy of comment. First, the critical wave number $m^{c}$ in the theory is over-estimated compared to the numerical values in table 2, roughly by a factor three. This confirms that local theories can only predict the order of magnitude of the wave number at the onset (e.g. Busse, 1970). On the marginal stability curve within the double-diffusive tongue, SINGE always find an $m=1$ mode.

Second, the reduced model recovers the non-rotating limit of finger convection. Indeed, the non-rotating limit (18), i.e. $Ra_{T}=-Ra_{C}$, is asymptomatically reached for large enough Rayleigh numbers. Note however that we found convective motion beyond this limit with SINGE for some boundary conditions (see §A).

Finally, double-diffusive effects are over-estimated in the reduced model for unstably stratified fluids (above the dashed-line in figure 21), predicting unstable regions where the system is in fact stable. In addition, in the reduced model, the unstable double-diffusive tongue widens without bound when increasing $L$, whereas it reaches a limit for $L\gtrsim 10^{3}$ in our numerical computations (see figure 4). Quantitatively, these discrepancies increase when $Ek$ decreases.